OK, it’s time to list and figure out the Subroutine structure of Model II.

I’m doing this on a new download onto the Devuan chip using userid pi.

This is so that the structure is not tied to my userid.

This is so that (if I’m lucky) I can offer clones of the chip to any Pi user. We’ll see…

First, we just find all instances of SUBROUTINE in the FORTRAN code. Then we’ll make a very reduced list (one line per subroutine defintion) Don’t let this large block of text scare you. I break it down to smaller bits below. This is just the raw data capture from the source code.

As I remove the comments down below, this lets you have a reverence to see a bit of the comments near the active code.

pi@raspberrypi:/Climate/GCM/modelII/trunk $ grep SUB *.f DB11pdC9.f: SUBROUTINE DIAGA (U,V,T,P,Q) DB11pdC9.f: SUBROUTINE DIAGB (U,V,T,P,Q) DB11pdC9.f:c ** ZERO OUT SUBSURFACE VERTICAL WINDS DB11pdC9.f: SUBROUTINE DIAGJ DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES AREA WEIGHTED STATISTICS OF DB11pdC9.f:c ** TITLES FOR SUBROUTINE DIAGJK DB11pdC9.f: SUBROUTINE DIAGJK DB11pdC9.f: SUBROUTINE JKMAP (NT,PM,AX,SCALE,SCALEJ,SCALEK,KMAX,JWT,J1) DB11pdC9.f:c ** TITLES FOR SUBROUTINE DIAGJL DB11pdC9.f: 2'SUBGRID SCALE TEMPERATURE DEVIATION (HUNDREDTHS OF DEGREES KEL.)' DB11pdC9.f: SUBROUTINE DIAGJL DB11pdC9.f:c ** SUBGRID SCALE TEMPERATURE DEVIATION DB11pdC9.f: SUBROUTINE JLMAP (NT,PL,AX,SCALE,SCALEJ,SCALEL,LMAX,JWT,J1) DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES LAYER BY LATITUDE TABLES ON THE LINE DB11pdC9.f:c ** TITLES FOR SUBROUTINE DIAGIL DB11pdC9.f: SUBROUTINE DIAGIL DB11pdC9.f: SUBROUTINE ILMAP (NT,PL,AX,SCALE,SCALEL,LMAX,JWT,ISHIFT) DB11pdC9.f:c ** TITLES FOR SUBROUTINE DIAG7 DB11pdC9.f: SUBROUTINE DIAG7A DB11pdC9.f:c ** THIS SUBROUTINE ACCUMULATES A TIME SEQUENCE FOR SELECTED DB11pdC9.f: SUBROUTINE MEM (SERIES,ITM,MMAX,NUAMAX,NUBMAX,POWER,FPE,VAR,PNU) DB11pdC9.f: SUBROUTINE DIAGIJ DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES LATITUDE BY LONGITUDE MAPS OF DB11pdC9.f: SUBROUTINE IJMAP (NT,ARRAY,BYIACC) DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES NUMERICAL LATITUDE BY LONGITUDE MAPS OF DB11pdC9.f:c ** TITLES FOR SUBROUTINE DIAG9 DB11pdC9.f: SUBROUTINE DIAG9A (M) DB11pdC9.f: SUBROUTINE DIAG5A (M5,NDT) DB11pdC9.f:c ** TITLES FOR SUBROUTINE DIAG6 DB11pdC9.f: SUBROUTINE DIAG6 DB11pdC9.f:c ** THIS SUBROUTINE PRINTS THE DIURNAL CYCLE OF SOME QUANTITIES DB11pdC9.f: SUBROUTINE DIAG4A DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES A TIME HISTORY OF ENERGIES DB11pdC9.f: SUBROUTINE DIAGKS DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES A SUMMARY OF KEY NUMBERS CALCULATED IN DB11pdC9.f:c ** OTHER DIAGNOSTIC SUBROUTINES DB11pdC9.f: SUBROUTINE DIAG10(IPFLAG) DB11pdC9.f:c ** THIS SUBROUTINE SAVES THE INSTANTANEOUS SEA LEVEL PRESSURES DB11pdC9.f: SUBROUTINE SIJMAP DB11pdC9.f:c ** THIS SUBROUTINE PRODUCES LATITUDE BY LONGITUDE MAPS OF: FFT36macDBL.f: SUBROUTINE FRTR0 (KM) FFT36macDBL.f: 901 FORMAT ('0THIS FOURT SUBROUTINE NOT SUITED FOR KM = ',I8) FFT36macDBL.f: SUBROUTINE FRTR (F) FFT36macDBL.f:c ** THIS SUBROUTINE PERFORMS A FOURIER ANALYSIS ON THE ONE DIMENSIONAL FFT36macDBL.f: SUBROUTINE GETAN (F,G) FORCINGSjalC9.f: SUBROUTINE SOLARFORCING FORCINGSjalC9.f: SUBROUTINE DATAREFYEAR(XGAS,YEAR,NGAS) FORCINGSjalC9.f: SUBROUTINE DATAFILETREND(XGAS,YEAR,NGAS) Mjal2cpdC9.f: SUBROUTINE GEOM 401. Mjal2cpdC9.f: SUBROUTINE GEOM_8x10 Mjal2cpdC9.f: SUBROUTINE INPUT 501. Mjal2cpdC9.f:C**** THIS SUBROUTINE SETS THE PARAMETERS IN THE C ARRAY, READS IN THE 503. Mjal2cpdC9.f:C IN SUBROUTINE SURFACE THE MORE PRECISE VALUE 99.66 IS USED 850.3 Mjal2cpdC9.f: * ' RESUBMIT JOB WITH AN EARLIER TAUE CARD' 882.2 Mjal2cpdC9.f: SUBROUTINE DYNAM Mjal2cpdC9.f: SUBROUTINE AFLUX (U,V,P) 1001. Mjal2cpdC9.f:C**** THIS SUBROUTINE CALCULATES THE HORIZONTAL AIR MASS FLUXES 1003. Mjal2cpdC9.f: SUBROUTINE ADVECM (P,PA,DT1) 1501. Mjal2cpdC9.f:C**** THIS SUBROUTINE CALCULATES UPDATED COLUMN PRESSURES AS 1503. Mjal2cpdC9.f: SUBROUTINE ADVECV (PA,UT,VT,PB,U,V,P,DT1) 2001. Mjal2cpdC9.f:C**** THIS SUBROUTINE ADVECTS MOMENTUM (INCLUDING THE CORIOLIS FORCE) 2003. Mjal2cpdC9.f: SUBROUTINE PGF (UT,VT,PB,U,V,T,P,DT1) 2501. Mjal2cpdC9.f:C**** THIS SUBROUTINE ADDS TO MOMENTUM THE TENDENCIES DETERMINED BY 2503. Mjal2cpdC9.f: SUBROUTINE ADVECT (PA,TT,PB,T,DT1) 3001. Mjal2cpdC9.f:C**** THIS SUBROUTINE ADVECTS POTENTIAL TEMPERATURE 3003. Mjal2cpdC9.f: SUBROUTINE ADVECQ (PA,QT,PB,Q,DT1) 3501. Mjal2cpdC9.f:C**** THIS SUBROUTINE ADVECTS HUMIDITY AS DETERMINED BY DT1 AND THE 3503. Mjal2cpdC9.f: SUBROUTINE AVRX (PU) 1801. Mjal2cpdC9.f:C**** THIS SUBROUTINE SMOOTHES THE ZONAL MASS FLUX AND GEOPOTENTIAL 1803. Mjal2cpdC9.f: SUBROUTINE SDRAG 7001. Mjal2cpdC9.f:C**** THIS SUBROUTINE PUTS A DRAG ON THE WINDS ON THE TOP LAYER OF 7003. Mjal2cpdC9.f: SUBROUTINE FILTER 7501. Mjal2cpdC9.f:C**** THIS SUBROUTINE PERFORMS AN 8-TH ORDER SHAPIRO FILTER ON 7503. Mjal2cpdC9.f: SUBROUTINE SHAP1D (NORDER) 7801. Mjal2cpdC9.f:C**** THIS SUBROUTINE SMOOTHES THE ARRAY X IN THE ZONAL DIRECTION 7803. Mjal2cpdC9.f: SUBROUTINE DAILY 8001. Mjal2cpdC9.f:C**** THIS SUBROUTINE PERFORMS THOSE FUNCTIONS OF THE PROGRAM WHICH 8003. Mjal2cpdC9.f: SUBROUTINE CHECKT (N) 9001. Mjal2cpdC9.f:C**** THIS SUBROUTINE CHECKS WHETHER THE TEMPERATURES ARE REASONABLE 9003. Pjal0C9.f: SUBROUTINE CONDSE 2001. Pjal0C9.f:C**** THIS SUBROUTINE MIXES AIR CAUSED BY MOIST CONVECTION, AND 2003. Pjal0C9.f:C**** SET UP VERTICAL ARRAYS, OMITTING THE J AND I SUBSCRIPTS 2170. Pjal0C9.f:C**** SUBSIDENCE AND MIXING 2306. Pjal0C9.f: SUBROUTINE PRECIP 3001. Pjal0C9.f:C**** THIS SUBROUTINE USES THE PRECIPITATION TO CALCULATE THE GROUND 3003. Pjal0C9.f: SUBROUTINE RADIA0 3501. Pjal0C9.f:C**** THIS SUBROUTINE SETS THE RADIATION CONTROL PARAMETERS AND 3503. Pjal0C9.f: SUBROUTINE RADIA 4001. Pjal0C9.f:C**** THIS SUBROUTINES ADDS THE RADIATION HEATING TO THE TEMPERATURES 4003. Pjal0C9.f: SUBROUTINE SURFCE 4801. Pjal0C9.f:C**** THIS SUBROUTINE CALCULATES THE SURFACE FLUXES WHICH INCLUDE 4803. Pjal0C9.f: SUBROUTINE GROUND 6001. Pjal0C9.f:C**** THIS SUBROUTINE USES THE SURFACE FLUXES TO PREDICT IN TIME THE 6003. Pjal0C9.f:C**** TG1 IS AT FREEZING POINT, SUBTRACT WATER THAT DIFFUSES OUT OF G1 6481. Pjal0C9.f: SUBROUTINE DRYCNV 6801. Pjal0C9.f:C**** THIS SUBROUTINE MIXES AIR CAUSED BY DRY CONVECTION. SINCE DRY 6803. Pjal0C9.f:C**** CONVECTION IN THE BOUNDARY LAYER IS DONE IN SUBROUTINE SURFCE, 6804. Pjal0C9.f: SUBROUTINE ORBIT (OBLIQ,ECCN,OMEGT,DAY,SDIST,SIND,COSD,LAMBDA) 8201. Pjal0C9.f: SUBROUTINE OSTRUC 8501. Pjal0C9.f:C**** THIS SUBROUTINE RESTRUCTURES THE OCEAN TEMPERATURE PROFILE 8503. Pjal0C9.f:C**** THE SUBROUTINE ALSO MELTS ICE WHEN TGO > 0 (C). 8506. Pjal0C9.f: SUBROUTINE ODIFS 8601. Pjal0C9.f:C**** THIS SUBROUTINE CALCULATES THE ANNUAL OCEAN TEMPERATURE AT THE 8603. Pjal0C9.f:C**** TEMPERATURE, CALLS SUBROUTINE DIFFUS, AND REDUCES THE UPPER 8605. Pjal0C9.f: SUBROUTINE DIFFUS (IM,JM,DT,ALPHA,ED,Z12O,PWATER,R) 8701. Pjal0C9.f:C**** THIS SUBROUTINE CALCULATES THE VERTICAL MIXING OF A TRACER, R, 8703. R83ZAmacDBL.f: SUBROUTINE RCOMP1(NFTTTR,NFTTSR,NFTFOR) 1. R83ZAmacDBL.f: SUBROUTINE SETALB 179. R83ZAmacDBL.f: SUBROUTINE SETGAS 587. R83ZAmacDBL.f: SUBROUTINE SETAER 1011. R83ZAmacDBL.f: SUBROUTINE TAUGAS 1236. R83ZAmacDBL.f: SUBROUTINE THERML 1515. R83ZAmacDBL.f: SUBROUTINE SOLAR 1958. R83ZAmacDBL.f: SUBROUTINE SETAO2(O2CMA,NL) 2609. R83ZAmacDBL.f: SUBROUTINE SETO3D 2877. R83ZAmacDBL.f: SUBROUTINE PHDATM(P,H,D,T,O,Q,S,OCM,WCM,NPHD,NATM) 4238. R83ZAmacDBL.f:C NATM=4 GIVES ATMOSPHERE DATA FOR SUBARCTIC SUMMER 4252. R83ZAmacDBL.f:C NATM=5 GIVES ATMOSPHERE DATA FOR SUBARCTIC WINTER 4253. R83ZAmacDBL.f:C4444 SUBARCTIC SUMMER MCCLATCHY (1972) ATMOSPHERE DATA VS HEIGHT4422. R83ZAmacDBL.f:C5555 SUBARCTIC WINTER MCCLATCHY (1972) ATMOSPHERE DATA VS HEIGHT4453. R83ZAmacDBL.f: SUBROUTINE WRITER(INDEX,KPAGE) 4813. R83ZAmacDBL.f: SUBROUTINE SOLARZ(NG,KWRITE) 5668. R83ZAmacDBL.f: SUBROUTINE GAUSST(NG,X1,X2,XP,WT) 5782. R83ZAmacDBL.f: SUBROUTINE SETATM 5833. R83ZAmacDBL.f: SUBROUTINE SETFOR(NFTFOR) 6108. R83ZAmacDBL.f: SUBROUTINE HGAER1(XMU,TAU,G,GG) 6301. R83ZAmacDBL.f: SUBROUTINE HGCLD1(XMU,TAU,G,GG) 6820. R83ZAmacDBL.f: SUBROUTINE ADDVOL 7472. RFRCmacDBL.f: SUBROUTINE FORSET(TREF,KTREND,KWRITE) 1. RFRCmacDBL.f: SUBROUTINE ATREND(XGAS,YEAR,NGAS) 1001. RFRCmacDBL.f: SUBROUTINE BTREND(XGAS,YEAR,NGAS) 2001. RFRCmacDBL.f: SUBROUTINE CTREND(XGAS,YEAR,NGAS) 3001. RFRCmacDBL.f: SUBROUTINE DTREND(XGAS,YEAR,NGAS) 4001. RFRCmacDBL.f: SUBROUTINE YTREND(XGAS,YEAR,NGAS) 5001. RFRCmacDBL.f: SUBROUTINE ZTREND(XGAS,YEAR,NGAS) 6001. RFRCmacDBL.f: SUBROUTINE DTDX1D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7101. RFRCmacDBL.f: SUBROUTINE DXDT1D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7201. RFRCmacDBL.f: SUBROUTINE DXDT3D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7301. RFRCmacDBL.f: SUBROUTINE DTDX3D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7401. RFRCmacDBL.f: SUBROUTINE DTXCFY(DTSUM,DTNR,XNOW,XREF,NGAS) 7501. RFRCmacDBL.f: SUBROUTINE DTXCFZ(DTSUM,DTNR,XNOW,XREF,NGAS) 7601. RFRCmacDBL.f: SUBROUTINE VMSSET(KDTAUA,KWRITE) 8001. setpath.f: SUBROUTINE SETPATH UTILmacDBL.f: SUBROUTINE CLOCKS(IHSC) UTILmacDBL.f: SUBROUTINE DREAD (IUNIT,AIN,LENGTH,AOUT) UTILmacDBL.f: SUBROUTINE MREAD (IUNIT,M,NSKIP,AIN,LENGTH,AOUT) UTILmacDBL.f: SUBROUTINE READT (IUNIT,NSKIP,AIN,LENGTH,AOUT,IPOS) UTILmacDBL.f: SUBROUTINE TIMER (MNOW,MINC,MSUM) UTILmacDBL.f:c ** MINC (.01 S) = TIME SINCE LAST CALL TO SUBROUTINE TIMER UTILwinDBL.f: SUBROUTINE CLOCKS(IHSC) UTILwinDBL.f: SUBROUTINE DREAD (IUNIT,AIN,LENGTH,AOUT) UTILwinDBL.f: SUBROUTINE MREAD (IUNIT,M,NSKIP,AIN,LENGTH,AOUT) UTILwinDBL.f: SUBROUTINE READT (IUNIT,NSKIP,AIN,LENGTH,AOUT,IPOS) UTILwinDBL.f: SUBROUTINE TIMER (MNOW,MINC,MSUM) UTILwinDBL.f:c ** MINC (.01 S) = TIME SINCE LAST CALL TO SUBROUTINE TIMER

### Breaking It Down by Program

OK, now to reduce it to “one per definition” (or ‘weed out the comments and crap’)

DB11pdC9.f: SUBROUTINE DIAGA (U,V,T,P,Q) DB11pdC9.f: SUBROUTINE DIAGB (U,V,T,P,Q) DB11pdC9.f: SUBROUTINE DIAGJ DB11pdC9.f: SUBROUTINE DIAGJK DB11pdC9.f: SUBROUTINE JKMAP (NT,PM,AX,SCALE,SCALEJ,SCALEK,KMAX,JWT,J1) DB11pdC9.f: SUBROUTINE DIAGJL DB11pdC9.f: SUBROUTINE JLMAP (NT,PL,AX,SCALE,SCALEJ,SCALEL,LMAX,JWT,J1) DB11pdC9.f: SUBROUTINE DIAGIL DB11pdC9.f: SUBROUTINE ILMAP (NT,PL,AX,SCALE,SCALEL,LMAX,JWT,ISHIFT) DB11pdC9.f: SUBROUTINE DIAG7A DB11pdC9.f: SUBROUTINE MEM (SERIES,ITM,MMAX,NUAMAX,NUBMAX,POWER,FPE,VAR,PNU) DB11pdC9.f: SUBROUTINE DIAGIJ

I note in passing that we already have confusingly similar names in DIAGJK DIAGJL DIAGIL DIAGIJ … That’s just very human hostile coding practice…

DB11pdC9.f: SUBROUTINE IJMAP (NT,ARRAY,BYIACC) DB11pdC9.f: SUBROUTINE DIAG9A (M) DB11pdC9.f: SUBROUTINE DIAG5A (M5,NDT) DB11pdC9.f: SUBROUTINE DIAG6 DB11pdC9.f: SUBROUTINE DIAG4A DB11pdC9.f: SUBROUTINE DIAGKS DB11pdC9.f: SUBROUTINE DIAG10(IPFLAG) DB11pdC9.f: SUBROUTINE SIJMAP

So at this point, we have 20 subroutines devoted to “Diagnostics”. I suspect these are not very important at this stage of the analysis, but we’ll see when we get more into it. The diagnostic part is not high on my investigation rank of what to look at next…

FFT36macDBL.f: SUBROUTINE FRTR (F) FFT36macDBL.f: SUBROUTINE GETAN (F,G)

Two subroutines in the Fast Fourier Transform module. OK, might be worth a look to assure it’s done right and efficiently, but most likely it is and it is fine. Wonder why they did their own Get Tangent function, though…

FORCINGSjalC9.f: SUBROUTINE SOLARFORCING FORCINGSjalC9.f: SUBROUTINE DATAREFYEAR(XGAS,YEAR,NGAS) FORCINGSjalC9.f: SUBROUTINE DATAFILETREND(XGAS,YEAR,NGAS)

Gee, “FORCINGS” that have a (IMHO too simple) Solar component of seasonal and orbital shift of TSI / latitude and then two that look at GAS as “FORCINGS”… Worth looking into the SOLARFOCING routine to see what makes it tick, to confirm my suspicions. (Or refute them…)

This next block seems to do the bulk of the world modeling, with air flux and geometry and such.

Mjal2cpdC9.f: SUBROUTINE GEOM 401. Mjal2cpdC9.f: SUBROUTINE GEOM_8x10 Mjal2cpdC9.f: SUBROUTINE INPUT 501. Mjal2cpdC9.f: SUBROUTINE DYNAM Mjal2cpdC9.f: SUBROUTINE AFLUX (U,V,P) 1001. Mjal2cpdC9.f: SUBROUTINE ADVECM (P,PA,DT1) 1501. Mjal2cpdC9.f: SUBROUTINE ADVECV (PA,UT,VT,PB,U,V,P,DT1) 2001. Mjal2cpdC9.f: SUBROUTINE PGF (UT,VT,PB,U,V,T,P,DT1) 2501. Mjal2cpdC9.f: SUBROUTINE ADVECT (PA,TT,PB,T,DT1) 3001. Mjal2cpdC9.f: SUBROUTINE ADVECQ (PA,QT,PB,Q,DT1) 3501.

This next one has the comment that it covers the advection of humidity. I need to look inside it to see if it covers vertical movement as well (though advection ought not…). IMHO it is likely one of the ones with “something wrong” in how moisture is handled.

Mjal2cpdC9.f:C**** THIS SUBROUTINE ADVECTS HUMIDITY AS DETERMINED BY DT1 AND THE 3503.

Mjal2cpdC9.f: SUBROUTINE AVRX (PU) 1801. Mjal2cpdC9.f: SUBROUTINE SDRAG 7001.

As I recall it, the Hansen paper said they needed to add this drag to get things to stabilize. An interesting “Dig Here!” would be the justification for the particular value used. Near as I could tell from the paper, it was a tuning point…

Mjal2cpdC9.f:C**** THIS SUBROUTINE PUTS A DRAG ON THE WINDS ON THE TOP LAYER OF 7003.

Continuing…

Mjal2cpdC9.f: SUBROUTINE FILTER 7501. Mjal2cpdC9.f: SUBROUTINE SHAP1D (NORDER) 7801. Mjal2cpdC9.f: SUBROUTINE DAILY 8001. Mjal2cpdC9.f: SUBROUTINE CHECKT (N) 9001. Mjal2cpdC9.f:C**** THIS SUBROUTINE CHECKS WHETHER THE TEMPERATURES ARE REASONABLE 9003.

Another interesting question: What is done when temperatures have gotten “unreasonable”? Is this just a lid crammed on instabilities in the model? What does it do, and why? With what justification other than bringing stability to an out of control runaway process?

OK, 16 subroutines in Mjal* to look at. Not too bad for the global driver part.

This next one (Pjal*) looks to have the physics in it.

Pjal0C9.f: SUBROUTINE CONDSE 2001. Pjal0C9.f:C**** THIS SUBROUTINE MIXES AIR CAUSED BY MOIST CONVECTION, AND 2003.

Another place where one ought to look for how humidity ends up at altitude in the model, but not in reality…

Pjal0C9.f: SUBROUTINE PRECIP 3001. Pjal0C9.f: SUBROUTINE RADIA0 3501. Pjal0C9.f: SUBROUTINE RADIA 4001. Pjal0C9.f: SUBROUTINE SURFCE 4801. Pjal0C9.f: SUBROUTINE GROUND 6001. Pjal0C9.f: SUBROUTINE DRYCNV 6801. Pjal0C9.f:C**** THIS SUBROUTINE MIXES AIR CAUSED BY DRY CONVECTION. SINCE DRY 6803. Pjal0C9.f:C**** CONVECTION IN THE BOUNDARY LAYER IS DONE IN SUBROUTINE SURFCE, 6804.

More clue on where to look for the wrong moisture at altitude part. Dry Convection likely needs to be higher, and mixing less? Or maybe just the precipitation routine needs to wring out more water…

Pjal0C9.f: SUBROUTINE ORBIT (OBLIQ,ECCN,OMEGT,DAY,SDIST,SIND,COSD,LAMBDA) 8201. Pjal0C9.f: SUBROUTINE OSTRUC 8501. Pjal0C9.f: SUBROUTINE ODIFS 8601. Pjal0C9.f: SUBROUTINE DIFFUS (IM,JM,DT,ALPHA,ED,Z12O,PWATER,R) 8701.

I make that 11 physics routines to work through. IMHO, it is likely where things need the most help. Then again, the next block does a lot of the actual stirring of the pot.

R83ZAmacDBL.f: SUBROUTINE RCOMP1(NFTTTR,NFTTSR,NFTFOR) 1. R83ZAmacDBL.f: SUBROUTINE SETALB 179. R83ZAmacDBL.f: SUBROUTINE SETGAS 587. R83ZAmacDBL.f: SUBROUTINE SETAER 1011. R83ZAmacDBL.f: SUBROUTINE TAUGAS 1236. R83ZAmacDBL.f: SUBROUTINE THERML 1515. R83ZAmacDBL.f: SUBROUTINE SOLAR 1958. R83ZAmacDBL.f: SUBROUTINE SETAO2(O2CMA,NL) 2609. R83ZAmacDBL.f: SUBROUTINE SETO3D 2877. R83ZAmacDBL.f: SUBROUTINE PHDATM(P,H,D,T,O,Q,S,OCM,WCM,NPHD,NATM) 4238. R83ZAmacDBL.f: SUBROUTINE WRITER(INDEX,KPAGE) 4813. R83ZAmacDBL.f: SUBROUTINE SOLARZ(NG,KWRITE) 5668. R83ZAmacDBL.f: SUBROUTINE GAUSST(NG,X1,X2,XP,WT) 5782. R83ZAmacDBL.f: SUBROUTINE SETATM 5833. R83ZAmacDBL.f: SUBROUTINE SETFOR(NFTFOR) 6108. R83ZAmacDBL.f: SUBROUTINE HGAER1(XMU,TAU,G,GG) 6301. R83ZAmacDBL.f: SUBROUTINE HGCLD1(XMU,TAU,G,GG) 6820. R83ZAmacDBL.f: SUBROUTINE ADDVOL 7472.

Looks like 18 that do the setup of various batches of input and such. Ought to look through it to see if this is where the input files get read in too.

RFRCmacDBL.f: SUBROUTINE FORSET(TREF,KTREND,KWRITE) 1. RFRCmacDBL.f: SUBROUTINE ATREND(XGAS,YEAR,NGAS) 1001. RFRCmacDBL.f: SUBROUTINE BTREND(XGAS,YEAR,NGAS) 2001. RFRCmacDBL.f: SUBROUTINE CTREND(XGAS,YEAR,NGAS) 3001. RFRCmacDBL.f: SUBROUTINE DTREND(XGAS,YEAR,NGAS) 4001. RFRCmacDBL.f: SUBROUTINE YTREND(XGAS,YEAR,NGAS) 5001. RFRCmacDBL.f: SUBROUTINE ZTREND(XGAS,YEAR,NGAS) 6001. RFRCmacDBL.f: SUBROUTINE DTDX1D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7101. RFRCmacDBL.f: SUBROUTINE DXDT1D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7201. RFRCmacDBL.f: SUBROUTINE DXDT3D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7301. RFRCmacDBL.f: SUBROUTINE DTDX3D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 7401. RFRCmacDBL.f: SUBROUTINE DTXCFY(DTSUM,DTNR,XNOW,XREF,NGAS) 7501. RFRCmacDBL.f: SUBROUTINE DTXCFZ(DTSUM,DTNR,XNOW,XREF,NGAS) 7601. RFRCmacDBL.f: SUBROUTINE VMSSET(KDTAUA,KWRITE) 8001.

Looks to me like 14 subroutines devoted to Radiative Forcing calculations. (RFR?) Gee, wonder where the ocean cycle forcing and the cosmic ray forcing and the solar UV vs IR forcing and the… where are they handled?…

setpath.f: SUBROUTINE SETPATH UTILmacDBL.f: SUBROUTINE CLOCKS(IHSC) UTILmacDBL.f: SUBROUTINE DREAD (IUNIT,AIN,LENGTH,AOUT) UTILmacDBL.f: SUBROUTINE MREAD (IUNIT,M,NSKIP,AIN,LENGTH,AOUT) UTILmacDBL.f: SUBROUTINE READT (IUNIT,NSKIP,AIN,LENGTH,AOUT,IPOS) UTILmacDBL.f: SUBROUTINE TIMER (MNOW,MINC,MSUM) UTILwinDBL.f: SUBROUTINE CLOCKS(IHSC) UTILwinDBL.f: SUBROUTINE DREAD (IUNIT,AIN,LENGTH,AOUT) UTILwinDBL.f: SUBROUTINE MREAD (IUNIT,M,NSKIP,AIN,LENGTH,AOUT) UTILwinDBL.f: SUBROUTINE READT (IUNIT,NSKIP,AIN,LENGTH,AOUT,IPOS) UTILwinDBL.f: SUBROUTINE TIMER (MNOW,MINC,MSUM)

Then 11 that do utility functions. Not likely to be ‘an issue’, so lower priority to look over.

So what’s that total make?

Diagnostic: 20 FFT: 2 FORCINGS 3 Model Outer 16 Physics 11 Setting Up 18 Radiative 14 Utilities 11

I make that a total of 95 subroutines. OK, that’s going to take some time.

I’m going to start with disposing of the 2 FFT and 3 FORCINGS, just because they are small.

Then look at the 11 physics as I think they are likely the meat of it.

I suspect the Diagnostics are for postmortem of runs and will leave them for last (or a cursory examination may make it ‘never’ if they don’t contribute anything to the actual model runs.)

That will be followed by the 16 “Model Outer” as I seems to do a lot of interesting stuff in terms of how it coordinates the physics and such. (And has a lot of modification comments so a high ‘dick with’ factor implying it may be where tuning happened…)

IFF, after that, there’s still too much mystery, I’ll move on into the Setting Up 18 (or if I need to look at it to get the thing to run…) and the 14 Radiative (though I think they have the radiative physics right, just messed up the actual physical physics and left out important processes…)

I doubt “Utilities” has much interest, but I’ll poke my nose into it long enough to assure that.

### In Conclusion

So that’s the “big lumps” on the block diagram of Subroutines.

I still need to search through the code for functions and such. Then make sure I understand the outer driving cycle of the MAIN body (when I find it ;-)

Still, as a first cruise through, I think this is likely coming together nicely. We have some math and science parts (FFT and Radiative and Physics) that can be vetted (both for correctness and completion). We have an apparent Setting Up part and some Utilities and Diagnostics that likly have more to do with operations than execution of science, Then we have a FORCINGS that provides an understanding of what the modeler decided are the only relevant FORCINGS. Finally, the model MAIN driving the process of connecting all the bits and the outer loops. (Time and Geography…)

Doesn’t look to big, or too hard, just a bit long. Even at “one subroutine a day”, we’re talking 95 days. With 20 work days a month, that’s 5 months. And that’s just for the subroutines…

But every journey begins with just one step, and this is about step 3 on MODEL II, so we’re well on our way.

Anyone wishing to download, unpack, and take a whack at one of these, feel free (Hint! Hint!) and post anything you find in a comment. The download is at: http://edgcm.columbia.edu/ModelII/modelII_source.zip as a simple zip file, so just about any OS can download and unpack it… (HINT! HINT!! HINT!!!)

With that, I’m off to cook dinner (We’re having pepperoni and kielbasa on Naan bread pizza … yum! With olives and lots of Mexican Mix cheese on it…) then I’ll be back to start looking at FFT and FORCINGS.

Hmm, rolled their own FFT. Somewhere in the recesses of my old head I recall some nice little collection of numerical methods. If I am not mistaken, it has been around for at least 30 years in one form or another, since Knuth. My own copy is from the 90s or 00s and is packed up from the last move and is still in the box at the local Adding Space.

PS, my local EdGCM archive is unzipped and is some 79MB in size. I don’t recall ever running it and it is obviously well past its free trial use-by date. There were some small files in folders that started with input.

The FORTRAN source is all of 2 MB (result of “du -ms .” in that directory).

Perhaps once I figure out the inputs you could see if you have those files to get it started?

This is what I’m missing ( 11 files) from:

http://edgcm.columbia.edu/ModelII/InputFiles.html

Ground Data

Unit: 7

Name: Ground Data

Ex: G8X10

Size: 48392 bytes

Description: This file is used for istart=5 runs, paleo runs where the ground data does not match the information in the atmospheric conditions file because the ground has moved.

[…]

Ocean

Unit: 15

Name: Ocean

Ex: O8X10, O8X10.100MLD

Size: 124,560 bytes (3 array), 166,032 bytes (4 array), 16,579,248 bytes (observed)

Description: There are three versions of this file. The smaller ocean file has just the surface temperatures for the ocean for every month. When the smaller file is used the ocean temperature is fixed. The larger 4 array ocean file adds in a mixed layer depth so that the water in the mixed layer can change temperature. For the mixed layer to work some additional files must be read in by the model. The third version of the file is observed mode. The ocean file contains a number of years of ocean data starting from February of the first year. TAUI is matched to the start of the first year of the ocean file so that the model can advanced the file to the correct month. In all cases the ocean files contain an averaged value for the middle of each month and days are linearly interpolated between the two nearest monthly values

Arrays: Two version of the ocean include 12 months of data with the following format for the 3 array ocean:

MON, O(IM,JM,3)

The 4 array ocean adds a mixed layer depth:

MON, O(IM,JM,3), OMLD(IM,JM)

Arrays: Observed oceans contain these arrays repeated from 02/01/1870-12/01/1998

The contents of the (IM,JM,3) array is

1. OCEAN SURFACE TEMPERATURE (DEGREES C)’

2. RATIO OF OCEAN ICE TO WATER (PERCENT)’,

3. OCEAN ICE AMOUNT OF LAYER 2 (10**2 KG/M**2)’,

[…]

Ocean Transports

Unit: 17

Name: Ocean transports spectral coefficient (OTSpec)

Ex: OTSPEC.RP041bC9.M250ZD

Size: 62,216 bytes

Description: fourth order spectral coefficients for ocean heat transports in the qflux and deep ocean. Gary developed a complex scheme where the annual temperature variations of an ocean grid cell’s mixed layer are encoded into spectral coefficients. The upshot is that the otspec file adjusts the temperature of the mixed layer to incorporate season and the energy that would be added or removed by ocean currents (if they existed).

Arrays:

OTA1(IM,JM),OTB1(IM,JM),OTC(IM,JM)

OTA2(IM,JM),OTB2(IM,JM),OTA3(IM,JM)

OTB3(IM,JM),OTA4(IM,JM),OTB4(IM,JM)

Drag Coefficient

Unit: 19

Name: Drag coefficient

Ex: CD8X10

Size: 3464 bytes

Description: Drag coefficient

Arrays: ROUGHL(IM,JM)

Radiation

Unit: 21

Name: Radiation RTAU

Ex: RTAU.G25L15

Size: 316 kb

Description: Radiation RTAU

Arrays: ?

Radiation

Unit: 22

Name: Radiation RPLK

Ex: RPLK25

Size: 28 kb

Description: Radiation RPLK

Arrays: ?

Vegetation

Unit 23

Name: Vegetation

Ex: V8X10

Size: 28040 bytes

Description: There are eight types of vegetation defined by the model. Eight arrays store the distributions of vegetation as a percentage of the land in each grid box. There is also a veg-table (vadata) that defines the properties of each of the kinds of vegetation. The four (!) sections of the veg-table are: bare earth visual albedos (seasonal), bare earth near ir albedos (seasonal), water field capacities in kg/m**2 (3 arrays), masking depths in m (1 array). Since the eight simple types are mappings from thirty complex types the distributions are difficult to change for the modern. See NASA tech memo 86096, June 1984, Elaine Matthews.

Arrays: V(IM,JM,8), VADATA(8,4,3)

The order of the vegetation in the file is:

1’RATIO OF DESERT TO EARTH (PERCENT)’,

2’RATIO OF TUNDRA TO EARTH (PERCENT)’,

3’RATIO OF GRASSLAND TO EARTH (PERCENT)’,

4’RATIO OF SHRUB/GRASSLAND TO EARTH (PERCENT)’,

5’RATIO OF TREE/GRASSLAND TO EARTH (PERCENT)’,

6’RATIO OF DECIDUOUS FOREST TO EARTH (PERCENT)’,

7’RATIO OF EVERGREEN FOREST TO EARTH (PERCENT)’,

8’RATIO OF RAINFOREST TO EARTH (PERCENT)’/

QFlux Max Mixed Layer Depth

Unit: 25

Name: QFLUX Max Mixed Layer Depth

Ex: Z1OMAX.8X10.250M

Size: 3464 bytes

Description: The mixed layer is the surface part of the ocean handled in qflux runs. The max mixed layer depth is simply the maximum thickness of the mixed layer during the year for every grid point. This array is trivial to calculate and later versions of the GCM do this internally instead of using this file. For not very good reasons this file usually starts with a Z however it is not connected to topography in any way.

Arrays: (IM,JM)

Topography

Unit: 26

Name: Topography

Ex: Z8X101

Size: 10376 bytes

Description: The topography file defines the three inorganic properties of land in Model II. It is important to realize that this file is not the same as in Model II’ or Model E. The three properties are: topographic height in meters scaled by gravity, land cover, and land ice. The topographic height is the area weighted average height of the land in a grid box. The wacky part is that topographic height is scaled by gravity. Land cover is just the percent of a grid box covered by land. Land ice is the percentage of land in a grid box covered by ice. Thus is a grid box is 12% land by all that land is ice covered then the land ice is 100%. Land ice is special unmelting ice that doesn’t interact with the ice created by repeated snowfall.

Arrays: PHIS(IM,JM),PLAND(IM,JM),RLICE(IM,JM)

23 Special Regions

Unit: 29

Name: 23 special regions

Ex: REG8X10

Size: 3728 bytes

Description: The 23 special regions in model have extra information collected for them. Unlike all the other input files there is a hidden special regions file in the model that is used if no file is read in. The file contains an 80 byte title and then a single array of the 23 special regions. The regions are numbered 1-23 with 24 as no special region. Since there is only a single array the regions cannot overlap, but they do not have to be contiguous. The last array contains the names of the 23 special regions as 8 character strings.

Arrays: TITLE, JREG(IM,JM), IXJX(23)

Deep Ocean

Unit: 62

Name: ED Deep Ocean

Ex: ED8X10

Size: 3464 bytes

Description: Eddie diffusions to the deep ocean from modern observations. How fast the bottom of the mixed layer is diffused into the deep ocean. This file is only used when kocean=2.

Arrays: EDO(IM,JM)

I am saddened that they still rely on ocean data being an input rather than trying to calculate currents and other mixing from first principals. (I suppose a lot of that is that their grid size is far too large to capture those processes, but still).

@Soronel:

Remember that this Model II code is from about 1980… (some of it even earlier…)

It may well be that other newer models have gone well beyond it.

With that said, there is a lot of “fitting” of the model to “observations” done Those input files do much of it. IMHO, that’s one of the biggest problems with this model. You fit it just too much. That can mask all sorts of errors or omisions.

So, say, the atmosphere is unstable and results in runaways. They adjust the input for atmospheric ‘drag’ until it makes more “normal” results. That can be valid (IF drag really is the problem and IF they get it right and IF they are not missing something else). It can also be very invalid. Say, for example, there were a mysterious “ballooning factor” that caused atmosphere thickness to change, but you didn’t have it in the model. You fiddle with the ‘drag’ parameters until it reflects reality nicely. Then run the model for your “conclusions”… BUT, you left out that other thickness factor during the tuning of the drag setting, so changes in thickness can now result in big changes you can’t get in your model. Worse, if you are correctly showing temperature changes based on drag of one value, but it really ought to be another when thickness is right, well, your temperature results are now wrong… and your sensitivity to other factors now wrong too.

Note that this example is very roughly modeled on something they actually did. Drag was wrong and they had to adjust it to prevent various failure modes… From Hansen 1983:

http://pubs.giss.nasa.gov/docs/1983/1983_Hansen_ha05900x.pdf

page 7 of the pdf, but page 615 as published

The paper has more of that kind of thing. ~”This value was yielding wrong results, so we changed it to get what looked better”. That’s a lot of tuning to justify. Then it clearly isn’t quite right, yet that mal-fit is ignored. (“unrealistically low altitude”)…

That paper is just a scan (no OCR) so hard to copy from. I’m looking for one a bit more amenable to copy / paste to do a review on it. Still, anyone wanting to ‘read ahead’, it’s the one that describes this Model II code, by the head author…

E.M.,

Based on past experience in writing stuff for the ‘metrology of circular objects’ (don’t ask) I’d say the get tangent routine was probably because tangents have that pesky habit of heading towards infinity at the most inopportune times?

I’d guess it’s just a wrapper to catch that and prevent crashes?

Most trigonometric (arc) routines back then would also only give results in 2 quadrants out of the 4, and needed deep and arcane magic and tweaking to extend the results to the full 360 degrees

Down in the subroutine declarations in the main body program, we have this interesting note:

Originally, the main loop called a uniform DYNAM subroutine, but in this version a comment says it was divided into subroutines of its own. This is one of them, handling a special case of advection (sideway movement).

So they need to give the poles special handling or it has “computational instability”…

There’s another interesting bit down in the ocean / ice subroutine. First we read in some data, but then if TGO is not defined, we force it to a ‘reasonable value’:

So what makes -10 “reasonable” in the absence of any data? …

Then it reduced the ocean ice / water ratio. Why? Don’t ask why… there was just too much ice, that’s all, and it had to go. /sarc;

Having reduced the ocean ice, if it qualifies as ‘zero’ we get rid of snow and stuff too.

In the absence of any physical explanation for ‘why’ it sure looks like a thumb on the ocean ice scale… I hope “why” is covered in one of the papers somewhere…

Interesting. Subroutine CHECKT seems to force having below zero (C) ground or water temperatures for ice or snow to exist on top of them. Now I can see that long term, but when snow first falls, it is often sitting on top of warm ground. Very slowly melting into it and increasing soil moisture. Similarly, ocean ice takes time to melt and depending on the size, can be a log way south of the arctic…

So you can’t model the first snow of spring? Or a small iceberg melting off of the St. Lawrence Seaway?

I wonder if they just force that snow to not exist to make an acceptable model run… but that’s part of the “fixing up the data” that happens in the input files that are missing…

Pingback: GCM Model II – DYNAM and the Dynamic Subroutines | Musings from the Chiefio

Hmm, I think that I do have some of the needed input files. I don’t know if I have them all. They’re small. Zipping the folders named Data and Input, I get five folders that are about 9MB total. [I just love the cryptic names ;P.]

EM, that comment about drag hits the nail on the head. Dr. Gerald Browning discusses this problem extensively in numerical weather prediction models, which do have the advantage of measurements to constrain the inherent numerical instabilities due to finite precision mathematics and the oversimplifications done to get these things to run in reasonable time frames. Dr. Browning also pans these errors in the worse-off GCMs that can’t be checked against real world data (with the warts that come with it).