The top most program in the Model II GCM (Global Circulation Model) is named “Mjal2cpdC9.f”
Now “names matter”, but sometimes more than others. The “.f” tells the compiler (thing that turns a high level language like FORTRAN into the binary stuff the computer actually runs) that this is an old school FORTRAN program and to act accordingly. (Those ending in “.f90” tell it to use the Fortran from 1990 and the newer rules and syntax). Typically “.f” means use the F77 (FORTRAN from 1977) language.
Having read a bunch of the model code now, it seems they have a high “dick with” factor on names. They like to encode metadata into the name itself. In database design, it is considered Very Bad Form to encode data into a name. Why? Because when something changes you get a maintenance task that can induce errors. Here, the encoding part is the ‘jal2cpdC9″. What does it mean? Who knows… There is evidence in the top COMMENT block for lots of changes over many years. I would speculate it is a kind of “version number”, but without the virtue of the typical version schemes that let you know their age order.
I think (speculate) that M stands for MAIN. In any case, it is an easy handle to use for this program.
The Core Of It
Often a program has a bunch of “stuff” that gets in the way of seeing the core function. This one is no different. There’s a couple of pages of comments on historical modifications (we get to that down below) and LOTS of diagnostic saving, along with a fair amount of timing, checkpoint and restart stuff. Nice to have, gets in the way of seeing the basic main loop.
So I’m going to put that main loop, the core of the actual operation of the program, here at the top. This is the guts of it doing the function for which it was designed. The rest is wrapper to set up, tear down, diagnose, restart, or implement complicated functions (usually as subroutines for easier understanding and maintenance). This part gets repeated down below in conclusions (where I first wrote it):
So there is a ‘Dynamic” part “CALL DYNAM” to investigate, then we do the non-DYNAM parts. DYNAM is a set of subroutines in this program, but will be examined separately. I’m cutting the list of non-DYNAM steps down to just the most pertinent lines:
C**** CONDENSTATION, SUPER SATURATION AND MOIST CONVECTION 140. CALL CONDSE 141. CALL PRECIP C**** RADIATION, SOLAR AND THERMAL 149. CALL RADIA 150. C**** SURFACE INTERACTION AND GROUND CALCULATION 157. 400 CALL SURFCE 158. CALL GROUND 159. CALL DRYCNV 160. C**** STRATOSPHERIC MOMENTUM DRAG 166. CALL SDRAG 167. C**** SEA LEVEL PRESSURE FILTER 177. CALL FILTER 181. c ** This code calls the ocean diffusion in the deep ocean mode CALL ODIFS C**** RESTRUCTURE THE OCEAN LAYERS AND ELIMINATE SMALL ICE BERGS 205. CALL OSTRUC 206.
I’ve cut out the other bits from around this core just to make it easy to see what happens in what order.
The Preamble
Starting at the top, it does a lot of ‘setup’. Things like setting the initial state of variables and finding out what time it is (for tracking how long it runs). We’re skipping most of that stuff here.
The default “unit numbers” in FORTRAN are Unit 5 for reading in, unit 6 for writing out. Here it writes out some information (initial time and condition) to the report it makes. But this part of the “setup” is interesting:
WRITE (6,'(A,13X,A,I6,A,F6.2,I6,A5,I27,I7,F7.1,A,F11.2)') 39. * '0CLIMATE MODEL STARTED UP','DAY',IDAY,', HR',TOFDAY, 39.1 * JDATE,JMONTH,MINC,MELSE,PERCNT,' TAU',TAU 39.2 DOPK=1. 40. MODD5K=1000 41. CALL CHECKT (1) 41.4 CF PALMER(36,1)=-1. 41.5 RUNON=-1. 41.6 IF(TAU+TAUERR.GE.TAUE) GO TO 630 42. CF PALMER(36,1)=1. 42.5 RUNON=1. 42.501 RUNON=1. 42.6
Note the CF lines? A “C” or “c” in newer Fortran indicates a comment or inactive line. “PALMER” has been deactivated. One wonders why, but that is for later. Also note it calls “CHECKT”. That is a program that checks temperaures. It is a subroutine of this Main program, down at the bottom. Here’s what it says in comments:
SUBROUTINE CHECKT (N) 9001. C**** 9002. C**** THIS SUBROUTINE CHECKS WHETHER THE TEMPERATURES ARE REASONABLE 9003. C**** FOR DEBUGGING PURPOSES. IT IS TURNED ON BY SETTING IDACC(11) 9004. C**** TO BE POSITIVE. REMEMBER TO SET IDACC(11) BACK TO ZERO AFTER 9005. C**** THE ERRORS ARE CORRECTED. 9006. C**** 9007. INCLUDE 'BA94jalC9.COM' COMMON U,V,T,P,Q 9008.2 IF(IDACC(11).LE.0) RETURN 9009. C**** 9010. C**** CHECK WHETHER GDATA ARE REASONABLE AND CONSISTENT OVER EARTH 9011.
This is a very interesting bit of clue. One wonders what bounds are set on “REASONABLE”. There’s a variable used to turn it on (so what happens if you leave it off and run a result? Is “unreasonable” sometimes reasonable enough to use in scaring the children?) Then, once “errors are corrected” it can be turned back off. So what errors and how does one “correct” them? OK, benefit of the doubt. Typically in all programs you ‘sanity check’ the input data. Assume until shown otherwise that’s all this is and it is looking to find where someone entered 220 instead of 22.0 (or whatever). Still, note to self, it needs checking to see just what it does, and why, and there is a mystery about just how much hand holding and manual intervention is done on model “runs”…
The Main Loop
At the far right they have line numbers. Those are NOT active in the program. Just let you put your deck of punched cards in proper order if you drop it. No, really! I’ve done it! ;-) That they are using fractional numbers is interesting. It shows where lines were added after earlier versions need changing. Here we see the original line 46 got expanded into 8 lines of code.
The actual “line numbers” used in the program are on the left. See that one “98 CONTINUE”? That 98 is a target line number that a “GOTO” can jump to for continued operation. The “CONTINUE” statement itself doesn’t do anything, it just lets you lay down a marker that “something might jump to here later”. It can also be used to emphasize visually that a new block of code starts here, much like the comment block saying “MAIN LOOP”.
At this point, this part of the MAIN LOOP does some more value setting and resetting, but this is very likely to be done many times (unlike the preamble that is typically only done once). It also writes out some “restart” information when it passes through this part of the code again. Nice touch. Checkpoint / Restart code lets you save intermediate points in a given model run and restart from them later. It also lets you take a powerfail and only lose the last cycle of a run. VERY useful. Nice touch.
C**** 43. C**** MAIN LOOP 44. C**** 45. 98 CONTINUE 45.5 IF(USET.LE.0.) GO TO 100 46. IF(.NOT.EVENT(USET)) GO TO 100 46.005 C**** ZERO OUT INTEGRATED QUANTITIES EVERY USET HOURS 46.01 DO 99 K=4,11 46.02 DO 99 J=1,JM 46.03 DO 99 I=1,IM 46.04 99 OA(I,J,K)=0. 46.05 100 IF(.NOT.EVENT(TAUT)) GO TO 200 46.06 C**** WRITE RESTART INFORMATION ONTO DISK 47. 120 CALL RFINAL (IRAND) 48. REWIND KDISK 49. WRITE (KDISK) TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50. c ** MFS (CHANGED) c ** write out deep ocean arrays as well as normal arrays c ** Reto says this is the correct write statment for istart=10 c * RQT,SRHR,TRHR,TSFREZ,(AJ(K,1),K=1,KACC),TDIURN,OA,TAU 51. * RQT,SRHR,TRHR,TSFREZ,TG3M,RTGO,STG3,DTG3, 51. * (AJ(K,1),K=1,KACC),TDIURN,TAU 51.1 c ** END (CHANGED) REWIND KDISK 52. CALL CLOCKS (MNOW) 53. MINC=MLAST-MNOW 54. MELSE=MELSE+MINC 55. PERCNT=100.*MELSE/(MSTART-MNOW+1.D-5) 56. WRITE (6,'(A,I3,55X,2I7,F7.1,A,F11.2)') ' OUTPUT RECORD WRITTEN ON 57. * UNIT',KDISK,MINC,MELSE,PERCNT,' TAU',TAU 57.1 KDISK=3-KDISK 58. MLAST=MNOW 59.
Note that the Reto comment shows change over time to the restart file. I’d guess when ocean diffusion was added, but we’ll find out another day.
Like all good looops, the program then checks if it ought to end now (and potentially does some error handling):
C**** TEST FOR TERMINATION OF RUN 60. 200 READ (3,END=210) LABSSW 61.1 210 REWIND 3 61.2 IF(LABSSW.EQ.LABEL1) KSS6=1 61.3 c ** MFS (ADDED) c ** Check if a file named SSW.STOPGCM exists and stop the model if c ** it does. Errors retuned are both MacOS style and Fortran c ** style. Basically if SSW.STOPGCM is held open then the model c ** will not stop, but if the file is closed then the model c ** will stop. In this version of the code send switch also works. OPEN(3,FILE='SSW.STOPGCM',FORM='unformatted', * STATUS='old',IOSTAT=IOErr) c 0 = NoErr IF(IOErr.EQ.0) THEN KSS6=1 CLOSE(3,IOSTAT=IOErr) END IF c ** END (ADDED) C IF(MBEGIN-MNOW.GT.2880000) KSS6=1 61.5 c ** MFS (CHANGED) c ** For some reason the qflux/sst code used IF(KSS6.EQ.1) GO TO 800. c ** I don't know why the MODRD.EQ.NRAD-NDYN was added. c ** the old code in this model was c IF(KSS6.EQ.1) GO TO 800 62. c ** but the old code in the deep ocean model was c IF(KSS6.EQ.1.AND.MODRD.EQ.NRAD-NDYN) GO TO 800 62. IF(KOCEAN.EQ.2) THEN IF(KSS6.EQ.1.AND.MODRD.EQ.NRAD-NDYN) GO TO 800 ELSE IF(KSS6.EQ.1) GO TO 800 62. END IF c ** END (CHANGED) IF(TAU+TAUERR.GE.TAUE) GO TO 810 63.
Now that block may look complicated, but all it is doing is asking the question “Stop now?” and if so, going to line 800. (see that “IF(KSS6.EQ.1) GO TO 800” line?) So what is the 800 line?
C**** 302. C**** END OF MAIN LOOP 303. C**** 304. C**** RUN TERMINATED BECAUSE SENSE SWITCH 6 WAS TURNED ON 305. 800 WRITE (6,'("0SENSE SWITCH 6 HAS BEEN TURNED ON.")') 306. 810 IF(EVENT(TAUT).OR.TAUT.GE.TAU+1000.) GO TO 820 307. CALL RFINAL (IRAND) 308. REWIND KDISK 309.
So you can see that we’ve got about 250 lines between the start of the actual loop and the place where it ends. (that 309 minus the 63 so far, modulo the .xxx additions). That includes a lot of comments, so not a lot of actually active lines. Similarly the “Go To 810” and references to 820 will also be in that exit handling block.
OK, if we don’t exit, what is done? Well, more resetting variables and arrays so they can be used in the actual work, again. Out with the old, in with the new (data that is…). This is the kind of thing we skipped in the preamble, but it gets done each turn of the MAIN LOOP, so we need to know it here. At line 82 we get a comment block about “prt” files. Near as I can tell, that is what they call “print out files”, or dumps of state to readable media. (Highly useful for further “tuning”…)
C**** IF TIME TO ZERO OUT DIAGNOSTIC ACCUMULATING ARRAYS, DO SO 64. C**** (ALWAYS AT THE BEGINNING OF THE RUN ...... 65. IF(TAU-TAUERR.LE.TAUI) GO TO 260 66. IF(.NOT.EVENT(HR24)) GO TO 300 67. C**** .... AND NORMALLY ALSO AT THE BEGINNING OF EACH MONTH) 67.5 DO 250 K=1,13 68. IF(JDAY.EQ.NDZERO(K)) GO TO 260 69. 250 CONTINUE 70. GO TO 290 71. 260 TAU0=TAU 72. IDAY0=IDAY 73. TOFDY0=TOFDAY 74. JDATE0=JDATE 75. JMNTH0=JMONTH 76. JYEAR0=JYEAR 77. DO 270 I=1,10 78. 270 IDACC(I)=0 79. DO 275 K=1,KACC 80. 275 AJ(K,1)=0. 81. DO 280 J=1,JM 81.1 DO 280 I=1,IM 81.2 AIJ(I,J,76)=-1000. 81.3 AIJ(I,J,77)=1000. 81.4 280 AIJ(I,J,78)=1000. 81.5 C**** INITIALIZE SOME ARRAYS AT THE BEGINNING OF SPECIFIED DAYS 82. c ** MFS (CHANGED) c ** Rename the prt files here to get the correct names and split the c ** prt into monthly files. Because not everything has been c ** initialized yet durring the first hour I have to init a few c ** variables. c 290 IF(JDAY.NE.32) GO TO 294 83. 290 CONTINUE c init variables that I'm ahead of LLAB1 = INDEX(LABEL1,'(') -1 c IF(LABEL1(4:4).EQ.' ') LLAB1=4 WRITE(CYEAR,'(I4)')JYEAR0 c Create the proper file name which is the month plus 1: c ** 9/X (SWITCH) c fName = ':prt:'//JMNTH0(1:3)//CYEAR//'.prt'//LABEL1(1:LLAB1) fName = './prt/'//JMNTH0(1:3)//CYEAR//'.prt'//LABEL1(1:LLAB1) c fName = './'//JMNTH0(1:3)//CYEAR//'.prt'//LABEL1(1:LLAB1) c Open a new file with that name: c OPEN(6,FILE=fName,FILETYPE='TEXT',CREATOR='R*ch') c OPEN(6,FILE=fName,FILETYPE='TEXT',CREATOR='SSPP') OPEN(6,FILE=fName) c ** END (SWITCH) c put back old code IF(JDAY.NE.32) GO TO 294 c ** END (CHANGED) JEQ=1+JM/2 83.01 DO 292 J=JEQ,JM 83.02 DO 292 I=1,IM 83.03 292 TSFREZ(I,J,1)=JDAY 83.04 JEQM1=JEQ-1 83.05 DO 293 J=1,JEQM1 83.06 DO 293 I=1,IM 83.07 293 TSFREZ(I,J,2)=JDAY 83.08 GO TO 296 83.09 294 IF(JDAY.NE.213) GO TO 296 83.1 JEQM1=JM/2 83.11 DO 295 J=1,JEQM1 83.12 DO 295 I=1,IM 83.13 295 TSFREZ(I,J,1)=JDAY 83.14 C**** INITIALIZE SOME ARRAYS AT THE BEGINNING OF EACH DAY 83.4 296 DO 297 J=1,JM 83.5 DO 297 I=1,IM 83.51 TDIURN(I,J,1)=1000. 83.511 TDIURN(I,J,2)=-1000. 83.512 TDIURN(I,J,3)=1000. 83.513 TDIURN(I,J,4)=-1000. 83.514 TDIURN(I,J,5)=0. 83.52 TDIURN(I,J,6)=-1000. 83.521 PEARTH=FDATA(I,J,2)*(1.-FDATA(I,J,3)) 83.53 IF(PEARTH.GT.0.) GO TO 297 83.54 TSFREZ(I,J,1)=365. 83.55 TSFREZ(I,J,2)=365. 83.56
Quite a lot of resetting variables, eh? Folks think programming is complicated. Actually lots of it is dead boring stuff like zero filling arrays…
Now we get into the actual order of functions being done IN the model portion proper. It starts with a target “CONTINUE’ statement so returns to this point are easy. Line 83.54 above jumps here (and skips setting TSFREZ) if PEARTH is set to bigger than zero. I think that may be a geometry setting, but it’s a ‘Dig Here’ for the future.
This first bit does some diagnostics saving, calls that CHECKT subroutine to make sure temperatures are sane, calls CLOCKS to do the time tracking (for how long things run) and then we get to the next block where it starts actually doing things.
Notice how much of the time is “setting the table” and how little is actually doing?… Programming is like that.
297 CONTINUE 83.57 C**** 84. C**** INTEGRATE DYNAMIC TERMS 85. C**** 86. 300 MODD5D=MOD(NSTEP,NDA5D) 87. IF(MODD5D.EQ.0) CALL DIAG5A (2,0) 88. IF(MODD5D.EQ.0) CALL DIAG9A (1) 88.5 CALL DYNAM DOPK=1. 123. CALL CHECKT (2) 123.5 CALL CLOCKS (MNOW) 124. MINC=MLAST-MNOW 125. MDYN=MDYN+MINC 126. MLAST=MNOW 127. PERCNT=100.*MDYN/(MSTART-MNOW) 128. C WRITE (6,'(A,I1,3X,A,I6,A,F6.2,I6,A5,2I7,F7.1,21X,A,F11.2)') 129. C * ' DYNAMIC TERMS INTEGRATED, MRCH=',MRCH,'DAY',IDAY, 129.1 C * ', HR',TOFDAY,JDATE,JMONTH,MINC,MDYN,PERCNT,'TAU',TAU 129.2 IF(MODD5D.EQ.0) CALL DIAG5A (7,NDYN) 130. IF(MODD5D.EQ.0) CALL DIAG9A (2) 131. IF(EVENT(HR12)) CALL DIAG7A 132.
And, at last, we get to the actual working core of the MAIN LOOP… AFTER it does some diagnostic subroutine calls itself… Sigh. I note in passing that line 140.5 has been commented out.
So what does it do? In order: Condensation, ANOTHER CHECKT, Precipitation, ANOTHER CHECKT, and time stamp. Then yet more diagnostics prep / dump. This thing has a LOT of TLC and hand holding needed / possible…
C**** 133. C**** INTEGRATE SOURCE TERMS 134. C**** 135. MODRD=MOD(NSTEP,NRAD) 136. MODD5S=MOD(NSTEP,NDA5S) 137. IF(MODD5S.EQ.0) IDACC(8)=IDACC(8)+1 138. IF(MODD5S.EQ.0.AND.MODD5D.NE.0) CALL DIAG5A (1,0) 139. IF(MODD5S.EQ.0.AND.MODD5D.NE.0) CALL DIAG9A (1) 139.5 C**** CONDENSTATION, SUPER SATURATION AND MOIST CONVECTION 140. C IF(MOD(NSTEP,NCNDS).NE.0) GO TO 400 140.5 CALL CONDSE 141. CALL CHECKT (3) 141.5 CALL PRECIP 142. CALL CHECKT (4) 142.5 CALL CLOCKS (MNOW) 143. MINC=MLAST-MNOW 144. MCNDS=MCNDS+MLAST-MNOW 145. MLAST=MNOW 146. IF(MODD5S.EQ.0) CALL DIAG5A (9,NCNDS) 147. IF(MODD5S.EQ.0) CALL DIAG9A (3) 148.
So that’s the basic structure of these parts of the MAIN LOOP. Diagnostics, couple of physics routines (with constant sanity checking on the temperatures), more diagnostics, and occasional clock watching.
Next up is the Radiation part:
C**** RADIATION, SOLAR AND THERMAL 149. CALL RADIA 150. CALL CHECKT (5) 150.5 CALL CLOCKS (MNOW) 151. MINC=MINC+MLAST-MNOW 152. MRAD=MRAD+MLAST-MNOW 153. MLAST=MNOW 154. IF(MODD5S.EQ.0) CALL DIAG5A (11,NCNDS) 155. IF(MODD5S.EQ.0) CALL DIAG9A (4) 156.
Pretty simple, isn’t it? Do the radiative model part, check you didn’t make insane temperatures, note the time, do diagnostics if that is set, move on.
Next is surfaces and the ground:
C**** SURFACE INTERACTION AND GROUND CALCULATION 157. 400 CALL SURFCE 158. CALL CHECKT (6) 158.5 CALL GROUND 159. CALL CHECKT (7) 159.5 CALL DRYCNV 160. CALL CHECKT (8) 160.5 CALL CLOCKS (MNOW) 161. MINC=MINC+MLAST-MNOW 162. MSURF=MSURF+MLAST-MNOW 163. MLAST=MNOW 164. IF(MODD5S.EQ.0) CALL DIAG9A (5) 165.
Note that DRYCNV means Dry Convection is handled as part of the ground layer.
Next we get the Stratosphere and that whole “momentum” and “drag” bit. There is also a commented out write statement that looks like a ‘diagnostic write’ to me. I use them a lot. At some point in the program, during debugging, you have it tell you what it is doing with some variables. Latter you comment it out for regular runs, but leave it there for future debugging if any is needed.:
C**** STRATOSPHERIC MOMENTUM DRAG 166. CALL SDRAG 167. CALL CHECKT (9) 167.5 CALL CLOCKS (MNOW) 168. MINC=MINC+MLAST-MNOW 169. MSURF=MSURF+MLAST-MNOW 170. MLAST=MNOW 171. IF(MODD5S.EQ.0) CALL DIAG5A (12,NCNDS) 172. IF(MODD5S.EQ.0) CALL DIAG9A (6) 173. MSRCE=MCNDS+MRAD+MSURF 174. PERCNT=100.*MSRCE/(MSTART-MNOW) 175. C WRITE (6,'(A,64X,2I7,F7.1)') 176. C * ' SOURCE TERMS INTEGRATED', MINC,MSRCE,PERCNT 176.1
Here we get a “filter”. I suspect they had “issues” with sea level excursions and needed to crowbar them flatter when extreme, but that is just speculation. One would need to look at the FILTER subroutine to see what it does, and from that infer why. (Or maybe it’s in the Hansen ’83 paper somewhere…) Note that the first line skips this ‘patch’ of sea level pressure if you have MFILTR set to below zero.. or are on a modulo of NSTEP of it… That jumps to line 500 just below where date and time are updated.
C**** SEA LEVEL PRESSURE FILTER 177. IF(MFILTR.LE.0.OR.MOD(NSTEP,NFILTR).NE.0) GO TO 500 178. IDACC(10)=IDACC(10)+1 179. IF(MODD5S.NE.0) CALL DIAG5A (1,0) 180. CALL DIAG9A (1) 180.5 CALL FILTER 181. CALL CHECKT (10) 181.5 CALL CLOCKS (MNOW) 182. MDYN=MDYN+MLAST-MNOW 183. MLAST=MNOW 184. CALL DIAG5A (14,NFILTR) 185. CALL DIAG9A (7) 186.
That’s the core of the guts of it. They did a ‘glue on’ of ocean deep diffusion that shows up down below
This next bit fiddles the time. Again we have some diagnostics calls if set, and some clock watching. It looks to me like the ocean deep diffusion was added at a long time scale, so outside this day cycle. I need to actually work through when date values change to figure that out for sure. It might also be that they do “land and air” then oceans just as two different spaces.
C**** 187. C**** UPDATE MODEL TIME AND CALL DAILY IF REQUIRED 188. C**** 189. 500 NSTEP=NSTEP+NDYN 190. ITAU=(NSTEP+NSTEP0)*IDTHR 191. TAU=DFLOAT(ITAU)/XINT 192. IDAY=1+ITAU/I24 193. TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 194. IF(.NOT.EVENT(HR24)) GO TO 590 195. CALL DIAG5A (1,0) 196. CALL DIAG9A (1) 196.5 CALL DAILY 197. CALL CLOCKS (MNOW) 198. MELSE=MELSE+(MLAST-MNOW) 199. MLAST=MNOW 200. NDAILY=SDAY/DT 201. CALL DIAG5A (16,NDAILY) 202. CALL DIAG9A (8) 203. DO 530 J=1,JM 203.21 IMAX=IM 203.22 IF((J.EQ.1).OR.(J.EQ.JM)) IMAX=1 203.23 DO 530 I=1,IMAX 203.24 TSAVG=TDIURN(I,J,5)/(HR24*NSURF) 203.31 IF(32.+1.8*TSAVG.LT.65.)AIJ(I,J,52)=AIJ(I,J,52)+(33.-1.8*TSAVG) 203.32 AIJ(I,J,54)=AIJ(I,J,54)+18.*((TDIURN(I,J,2)-TDIURN(I,J,1)) 203.41 * /(TDIURN(I,J,4)-TDIURN(I,J,3)+1.D-20)-1.) 203.42 AIJ(I,J,30)=AIJ(I,J,30)+(TDIURN(I,J,4)-TDIURN(I,J,3)) 203.425 AIJ(I,J,80)=AIJ(I,J,80)+(TDIURN(I,J,4)-273.16) 203.426 IF (TDIURN(I,J,6).LT.AIJ(I,J,78)) AIJ(I,J,78)=TDIURN(I,J,6) 203.48 530 CONTINUE 203.51 c ** MFS (CHANGED) c ** Old code assumed that kocean 1 was kocean = 0, which is no longer c ** true. c IF(KOCEAN.NE.1) GO TO 590 204. IF(KOCEAN.EQ.0) GO TO 590 c ** END (CHANDED) DO 540 J=1,JM 204.3 DO 540 I=1,IM 204.4 AIJ(I,J,59)=AIJ(I,J,59)+ODATA(I,J,4) 204.5 540 AIJ(I,J,60)=AIJ(I,J,60)+ODATA(I,J,5) 204.6
NOW we do the oceans. First ocean diffusion, then ‘restructuring’ and for some reason they eliminate small icebergs. One would think icebergs matter, but who knows… The comments indicate this part was a ‘glue on’ to the basic model at some later time. KOcean lets you choose to use it or not. Set to “2” one calls ODIFS, if not, you skip it. CALL OSTRUCT ‘restructures’ the ocean, then again a sanity check on temps and clock watching.
c ** MFS (ADDED) c ** This code calls the ocean diffusion in the deep ocean mode ( c ** kOcean = 2). The diffusion code lets heat move downward from c ** the bottom of the mixed layer into the rest of the ocean. c ** in the deep ocean code the old line was c CALL ODIFS 204.61 IF(KOCEAN.EQ.2)THEN CALL ODIFS END IF c ** END (ADDED) C**** RESTRUCTURE THE OCEAN LAYERS AND ELIMINATE SMALL ICE BERGS 205. CALL OSTRUC 206. CALL CHECKT (11) 206.5 CALL CLOCKS (MNOW) 207. MSURF=MSURF+(MLAST-MNOW) 208. MLAST=MNOW 209.
With that, it is basically done. Just write out the information every USET hours of cycle time and “wash and repeat”. One is left wondering why “Reto said that tauo never worked”, and what tauo is, and why you needed it, and what it did wrong when it never worked… but it is fixed now, so Nevermind…
C**** 210. C**** WRITE INFORMATION ONTO A TAPE EVERY USET HOURS 211. C**** 212. 590 CONTINUE 212.5 c ** RAR (CHANGED) c ** Reto said that tauo never worked because the if was missing here. c ** "I guess we saved an if statement." This has now been fixed so c ** that at the end of the first day after tauo ocean data is c ** collected. Code works now, collection starts on the first day of c ** TAUO (the +.5 is to prevent the previous day from being saved). IF(USET.LE.0.) GO TO 600 213. IF((TAUO+.5).GT.TAU) GO TO 600 c ** END (CHANDED) IF(.NOT.EVENT(USET/2.)) GO TO 600 214. IF(EVENT(USET)) GO TO 552 214.5 C**** SAVE INSTANTANEOUS QUANTITIES 215. DO 551 J=1,JM 215.1 DO 551 I=1,IM 215.2 OA(I,J,1)=.1*916.6+GDATA(I,J,1) 215.3 OA(I,J,2)=GDATA(I,J,3) 215.4 551 OA(I,J,3)=GDATA(I,J,7) 215.5 GO TO 600 215.6 552 WRITE (20) TAU,OA 216. c ** RAR (CHANGED) c ** old code to force an immediate write to tape, this doesn't work c ** anymore. c ENDFILE 20 216.5 c BACKSPACE 20 216.6 c ** END (CHANDED) CALL CLOCKS (MNOW) 217. MINC=MLAST-MNOW 218. MELSE=MELSE+MINC 219. PERCNT=100.*MELSE/(MSTART-MNOW) 220. COD WRITE (6,'(A,78X,A,F11.2)') ' INFORMATION WRITTEN ON UNIT 20', 221. COD * ' TAU',TAU 221.1
Having finished that cycle, a whole LOAD of Diagnostic data get dumped. This just shouts “unstable model” and “lots of manual tuning and intervention” to me, but that will have to wait for confirmation in the running…
C**** 222. C**** CALL DIAGNOSTIC ROUTINES 223. C**** 224. 600 IF(MOD(NSTEP+NDA4,NDA4).EQ.0) CALL DIAG4A 225. CF RUNON=PALMER(36,1) 225.1 IF(USESLP.EQ.0.) GO TO 603 225.7 IF(USESLP.LT.0..AND.EVENT(-USESLP)) USESLP=-USESLP 225.8 IF(USESLP.GT.0..AND.EVENT(USESLP)) CALL DIAG10(0) 225.9 603 CONTINUE 226. IF(NDPRNT(1).GE.0) GO TO 610 227. C**** PRINT CURRENT DIAGNOSTICS (INCLUDING THE INITIAL CONDITIONS) 228. IF(KDIAG(1).LT.9) CALL DIAGJ 229. IF(KDIAG(2).LT.9) CALL DIAGJK 229.5 IF(KDIAG(2).LT.9) CALL DIAGJL 230. IF(KDIAG(7).LT.9) CALL DIAG7P 231. IF(KDIAG(3).LT.9) CALL DIAGIJ 232. IF(KDIAG(9).LT.9) CALL DIAG9P 233. IF(KDIAG(5).LT.9) CALL DIAG5P 234. IF(KDIAG(4).LT.9) CALL DIAG4 235. IF(TAU.LE.TAUI+DTHR*(NDYN+.5)) CALL DIAGKN 236. NDPRNT(1)=NDPRNT(1)+1 237. IF(TAU.LE.TAUI+DTHR*(NDYN+.5)) GO TO 610 237.1 C**** RESET THE UNUSED KEYNUMBERS TO ZERO 237.4 DO 605 I=1,42 237.5 605 KEYNR(I,KEYCT)=0 237.6 610 IF(.NOT.EVENT(HR24)) GO TO 690 238. C**** PRINT DIAGNOSTIC TIME AVERAGED QUANTITIES ON NDPRNT-TH DAY OF RUN 239. DO 620 K=1,13 240. IF(JDAY.EQ.NDPRNT(K)) GO TO 630 241. 620 CONTINUE 242. GO TO 640 243. 630 WRITE (6,'("1"/64(1X/))') 244. IF(KDIAG(1).LT.9) CALL DIAGJ 245. IF(KDIAG(2).LT.9) CALL DIAGJK 245.5 IF(KDIAG(2).LT.9) CALL DIAGJL 246. IF(KDIAG(7).LT.9) CALL DIAG7P 247. IF(KDIAG(3).LT.9) CALL DIAGIJ 248. IF(KDIAG(9).LT.9) CALL DIAG9P 249. IF(KDIAG(5).LT.9) CALL DIAG5P 250. IF(KDIAG(6).LT.9) CALL DIAG6 251. IF(KDIAG(4).LT.9) CALL DIAG4 252. C**** THINGS TO DO BEFORE ZEROING OUT THE ACCUMULATING ARRAYS 254. C**** (NORMALLY DONE AT THE END OF A MONTH) 255. 640 DO 650 K=1,13 256. IF(JDAY.EQ.NDZERO(K)) GO TO 660 257. 650 CONTINUE 258. GO TO 690 259. C**** PRINT THE KEY DIAGNOSTICS 260. 660 CALL DIAGKN 261. IF(RUNON.EQ.-1.) STOP 13 261.5 C**** PRINT AND ZERO OUT THE TIMING NUMBERS 262. CALL CLOCKS (MNOW) 263.
The End of the Core Of Main
I’m not going to put the rest of that Mjal2cpdC9.f program in here. It is substantially some setting up / zeroing data items at the start, and subroutines at the bottom. Those subroutines are best handled separately.
This program is the main flow of control program wrapped around the others. We see how the model runs the big lumps, and what is considered in what order. The rest is largely physics, math, and details. Those do matter, but best handled as their own bits.
Here we got confirmation that the DIAGxx subroutines are not part of the model proper, but gather and handle diagnostic data. Useful to look at “someday”, but not needed for a first pass understanding of the model itself. We also see that CLOCKS is just watching how much time is used, so nice housekeeping but not the model itself.
For each of the main actions, we now know the order of execution, and can look at each one does a bit more deeply. We can also look at what is NOT done…
So there is a ‘Dynamic” part “CALL DYNAM” to investigate, then we do the non-DYNAM parts. DYNAM is a set of subroutines in this program, but will be examined separately. I’m cutting the list of non-DYNAM steps down to just the most pertinent lines:
C**** CONDENSTATION, SUPER SATURATION AND MOIST CONVECTION 140. CALL CONDSE 141. CALL PRECIP C**** RADIATION, SOLAR AND THERMAL 149. CALL RADIA 150. C**** SURFACE INTERACTION AND GROUND CALCULATION 157. 400 CALL SURFCE 158. CALL GROUND 159. CALL DRYCNV 160. C**** STRATOSPHERIC MOMENTUM DRAG 166. CALL SDRAG 167. C**** SEA LEVEL PRESSURE FILTER 177. CALL FILTER 181. c ** This code calls the ocean diffusion in the deep ocean mode CALL ODIFS C**** RESTRUCTURE THE OCEAN LAYERS AND ELIMINATE SMALL ICE BERGS 205. CALL OSTRUC 206.
That’s the real core of it. So to understand what it does in detail, those 9 subroutines are the places to explore.
This level is mostly wrapper around that sequence. One wonders why they start with condensation, when IMHO it starts with radiation from the sun driving evaporation, but in the cycle of things, that may not be too important.
Were I designing it, I’d follow the flow of energy from the sun, to ground and oceans, and then into evaporation / condensation / precipitation. But that’s just me…
They seem to conflate solar input with “thermal” (i.e. IR ‘backradiation’) stuff into one subroutine. IMHO they ought to be treated as distinct, but again, not my design.
Then stratosphere, drag, and ocean diffusion glued on and with sea level pressures optionally patched over.
Overall, it seems “unphysical” in the approach, but we’ll have to see if the execution saves it.
With that, I’m taking a lunch break…
Next step will be to “Dig Here!” into some of those subroutines…
Epilog
Here is a straight copy of the comments from the maintenance folks over the years. It is interesting to read just for the insight into their approach and thinking. Presented without further commentary, but just noting in passing that the code seems to use a 1958 base year atmosphere and has had a lot of changes over the years. Note that this detailed commentary starts LONG after the 1983 origination date, so much change is not recorded…
cat Mjal2cpdC9.f C DYNAM broken into subroutines c ********************************************************************* c ********************************************************************* c ** c ** Model IImac c ** Based on GCMII code for IBM RS/6000 computers created at GISS c ** Modified to compile under Absoft Pro Fortran 8.0 for MacOS 9/X and c ** Pro Fortran 8.2 for MacOS X. Also compiles under IRIX & Linux. c ** Based on MP030bC9 c ** c ** CHANGE HISTORY: c ** c ** 09/26/01 Suki identified why we want to use MP030bC9 (SXS) c ** 09/26/01 minimal changes to get it to compile (MFS) c ** 09/28/01 put in some changes to make IO work (MFS) c ** 09/29/01 fixed critical bug with namelist (MFS) c ** 10/01/01 works (MFS) c ** 11/26/01 put back better stop messages, fixed date printout (MFS) c ** 11/27/01 deep ocean, kocean=2 (MFS/GLR) c ** 01/11/02 failed attempt to fix things (reverted) (MFS) c ** 01/14/02 fixed error message (MFS) c ** 01/15/02 fixed restart (MFS/RAR) c ** 02/05/02 restart bug fixed (RAR) c ** 02/07/02 :oda: -> :ocean: (MFS) c ** 02/13/02 OHT multiplier (deleted) (MFS/GLR) c ** 02/14/02 trend code, deep ocean fixes (MFS) c ** 02/19/02 trend code with debug, fixes (MFS) c ** 02/27/02 new trend code, data file trends only (MFS) c ** 05/27/02 try old S0X & old trends (MSC) c ** 06/03/02 put back mostly old trend code (MSC/MFS) c ** 06/07/02 add ktrend to the namelist (MSC/MFS) c ** 06/12/02 more trend fixes, path fixes (MFS) c ** 06/14/02 found year bug that killed trends (MFS) c ** 06/17/02 Andy says keep 1958 as the default for trends (MFS/AAL) c ** 06/27/02 trend still don't work, remove (MFS) c ** 06/28/02 older radiation, 20 char run numbers (MFS) c ** 07/01/02 pass ktrend again to keep it in scope (MFS) c ** 07/03/02 allow you to bypass trend code (MFS) c ** 07/23/02 Xserve test version (MFS) c ** 07/24/02 call to writer when trends change (MFS/RAR) c ** 09/18/02 disable calls to writer (MFS) c ** 09/26/02 command line version for MacOS X/IRIX (MFS) c ** 09/28/02 put in Jeff's RANVAX to replace Power series RANDU (MFS) c ** 04/08/04 2nd order version, modular dynamics, diag fixes, c ** observed ocean mode, drag preserves wind direction (JAL) c ** 04/09/04 merged in changes from MP030b (MFS/JAL) c ** 05/06/04 setpath for PC support (MFS) c ** 05/07/04 options experiment, works (MFS) c ** 05/19/04 climatological change (JAL) c ** 05/20/04 eccentricity, axial tilt, and omegat (JAL) c ** 05/24/04 removed tolerance for dates (MAC) c ** 06/22/04 try xlf, remove Absoft io options, iint to int (MFS) c ** 06/24/04 more xlf options (MFS) c ** 07/15/04 mrwe_clear to fix PC buffer overflow (MFS) c ** 01/27/05 work on German crash and CLOCKS to MCLOCK (MFS) c ** 04/07/05 solar trend fix and sync with PC for EdGCM 2.3.6 (MFS) c ** 04/26/05 use DLAT to fix bug with 8x10 in RADIA0 (GLR/MAC) c ** 07/22/05 fixes for Pro Fortran 9.2 & command line (MFS) c ** 07/28/05 fix for surface T(I,J,1)=TH1 bug (JAJ) c ** 10/17/05 fix added SHW and SHI from Previp (JAJ) c ** 10/27/05 test no stop ssw (MFS/MAC/JAJ) c ** 12/19/05 test new lunar frontend (MFS) c ** 01/17/06 reverted CLOCKS to support Windows, reverted surface (MFS) c ** 03/23/06 fix print to behave on Windows and readme2 (MFS) c ** c ** NOTES: c ** Suki redeemed, she gave us MP030bC9 originally because of fixes c ** not present in MP008C9. Since MP030bC9 is stable and MP008C9 is c ** no I have to go in and add in the deep ocean mode to MP030bC9. c ** c ** old MP030bC9 c ** 07/29/99 Applied SGI patch for Read(5) (MFS) c ** 08/01/99 IO changes for VAX compadiblity (MAC) c ** 08/04/99 redirected standard out to a file for speed (MFS) c ** 08/06/99 separate and correctly name montly prt files (MFS) c ** 08/06/99 namelist=unit 110, restart=unit 111 (MFS) c ** 08/13/99 finally fixed the prt files (MFS) c ** 01/08/00 fixed the rsf problem, optimized (MFS) c ** 01/09/00 fix to overwrite old output files, need for restart (MFS) c ** 03/06/00 allowed model to start with old startup.prt file (MFS) c ** 03/06/00 fixed bad read 9 to 111 (MFS) c ** 06/05/00 changed to new folder structure (MFS) c ** 08/08/00 changed code to use 'ssw.stopgcm' instead of unit 3 (MFS) c ** 10/04/00 fixed uset and useslp to output into the oda folder (MFS) c ** 10/05/00 changed name of file and some output (MFS) c ** c ** old MP008macC9 c ** 12/20/00 merged MP030aC9, MP008C9 and MA94CC9 (MAC) c ** 12/21/00 start to put together new version (MFS) c ** 12/26/00 put back the modifications to MP030bC9 that c ** David and I made from 07/27/99-10/04/00 (MFS/DCH/MAC) c ** 12/28/00 finish IO modifications (MFS/MAC) c ** 12/28/00 updated the error messages (MFS) c ** 12/29/00 try to read in old restart files (MFS) c ** 01/02/01 hacked istart 300 to read in old files (MFS) c ** 01/03/01 hacked acc files to work with old post-processors (MFS) c ** 01/03/01 started changes for kocean = 2 (MFS) c ** 01/05/01 inserted changes based on differcing MP030bC9 (MFS) c ** 01/23/01 changed unit 111 to 109 (originally 9) (MFS) c ** 01/25/01 fixed unit 20 to use unformatted (MFS) c ** 02/10/01 changed istart=5 to work for paleo runs (gdata) c ** also changes to make deep ocean work (MFS) c ** 02/15/01 version given to Jeff (MFS) c ** 02/20/01 changed TAUT to 240 (MFS) c ** 02/21/01 all istarts distingish kocean = 2 (MFS) c ** 03/02/01 fix flux collection (MFS) c ** 03/12/01 impliment tauo (RAR) c ** 03/13/01 write out tau along with daily restart file (MFS) c ** 04/06/01 debug info for ktrend, volcanos are disabled (MFS) c ** 04/08/01 pass ktrend explicitly (MFS) c ** 04/09/01 try use JYEAR as reference instead of 1958 (MFS) c ** 04/12/01 added INPUTFORCINGS, removed stuff from INPUTZ (MFS) c ** 04/13/01 use new com file BRFRCmac.com (MFS) c ** 04/16/01 eliminate ktrend (MFS) c ** 04/16/01 trends finally work again, fixed comments (MFS) c ** 04/26/01 variable orbit parameters, eg different length days (MFS) c ** 04/27/01 take Reto's advice and change length of DT (MFS) c ** 04/30/01 use the new and improved ORBIT subroutine (GLR) c ** 04/30/01 modify ORBIT to handle variable orbits (MFS) c ** 05/01/01 added distance of orbit in solar units (MFS) c ** 05/02/01 fixed namelist to match orbitial parameters (MFS) c ** 05/07/01 fixed a few more date bugs (MFS) c ** 05/08/01 allow misalligned restarts (MFS) c ** 05/09/01 debug orbitial param (MFS) c ** 05/10/01 fix ocean bug with variable length years (MFS) c ** 05/11/01 more date corrections for JDAY (MFS) c ** 05/14/01 pass copies of orbitial param to ORBIT (MFS) c ** 05/21/01 problems with S0X (MFS) c ** 05/22/01 changed S0X calc to always use real units (MFS) c ** 05/24/01 change OHT strength via namelist (GLR) c ** 05/25/01 fixed some S0X bugs (MFS) c ** 06/15/01 split into Classic and MacOSX versions, this is X (MFS) c ** 06/15/01 altered paths to work with Unix (MFS) c ** 06/21/01 use randvax to use different random numbers, problems (MFS) c ** 06/22/01 fix special regions, don't read in unit 29 (MFS) c ** 07/03/01 revered to non-modular dynamics, readded tauo fix (MFS) c ** 07/10/01 put back old random numbers (MFS) c ** 07/11/01 try using dt multiplier of 4 (MFS) c ** 07/13/01 orbital parameter fixes (MFS) c ** 07/16/01 updates for Carbon version (MFS) c ** 07/23/01 put back Reto's stuff (MFS) c ** 07/24/01 IDACC checks for all debug write statements (MFS) c ** 07/24/01 put in 9/X switches to quickly change versions (MFS) c ** 07/24/01 moved inittrend to be after DAILY0 so JYEAR is right (MFS) c ** 07/27/01 updated daily message (MFS) c ** 07/31/01 this branch does not have the AIX revisions (MFS) c ** 08/02/01 put in Reto's revisions to checkt and other changes (RAR) c ** 08/10/01 removed debug write to radiation that filled up disk (MFS) c ** 08/21/01 fixed bug with long labels and file names (MFS) c ** 08/28/01 removed buggy oht multiplier (MFS) c ** c ** MP008macC9 in Private Beta 1 c ** 12/20/00 merged MP030aC9, MP008C9 and MA94CC9 (MAC) c ** 12/21/00 start to put together new version (MFS) c ** 12/26/00 put back the modifications to MP030bC9 that c ** David and I made from 07/27/99-10/04/00 (MFS/DCH/MAC) c ** 12/28/00 finish IO modifications (MFS/MAC) c ** 12/28/00 updated the error messages (MFS) c ** 12/29/00 try to read in old restart files (MFS) c ** 01/02/01 hacked istart 300 to read in old files (MFS) c ** 01/03/01 hacked acc files to work with old post-processors (MFS) c ** 01/03/01 started changes for kocean = 2 (MFS) c ** 01/05/01 inserted changes based on differcing MP030bC9 (MFS) c ** 01/23/01 changed unit 111 to 109 (originally 9) (MFS) c ** 09/17/01 few more changes to file names via SWITCH (MFS) c ** 09/24/01 file for EdGCM Private Beta 1 (MFS) c ** c ********************************************************************* c ********************************************************************* C**** MP030bC9 BA94jalC9 MP000C9 1/27/92 0.1 C**** basic model II with 1958 Atmosphere and mean strat aerosols (.12) 0.21 C**** This update includes the correct orbit as well as 5 harmonics 0.22 C**** for the ocean heat transport. 0.23 C**** 0.24 C**** Thinner ocean ice with smaller leads: minimal thickness Z=.5m 0.3 C**** rather than 1m, minimal leads = .06*(1/Z-1/5.) rather than .1/Z. 0.4 C**** If warm ocean melts ice, it reduces the ice depth. 0.5 C**** The mixed layer ocean temperature is not affected by the sensible 0.6 C**** heat of precipitation and evaporation (mixed layer conserves 0.71 C**** energy). All computations are done in double precision. 0.8 C**** subroutines in MB08M9 : PRECIP,GROUND,DAILY,OSTRUC 0.85 C**** 0.9 C**** MP005C9 BA94jalC9 MA94M9 03/06/90 0.1 C**** 0.3 C**** basic model II (OA,PALMER,AIJL omitted) .5 box longitude shift 0.4 C**** subroutines in MP005M9: all model routines 0.5 ***** MA94M9(NewModel) BA94M9(CommonBlock) MA94M9(OldSource) 07/19/91 0.1 ***** OPT(3) 0.2 ***** S35V MOD modified for double precision on work station 2/11/92 0.3 ***** MODEL II but with 1958 Freons and mean strat aerosols (.12) 0.4 ***** corrected line 6361.5 (diag) on 2/2/93 0.5 INCLUDE 'BA94jalC9.COM' 1. c ** MFS (ADDED) INCLUDE 'FORCINGSmac.COM' c ** extra ocean arrays for deep ocean COMMON/OCN/TG3M(IM,JM,12),RTGO(IM,JM,LM),STG3(IM,JM), 1.1 * DTG3(IM,JM) 1.2 c ** END (ADDED) COMMON U,V,T,P,Q 2. COMMON/WORKO/OA(IM,JM,11) C**** 7.12 C**** DATA SAVED IN ORDER TO CALCULATE OCEAN TRANSPORTS 7.13 C**** 7.14 C**** 1 ACE1I+SNOWOI (INSTANTANEOUS AT NOON GMT) 7.15 C**** 2 TG1OI (INSTANTANEOUS AT NOON GMT) 7.16 C**** 3 TG2OI (INSTANTANEOUS AT NOON GMT) 7.17 C**** 4 ENRGP (INTEGRATED OVER THE DAY) 7.18 C**** 5 SRHDT (INTEGRATED OVER THE DAY) 7.19 C**** 6 TRHDT (FOR OCEAN, INTEGRATED OVER THE DAY) 7.2 C**** 7 SHDT (FOR OCEAN, INTEGRATED OVER THE DAY) 7.21 C**** 8 EVHDT (FOR OCEAN, INTEGRATED OVER THE DAY) 7.22 C**** 9 TRHDT (FOR OCEAN ICE, INTEGRATED OVER THE DAY) 7.23 C**** 10 SHDT (FOR OCEAN ICE, INTEGRATED OVER THE DAY) 7.24 C**** 11 EVHDT (FOR OCEAN ICE, INTEGRATED OVER THE DAY) 7.25 C****
Very interesting observations about the code.
As you say some hints about where the knobs and levers are that can be manipulated to “tweak” the model to get good fit with historical data.
I find the change comments a bit more terse than what most programmers I know do.
“Fixed broken do-hicky”
Is not as useful, as something like :
Corrected broken For loop in routine do-hicky to properly index the counter ( one off error in counter)
If you are going to bother doing change comments at least make them verbose enough that they give a clue what problem was fixed without having inside historical knowledge of the coding development.
Yeah… nothing like “Fixed bugs” from the other guy to make your job easier when you need to back out that fix…
Pingback: GCM Model II – DYNAM and the Dynamic Subroutines | Musings from the Chiefio
Heh, dropping the card deck. I’ve done it, too; but when I learned keypunching way back when, I learned very early to bind the deck. That way, when I dropped it, at least the cards wouldn’t fly everywhere and were less likely to be damaged, necessitating punching a replacement.