----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : E Dated : June 13, 2008 ----------------------------------------------------------------------- The following code modifications address problems identified in the Version 5 CALPUFF modeling system currently approved for regulatory applications. Bulletin MCB-E(080613) outlines these changes and together with those in MCB-A(040716), MCB-B(051216), MCB-C(060804), and MCB-D(070622) represent all of the corrections made to these programs since EPA first approved CALPUFF as a Guideline model in 2003. Also, several changes are noted that do not alter model results. ----------------------------------------------------------------------- OVERVIEW ----------------------------------------------------------------------- CALUTILS.FOR (to Version 2.56, Level 080407) COORDLIB.FOR (to Version 1.99, Level 070921) -------------------------------------------- Components used in the CALPUFF system for coordinate conversions are packaged as the FORTRAN include-files: COORDLIB.FOR BLOCKDAT.CRD NIMA.CRD Furthermore, many subroutines used throughout the system are packaged as the FORTRAN include-file: CALUTILS.FOR This MCB includes changes to both COORDLIB and CALUTILS components. The following model versions incorporate these updates: CALMET (Version 5.812, Level 080512) CALUTILS & COORDLIB CALPUFF (Version 5.831, Level 080520) CALUTILS & COORDLIB CALPOST (Version 5.6396, Level 080407) CALUTILS & COORDLIB Several processors within the system also incorporate these components: BUOY (Version 1.247, Level 080407) CALUTILS & COORDLIB CTGCOMP (Version 2.251, Level 080407) CALUTILS CTGPROC (Version 2.683, Level 080407) CALUTILS & COORDLIB MAKEGEO (Version 2.291, Level 080407) CALUTILS PMERGE (Version 5.321, Level 080407) CALUTILS PXTRACT (Version 4.251, Level 080407) CALUTILS POSTUTIL(Version 1.58, Level 080407) CALUTILS PRTMET (Version 4.341, Level 080407) CALUTILS READ62 (Version 5.541, Level 080407) CALUTILS SMERGE (Version 5.571, Level 080407) CALUTILS TERREL (Version 3.686, Level 080407) CALUTILS & COORDLIB CALUTILS.FOR ------------ --- Problem Area 1 --- Exponential notation processing in ALTONU did not properly interpret a control-file entry without a decimal point. For example, a real number entered as 2E+03 was interpreted as 0.2E+03. This only applies to exponential (scientific) notation. Control files prepared with the GUI use the decimal point for all real numbers, and so are not affected by this bug. Modified: ALTONU COORDLIB.FOR ------------ -- Problem Area 1 -- A non-zero False Northing appropriate for UTM-S was used when converting S. hemisphere locations to UTM-N coordinates. This causes the coordinates returned to be in UTM-S. Modified: COORDS, PJINIT -- Problem Area 2 -- Three work-array variables and one output variable were not completely initialized, causing the program to halt when extra compiler-checking options are used. These do not affect results. Modified: COORDS CALPUFF (to Version 5.831, Level 080520) ---------------------------------------- -- Problem Area 1 -- When performing cavity sampling for PRIME downwash, restrict primary source calculations to receptors downwind of primary source and add screen for receptors located far to the side (no impact). Without this restriction, the model may halt with an attempted division by zero. Receptors upwind of the source are processed for cavity impacts starting with Version 5.8, Level 070623. Modified: CAV_SAMP -- Problem Area 2 -- Fix bug in wet flux calculation for sampling puffs (not slugs). Horizontal sampling factors were not recalculated if puff mass did not diffuse to the surface. But these factors are needed for the wet fluxes due to elevated puffs. The result is that wet fluxes contributed by elevated puffs do not have the correct horizontal puff-receptor distance relationship. Wet removal mass depletion calculations that influence concentrations are correct, as are wet fluxes from puffs in contact with the ground. Modified: CALCPF -- Problem Area 3 -- Add cap on sigma-z to avoid a floating-point error when computing virtuals, which halts a run. Default cap is set to 5000km, which is expected to have no influence on concentrations. Modified: COMPARM.PUF BLOCK DATA, READCF, COMP -- Problem Area 4 -- Add check for ATAN2(0.0,0.0) in FOGREC, as function will halt execution if both arguments are zero. Set flow direction to zero (north) in this case, and skip ATAN2 call. Modified: FOGREC -- Problem Area 5 -- Assign several undefined variables that do not affect results: - PROBLEM variable initialized to false in LN2FILL - Calm logical is not set in call to VCBAR from DRY, so a frozen puff treatment for dry depletion in calms is never triggered. This would have made a single calculation of the vertical distribution at the middle of the step in place of 3 calculations during a step. Logical is explicitly set to false to ensure current 3-step calculation is always done. - Initialize ISTA=ILQ=IPRECIP=0 in WET for debug output (they may not be computed) - Change J to I in debug output in VOLS, AREAS1/2 - Initialize INDEX=0 in POINTS1/2 for debug output - Initialize DUMMY=0.0 in RDTIEM3 Modified: LN2FILL, RDTIEM3 DRY, VCBAR, WET, VOLS, POINTS1/2, AREAS1/2 -- Problem Area 6 -- Relax requirement that the input restart file be from the same version and level of the code, and report a WARNING. Modified: RESTARTQ -- Problem Area 7 -- Refine mixing height adjustment to the extent of the layer that spans a puff (used for obtaining the transport wind). The layer is nominally the puff height +/- one sigma-z, but the top should not exceed the reflecting lid height for the mass distribution. The previous mixing height was selected as the top when the puff center exceeds the current mixing height, but this height might also be less than the puff center when there is a 2-layer mass distribution for partial penetration or fumigation. Use the maximum mixing height if this occurs, and add checks for a layer top less than a layer bottom to routines that define or use layer interfaces to obtain mean winds. Modified: PUFFDZ, ADVECT -- Problem Area 8 -- Treat case of falling puff in the procedure that determines the average transport during a step that includes gradual rise near a point source. The layer for averaging the wind should extend from the bottom of the puff at the start of the step to the top of the puff at the end of the step, regardless of whether the puff is rising or falling. Also, the gradual rise for point sources with Schulman-Scire downwash active must be explicitly addressed in RISEWIND because the GRISE calls do not include it. Modified: RISEWIND -- Problem Area 9 -- Procedure that determines the average transport during a step that includes gradual rise near a point source does not include the stack-tip downwash adjustment to the puff height. This adjustment applies to point sources that are not subject to PRIME or Schulman-Scire building downwash, when the MTIP option is used. The tip downwash adjustment is now applied within the gradual rise subroutine and not elsewhere. Modified: GRISE, PUFRECS, SLGRECS, PLGRECS, ZTRACE -- Problem Area 10 -- Replace source-based storage arrays related to plume rise and PRIME downwash tabulations with puff-based storage arrays, implemented via a direct access (DA) file. This allows all time-of-release tabulations to be accessed for meteorological periods that precede the current period. Previous implementation only provided the previous period rise tables for buoyant line sources, and otherwise used final rise properties for puffs from the previous period even when still within the gradual rise distance from the source. To RESTART file, add source tabulations stored for previous met periods and introduce dataset name and version record. Without the table for a previous met period in a restart file, table values of zero at start-up may be accessed and halt the run with a divide-by-zero. Additional checks on the index used to identify the met period are added. New: SRCTAB.PUF SWAPTAB, SRCTABIN, SRCTABOUT Modified: PARAMS.PUF, PT1.PUF, PT2.PUF, AR2.PUF, LN1.PUF, LN2.PUF POINTS1, POINTS2, AREAS2, LINES1, LINES2, AREAS1, VOLS, BCS1, INITPUF, SETPUF, PUFRECS, OPENOT, SWAP, COMP, GRISE, WAKE_TAB, RECSPEC0, RESTARTQ, RESTARTO, RESTARTI -- Problem Area 11 -- Fix logic to implement QA on MTILT. An errant elseif condition caused checks for NSPEC, SG, and MCTADJ to be skipped. MTILT Restrictions are: - MDRY = 1 - NSPEC = 1 (must be particle species as well) - sg = 0 GEOMETRIC STANDARD DEVIATION in Group 8 is set to zero for a single particle diameter - MCTADJ= 0,1,3 (Interaction with option 2 not supported) Modified: QAINP CALMET (to Version 5.812, Level 080512) --------------------------------------- -- Problem Area 1 -- The date test performed in the case of a simulation ending on the first hour of a year following a leap year was not done correctly and neither was the computation of dtinc for a run ending on the first hour of a new year. Modified: RDMM5 -- Problem Area 2 -- When no upper air observations are used in the simulation, CALMET uses the wrong coordinates for the precip stations when the coordinates are included in the control (.inp) file, as opposed to included in the precip.dat file. In that case CALMET does not convert the coordinates to the internal relative units used for precipitation interpolation. Modified: READCF -- Problem Area 3 -- CALMET stops searching for a first valid MM5 record after it has read through the first MM5 file listed as input data, even if there subsequent MM5 data files ,possibly issuing a "ran out of MM5 data before start" error message too soon. Modified: RDMM5 --- Problem Area 4 --- For overwater gridpoints, as defined by JWAT1-JWAT2, and for cases when there is no SEA.DAT file, the RH overwater default was setto 80% instead of 100%. This was inconsistent with the RH value overwater in the rest of the code and the default when there is a SEA.DAT file with missing RH observations. Modified: SURFVAR --- Problem Area 5 --- In NOOBS mode using MM4 data, when the simulation started during convective conditions (e.g.in high latitudes or when the base time zone is different from local time zone), the prognostic sounding used to compute the convective mixing height growth was not yet initialized, causing the simulation to stop. Modified : RDMM4 --- Problem Area 6 --- In NOOBs-temperature mode (ITPROG=2),at the first timestep, the surface temperatures are initialized with not-yet defined variables, which could cause the simulation to stop with some compilers. However results are not affected because the surface temperatures are correctly updated later on. --- Problem Area 7 --- In NOOBS mode, there is no check made on the values of ISURFT or IUPT (because they do not matter). However using a MOD5 CALMET.INP file that has negative values with MOD6 can trigger a non-consistent simulation (in particular, ISURFT=-1 and IUPT=-1 trigger 2D surface temperatures and lapse rates instead of domain-average ones in MOD6) -- Problem Area 8 -- When the coordinates of the surface stations are provided in the SURF.DAT file and CALMET also uses precipitation observations, CALMET checks for precipitation IDs in the PRECIP.DAT file even when those coordinates are only listed in the control (.inp) file Modified: READHD -- Problem Area 9 -- When cloud ceiling heights are computed based on prognostic cloud mixing ratios (ICLOUD=3), prognostic data is of MM5 type (IPROG>5), prognostic timestep is longer than CALMET timestep (ISTEPPG>1), and MM5 cloudy skies are replaced by clear skies at the next MM5 timestep,the ceiling heights at intermediate hours were underestimated. Modified: RDMM5 -- Problem Area 10 -- Prognostic ceiling heights (ICLOUD=3) were not correctly interpolated to the CALMET gridpoints Modified: RDMM5 -- Problem Area 11 -- No vertical extrapolation of temperature from lowest 3D.DAT level to lower CALMET levels was performed overwater when ITWPROG=2. Constant T profile was assumed below lowest 3D.DAT level instead of the intented interpolation between SST and lowest 3D.DAT level temperature. Only applies to overwater gridpoints when prognostic temperatures are used. Modified: RDMM5 ----------------------------------------------------------------------- ----------------------------------------------------------------------- 1. CALPUFF Changes (Version 5.831, Level 080520) ----------------------------------------------------------------------- *********************************************************************** --- Problem Area 1 -- - When performing cavity sampling for PRIME downwash, restrict primary source calculations to receptors downwind of primary source. *********************************************************************** ------------------------------------ Subroutine CAV_SAMP (Problem Area 1) ------------------------------------ BEFORE change: (a) -------------- data rt2pi/2.5066283/, zero/0.0/, pi/3.1415927/ AFTER change: (a) -------------- data rt2pi/2.5066283/, zero/0.0/, pi/3.1415927/ data sigcut/12.0/ BEFORE change: (b) -------------- c --- Primary source c ------------------ if(ipositn.GE.3) then c --- Primary source contributes outside of cavity c --- Set primary source height call INTERTAB(xr,ntr,xtr,i1,i2,factor) ht=ztr(i1)+factor*(ztr(i2)-ztr(i1)) fy0=(EXP(-0.5*(yr/sy)**2))/(rt2pi*sy) cbyq0=fy0*VCOUP(icode,zr,ht,sz,hlid)/Urh endif AFTER change: (b) -------------- c --- Primary source c ------------------ c --- Primary source does not impact upwind if(ipositn.GE.3 .AND. xr.GT.0.0) then c --- Drop calculation for receptors too far to the side if(ABS(yr) .LE. (sigcut*sy)) then c --- Primary source contributes outside of cavity c --- Set primary source height call INTERTAB(xr,ntr,xtr,i1,i2,factor) ht=ztr(i1)+factor*(ztr(i2)-ztr(i1)) fy0=(EXP(-0.5*(yr/sy)**2))/(rt2pi*sy) cbyq0=fy0*VCOUP(icode,zr,ht,sz,hlid)/Urh endif endif *********************************************************************** --- Problem Area 2 --- Horizontal sampling factors for puffs were not recalculated for wet flux calculation if puff mass did not diffuse to the surface. *********************************************************************** ---------------------------------- Subroutine CALCPF (Problem Area 2) ---------------------------------- BEFORE change: (a) -------------- if(f1.GT.0.) then c --- Sample horizontal distribution call pfsamp(samxo,samyo,samxn,samyn,xr,yr,syrgu,samd2, & xi1,xi2) c c --- Normalize horizontal sampling result by 1/(2*pi*syrec^2) c --- and scale by the Area Source Distribution Factor c --- (nominally 1.0); syrec is in meters; samdm=SQRT(samd2)*delsam dinv=ASDF(ldbhr,samdm,syrec)/(6.2831853*syrec**2) t2=xi2*dinv t1=xi1*dinv-t2 c --- (QOLD*T1+QNEW*T2 IS THE HORIZONTAL TERM * SOURCE TERM OF c --- THE GAUSSIAN EQN.) c c --- Concentration and Dry Fluxes xtemp=f1*tfract do ipol=1,nspec conc=vdpvd(ipol)*xtemp*(qold(ipol)*t1+qnew(ipol)*t2) chisam(isamp,jsamp,ipol)=chisam(isamp,jsamp,ipol)+conc if(mdry.EQ.1) then dfsam(isamp,jsamp,ipol)=dfsam(isamp,jsamp,ipol)+ & (q01dry(ipol,1)*t1+q01dry(ipol,2)*t2)* & vd(ipol)*vdpvd(ipol)*xtemp endif enddo endif AFTER change: (a) -------------- c --- Sample horizontal distribution if puff impact reaches the c --- surface due to diffusion or washout if(f1.GT.0. .OR. mwet.EQ.1) then call pfsamp(samxo,samyo,samxn,samyn,xr,yr,syrgu,samd2, & xi1,xi2) c --- Normalize horizontal sampling result by 1/(2*pi*syrec^2) c --- and scale by the Area Source Distribution Factor c --- (nominally 1.0); syrec is in meters; samdm=SQRT(samd2)*delsam dinv=ASDF(ldbhr,samdm,syrec)/(6.2831853*syrec**2) t2=xi2*dinv t1=xi1*dinv-t2 c --- (QOLD*T1+QNEW*T2 IS THE HORIZONTAL TERM * SOURCE TERM OF c --- THE GAUSSIAN EQN.) endif c --- Concentration and Dry Fluxes if(f1.GT.0.) then xtemp=f1*tfract do ipol=1,nspec conc=vdpvd(ipol)*xtemp*(qold(ipol)*t1+qnew(ipol)*t2) chisam(isamp,jsamp,ipol)=chisam(isamp,jsamp,ipol)+conc if(mdry.EQ.1) then dfsam(isamp,jsamp,ipol)=dfsam(isamp,jsamp,ipol)+ & (q01dry(ipol,1)*t1+q01dry(ipol,2)*t2)* & vd(ipol)*vdpvd(ipol)*xtemp endif enddo endif BEFORE change: (b) -------------- if(f1.GT.0.) then c --- Sample horizontal distribution call pfsamp(xold,yold,xnew,ynew,xng(i),yng(i),syrgu,d2, & xi1,xi2) c c --- Normalize horizontal sampling result by 1/(2*pi*syrec^2) c --- and scale by the Area Source Distribution Factor c --- (nominally 1.0); syrec is in meters; dm=SQRT(d2)*dgrid dinv=ASDF(ldbhr,dm,syrec)/(6.2831853*syrec**2) t2=xi2*dinv t1=xi1*dinv-t2 c c --- Compute concentration xtemp=f1*tfract do ipol=1,nspec conc=vdpvd(ipol)*xtemp*(qold(ipol)*t1+qnew(ipol)*t2) chirec(i,ipol)=chirec(i,ipol)+conc if(ipol.EQ.2 .AND. LADTFOG) then c --- FOG: SUM the T excess (ipol=2) at receptor but c --- do not allow sum to exceed excess at release chirec(i,ipol)=AMIN1(txsmxfog,chirec(i,ipol)) endif enddo endif AFTER change: (b) -------------- c --- Sample horizontal distribution if puff impact reaches the c --- surface due to diffusion or washout if(f1.GT.0. .OR. mwet.EQ.1) then call pfsamp(xold,yold,xnew,ynew,xng(i),yng(i),syrgu,d2, & xi1,xi2) c --- Normalize horizontal sampling result by 1/(2*pi*syrec^2) c --- and scale by the Area Source Distribution Factor c --- (nominally 1.0); syrec is in meters; dm=SQRT(d2)*dgrid dinv=ASDF(ldbhr,dm,syrec)/(6.2831853*syrec**2) t2=xi2*dinv t1=xi1*dinv-t2 endif if(f1.GT.0.) then c --- Compute concentration xtemp=f1*tfract do ipol=1,nspec conc=vdpvd(ipol)*xtemp*(qold(ipol)*t1+qnew(ipol)*t2) chirec(i,ipol)=chirec(i,ipol)+conc if(ipol.EQ.2 .AND. LADTFOG) then c --- FOG: SUM the T excess (ipol=2) at receptor but c --- do not allow sum to exceed excess at release chirec(i,ipol)=AMIN1(txsmxfog,chirec(i,ipol)) endif enddo endif BEFORE change: (c) -------------- if(f1.GT.0.) then c --- Sample horizontal distribution call pfsamp(xold,yold,xnew,ynew,xrct(i),yrct(i),syrgu,d2, & xi1,xi2) c c --- Normalize horizontal sampling result by 1/(2*pi*syrec^2) c --- and scale by the Area Source Distribution Factor c --- (nominally 1.0); syrec is in meters; dm=SQRT(d2)*dgrid dinv=ASDF(ldbhr,dm,syrec)/(6.2831853*syrec**2) t2=xi2*dinv t1=xi1*dinv-t2 c c --- Compute concentration xtemp=f1*tfract do ipol=1,nspec conc=vdpvd(ipol)*xtemp*(qold(ipol)*t1+qnew(ipol)*t2) chict(i,ipol)=chict(i,ipol)+conc enddo endif AFTER change: (c) -------------- c --- Sample horizontal distribution if puff impact reaches the c --- surface due to diffusion (no Wet fluxes at CTSG receptors) if(f1.GT.0.) then call pfsamp(xold,yold,xnew,ynew,xrct(i),yrct(i),syrgu,d2, & xi1,xi2) c --- Normalize horizontal sampling result by 1/(2*pi*syrec^2) c --- and scale by the Area Source Distribution Factor c --- (nominally 1.0); syrec is in meters; dm=SQRT(d2)*dgrid dinv=ASDF(ldbhr,dm,syrec)/(6.2831853*syrec**2) t2=xi2*dinv t1=xi1*dinv-t2 endif if(f1.GT.0.) then c --- Compute concentration xtemp=f1*tfract do ipol=1,nspec conc=vdpvd(ipol)*xtemp*(qold(ipol)*t1+qnew(ipol)*t2) chict(i,ipol)=chict(i,ipol)+conc enddo endif *********************************************************************** -- Problem Area 3 -- Add cap on sigma-z to avoid a floating-point error when computing virtuals, which halts a run. Default cap is set to5000km. *********************************************************************** -------------------------------------- Subroutine BLOCK DATA (Problem Area 3) -------------------------------------- BEFORE change: -------------- c --- Initialize the index pointer for the variable emissions scaling c --- factor array (VQFAC) data iqnext/1/ AFTER change: -------------- c --- Initialize the index pointer for the variable emissions scaling c --- factor array (VQFAC) data iqnext/1/ c --- Default for maximum sigma-z allowed (sigma-z cap) data szcap_m/5.0e06/ ---------------------------------- Subroutine READCF (Problem Area 3) ---------------------------------- BEFORE change: (a) -------------- c Common block /COMPARM/: c XMXLEN, MXNEW, XSAMLEN, MXSAM, SL2PF, CDIV, c WSCALM,WSCAT(5),TCCAT(11),VQFAC(12,mxq12),IQNUM(5),IQNEXT, c SVMIN(6,2), SWMIN(6,2), SYMIN, SZMIN, XMINZI, XMAXZI, AFTER change: (a) -------------- c Common block /COMPARM/: c XMXLEN, MXNEW, XSAMLEN, MXSAM, SL2PF, CDIV, c WSCALM,WSCAT(5),TCCAT(11),VQFAC(12,mxq12),IQNUM(5),IQNEXT, c SVMIN(6,2),SWMIN(6,2),SYMIN,SZMIN,SZCAP_M,XMINZI,XMAXZI, BEFORE change: (b) -------------- l 'MDEPBC','HTMINBC','RSAMPBC','MXAGEHR', 5*' ', AFTER change: (b) -------------- l 'MDEPBC','HTMINBC','RSAMPBC','MXAGEHR','SZCAP_M', 4*' ', BEFORE change: (c) -------------- data ivleng2/ i 13*1, 47*0, j 60*0, k 1,2*12,4*1,4*12, 49*0, l 16*1,2*12,1,5,6,2,6,1,1,24,13*1,2,2*1,2*mxss, l 3,4*1,mxspec,5*1, 5*0, AFTER change: (c) -------------- data ivleng2/ i 13*1, 47*0, j 60*0, k 1,2*12,4*1,4*12, 49*0, l 16*1,2*12,1,5,6,2,6,1,1,24,13*1,2,2*1,2*mxss, l 3,4*1,mxspec,6*1, 4*0, BEFORE change: (d) -------------- data ivtype2/ i 11*1, 2*2, 47*0, j 60*0, k 2, 5*1, 2, 4*1, 49*0, l 1,2*2,2*1,2*2,1,2,3*1,2,11*1,2*2,3*1,3*2,8*1, l 2,2,4*1,2,4*1,2,2*1,2, 5*0, AFTER change: (d) -------------- data ivtype2/ i 11*1, 2*2, 47*0, j 60*0, k 2, 5*1, 2, 4*1, 49*0, l 1,2*2,2*1,2*2,1,2,3*1,2,11*1,2*2,3*1,3*2,8*1, l 2,2,4*1,2,4*1,2,2*1,2,1, 4*0, BEFORE change: (e) -------------- 7 MDEPBC,HTMINBC,RSAMPBC,MXAGEHR, 8 idum,idum,idum,idum,idum) c c --- Reset calm wind speed threshold by "small" amount wscalm=0.99999*wscalm AFTER change: (e) -------------- 7 MDEPBC,HTMINBC,RSAMPBC,MXAGEHR,SZCAP_M, 8 idum,idum,idum,idum) c c --- Reset calm wind speed threshold by "small" amount wscalm=0.99999*wscalm BEFORE change: (f) -------------- write(io6,*)'symin = ',symin write(io6,*)'szmin = ',szmin write(io6,*)'xminzi = ',xminzi write(io6,*)'xmaxzi = ',xmaxzi AFTER change: (f) -------------- write(io6,*)'symin = ',symin write(io6,*)'szmin = ',szmin write(io6,*)'szcap_m = ',szcap_m write(io6,*)'xminzi = ',xminzi write(io6,*)'xmaxzi = ',xmaxzi -------------------------------- Subroutine COMP (Problem Area 3) -------------------------------- BEFORE change: (a) -------------- logical ltest integer nsource(9) AFTER change: (a) -------------- logical ltest logical lszcap integer nsource(9) BEFORE change: (b) -------------- c --- Set logical for reading PROFILE.DAT as a secondary met file c --- (e.g. for turbulence data or inversion strength) AFTER change: (b) -------------- c --- Set logical for applying cap to sigma-z to protect virtuals if(szcap_m.GT.0.0) then lszcap=.TRUE. else lszcap=.FALSE. endif c --- Set logical for reading PROFILE.DAT as a secondary met file c --- (e.g. for turbulence data or inversion strength) BEFORE change: (c) -------------- c c --- Transfer ozone concentration for chemistry ldb2=.false. if(mchem.EQ.1 .OR. mchem.EQ.3 .OR. mchem.EQ.4) & call GETOZ(ixmid,iymid,dgrid,ldb2,chioz) AFTER change: (c) -------------- c --- CAP c --- Place a cap on the size of the sigma-z at the start of the step c --- to keep it from growing (much) beyond a very large value c --- This keeps sigma within the range of the host curves if(lszcap) then sigzb(ii)=AMIN1(sigzb(ii),szcap_m) sigze(ii)=AMIN1(sigze(ii),szcap_m) if(ldbhr) then write(io6,*) ' ----- SIGMA-Z CAP (m): ',szcap_m write(io6,*) ' ----- sigzb,sigze = ', & sigzb(ii),sigze(ii) endif endif c c --- Transfer ozone concentration for chemistry ldb2=.false. if(mchem.EQ.1 .OR. mchem.EQ.3 .OR. mchem.EQ.4) & call GETOZ(ixmid,iymid,dgrid,ldb2,chioz) *********************************************************************** --- Problem Area 4 --- Add check for ATAN2(0.0,0.0) in FOGREC, as function will halt execution if both arguments are zero. *********************************************************************** ---------------------------------- Subroutine FOGREC (Problem Area 4) ---------------------------------- BEFORE change: -------------- c --- Convert vector components to a flow direction flow0=90.-(ATAN2(vadv,uadv))/d2rad flow0=AMOD(flow0,360.) AFTER change: (a) -------------- c --- Convert vector components to a flow direction if(vadv.EQ.0.0 .AND. vadv.EQ.0.0) then c --- Point calms to the north flow0=0.0 else flow0=90.-(ATAN2(vadv,uadv))/d2rad flow0=AMOD(flow0,360.) endif *********************************************************************** --- Problem Area 5 --- Assign several undefined variables that do not affect results *********************************************************************** ----------------------------------- Subroutine LN2FILL (Problem Area 5) ----------------------------------- BEFORE change: -------------- logical problem,lqneg AFTER change: -------------- logical problem,lqneg c --- Initialize PROBLEM problem=.FALSE. ----------------------------------- Subroutine RDTIEM3 (Problem Area 5) ----------------------------------- BEFORE change: -------------- data iv54/54/ AFTER change: -------------- data iv54/54/ c --- Initialize dummy variable dummy=0.0 ------------------------------- Subroutine DRY (Problem Area 5) ------------------------------- BEFORE change: (a) -------------- logical ldbhr,ldb2 AFTER change: (a) -------------- logical ldbhr,ldb2 logical lcalmf c --- Explicitly disable single-point VCBAR sampling by setting local c --- variable to false ("calm" sampling was never fully implemented) lcalmf=.false. BEFORE change: (b) -------------- c --- Puff is Gaussian in the vertical mfact0=icode-iicode f1=vcbar(jp,mfact0,hlid,ppcoef,lcalm) const=tsamp*f1 endif AFTER change: (b) -------------- c --- Puff is Gaussian in the vertical mfact0=icode-iicode f1=vcbar(jp,mfact0,hlid,ppcoef,lcalmf) const=tsamp*f1 endif ------------------------------- Function VCBAR (Problem Area 5) ------------------------------- BEFORE change: -------------- c --- Set list file output switch for testing ldbg=.FALSE. AFTER change: -------------- c --- Set list file output switch for testing ldbg=.FALSE. c --- Require lcalm = false if(LCALM) then write(io6,*)'VCBAR: LCALM must be FALSE, LCALM = ',lcalm stop 'Halted in VCBAR -- See list file.' endif ------------------------------- Subroutine WET (Problem Area 5) ------------------------------- BEFORE change: -------------- c --- Missing value indicators for integer, real variables data imiss/9999/ AFTER change: -------------- c --- Missing value indicators for integer, real variables data imiss/9999/ c --- Initialize variables written as debug output that may not c --- be assigned ista=0 ilq=0 iprecip=0 -------------------------------- Subroutine VOLS (Problem Area 5) -------------------------------- BEFORE change: -------------- write(io6,208)j,np,iru,dt,sqrts 208 format(5x,'J: ',i5,2x,'NP: ',i5,2x, AFTER change: -------------- write(io6,208)i,np,iru,dt,sqrts 208 format(5x,'I: ',i5,2x,'NP: ',i5,2x, ----------------------------------- Subroutine POINTS1 (Problem Area 5) ----------------------------------- BEFORE change: -------------- c --- Initialize local ISC building downwash variables hb=0.0 AFTER change: -------------- c --- Initialize local ISC building downwash variables index=0 hb=0.0 ----------------------------------- Subroutine POINTS2 (Problem Area 5) ----------------------------------- BEFORE change: -------------- c --- Initialize local ISC building downwash variables hb=0.0 AFTER change: -------------- c --- Initialize local ISC building downwash variables index=0 hb=0.0 ---------------------------------- Subroutine AREAS1 (Problem Area 5) ---------------------------------- BEFORE change: -------------- write(io6,208)j,np,iru,dt,sqrts 208 format(5x,'J: ',i5,2x,'NP: ',i5,2x, AFTER change: -------------- write(io6,208)i,np,iru,dt,sqrts 208 format(5x,'I: ',i5,2x,'NP: ',i5,2x, ---------------------------------- Subroutine AREAS2 (Problem Area 5) ---------------------------------- BEFORE change: -------------- write(io6,208)j,np,iru,dt,sqrts 208 format(5x,'J: ',i5,2x,'NP: ',i5,2x, AFTER change: -------------- write(io6,208)i,np,iru,dt,sqrts 208 format(5x,'I: ',i5,2x,'NP: ',i5,2x, *********************************************************************** --- Problem Area 6 --- Relax requirement that the input restart file be from the same version and level of the code, and report a WARNING. *********************************************************************** ------------------------------------ Subroutine RESTARTQ (Problem Area 6) ------------------------------------ BEFORE change: (a) -------------- logical problem problem=.FALSE. c --- File is already open; read first record into local variables read(io3) ver0,level0,npuffs0,nspec0,ndathr0,nar10,nar20, & nln10,nln20,npt10,npt20,nvl10,nvl20 c --- Check against variables for current run if(ver.NE.ver0 .OR. level.NE.level0) problem=.TRUE. AFTER change: (a) -------------- logical problem, warning problem=.FALSE. warning=.FALSE. c --- File is already open; read first record into local variables read(io3) ver0,level0,npuffs0,nspec0,ndathr0,nar10,nar20, & nln10,nln20,npt10,npt20,nvl10,nvl20 c --- Check against variables for current run if(ver.NE.ver0 .OR. level.NE.level0) warning=.TRUE. BEFORE change: (b) -------------- if(problem) then c --- Report data and quit write(io6,*)'FATAL Problem in Restart File!' write(io6,*)' --- Run Data / Restart Data' write(io6,*)'Version: ',ver,ver0 write(io6,*)'Level : ',level,level0 write(io6,*)'Species: ',nspec,nspec0 AFTER change: (b) -------------- if(warning) then c --- Report data and continue write(io6,*) write(io6,*)'WARNING: Header difference in Restart File!' write(io6,*)' --- Run Data / Restart Data' write(io6,*)'Version: ',ver,ver0 write(io6,*)'Level : ',level,level0 write(io6,*) write(*,*) write(*,*) 'WARNING reported in RESTARTQ -- see list file.' endif if(problem) then c --- Report data and quit write(io6,*) write(io6,*)'FATAL Problem in Restart File!' write(io6,*)' --- Run Data / Restart Data' write(io6,*)'Species: ',nspec,nspec0 *********************************************************************** -- Problem Area 7 -- Refine mixing height adjustment to the extent of the layer that spans a puff (used for obtaining the transport wind). *********************************************************************** ------------------------------------ Subroutine PUFFDZ (Problem Area 7) ------------------------------------ BEFORE change: -------------- zbot=htmet-sigmaz if(zbot.LT.0.) zbot=0.0 ztop=htmet+sigmaz if(icode.NE.3 .and. icode.NE.13 .AND. istab.LE.4) then if(dpbl.GT.htmet) then ztop=amin1(ztop,dpbl) else ztop=amin1(ztop,ziold(ii)) endif endif AFTER change: -------------- zbot=htmet-sigmaz if(zbot.LT.0.) zbot=0.0 ztop=htmet+sigmaz c --- For non-stable surface layer, check for lid cap on ztop if(icode.NE.3 .and. icode.NE.13 .AND. istab.LE.4) then if(dpbl.GT.htmet) then ztop=amin1(ztop,dpbl) elseif(ziold(ii).GT.htmet) then ztop=amin1(ztop,ziold(ii)) elseif(zimax(ii).GT.htmet) then ztop=amin1(ztop,zimax(ii)) endif endif ------------------------------------ Subroutine ADVECT (Problem Area 7) ------------------------------------ BEFORE change: -------------- if(ldbg) then write(io6,*) write(io6,*)'ADVECT: metfm = ',metfm write(io6,*)'zbot,ztop,ht(release) ',zbot,ztop,ht endif AFTER change: -------------- if(ldbg) then write(io6,*) write(io6,*)'ADVECT: metfm = ',metfm write(io6,*)'zbot,ztop,ht(release) ',zbot,ztop,ht endif c --- Stop if ztop 1.001 if (dtinc.gt.1.001) dtinc = 1.*(ndathr-100+24-it1)/isteppg AFTER change: (b) ------------- c --- Interpolate in time c --- dtinc: 0.<=real<1 during run; dtinc possibly =1 in last step of run) call deltsec(it1,0,ndathr,0,incsec) dtinc=incsec/(isteppg*3600.) c --- Enforce bounds if (dtinc.lt.0.or.dtinc.gt.1.) then write(6,*)'STOP in RDMM5- DTINC out of bounds' write(io6,*)'STOP in RDMM5- DTINC out of bounds' write(io6,*)'Current Time (YYYYJJJHH) (LST): ',ndathr write(io6,*)'Current M3D times in memory (LST): ',it1,it2 write(io6,*)'Value of dtinc: ',dtinc STOP endif *********************************************************************** --- Problem Area 2 --- When no upper air observations are used in the simulation, CALMET use the wrong coordinates for the precip stations when the coordinates are included in the control (.inp) file, as opposed to included in the precip.dat file. In that case CALMET does not convert the coordinates to the internal relative units used for precipitation interpolation by CALMET *********************************************************************** ---------------------------------- Subroutine READCF (Problem Area 2) ---------------------------------- Change the logical variable triggering the conversion to internal relative coordinates of precip stations from LCFUPR to LCFPRC BEFORE change: -------------- c if(LCFUPR) then do 1610 i=1,npsta c c --- convert from km to m, and into relative coordinates AFTER change: -------------- if(LCFPRC) then do 1610 i=1,npsta c c --- convert from km to m, and into relative coordinates *********************************************************************** --- Problem Area 3 --- CALMET stops searching for a first valid MM5 record after it has read through the first MM5 file listed as input data, even if there subsequent MM5 data files, possibly issuing a "ran out of MM5 data before start" error message too soon. *********************************************************************** --------------------------------- Subroutine RDMM5 (Problem Area 3) --------------------------------- Continue searching for the first valid MM5 record throughout the multiple MM5 files BEFORE change: -------------- c --- ran out of data, at least in the current file 999 if (iskip .eq. 1) then write(io6,*) ' ran out of MM5 data before start!' stop c --- Multiple MM5.DAT: open next MM5.DAT file (if available) c --- and resume loop (at point where overlapping hours can be skipped) else if (nfm3d.lt.nm3d) then nfm3d=nfm3d+1 if (imm53d.eq.0) call rdhd51(nfm3d) if (imm53d.eq.1) call rdhd52(nfm3d) if (imm53d.eq.2) call rdhd53(nfm3d,itwprog) goto 1 else write(io6,*) 'ran out of MM5 data during run...' write(io6,*) 'On (yyyyjulhr): ',ndathr stop end if AFTER change: -------------- c --- ran out of data, at least in the current file 999 continue c --- Multiple MM5.DAT: open next MM5.DAT file (if available) c --- and resume loop (at point where overlapping hours can be skipped) if (nfm3d.lt.nm3d) then nfm3d=nfm3d+1 if (imm53d.eq.0) call rdhd51(nfm3d) if (imm53d.eq.1) call rdhd52(nfm3d) if (imm53d.eq.2) call rdhd53(nfm3d,itwprog) goto 1 else if (iskip .eq. 1) then write(io6,*) ' ran out of MM5 data before start!' else write(io6,*) 'ran out of MM5 data during run...' write(io6,*) 'On (yyyyjulhr): ',ndathr endif stop end if ***************************************************************** -- Problem Area 4 -- For overwater gridpoints, as defined by JWAT1-JWAT2, and for cases when there is no SEA.DAT file, the RH overwater default was set at 80% instead of 100%. This was inconsistent with the default RH value overwater in the rest of the code and when there is a SEA.DAT file with missing RH observations ***************************************************************** ------------------------------------ Subroutine SURFVAR (Problem Area 4) ------------------------------------ Set RH to 100% overwater when there are no overwater stations (i.e. no SEA.DAT file) BEFORE change: -------------- c --- if no valid site, assume 80% RH (overland only -Overwater:100%) if (numsta .eq. 0) then rhsurf=80. goto 9000 end if AFTER change: ------------- c --- if no valid site, assume 80% RH (overland only -Overwater:100%) if (numsta .eq. 0) then c --- 80% overland and 100% overwater if(lndwat.eq.1) then c numsta=0 can happen for overwater pts if there are c no SEA.DAT rhsurf=100. else rhsurf=80. endif goto 9000 end if ***************************************************************** -- Problem Area 5 -- In NOOBS mode using MM4 data, when the simulation started during convective conditions (e.g.in high latitudes or when the base time zone is different from local time zone), the prognostic sounding used to compute the convective mixing height growth was not yet initialized, causing the simulation to stop. ***************************************************************** --------------------------------- Subroutine RDMM4 (Problem Area 5) --------------------------------- Use the current hour MM4 soundings to compute the convective mixing height growth at the first timestep (subsequent steps keep using the previous hour soundings) BEFORE change: -------------- zl(i,j,1) = 0. tz(i,j,1)=tsurf(igrab(i,j,1),jgrab(i,j,1)) 125 continue AFTER change: ------------- zl(i,j,1) = 0. tz(i,j,1)=tsurf(igrab(i,j,1),jgrab(i,j,1)) c --- Fill in temp soundings with current value for the first valid timestep c --- to avoid all zeroes in MIXDT2 if first timestep is convective if (ifirstpg.eq.1) then nlevag0=nlevag1 if(i.eq.nx.and.j.eq.ny) ifirstpg=2 do 6622 k = 1,nlevag0 tz0(i,j,k)=tz(i,j,k) zl0(i,j,k)=zl(i,j,k) 6622 continue endif 125 continue ***************************************************************** --- Problem Area 6 --- In NOOBs-temperature mode (ITPROG=2),at the first timestep, the surface temperatures are initialized with not-yet defined variables, which could cause the simulation to stop with some compilers. However results are not affected because the surface temperatures are correctly updated later on. ***************************************************************** --------------------------------- Subroutine SURFVAR (Problem Area 6) --------------------------------- Initialize surface temperatures with default temperatures at first call to SURFVAR BEFORE change: -------------- if ( itprog .eq. 2 ) then temp2d(i,j) = tprog(i,j,1) goto 50 endif AFTER change: ------------- if ( itprog .eq. 2 ) then if(iall.ne.2)temp2d(i,j) = tprog(i,j,1) goto 50 endif ***************************************************************** --- Problem Area 7 --- n NOOBS mode, there is no check made on the values of ISURFT or IUPT (because they do not matter). However using a MOD5 CALMET.INP file that has negative values with MOD6 can trigger a non-consistent simulation (in particular, ISURFT=-1 and IUPT=-1 trigger 2D surface temperatures and lapse rates instead of domain-average ones in MOD6) ***************************************************************** --------------------------------- Subroutine READCF (Problem Area 7) --------------------------------- Do not allow negative values of ISURFT or IUPT in noobs mode BEFORE change: -------------- if(noobs.ne.2 .AND. idiopt(1).eq.0 .AND. & (isurft.lt.1.or.isurft.gt.nssta)) then write(io6,*) write(io6,*) 'READCF: Error in Input Group 5' write(io6,*) 'Invalid sfc station ISURFT =',isurft write(io6,*) 'Number must be 1 to ',nssta lerrcf=.TRUE. endif c c frr (09/01) perform QA checks only if matters (i.e. not in noobs mode) if(noobs .eq. 0) then if(idiopt(2).eq.0 .AND. (iupt.lt.1.or.iupt.gt.nusta))then write(io6,*) write(io6,*) 'READCF: Error in Input Group 5' write(io6,*) 'Invalid upper station IUPT =',iupt write(io6,*) 'Number must be 1 to ',nusta lerrcf=.TRUE. endif c if(iextrp .eq. 0 .or. iextrp .gt. 4 .or. iextrp .lt. -4) then AFTER change: ------------- c --- Test should be based on ITPROG not NOOBS c if(noobs.ne.2 .AND. idiopt(1).eq.0 .AND. if(itprog.ne.2 .AND. idiopt(1).eq.0 .AND. & (isurft.lt.1.or.isurft.gt.nssta)) then write(io6,*) write(io6,*) 'READCF: Error in Input Group 5' write(io6,*) 'Invalid sfc station ISURFT =',isurft write(io6,*) 'Number must be 1 to ',nssta lerrcf=.TRUE. endif c --- Do not accept negative values for ISURFT because they have c --- a meaning in MOD6 (if not in MOD5) if (itprog.eq.2.and.isurft.lt.0) then write(io6,*) write(io6,*) 'READCF: Error in Input Group 5' write(io6,*) 'Invalid ISURFT =',isurft write(io6,*) 'ISURFT cannot be negative' lerrcf=.TRUE. endif c c Value check of IUPT should be based on ITPROG not NOOBs (071204) c Do not allow negative IUPT (071204) if(idiopt(2).eq.0.and.itprog.eq.0 .AND. : (iupt.lt.1.or.iupt.gt.nusta))then write(io6,*) write(io6,*) 'READCF: Error in Input Group 5' write(io6,*) 'Invalid upper station IUPT =',iupt write(io6,*) 'Number must be 1 to ',nusta lerrcf=.TRUE. endif if (itprog.gt.0.and.iupt.lt.0) then write(io6,*) write(io6,*) 'READCF: Error in Input Group 5' write(io6,*) 'Invalid IUPT =',iupt write(io6,*) 'IUPT cannot be negative' lerrcf=.TRUE. endif c if(noobs .eq. 0) then if(iextrp .eq. 0 .or. iextrp .gt. 4 .or. iextrp .lt. -4) then *************************************************************************** -- Problem Area 8 -- When the coordinates of the surface stations are provided in the SURF.DAT file and CALMET also uses precipitation observations, CALMET checks for precipitation IDs in the PRECIP.DAT file even when those coordinates are only listed in the control (.inp) file *************************************************************************** ---------------------------------- Subroutine READHD (Problem Area 8) ---------------------------------- Change the logical variable triggering the check on precipitation station IDs from LCFSFC to LCFPRC BEFORE change: -------------- if(LCFSFC) then c --- precip. station IDs from precip. file must match those in c --- control file nmatch=0 AFTER change ------------- if(LCFPRC) then c --- precip. station IDs from precip. file must match those in c --- control file nmatch=0 *************************************************************************** -- Problem Area 9 -- When cloud ceiling heights are computed based on prognostic cloud mixing ratios (ICLOUD=3), prognostic data is of MM5 type (IPROG>5), prognostic timestep is longer than CALMET timestep (ISTEPPG>1), and MM5 cloudy skies are replaced by clear skies at the next MM5 timestep, the ceiling heights at intermediate hours were underestimated. *************************************************************************** ---------------------------------- Subroutine RDMM5 (Problem Area 9) ---------------------------------- Check values of prognostic ceiling heights for clear sky conditions at both t1 and t2 before performing a time interpolation, BEFORE change: -------------- do 137 i=1,nx c --- linearly interpolate ceiling height if non zero values at t1 and t2 c --- otherwise pick the non zero value if ( iceilg1(i,j).ne.0 .and. iceilg1(i,j).ne.0 )then AFTER change ------------- do 137 i=1,nx c --- linearly interpolate ceiling height if non zero values at t1 and t2 c --- otherwise pick the non zero value if ( iceilg1(i,j).ne.0 .and. iceilg2(i,j).ne.0 )then *********************************************************************** -- Problem Area 10 -- Prognostic ceiling heights (ICLOUD=3) were not correctly interpolated to the CALMET gridpoints *************************************************************************** ---------------------------------- Subroutine RDMM5 (Problem Area 10) ---------------------------------- Update CALMET gridpoint coordinates (x,y) before interpolating the prognostic ceiling heights onto the CALMET grid BEFORE change: -------------- do 127 j=1,ny do 127 i=1,nx AFTER change ------------- do 127 j=1,ny c --- Compute Y of cell center y = yorigcc + (j - 1) * dyk do 127 i=1,nx c --- Compute X of cell center x = xorigcc + (i - 1) * dxk *************************************************************** -- Problem Area 11 -- No vertical extrapolation of temperature from lowest 3D.DAT level to lower CALMEt levels was performed overwater when ITWPROG=2 (constant T profile was assumed below lowest 3D.DAT level instead of the intented interpolation between SST and lowest 3D.DAT level temperature) *************************************************************************** ---------------------------------- Subroutine RDMM5 (Problem Area 11) ---------------------------------- Apply temperature interpolation overwater at CALMEt levels below lowest M3D level to actual temperature array BEFORE change: -------------- else c daytime: assume "dry" adiab (constant Tetha virt) c except overwater if SSTP is available (060322-frr) if (itwprog.eq.2.and.ilu4(i,j).eq.iluoc3d) then zx=cellzc(k)/zp(i,j,index) ptmm4=sstp2(i,j)*(1.-zx) : +pt(index) *zx else vptdat(i,j,k)=vptdat(i,j,index) TMM4(I,J,K)= tp(i,j,index)+ : 0.0098*(zp(i,j,index)-cellzc(k)) endif endif AFTER change ------------- else c daytime: assume constant Tetha virt. and "dry" adiabat c except overwater if SSTP is available (060315-frr) if (itwprog.eq.2.and.ilu4(i,j).eq.iluoc3d) then c Interpolate Temp. between SST and lowest M3D level (080512) zx=cellzc(k)/zp(i,j,index) tmm4(i,j,k)=sstp2(i,j)*(1.-zx) : +tp(i,j,index) *zx else TMM4(I,J,K)= tp(i,j,index)+ : 0.0098*(zp(i,j,index)-cellzc(k)) endif vptdat(i,j,k)=vptdat(i,j,index) endif *********************************************************************** ----------------------------------------------------------------------- 3. COORDLIB Changes (Version 5.82, Level 071109) ----------------------------------------------------------------------- *********************************************************************** --- Problem Area 1 --- A non-zero False Northing appropriate for UTM-S was used when converting S. hemisphere locations to UTM-N coordinates. This causes the coordinates returned to be in UTM-S. *********************************************************************** ---------------------------------- Subroutine PJINIT (Problem Area 1) ---------------------------------- BEFORE change: -------------- IF (DATA(2) .LT. ZERO) BUFFL(8) = 10000000.0D0 AFTER change: -------------- c --- COORDS c --- Use just the ZONE provided when setting the false Northing c IF (DATA(2) .LT. ZERO) BUFFL(8) = 10000000.0D0 ---------------------------------- Subroutine COORDS (Problem Area 1) ---------------------------------- BEFORE change: (a) -------------- C --- STEP 2 PROJECTION FROM LAT-LON CALL GTPZ0(CRDIN,0,INZONE,TPARIN,INUNIT,IOSPH,IPR,JPR, . LEMSG,LPARM,CRDIO,IOSYS,IOZONE,TPARIO,IOUNIT,LN27,LN83, . FN27,FN83,LENGTH,IFLG) CALL ERRPRT(IFLG,LPR,1) C C --- SPECIAL FIX FOR NH CROSS OVER OF ZONE [IOZONE > 0 crdin(2) <0.0] IF(INSYS.EQ.0.AND.IOSYS.EQ.1.AND.IOZONE.GT.0.AND.CRDIN(2). 1 LT.0.0)THEN CRDIO(2) = CRDIO(2)-10000000.0 ENDIF AFTER change: (a) -------------- C --- STEP 2 PROJECTION FROM LAT-LON CALL GTPZ0(CRDIN,0,INZONE,TPARIN,INUNIT,IOSPH,IPR,JPR, . LEMSG,LPARM,CRDIO,IOSYS,IOZONE,TPARIO,IOUNIT,LN27,LN83, . FN27,FN83,LENGTH,IFLG) CALL ERRPRT(IFLG,LPR,1) C c --- Fix moved into PJINIT (070921) cC --- SPECIAL FIX FOR NH CROSS OVER OF ZONE [IOZONE > 0 crdin(2) <0.0] c IF(INSYS.EQ.0.AND.IOSYS.EQ.1.AND.IOZONE.GT.0.AND.CRDIN(2). c 1 LT.0.0)THEN c CRDIO(2) = CRDIO(2)-10000000.0 c ENDIF BEFORE change: (b) -------------- C --- REGULAR CASES IF(INSYS.NE.IOSYS)THEN CALL GTPZ0(CRDIN,INSYS,INZONE,TPARIN,INUNIT,INSPH,IPR,JPR, . LEMSG,LPARM,CRDIO,IOSYS,IOZONE,TPARIO,IOUNIT,LN27,LN83, . FN27,FN83,LENGTH,IFLG) CALL ERRPRT(IFLG,LPR,1) C C --- SPECIAL FIX FOR NH CROSS OVER OF ZONE [IOZONE > 0 crdin(2) <0.0] IF(INSYS.EQ.0.AND.IOSYS.EQ.1.AND.IOZONE.GT.0.AND.CRDIN(2). 1 LT.0.0)THEN CRDIO(2) = CRDIO(2)-10000000.0 ENDIF AFTER change: (b) -------------- C --- REGULAR CASES IF(INSYS.NE.IOSYS)THEN CALL GTPZ0(CRDIN,INSYS,INZONE,TPARIN,INUNIT,INSPH,IPR,JPR, . LEMSG,LPARM,CRDIO,IOSYS,IOZONE,TPARIO,IOUNIT,LN27,LN83, . FN27,FN83,LENGTH,IFLG) CALL ERRPRT(IFLG,LPR,1) C c --- Fix moved into PJINIT (070921) cC --- SPECIAL FIX FOR NH CROSS OVER OF ZONE [IOZONE > 0 crdin(2) <0.0] c IF(INSYS.EQ.0.AND.IOSYS.EQ.1.AND.IOZONE.GT.0.AND.CRDIN(2). c 1 LT.0.0)THEN c CRDIO(2) = CRDIO(2)-10000000.0 c ENDIF *********************************************************************** --- Problem Area 2 --- Three work-array variables and one output variable were not completely initialized, causing the program to halt when extra compiler-checking options are used. *********************************************************************** ---------------------------------- Subroutine COORDS (Problem Area 2) ---------------------------------- BEFORE change: -------------- C C --- Now converts the TPARIN, TPARIO FROM DD to DDDMMMSSS.SS C --- UTM's AFTER change: -------------- c --- Initialize full work arrays do k = 1,15 dwrk(k) = 0.0D+00 dwrk2(k) = 0.0D+00 tdum(k) = 0.0D+00 enddo c --- Initialize output variable UTMOUT utmout = 0.0D+00 C C --- Now converts the TPARIN, TPARIO FROM DD to DDDMMMSSS.SS C --- UTM's ----------------------------------------------------------------------- 4. CALUTILS Changes (Version 2.56, Level 080407) ----------------------------------------------------------------------- *********************************************************************** --- Problem Area 1 -- Exponential notation processing in ALTONU did not properly interpret an entry without a decimal point. *********************************************************************** ---------------------------------- Subroutine ALTONU (Problem Area 1) ---------------------------------- Treat case in which exponential notation is used without a decimal point. Trap case where no number appears in front the E or D in exponential notation. BEFORE change: -------------- 111 continue c c --- Convert integer numerics to real number rno=0.0 AFTER change: -------------- 111 continue c --- 080407 Update: c --- Correct for missing decimal point before decoding if(idec.EQ.0) idec=istop+1 c --- Trap missing number in front of E,D if(istop.LT.1 .OR. istart.GT.istop) then write(ioout,120)(alp(n),n=1,ncar) write(*,*) write(*,*)'Missing number!' stop 'Halted in ALTONU -- see list file.' endif c c --- Convert integer numerics to real number rno=0.0