Pattern Matching in GISS/NASA ModelE Coding
In a previous post I gave an illustration of how GISS/NASA employees have implemented new and innovative ways to produce inactive code using the capabilities provided by F90/95. I had run across the following statements in routine DIAG.f:
EWATER(J)=EWATER(J)+EL !+W*(SHV*T(I,J,L)*PK(L,I,J)+GRAV
! * *HSCALE*LOG(P(I,J)/PMID(L,I,J)))
The ‘!’ in the first line is going to be very difficult to remember it exists and correctly maintain. Someone might come along and say, “I wonder what that’s doing in the middle of an executable statement.” and promptly un-do the comment. Or un-do the comment of the second line while overlooking the comment in the first line. That would make a screw up on several levels.
Today I have found many more examples of innovative coding by employees of GISS/NASA. It is clear that the NASA Software Quality Assurance procedures are ignored by GISS/NASA. It is equally clear that there are no Software Quality Assurance procedures being applied to the GISS/NASA ModelE code. None.
Update November 2, 2008 down near the end.
Black box computer codes are universally acknowledged to have great potential for harboring inordinate numbers of bugs. But I had never seen anything like the above example of bad coding in my entire career. And believe me I have seen lots of bad coding. In my industry, legacy coding rules.
Recently I was thinking about that coding and said to myself, Hmmm … I wonder how many more applications of that particularly nasty application of code are in ModelE. I was thinking that I might have to try many arbitrary search patterns to run across even another single example. Imagine my surprise when the very first one I tried, ‘! +’ gave me 18 hits. Now some of these are in fact inline comments and do not affect the calculations. And in this case, inline means actually on the line of coding that performs calculations. Such comments are generally very much less than useful.
Here is an example of what that first EWAG search found, in routine SEAICE.f
C**** reconstitute upper layers
MSI1 = ACE1I + SNOW
IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
HSIL(1) = (HICE-FHSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I + HSNOW
SSIL(1) = (SICE-FSSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I ! + SSNOW
Those inline comments that modify the calculations are hard to spot, I think.
Well, given that this was much easier than I ever expected, I gave it another shot with ‘ !*’ and got 39 hits. Again some are actual comments, but the majority are not comments but instead modify the calculations. This search gave hits in the several versions of the CLOUDS_S12000 routines
C**** COMPUTE EVAPORATION OF RAIN WATER, ER
RHN=MIN(RH(L),RHF(L))
IF(WMX(L).GT.0.) THEN
ER(L)=(1.-RHN)*LHX*PREBAR(L+1)*GRAV*BYAM(L) !*GbyAIRM0 ! GRAV/AIRM0
and here’s an example from GHY_DRV.f
#ifdef TRACERS_WATER
ccc accumulate tracer evaporation and runoff
do nx=1,ntg
n=ntixw(nx)
trevapor(n,itype,i,j) = trevapor(n,itype,i,j) + atr_evap(nx)
!trevapor(n,itype,i,j) = trevapor(n,itype,i,j) + aevap !*rhow
Actual blocks of comments that explain the modifications are not present in the code, but it seems to me that this change had a great potential to change the units for the aevap contribution. Apparently the coders were not so sure about this one either, because the entire line is also commented out. Makes you wonder the order of the modifications.
This was way too easy. So next I tried ‘! *’ and I got 20 hits. Here’s one from CLOUDS_S12000,f
ELSE
RCLD=25.0*(WTEM/4.2d-3)**BY3 ! * 500./(pl(l)+500.)
and from GHY.f
c atr_g(:ntg) = ! instantaneous ! atr_g(:ntg) +
c & ( fb*( tr_w(:ntg,1,1)*(1.d0-fr_snow(1))
c & + tr_wsn(:ntg,1,1) ) !*fr_snow(1)
c & + fv*( tr_w(:ntg,0,2)*(1.d0-fm*fr_snow(2)) ! *fw0 !! remove fw?
c & + tr_wsn(:ntg,1,2)*fm ) ) / !*fr_snow(2)
c & (rhow * tot_w1) ! * dts
I love that one. Note that in the case of the line containing !*fw0 !!remove fw? there was a display of We don’t know what we’re doing. And in the last line the multiplication by the time step, dts, has been commented out. To be fair, note that all these lines have been commented out by the use of ‘c’ in the first column. Given the clear lack of understanding about the coding, that is very likely a good thing.
There is also an application in SUBROUTINE SPLN44
XEXM=XE**2
!=1 CUSPWT=(1.D0-XEXM)*CUSPWM+XEXM*CUSPWE
!nu FFLINR=A+XF*F32
QK(K)=FFCUSP ! *CUSPWT+FFLINR*(1.D0-CUSPWT)
Next I tried ‘!+’ and got three hits that modify the calculations. One being the EWATER(j) case given at the top of this post. And this one also from DIAG.f
710 AJK(J,K-1,JK_TOTVTAM)=AJK(J,K-1,JK_TOTVTAM)+WU4I !+W4I*UEARTH
It looks like an error was discovered and fixed by this application. Again it seems that very likely the units could have been incorrect in the original coding. Or they may be incorrect in the modified coding.
And finally this one from GHY.f
! evaporation from the canopy
if ( process_vege ) then
ibv = 2
evap_max_wet(ibv) = w(0,2)/dt !+ pr ! pr doesn''t work for snow
So next I tried ‘! -‘ and got several hits that are difficult to figure out. For example, in ATMDYN_S12000.f we see:
do jk=iplimb+1,iplimt ! -1 ?????
zpmk(jk)=0.5*(pk(jk-1)+pk(jk))
Is that a correction to previous coding, or is it a question on the part of the person writing the code. Maybe they are hoping that the compiler will correctly figure it out.
And from GHY_DRV.f
c****
c**** outside loop over time steps, executed nisurf times every hour
c****
timez=jday+(mod(itime,nday)+(ns-1.)/nisurf)/nday ! -1 ??
if(jday.le.31) timez=timez+365.
yet additional sterling examples of certainty on the part of the code writers.
Well I could go on, but this exercise is getting to be just way too distressing to continue. Maybe you’ll fire up your favorite editor and look around in ModelE, too. You can get your personal copy from here.
November 2, 2008 update. I did two quick searches using ‘!/’ and ‘! /’ and got a few hits that modify the calculating. I’ll not give the details here.
There is probably a way to write a script that will find only the modifications that affect the calculations. Maybe I’ll give it a shot. My favorite scripting language for parsing strings is FORTRAN :-).
Additionally, I did several quick searches in the UCAR CAM 3.1.p2 source and did not get any hits. While looking through the coding, this code seems to be in far better shape than the ModelE coding.
I think I have been reasonable in all the posts on this blog. But I think these examples clearly illustrate an alarming lack of the most basic proper attention to the source of the numbers that are used in Climate Crisis analyses by employees of GISS/NASA. Actually I think this is bordering on being a national disgrace.
November 2, 2008 update. On second thoughts, it is a national disgrace. And, GISS/NASA doesn’t need “an extra thousand code checkers”, they need individual code writers that can follow the most simple concepts of writing reasonable code.
Models Methods Software 
Welcome to the U.S. Government’s way of doing things!
It is frightening that decisions are being made on this crap that will involve the expenditure of trillions of dollars.
A few questions in the GISS/NASA ModelE coding. There were many others. Almost all of the below have the potential to affect the calculations. It got to be too tedious trying to understand the coding.
conv_strat_rsf.f:32: REAL*8 TAUX,X ,rc(161) !(169)??
ocnIC.E001.f:73: ccc the above line could substitute for next 3 lines: same results ?
otspec.E001.f:144: ccc the above line could substitute for next 3 lines: same results ?
ATMDYN_SI2000.f:937: AKAP=KAPA-.205d0 ! what is this number?
ATMDYN_SI2000.f:1431: C??? currently only below 150mb – delete this line when this changes
ATMDYN_SI2000.f:1437: do l = 1,ls1-1 ! soon: 1,lm ???
ATMDYN_SI2000.f:1443: do l = 1,ls1-1 ! soon: 1,lm ???
ATMDYN_SI2000.f:1607: do jk=iplimb+1,iplimt ! -1 ?????
PBL_SI2000.f:299: call tfix(t,z,ttop,tgrnd,lmonin,n) !why should we do this?
RAD_DRV_rnd.f:897: AGESN(1)=SNOAGE(3,I,J) ! land ! ? why are these numbers
RAD_DRV_rnd.f:898: AGESN(2)=SNOAGE(1,I,J) ! ocean ice so confusing ?
ATMDYN.f:29: INTEGER NS, NSOLD,MODDA !? ,NIdynO
ATMDYN.f:43: !? NIdynO=MOD(NIdyn,2) ! NIdyn odd is currently not an option
ATMDYN.f:1162: AKAP=KAPA-.205d0 ! what is this number?
ATMDYN.f:2038: do jk=iplimb+1,iplimt ! -1 ?????
ATURB.f:145: cc tmomij(:,l)=tmom(:,i,j,l) ! vert. grad. should virtual ?
ATURB.f:292: C**** Does this ever happen for q? (put this in just in case)
DIAG.f:1353: C**** Warning: This diagnostic has 3 flaws (?)
DIAG.f:1354: 1 – It assumes that DTsrc=1hr, (DTsrc=3600.)
DIAG.f:1355: 2 – since DTdaa-Ndaa*DTsrc=2*DTdyn rather than 0,
DIAG.f:1356: some hours are skipped once in a while
DIAG.f:1357: * 3 – Some of the first Ndaa hours are skipped at the
DIAG.f:1358: beginning of a month and overcounted at the end;
DIAG.f:1359: this happens to balance out, if and only if
DIAG.f:1360: mod(days_in_month,ndaa)=0 (i.e. February if Ndaa=7)
DIAG.f:1696: PLO(L)=SP*SIG(L)+PTOP ! PL or PLO ??
DIAG.f:1697: 2025 PL(L)=SP*SIGE(L)+PTOP ! PLE or PL ??
DIAG.f:1699: IF (PM(K+1).GE.PLO(1)) GO TO 2090 ! really ?? not PL?
DIAG.f:1702: IF (PM(K).GE.PLO(1)) GO TO 2040 ! really ?? not PL?
DIAG_COM.f:46: !@var ASJL latitude/height supplementary diagnostics (merge with AJL?)
DIAG_COM.f:766: GO TO 10 ! or should that be just a warning ??
DIAG_COM.f:842: !**** beg=F?L where F/L=letter 1 of First/Last month for 2-11 mo.periods
DIAG_PRT.f:1088: COSBYPV(J)=0. ! are these values used?
DIAG_PRT.f:2297: MLAT(J)=NINT(MIN(1d5,MAX(-1d5,FLAT(J)))) ! prevent too large int?
DIAG_PRT.f:2299: 130 AGLOB=AGLOB+AHEML(JHEMI) !/JWT no longer needed?
DIAG_PRT.f:2489: if(sname.eq.’dudt_mtndrg’) then ! make sumfac an argument … ???
DIAG_PRT.f:2490: SUMFAC=10. ! … to avoid this if-block ???
DIAG_PRT.f:3421: off = TF ! should do in accum-phase ??
DIAG_PRT.f:3632: c**** length of growing season (not quite right ???)
GHY.f:157: !@var zw water table (not computed ?)
GHY.f:310: !@var trpr flux of tracers in precipitation (?/m^2 s)
GHY.f:311: !@var tr_w amount of tracers in the soil (?)
GHY.f:312: !@var tr_wsn amount of tracers in the snow (?)
GHY.f:316: !@var tr_evap flux of tracers to atm due to evaporation (?/m^2 s)
GHY.f:317: !@var tr_evap flux of tracers due to runoff (?/m^2 s)
GHY.f:853: ! need this ? if(igcm.ge.0 .and. igcm.le.3) xl=eddy/(z1-zs)
GHY.f:1601: ccc do we need this check ?
GHY.f:1854: aedifs=aedifs-dts*shw*dedifs*fv ! not used ?
GHY.f:1855: af0dt=af0dt-dts*(fb*fh(1,1)+fv*fch(0)+htpr) ! E0 excludes htpr?
GHY.f:1866: c & + fv*( w(0,2)*(1.d0-fm*fr_snow(2)) ! *fw0 !! remove fw ?
GHY.f:1872: c & + fv*( tr_w(:ntg,0,2)*(1.d0-fm*fr_snow(2)) ! *fw0 !! remove fw?
GHY.f:1927: abetad=0.d0 ! not accumulated : do we need it? YES
GHY.f:1928: abetav=0.d0 ! not accumulated : do we need it?
GHY.f:1929: abetat=0.d0 ! not accumulated : do we need it?
GHY.f:1929: abetap=0.d0 ! not sure how it is computed : probably wrong
GHY.f:1931: abetab=0.d0 ! not accumulated : do we need it?
GHY.f:1932: abeta=0.d0 ! not accumulated : do we need it?
GHY.f:1933: acna=0.d0 ! not accumulated : do we need it?
GHY.f:1934: acnc=0.d0 ! not accumulated : do we need it?
GHY.f:2445: !@var tr_wc ratio tracers/water in soil layers (?/m)
GHY.f:2446: !@var tr_wcc ratio tracers/water in canopy (?/m)
GHY.f:2449: !@var tr_wsnc ratio tracers/water in snow layers (?/m)
GHY.f:2473: ! tr_surf(:) = 1000.d0 !!!???
GHY.f:2600: wsni(i-1,ibv) = wsni(i-1,ibv) + flux_dt ! extra?
GHY.f:2617: wsni(i+1,ibv) = wsni(i+1,ibv) – flux_dt ! extra?
GHY_DRV.f:212: timez=jday+(mod(itime,nday)+(ns-1.)/nisurf)/nday ! -1 ??
GHY_DRV.f:1101: ccc ??? remove next 5 lines? -check the old version
GHY_DRV.f:1239: ! should also restrict to TYPE=nWATER ?
ICEDYN.f:59: !@var UIB,VIB velocity arrays (m/s) (????)
ICEDYN_DRV.f:734: c USE ICEGEOM, only : dxyp,dyp,dxp,dxv,bydxyp ?????
LAKES.f:302: C**** diffuse tracers using same KV as for heat?
LAKES.f:1114: C**** remove lakes that are too small ??
LAKES.f:1401: TAIJN(I,J,tij_lk1,:)=TAIJN(I,J,tij_lk1,:)+TRLAKEL(:,1) !*PLKICE?
LAKES.f:1402: TAIJN(I,J,tij_lk2,:)=TAIJN(I,J,tij_lk2,:)+TRLAKEL(:,2) !*PLKICE?
LANDICE_DRV.f:103: C**** Possibly this should be a function of in-situ freezing temp?
MODEL_COM.f:154: !@var ZATMO,HLAKE Topography arrays: elevation (m), lake depth (m) ???
MODEL_COM.f:158: !@var WFCS water field capacity of first ground layer (kg/m2) ???
MOMEN2ND.f:381: j = J_0S!STG ! why not staggered grid index?
MOMEN2ND.f:385: j = J_1!STG ! why not staggered grid index?
NCACC.f:31: c refresh some global parameters that change over the course of a run
NCACC.f:32: !!! should the following lines will be really removed?
NCACC.f:33: !!! there are no dparm_defs, iparm_defs in DEFACC.f any more.
NCACC.f:34: !!! Those who familiar with this code please check. I.A.
NCACC.f:35: !!! call iparm_defs
NCACC.f:36: !!! call dparm_defs
OCEAN.f:294: IF(MSI(I,J).eq.0) MSI(I,J)=MSINEW ! does this happen?
OCN_TRACER.f:98: t0m4(:,:,:)=(t0m4(:,:,:)*1d-3+1.)*trw0(n)
OCN_TRACER.f:99: tzm4(:,:,:)=0. ! tzm4(:,:,:)*1d-3*trw0(n) corrupted?
OCNGM.f:102: C**** ASY1(I,JM-1,L), ASY0(I,JM-1,LM1) and ASY3(I,2,L), ASY2(I,2,LM1)??
OCNKPP.f:1062: C**** Is this still necessary now that fluxes are saved?
PBL.f:346: c**** cannot update wsh without taking care that wsh used for tracers is
PBL.f:347: c**** the same as that used for q
PBL.f:348: c wsh = sqrt((u(1)-uocean)**2+(v(1)-vocean)**2+wstar2h)
PBL.f:349: wsm = wsh
PBL.f:355: C**** and GHY_DRV. This may need to be tracer dependent?
PBL.f:360: csgs wmin = wt ! minimum wind velocity (usually wt?)
PBL_SIMPLE_DRV.f:139: betas = rgas*tvs/pg ! ??? shouldn’t be PS ?
RAD_DRV.f:1114: AGESN(1)=SNOAGE(3,I,J) ! land ! ? why are these numbers
RAD_DRV.f:1115: AGESN(2)=SNOAGE(1,I,J) ! ocean ice so confusing ?
RADIATION.f:390: !sl K ,FTAUSL(33),TAUSL(33) ! input rather than output ?
RADIATION.f:407: !sl K ,FTAUSL,TAUSL ! input rather than output ?
RADIATION.f:1037: !? IF(LASTVC.GT.0) NRFUN=NRFN0
RADIATION.f:5866: DPF=0. ! initialise?
RADIATION.f:5897: DPF=0. ! initialise?
SNOW.f:374: k_ground = 3.4d0 !/* W K-1 m */ /* –??? */
SNOW.f:375: c_ground = 1.d5 !/* J m-3 K-1 */ /* — ??? */
SNOW.f:501: water_down = – evap_dt * delta_tsn_impl * dt ! ??
STRAT_DIAG.f:33: USE DYNAMICS, only : w=>conv ! I think this is right….?
STRAT_DIAG.f:42: C**** NOTE: AEP was a separate array but is now saved in AJL (pointer?)
STRAT_DIAG.f:290: C? FER1 does not include COSP, yet ER1 = 1/COSV*(FER1(J)-FER1(J-1))
STRAT_DIAG.f:585: C**** NOTE: AEP was a separate array but is now saved in AJL (pointer?)
STRAT_DIAG.f:589: c COMMON /PROGCB/ U,V,T,SX,SY,SZ,P,Q !not used?
STRAT_DIAG.f:618: C**** common needed to pass from XEP to individual arrays (better way??)
STRAT_DIAG.f:619: C**** Not clear how many of these are actually used.
STRAT_DIAG.f:620: cBMP COMMON /EPCOM/ FMY, FEY, FMZ, FEZ, FMYR, FEYR, FMZR, FEZR, COR,
STRAT_DIAG.f:621: cBMP * CORR, FER1, ER21, ER22, VR, WR, RX, UI, VI, WI, DUT, DUD
STRATDYN.f:57: !@var PK local P**Kapa array – should be done by DYNAMICS?
STRATDYN.f:247: C**** Fill in USURF,VSURF at poles (Shouldn’t this be done already?)
STRATDYN.f:276: C**** Fill in T,RHO at poles (again shouldn’t this be done already?)
STRATDYN.f:507: C**** Is this calculated in PBL?
SURFACE.f:166: MODDD=MOD(1+ITime/NDAY+NS,NIsurf) ! 1+ not really needed ??
TRACER_PRT.f:487: * ,jm,lm,0,lat_dg)
TRACER_PRT.f:487: !? * ,jm,lm,lm_req,lat_dg)
TRACER_SPECIAL_Lerner.f:1333: C**** We want to interpolate from an irregular to an irregular grid.
TRACER_SPECIAL_Lerner.f:1334: C**** This is a little sloppy because value at ltopx is not accurate
TRACER_SPECIAL_Lerner.f:1335: C**** But what WOULD be correct???
TRACERS_DRV.f:32: USE LANDICE_COM, only : trli0 ! should these be in tracer_com?
TRACERS_DRV.f:339: ntsurfsrc(n) = 1 ! (vegetation?)
TRACERS_DRV.f:3402: C**** for gradients defined on water mass (should be an option?)
TRACERS_DRV.f:4059: !@var SSFAC dummy variable (assumed units= kg water?)
TRACERS_DRV.f:4064: !@var WMXTR mixing ratio of water available for tracer condensation?
TRACERS_DRV.f:4172: !@var BELOW_CLOUD logical- is the current level below cloud?
TRDIAG_COM.f:420: go to 10 ! or should this be just a warning ??
VEG_DRV.f:132: do k=1,10 ! 11 ????
VEGETATION.f:148: !@var sigma Leaf scattering coefficient (?unitless).
VEGETATION.f:155: !@var k Canopy nitrogen extinction coefficient (?unitless).
VEGETATION.f:183: !@var rho Canopy reflectivity (?unitless).
VEGETATION.f:199: !@var m1 CO2 dependence of light-limited photosynthesis (?units).
VEGETATION.f:201: !@var m2 CO2 dependence of light-saturated photosynthesis (?units).
VEGETATION.f:203: !@var msat Nitrogen dependence of photosynthesis (?units).
VEGETATION.f:286: ! CO2 dependence of light-limited photosynthesis (?units).
VEGETATION.f:288: ! CO2 dependence of light-saturated photosynthesis (?units).
VEGETATION.f:290: ! Nitrogen dependence of photosynthesis (?units).
VEGETATION.f:329: ! CO2 dependence of light-limited photosynthesis (?units).
VEGETATION.f:331: ! CO2 dependence of light-saturated photosynthesis (?units).
VEGETATION.f:333: ! Nitrogen dependence of photosynthesis (?units).
VEGETATION.f:385: GPP=0.012D-6*Anet ! should be dependent on conductance ??
VEGETATION.f:434: ! Rule (Press et al., 19??).
VEGETATION.f:479: !@var alpha Intrinsic quantum efficiency (?units).
VEGETATION.f:481: !@var ka Chlorophyll extinction coefficient (?units).
VEGETATION.f:483: !@var n3 Ratio of foliage chlorophyll to N (?units).
VEGETATION.f:514: ! Check following OK. n3 is ratio of chlorophyll to foliage N (units?).
prtdag/jldag.f:82: xwon=twopi/(dlon*dble(im)) ! for wonderland model?
This use of “!” is very poor practice, seeing as N! means means N “factorial”.
Dan,
Do you have a summary of your impression of NASA/GISS code practices? One written so that it would be understandable to high school students? Something that highlights what you think the major problems are and a conclusion that you draw from them re: climate models, AGW, etc.
re:#5
stan,
I don’t have summary written about the GISS/NASA ModelE code. And, if I did it would not address all climate models or AGW. It would address only the coding of ModelE and not the model equations themselves. There are additional problems with the latter also. I have a couple of other posts here about the ModelE coding here and here.
I have a general document about some good software practices here.
The concept of a Design Review Team is discussed in that document. As mentioned there, because it would not be feasible to initially review every line of code, the team would be assigned to look at the coding for a few critical routines. Based on the results of that review, a decision regarding the code as a whole would be made. Based on what I’ve seen in the ModelE coding, my recommendation would be that the entire code would need to be reviewed in more detail.
[…] Pattern Matching in NASA/GISS ModelE Coding […]