----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : F Dated : October 25, 2010 ----------------------------------------------------------------------- The following code modifications address problems identified in the Version 5 CALPUFF modeling system currently approved for regulatory applications. Bulletin MCB-F(101025) outlines these changes and together with those in MCB-A(040716), MCB-B(051216), MCB-C(060804), MCB-D(070622), and MCB-E(080613) 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.57, Level 090202) -------------------------------------------- Many subroutines used throughout the CALPUFF System are packaged as the FORTRAN include-file: CALUTILS.FOR This MCB includes changes to CALUTILS components that have no influence on model results, but enhance the use of the system. --- Increase control file line length to 200 characters Modified: PARAMS.CAL, READIN --- Activate CPU clock using F95 system routine Modified: DATETM CALPUFF (to Version 5.835, Level 101025) ---------------------------------------- -- Problem Area 1 -- Relax test for small negative travel in SIGTZ following the approach used in SIGTY. Small negative travel should be interpreted as zero, but larger negative travel indicates a potential problem and the code should halt as it does now. Configuration for both SIGTY and SIGTZ: 0.00s >dt> -.01s dt reset to zero (incident not reported) -.01s >dt> -1.0s dt reset to zero (incident counted) -1.0s >dt CALPUFF Halted with error message 0.00m >dx> -.01m dx reset to zero (incident not reported) -.01m >dx> -1.0m dx reset to zero (incident counted) -1.0m >dx CALPUFF Halted with error message Incidents counted are reported to the list file at the end of the run with the most negative travel encountered. (No impact on results) Modified: FIN, SIGTY, SIGTZ, WARN -- Problem Area 2 -- Align specific details of PRIME implementation with those used in ISC-PRIME to more closely match that model when using the same meteorology and dispersion. a. Fix cavity concentrations (missing reflection from the ground, and impact from the end of the cavity to the end of the cavity transition was not included). Concentrations further downwind are not affected. b. Computing sigmas in building wake: Include BID in sigmas at the entry to the wake (xi) and allow these enhanced sigmas to grow within wake PDF region. Tabulated sigmas for PRIME source will now include the BID so that BID in the puff arrays must be set to zero. Previous method grew sigmas in the wake without BID, and tabulated sigmas from the source to the end of the wake without BID. Non-zero BID was then added whenever a sigma was extracted from the table. c. Compute receptor-specific sigmas downwind of the end of the PRIME wake (e.g., beyond the end of the wake table) using time or distance measured from this point rather than interpolating sigmas via virtual time (distance) at the start and end of the sampling step. This assures that the sigmas are continuous in this region. This applies to puff steps that start within the wake. Subsequent steps interpolate the virtuals. d. Align calculation of plume drdx in wake region with ISC-PRIME. The drdx includes extra dilution in the wake and this is then used in the numerical rise, so the resolution of tabulated values can alter rise and mass captured in the cavity. Implement by creating XWAKISC() and DRWAKISC() arrays. Since these are only used in NUMRISE for the current source, they are not saved to the SRCTAB file. e. Trap negative virtuals when computing source sigma that grows to wake-modified sigma at end of wake, and set to zero f. Use ISC-PRIME method of computing distances for numerical rise and wake arrays. g. Disable comparison of the wake-influenced sigma growth with the current ambient sigma growth (taking the larger of the 2) when filling the wake arrays (now conforms to the treatment used in ISC-PRIME). h. Remove approximate calculation of the BID-enhanced sigmas that was used when the BID term is small compared to the sigma. New: WAKE_ADDBID Modified: WAKEDAT.PUF CAV_CONC, CAV_SAMP, NUMPR1, NUMRISE, POINTS1, POINTS2, PUFRECS, SETPUF, SETSLG WAKE_DFSN, WAKE_DRDX, WAKE_FIN, WAKE_TAB, WAKE_XSIG -- Problem Area 3 -- Rename station ID returned from FINDR and FINDI (used when data are missing at nearest station) so that it does not replace the station ID obtained from the NEARS array. Data assigned to 2D or 3D arrays in CALMET are not affected. Modified: EXMET -- Problem Area 4 -- Roughness adjustment to PG sigmas is computed for 1 roughness length and should be used with non-gridded meteorological data, so restrict MROUGH = 1 to METFM = 2, 3, 4, or 5 to avoid unintended application with gridded data. Modified: QAINP -- Problem Area 5 -- Calculate mid-pt of sampling step using half the step rather than as the average of the 2 end-pts to improve precision. Modified: COMP -- Problem Area 6 -- Add concentration calculations at elevated receptors that are located above the mixing height when there is mass above. The impact at these receptors had been set to zero. Also, restrict deposition fluxes to receptors that are located ON the ground (currently only discrete receptors may be elevated). Modified: AREAINT, CALCBC, CALCPF, CALCSL, CAV_CONC, CAV_SAMP, COMP, CPQROMB, CPTRAP, LSSLINT, PLMFOG, POINTS1, POINTS2, ROMBTIN, SLUGAVE, SLUGSNP, SLUGINT, TRAPSL, VCBAR, VCOUP -- Problem Area 7 -- Trap case in which the shear adjustment to Briggs plume rise is applied in near-calm conditions. A distance to final rise of zero (calm) can produce an incorrect shear-modified rise that is zero. Also increase the number of digits used to write TSTAK and DTEMP to the list file (debug output). Contributed by Brian Zelt Modified: POINTS1, POINTS2 -- Problem Area 8 -- Assign power-law exponent to default value for met file types that do not use power-law profiles (CALMET, CTDM PROFILE, and AERMET PROFILE) so that exponent is available for Briggs plume rise wind shear modification option (MSHEAR). Modified: BLOCK DATA, ADVECT, AREAS1, AREAS2, LINES1, LINES2, NUMMET, POINTS1, POINTS2, POWLAW, PRFINSH, PRSHEAR, SETLINE, VOLS, WINDSET -- Problem Area 9 -- Add call to SRCTABOUT to update the source number of the current puff in the DA file when line-source slugs are processed during the first step. The source number given to slugs generated from line segments 2+ has '1000' added to it, and this '1000' is removed in the first model step. The model had halted in the second timestep. This only affects buoyant line sources when modeled with slugs. The DA file of source-related tabulations was introduced in MCB-E. Modified: COMP -- Problem Area 10 -- Use local source-table arrays in the subroutine used to move the puff index when inactive puffs are removed from DA file of source-related tabulations. The current contents in the arrays had been overwritten before being used because the /SRCTAB/ common had been included, and the roll-down process happens after the tables are created. This affects tabulated rise information whenever inactive puffs are removed when generating puffs/slugs. Roll-downs made when splitting puffs or when generating a restart file do not overwrite arrays that are currently in use. The DA file of source-related tabulations was introduced in MCB-E. Modified: SWAPTAB -- Problem Area 11 -- Do not allow the receptor-specific sigmas for a line source to be zero. Add a check for sy0 and sz0 when puffs/slugs are generated to assure that SYMIN and SZMIN are the smallest values used (avoids halt due to divide-by-zero). Modified: LINES1, LINES2 -- Problem Area 12 -- Add check for zero final rise from a buoyant line source when computing receptor-specific sigmas for the case of an emitting slug (attached to source). The zero rise can result in a divide-by-zero that halts the model. Modified: RECSPEC0 -- Problem Area 13 -- A number of variables were not properly defined or initialized. Results are not affected. Modified: COMP, FOGREC, LINES1, MET2, MET3, MET4, POINTS1, POINTS2, QAPLOT1, RDEMBC, RDHDBC, RDHDBC2, RDMET4, RDMET5, SLGRECS, SLUGINT, TFERCF, WAKE_DBG ----------------------------------------------------------------------- ----------------------------------------------------------------------- 1. CALPUFF Changes (Version 5.831, Level 080520) ----------------------------------------------------------------------- *********************************************************************** -- Problem Area 1 -- Relax test for small negative travel in SIGTZ following the approach used in SIGTY. Small negative travel should be interpreted as zero, but larger negative travel indicates a potential problem and the code should halt as it does now. *********************************************************************** ------------------------------------ Subroutine FIN (Problem Area 1) ------------------------------------ BEFORE change: -------------- c --- Process any warnings collected during the run call WARN('FIN ',0.) AFTER change: -------------- c --- Process any warnings collected during the run call WARN('FIN ',0.,0.) ------------------------------------ Subroutine SIGTY (Problem Area 1) ------------------------------------ BEFORE change: (a) -------------- elseif(dt.ge.zerosec(2)) then call WARN('SIGTY-t ',dt) AFTER change: (a) -------------- elseif(dt.ge.zerosec(2)) then call WARN('SIGTY-t ',dt,zerosec(1)) BEFORE change: (b) -------------- elseif(dxkm.ge.zerokm(2)) then call WARN('SIGTY-x ',dxkm) AFTER change: (b) -------------- elseif(dxkm.ge.zerokm(2)) then call WARN('SIGTY-x ',dxkm,zerokm(1)) ------------------------------------ Subroutine SIGTZ (Problem Area 1) ------------------------------------ BEFORE change: (a) -------------- c data crit/1.0e-4/,zero/0.0/ c c --- Set values to test for negative virtuals that are treated as zero data zerosec/-.01/, zerokm/-0.00001/ AFTER change: (a) -------------- c --- Declare variables to test negative virtuals c --- (threshold for warning, threshold for HALT) real zerosec(2), zerokm(2) c data crit/1.0e-4/,zero/0.0/ c c --- Set values to test for negative virtuals that are treated as zero c --- (threshold for warning, threshold for HALT) data zerosec/-.01,-1.0/, zerokm/-0.00001,-.001/ BEFORE change: (b) -------------- c --- Make sure travel time is positive. if(dt.lt.zero) then if(dt.ge.zerosec) then dt=zero else write(io6,*)'SR. SIGTZ: FATAL ERROR: dt = ',dt write(*,*) stop 'Halted in SIGTZ -- see list file.' endif endif AFTER change: (b) -------------- c --- Make sure travel time is positive. if(dt.lt.zero) then if(dt.ge.zerosec(1)) then dt=zero elseif(dt.ge.zerosec(2)) then call WARN('SIGTZ-t ',dt,zerosec(1)) dt=zero else write(io6,*)'ERROR in SUBR. SIGTZ -- Incremental travel ', 1 'time < 0 -- DT = ',dt write(*,*) stop 'Halted in SIGTZ -- see list file.' endif endif BEFORE change: (c) -------------- c --- Make sure downwind distance is positive. if(dxkm.lt.zero) then if(dxkm.ge.zerokm) then dxkm=zero else write(io6,*)'SR. SIGTZ: FATAL ERROR: dxkm = ',dxkm write(*,*) stop 'Halted in SIGTZ -- see list file.' endif endif AFTER change: (c) -------------- c --- Make sure downwind distance is positive. if(dxkm.lt.zero) then if(dxkm.ge.zerokm(1)) then dxkm=zero elseif(dxkm.ge.zerokm(2)) then call WARN('SIGTZ-x ',dxkm,zerokm(1)) dxkm=zero else write(io6,*)'ERROR in SUBR. SIGTZ -- Incremental travel ', 1 'distance < 0 -- DXKM = ',dxkm write(*,*) stop 'Halted in SIGTZ -- see list file.' endif endif ------------------------------------ Subroutine WARN (Problem Area 1) ------------------------------------ BEFORE change: (a) -------------- subroutine warn(astring,value) AFTER change: (a) -------------- subroutine warn(astring,value,zeroval) BEFORE change: (b) -------------- c VALUE - real - value to process AFTER change: (b) -------------- c VALUE - real - value to process c ZEROVAL - real - threshold travel time (sec) or distance (km) c for counting (calling this routine) BEFORE change: (c) -------------- logical LWARN character*12 astring c --- Initialize via data statements data lwarn/.FALSE./ data sigtysec/0.0/, ntsigty/0/ data sigtykm/0.0/, nxsigty/0/ AFTER change: (c) -------------- logical LWARNY,LWARNZ character*12 astring c --- Initialize via data statements data lwarny/.FALSE./,lwarnz/.FALSE./ data sigtysec/0.0/, ntsigty/0/ data sigtykm/0.0/, nxsigty/0/ data sigtzsec/0.0/, ntsigtz/0/ data sigtzkm/0.0/, nxsigtz/0/ BEFORE change: (d) -------------- if(astring.EQ.'SIGTY-t ') then c --- Process negative delta-t from SIGTY if(value.LT.sigtysec) sigtysec=value ntsigty=ntsigty+1 ntsigty=MIN(ntsigty,1000001) lwarn=.TRUE. elseif(astring.EQ.'SIGTY-x ') then c --- Process negative delta-x from SIGTY if(value.LT.sigtykm) sigtykm=value nxsigty=nxsigty+1 nxsigty=MIN(nxsigty,1000001) lwarn=.TRUE. AFTER change: (d) -------------- if(astring.EQ.'SIGTY-t ') then c --- Process negative delta-t from SIGTY if(value.LT.sigtysec) sigtysec=value ntsigty=ntsigty+1 ntsigty=MIN(ntsigty,1000001) lwarny=.TRUE. zerosecy=zeroval elseif(astring.EQ.'SIGTY-x ') then c --- Process negative delta-x from SIGTY if(value.LT.sigtykm) sigtykm=value nxsigty=nxsigty+1 nxsigty=MIN(nxsigty,1000001) lwarny=.TRUE. zeromy=zeroval*1000. elseif(astring.EQ.'SIGTZ-t ') then c --- Process negative delta-t from SIGTZ if(value.LT.sigtzsec) sigtzsec=value ntsigtz=ntsigtz+1 ntsigtz=MIN(ntsigtz,1000001) lwarnz=.TRUE. zerosecz=zeroval elseif(astring.EQ.'SIGTZ-x ') then c --- Process negative delta-x from SIGTZ if(value.LT.sigtzkm) sigtzkm=value nxsigtz=nxsigtz+1 nxsigtz=MIN(nxsigtz,1000001) lwarnz=.TRUE. zeromz=zeroval*1000. BEFORE change: (e) -------------- c --- Report diagnostic if(LWARN) then write(io6,*) write(io6,*) write(io6,*) write(io6,*)'Diagnostic for SIGTY:' write(io6,*)' Minimum dt(sec) = ',sigtysec if(ntsigty.GT.100000) then write(io6,*)' Number less than -.01s > ',ntsigty else write(io6,*)' Number less than -.01s = ',ntsigty endif write(io6,*)' Minimum dx(km) = ',sigtykm if(nxsigty.GT.100000) then write(io6,*)' Number less than -.01m > ',nxsigty else write(io6,*)' Number less than -.01m = ',nxsigty endif write(io6,*) write(io6,*) endif AFTER change: (e) -------------- if(LWARNY) then write(io6,*) write(io6,*) write(io6,*) write(io6,*)'Diagnostic for SIGTY:' write(io6,'(a25,f8.3)')' Most Negative dt(sec) = ', & sigtysec if(ntsigty.GT.100000) then write(io6,3)' Number less than ',zerosecy,'s > ', & ntsigty else write(io6,3)' Number less than ',zerosecy,'s = ', & ntsigty endif write(io6,'(a25,f8.3)')' Most Negative dx(m) = ', & (1000.*sigtykm) if(nxsigty.GT.100000) then write(io6,3)' Number less than ',zeromy,'m > ', & nxsigty else write(io6,3)' Number less than ',zeromy,'m = ', & nxsigty endif write(io6,*) write(io6,*) endif if(LWARNZ) then write(io6,*) write(io6,*) write(io6,*) write(io6,*)'Diagnostic for SIGTZ:' write(io6,'(a25,f8.3)')' Most Negative dt(sec) = ', & sigtzsec if(ntsigtz.GT.100000) then write(io6,3)' Number less than ',zerosecz,'s > ', & ntsigtz else write(io6,3)' Number less than ',zerosecz,'s = ', & ntsigtz endif write(io6,'(a25,f8.3)')' Most Negative dx(m) = ', & (1000.*sigtzkm) if(nxsigtz.GT.100000) then write(io6,3)' Number less than ',zeromz,'m > ', & nxsigtz else write(io6,3)' Number less than ',zeromz,'m = ', & nxsigtz endif write(io6,*) write(io6,*) endif BEFORE change: (f) -------------- return end AFTER change: (f) -------------- return 3 format(a19,f8.3,a4,i20) end *********************************************************************** --- Problem Area 2 --- Align specific details of PRIME implementation with those used in ISC-PRIME to more closely match that model when using the same meteorology and dispersion. *********************************************************************** ---------------------------------------- Subroutine CAV_CONC (Problem Areas 2, 6) ---------------------------------------- BEFORE change: (a) -------------- subroutine cav_conc(ldb,nspec,q,qcav,xsgrd,ysgrd,ntr,xtr,ztr, & flow,dpbl,xcgrd,ycgrd,x115,x115c) AFTER change: (a) -------------- subroutine cav_conc(ldb,nspec,q,qcav,xsgrd,ysgrd,ntr,xtr,ztr,htr, & flow,dpbl,hlidmax,xcgrd,ycgrd,x115,x115c) BEFORE change: (b) -------------- c ZTR(mxrise) - real - Height above ground (m) in rise table c FLOW - real - Flow vector (radians CW from N) c DPBL - real - Boundary layer depth (m) AFTER change: (b) -------------- c ZTR(mxrise) - real - Height above ground (m) in rise table c HTR(mxrise) - real - Rise (m) without streamline effects c FLOW - real - Flow vector (radians CW from N) c DPBL - real - Boundary layer depth (m) c HLIDMAX - real - Maximum lid height (m) for mass above HLID BEFORE change: (c) -------------- real xtr(mxrise),ztr(mxrise) AFTER change: (c) -------------- real xtr(mxrise),ztr(mxrise),htr(mxrise) BEFORE change: (d) -------------- c --- Get coupling call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,dpbl,cbyq0,cbyqc) AFTER change: (d) -------------- c --- Get coupling call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,htr,dpbl,hlidmax,cbyq0,cbyqc) BEFORE change: (e) -------------- c --- Get coupling call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,dpbl,cbyq0,cbyqc) AFTER change: (e) -------------- c --- Get coupling call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,htr,dpbl,hlidmax,cbyq0,cbyqc) BEFORE change: (f) -------------- c --- Get coupling call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,dpbl,cbyq0,cbyqc) AFTER change: (f) -------------- c --- Get coupling call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,htr,dpbl,hlidmax,cbyq0,cbyqc) ---------------------------------------- Subroutine CAV_SAMP (Problem Areas 2, 6) ---------------------------------------- BEFORE change: (a) -------------- subroutine cav_samp(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,hlid,cbyq0,cbyqc) AFTER change: (a) -------------- subroutine cav_samp(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,htr,hlid,hlidmax,cbyq0,cbyqc) BEFORE change: (b) -------------- c ZTR(mxrise) - real - Height above ground (m) in rise table c HLID - real - Boundary layer depth (m) AFTER change: (b) -------------- c ZTR(mxrise) - real - Height above ground (m) in rise table c HTR(mxrise) - real - Rise table (m) without streamline effects c HLID - real - Boundary layer depth (m) c HLIDMAX - real - Maximum lid height (m) for mass above HLID BEFORE change: (c) -------------- real xtr(mxrise),ztr(mxrise) logical ldb, nobidprm, nobidcav AFTER change: (c) -------------- real xtr(mxrise),ztr(mxrise),htr(mxrise) logical ldb, nobidprm BEFORE change: (d) -------------- c --- Add BID to sigmas for the primary source data nobidprm/.FALSE./ data nobidcav/.TRUE./ AFTER change: (d) -------------- c --- Add BID to sigmas for the primary source c --- (if not already in tabulated sigmas) data nobidprm/.FALSE./ BEFORE change: (e) -------------- c --- Set sigmas for primary (sy,sz) and cavity (syc,szc) sources call WAKE_XSIG(xr,timep,zero,nobidprm,sz,sy,szcdum,sycdum) call WAKE_XSIG(xr,timec,zero,nobidcav,szdum,sydum,szc,syc) AFTER change: (e) -------------- c --- Set primary source height and rise call INTERTAB(xr,ntr,xtr,i1,i2,factor) ht=ztr(i1)+factor*(ztr(i2)-ztr(i1)) rise=htr(i1)+factor*(htr(i2)-htr(i1)) c --- Set sigmas for primary (sy,sz) and cavity (syc,szc) sources call WAKE_XSIG(xr,timep,rise,nobidprm,sz,sy,szcdum,sycdum) call WAKE_XSIG(xr,timec,zero,nobidprm,szdum,sydum,szc,syc) BEFORE change: (f) -------------- 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 AFTER change: (f) -------------- c --- Primary source contributes outside of cavity fy0=(EXP(-0.5*(yr/sy)**2))/(rt2pi*sy) cbyq0=fy0*VCOUP(icode,zr,ht,sz,hlid,hlidmax)/Urh BEFORE change: (g) -------------- c --- Compute concentration due to cavity c --- Initialize receptor ht to value passed into subroutine zrc=zr if(ipositn.EQ.3) then c --- Sample vertical distribution when outside cavity, and c --- turn off the 'inside' cavity concentration cavnear=zero AFTER change: (g) -------------- c --- Compute concentration due to cavity material c --- Initialize receptor ht to value passed into subroutine zrc=zr if(ipositn.EQ.3) then c --- Sample vertical distribution when outside cavity BEFORE change: (h) -------------- else zrc=zero c --- 'Inside' cavity cbyq before transition region fycn=(EXP(-0.5*((yrb-ycb)/sycav(1))**2))/(rt2pi*sycav(1)) cavnear=fycn/(rt2pi*szcav(1)*Ub) endif c --- 'Outside' cavity cbyq fycf=(EXP(-0.5*((yrb-ycb)/syc)**2))/(rt2pi*syc) cavfar=fycf*VCOUP(icode,zrc,zero,szc,hlid)/Ub AFTER change: (h) -------------- else c --- Position 1 is within structure, and 2 is within cavity c --- Place receptor on the 'ground' to get inside cavity conc. zrc=zero endif c --- 'Inside' cavity cbyq before transition factor fycn=(EXP(-0.5*((yrb-ycb)/sycav(1))**2))/(rt2pi*sycav(1)) cavnear=fycn*VCOUP(icode,zrc,zero,szcav(1),hlid,hlidmax)/Ub c --- 'Outside' cavity cbyq before transition factor fycf=(EXP(-0.5*((yrb-ycb)/syc)**2))/(rt2pi*syc) cavfar=fycf*VCOUP(icode,zrc,zero,szc,hlid,hlidmax)/Ub ---------------------------------------- Subroutine NUMPR1 (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- nstep=1+(slast/(ds0*FLOAT(mxnw))) AFTER change: (a) -------------- nstep=1+((slast-0.5*ds0)/(ds0*FLOAT(mxnw))) BEFORE change: (b) -------------- c --- Define the meteorological grid c --- Set cell face heights at every 10 m from 10-200 m dz=10. nn=1 zfacea(nn)=0.0 do i=1,20 nn=nn+1 zfacea(nn)=zfacea(nn-1)+dz enddo c --- Set cell face heights every 50 m from 250-500 m dz=50. do i=21,26 nn=nn+1 zfacea(nn)=zfacea(nn-1)+dz enddo c --- Set cell face heights every 100 m from 600-2000 m dz=100. do i=27,41 nn=nn+1 zfacea(nn)=zfacea(nn-1)+dz enddo c --- Set cell face heights every 500 m from 2500-4000 m dz=500. do i=42,45 nn=nn+1 zfacea(nn)=zfacea(nn-1)+dz enddo if(nn.NE.nzap1) stop 'NUMPR1: nn.NE.nzap1' c --- Compute grid point heights do i=1,nza zgpta(i)=0.5*(zfacea(i+1)+zfacea(i)) enddo AFTER change: (b) -------------- c --- Define the meteorological grid c --- Set grid points every 10 m from 10-200 m dz=10. nn=1 zgpta(nn)=dz do i=2,20 nn=nn+1 zgpta(nn)=zgpta(nn-1)+dz enddo c --- Set grid points every 50 m from 250-500 m dz=50. do i=21,26 nn=nn+1 zgpta(nn)=zgpta(nn-1)+dz enddo c --- Set grid points every 100 m from 600-2000 m dz=100. do i=27,41 nn=nn+1 zgpta(nn)=zgpta(nn-1)+dz enddo c --- Set grid points every 500 m from 2500-4000 m dz=500. do i=42,45 nn=nn+1 zgpta(nn)=zgpta(nn-1)+dz enddo c --- Compute the cell face heights from the grid point values zfacea(1)=0.0 do i=2,nza zfacea(i)=0.5*(zgpta(i)+zgpta(i-1)) enddo zfacea(nzap1)=zgpta(nza)+0.5*(zgpta(nza)-zgpta(nza-1)) ---------------------------------------- Subroutine NUMRISE (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- c --- NUMRISE calls: ZMET,LUMP,RATE,MARCHING,UNLUMP c ZSTREAM, POSITION, WAKE_U c WAKE_DRDX, WAKE_DFSN, WAKE_FQC AFTER change: (a) -------------- c --- NUMRISE calls: ZMET,LUMP,RATE,MARCHING,UNLUMP c ZSTREAM, POSITION, WAKE_U c WAKE_DRDX, WAKE_DFSN, WAKE_FQC, WAKE_ADDBID BEFORE change: (b) -------------- c --- Restore step size (may have changed) ds=ds0 if(mprm.EQ.1) then c --- PRIME module c --- Determine maximum fraction of plume captured in cavity AFTER change: (b) -------------- c --- Restore step size (may have changed) ds=ds0 if(mprm.EQ.1) then c --- PRIME module c --- Use final form of tabulated rise array to add BID to c --- sigmas in the wake table upwind of xi call WAKE_ADDBID(ldbhr,ntr,xtr,htr) c --- Determine maximum fraction of plume captured in cavity ---------------------------------------- Subroutine POINTS1 (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- c --- Pass final wake structure to source arrays (without BID) call WAKE_FIN(ldbhr,nwk,xwk,szwk,sywk,drwk,szw0,syw0, & ncv,xcv,szcv,sycv,szc0,syc0,fqcvpt1(i)) AFTER change: (a) -------------- c --- Pass final wake structure to source arrays (with BID) call WAKE_FIN(ldbhr,nwk,xwk,szwk,sywk,drwk,szw0,syw0, & ncv,xcv,szcv,sycv,szc0,syc0,fqcvpt1(i)) BEFORE change: (b) -------------- c --- Reset BID and tip downwash tipdw=0.0 bidsq=(rise/3.5)**2 AFTER change: (b) -------------- c --- Reset BID and tip downwash tipdw=0.0 c --- The BID has been moved into the wake tables bidsq=0.0 ---------------------------------------- Subroutine POINTS2 (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- c --- Pass final wake structure to source arrays without BID call WAKE_FIN(ldbhr,nwk,xwk,szwk,sywk,drwk,szw0,syw0, & ncv,xcv,szcv,sycv,szc0,syc0,fqcvpt2(i)) AFTER change: (a) -------------- c --- Pass final wake structure to source arrays with BID call WAKE_FIN(ldbhr,nwk,xwk,szwk,sywk,drwk,szw0,syw0, & ncv,xcv,szcv,sycv,szc0,syc0,fqcvpt2(i)) BEFORE change: (b) -------------- c --- Reset BID and tip downwash tipdw=0.0 bidsq=(rise/3.5)**2 AFTER change: (b) -------------- c --- Reset BID and tip downwash tipdw=0.0 c --- The BID has been moved into the wake tables bidsq=0.0 ---------------------------------------- Subroutine PUFRECS (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- c ------------------------------------------------------------------- c --- Special case of PRIME building downwash --- c ------------------------------------------------------------------- c Plume rise is tabulated, and the sigmas are tabulated from c the source to the point where the turbulence reaches ambient. c Use WAKE_TAB to interpolate rise and sigmas within the c tabulated range of the sigmas, using the appropriate source c arrays. Beyond the range, do nothing here and allow the c normal processing to take place (idone=0). AFTER change: (a) -------------- c ------------------------------------------------------------------- c --- Special case of PRIME building downwash --- c ------------------------------------------------------------------- c Plume rise and sigmas are tabulated from the source to the c point where the turbulence reaches ambient. c 1. Receptors within the tabulated range of the sigmas use c WAKE_TAB to interpolate rise and sigmas c 2. Puffs starting within the tabulated range grow the sigmas c from the end of the table for receptors further downwind c 3. Puffs and receptors beyond the tabulated range are modeled c using normal (non-PRIME) processing later BEFORE change: (b) -------------- c --- Finish up and skip to terrain adjustments if sigmas c --- are extracted if(idone.EQ.1) then rfacsq=0.0 if(rise.GT.0.0) then c --- Add BID contribution to sigmas sigbidsq=(rise/3.5)**2 syr=sqrt(syr**2+sigbidsq) szr=sqrt(szr**2+sigbidsq) c --- Compute the corresponding BIDSQ adjustment factor rfacsq=sigbidsq/bidsq endif go to 500 endif AFTER change: (b) -------------- c --- Finish up and skip to terrain adjustments if sigmas c --- are extracted from wake table or grown from wake table if(idone.EQ.1) then rfacsq=0.0 if(rise.GT.0.0 .AND. bidsq.GT.0.0) then c --- Add BID contribution to sigmas sigbidsq=(rise/3.5)**2 syr=sqrt(syr**2+sigbidsq) szr=sqrt(szr**2+sigbidsq) c --- Compute the corresponding BIDSQ adjustment factor rfacsq=sigbidsq/bidsq endif go to 500 elseif(xttb1.LE.xlast) then c --- Puff starts in tabulated wake region and receptor is c --- beyond the wake so grow sigmas from end of table c --- Properties at XLAST: call WAKE_TAB(istype,isnum,idw0(ipnum),xlast, & sylast,szlast,hgr,rise,xlast2,idone) dxwakkm=(xtt1-xlast)*.001 dtwaksec=(xtt1-xlast)/uavg c --- Set selected data in /CSIGMA/ for sigma calls call SETCSIG(idopty,idoptz,iru,uavg,istab,el,bvf, & sigv,sigw,symin,szmin,ze1,dpbl) c --- Compute the receptor-specific sigmas via 'forward' c --- calls to sigma routines call SIGTZ(szlast,dxwakkm,dtwaksec,zb1,szr,dum1,dum2) call SIGTY(sylast,dxwakkm,dtwaksec,syr,dum1,dum2) c --- Height and rise at receptor call grise(xtt1,hgr,risefac) rfacsq=risefac**2 if(rise.GT.0.0 .AND. bidsq.GT.0.0) then c --- Add BID contribution to sigmas sigbidsq=(rise/3.5)**2 syr=sqrt(syr**2+sigbidsq) szr=sqrt(szr**2+sigbidsq) c --- Compute the corresponding BIDSQ adjustment factor rfacsq=sigbidsq/bidsq endif go to 500 endif ---------------------------------------- Subroutine SETPUF (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- c --- Process puff at start (within tabulated region) call SIGTY(syb1,zero,zero,sydum,vtyb1,vdyb1) call SIGTZ(szb1,zero,zero,zpb(ii),szdum,vtzb1,vdzb1) syb1=AMAX1(syb1,symin) szb1=AMAX1(szb1,szmin) c --- Process midpoint AFTER change: (a) -------------- c --- Process puff at start (within tabulated region) syb1=AMAX1(syb1,symin) szb1=AMAX1(szb1,szmin) call SIGTY(syb1,zero,zero,sydum,vtyb1,vdyb1) call SIGTZ(szb1,zero,zero,zpb(ii),szdum,vtzb1,vdzb1) c --- Process midpoint BEFORE change: (b) -------------- c --- Point within tabulated region call SIGTY(sym1,zero,zero,sydum,vtym1,vdym1) call SIGTZ(szm1,zero,zero,zpb(ii),szdum,vtzm1,vdzm1) AFTER change: (b) -------------- c --- Point within tabulated region sym1=AMAX1(sym1,symin) szm1=AMAX1(szm1,szmin) call SIGTY(sym1,zero,zero,sydum,vtym1,vdym1) call SIGTZ(szm1,zero,zero,zpb(ii),szdum,vtzm1,vdzm1) BEFORE change: (c) -------------- call SIGTZ(zero,dxkm,dts,zpb(ii),szm1,vtzm1,vdzm1) endif sym1=AMAX1(sym1,symin) szm1=AMAX1(szm1,szmin) AFTER change: (c) -------------- call SIGTZ(zero,dxkm,dts,zpb(ii),szm1,vtzm1,vdzm1) sym1=AMAX1(sym1,symin) szm1=AMAX1(szm1,szmin) endif BEFORE change: (d) -------------- c --- Point within tabulated region call SIGTY(sye1,zero,zero,sydum,vtye1,vdye1) call SIGTZ(sze1,zero,zero,zpb(ii),szdum,vtze1,vdze1) AFTER change: (d) -------------- c --- Point within tabulated region sye1=AMAX1(sye1,symin) sze1=AMAX1(sze1,szmin) call SIGTY(sye1,zero,zero,sydum,vtye1,vdye1) call SIGTZ(sze1,zero,zero,zpb(ii),szdum,vtze1,vdze1) BEFORE change: (e) -------------- call SIGTZ(zero,dxkm,dts,zpb(ii),sze1,vtze1,vdze1) endif sye1=AMAX1(sye1,symin) sze1=AMAX1(sze1,szmin) AFTER change: (e) -------------- call SIGTZ(zero,dxkm,dts,zpb(ii),sze1,vtze1,vdze1) sye1=AMAX1(sye1,symin) sze1=AMAX1(sze1,szmin) endif BEFORE change: (f) -------------- c --- Add quadrature contribution to sigmas in /CURRENT/ sysq=bidfnl(ii)+sy0sq if(sysq.GT.zero) then ratio1=sysq/(syb1*syb1) ratio2=sysq/(sym1*sym1) ratio3=sysq/(sye1*sye1) if(ratio1.LE.0.2) then c --- Approximate square root syb1=syb1*(one+half*ratio1) sym1=sym1*(one+half*ratio2) sye1=sye1*(one+half*ratio3) else syb1=syb1*sqrt(one+ratio1) sym1=sym1*sqrt(one+ratio2) sye1=sye1*sqrt(one+ratio3) endif endif AFTER change: (f) -------------- c --- Add quadrature contribution to sigmas in /CURRENT/ sysq=bidfnl(ii)+sy0sq if(sysq.GT.zero) then syb1=SQRT(syb1*syb1+sysq) sym1=SQRT(sym1*sym1+sysq) sye1=SQRT(sye1*sye1+sysq) endif BEFORE change: (g) -------------- c --- Add contribution of buoyancy-enhancement to sigmas in /CURRENT/ if(bidfnl(ii).GT.zero) then ratio1=bidfnl(ii)/(szb1*szb1) ratio2=bidfnl(ii)/(szm1*szm1) ratio3=bidfnl(ii)/(sze1*sze1) if(ratio1.LT.0.2) then c --- Approximate square root szb1=szb1*(one+half*ratio1) szm1=szm1*(one+half*ratio2) sze1=sze1*(one+half*ratio3) else szb1=szb1*sqrt(one+ratio1) szm1=szm1*sqrt(one+ratio2) sze1=sze1*sqrt(one+ratio3) endif endif AFTER change: (g) -------------- c --- Add contribution of buoyancy-enhancement to sigmas in /CURRENT/ if(bidfnl(ii).GT.zero) then szb1=SQRT(szb1*szb1+bidfnl(ii)) szm1=SQRT(szm1*szm1+bidfnl(ii)) sze1=SQRT(sze1*sze1+bidfnl(ii)) endif ---------------------------------------- Subroutine SETSLG (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- c --- Add quadrature contribution to new sigmas sye1=syae1 sze1=szae1 if(sysq.GT.zero .AND. iold.EQ.1) then ratio=sysq/(syae1*syae1) if(ratio.LE.0.2) then c --- Approximate square root sye1=syae1*(one+half*ratio) else sye1=syae1*sqrt(one+ratio) endif endif if(bidfnl(ii).GT.zero) then ratio=bidfnl(ii)/(szae1*szae1) if(ratio.LE.0.2) then c --- Approximate square root sze1=szae1*(one+half*ratio) else sze1=szae1*sqrt(one+ratio) endif endif AFTER change: (a) -------------- c --- Add quadrature contribution to new sigmas sye1=syae1 sze1=szae1 if(sysq.GT.zero .AND. iold.EQ.1) sye1=SQRT(syae1*syae1+sysq) if(bidfnl(ii).GT.zero) sze1=SQRT(szae1*szae1+bidfnl(ii)) BEFORE change: (b) -------------- c --- Add quadrature contribution to new sigmas syb2=syab2 szb2=szab2 if(sysq.GT.zero) then ratio=sysq/(syab2*syab2) if(ratio.LE.0.2) then c --- Approximate square root syb2=syab2*(one+half*ratio) else syb2=syab2*sqrt(one+ratio) endif endif if(bidfnl(ii).GT.zero) then ratio=bidfnl(ii)/(szab2*szab2) if(ratio.LE.0.2) then c --- Approximate square root szb2=szab2*(one+half*ratio) else szb2=szab2*sqrt(one+ratio) endif endif AFTER change: (b) -------------- c --- Add quadrature contribution to new sigmas syb2=syab2 szb2=szab2 if(sysq.GT.zero) syb2=SQRT(syab2*syab2+sysq) if(bidfnl(ii).GT.zero) szb2=SQRT(szab2*szab2+bidfnl(ii)) BEFORE change: (c) -------------- c --- Add contribution of buoyancy-enhancement sye2=syae2 sze2=szae2 if(sysq.GT.zero) then ratio=sysq/(syae2*syae2) if(ratio.LE.0.2) then c --- Approximate square root sye2=syae2*(one+half*ratio) else sye2=syae2*sqrt(one+ratio) endif endif if(bidfnl(ii).GT.zero) then ratio=bidfnl(ii)/(szae2*szae2) if(ratio.LE.0.2) then c --- Approximate square root sze2=szae2*(one+half*ratio) else sze2=szae2*sqrt(one+ratio) endif endif AFTER change: (c) -------------- c --- Add contribution of buoyancy-enhancement sye2=syae2 sze2=szae2 if(sysq.GT.zero) sye2=SQRT(syae2*syae2+sysq) if(bidfnl(ii).GT.zero) sze2=SQRT(szae2*szae2+bidfnl(ii)) ---------------------------------------- Subroutine WAKE_DFSN (Problem Area 2) ---------------------------------------- BEFORE change: (a) -------------- logical ldbhr,lcav,lwak,lrevcav AFTER change: (a) -------------- logical ldbhr,lcav,lwak,lrevcav c --- Local control for imposing ISC-PRIME constructs logical L_ISCPRIME BEFORE change: (b) -------------- c --- Target stepping length (m) data dx0/2.0/ AFTER change: (b) -------------- c --- Target stepping length (m) data dx0/2.0/ c --- ISC-PRIME mode data L_ISCPRIME/.TRUE./ BEFORE change: (c) -------------- c --- Initialize CAVITY parameters AFTER change: (c) -------------- c --- Add BID to entry sigmas szibid=SQRT(szi*szi+bidsq) syibid=SQRT(syi*syi+bidsq) c --- Recompute virtuals for these sigmas (km) call SIGTZ(szibid,zero,zero,Hb,szdum,tzi,xzikm) call SIGTY(syibid,zero,zero,sydum,tyi,xyikm) if(LDBHR) then write(io6,*) write(io6,*)'----- WAKE_DFSN: Wake Entry Sigmas' write(io6,*)'xikm,szi_amb,syi_amb = ',xikm,szi,syi write(io6,*)' szi_bid,syi_bid = ',szibid,syibid endif c --- Initialize CAVITY parameters BEFORE change: (c) -------------- c --- Set initial sigmas for cavity source using sigma-y at xi, c --- including BID syibid=SQRT(syi*syi+bidsq) call WAKE_CAV0(syibid,szcav0,sycav0) AFTER change: (c) -------------- c --- Set initial sigmas for cavity source using sigma-y at xi, c --- including BID call WAKE_CAV0(syibid,szcav0,sycav0) BEFORE change: (d) -------------- else lwak=.TRUE. xwak(1)=xi+xbadj szwak(1)=szi sywak(1)=syi drwak(1)=zero c --- Compute corresponding virtual distances AFTER change: (d) -------------- else lwak=.TRUE. xwak(1)=xi+xbadj szwak(1)=szibid sywak(1)=syibid drwak(1)=zero c --- Set corresponding virtual distances BEFORE change: (e) -------------- c --- Fill first element of marching arrays using values at xlow c --- (dist array measures distance downwind from primary source) dist(1)=xlow+xbadj if(lwak) then asigz(1)=szi asigy(1)=syi AFTER change: (e) -------------- c --- Fill first element of marching arrays using values at xlow c --- (dist array measures distance downwind from primary source) dist(1)=xlow+xbadj if(lwak) then asigz(1)=szibid asigy(1)=syibid BEFORE change: (f) -------------- c --- Check to see if cavity data should be revised based on c --- wake sigma from previous step, with BID added c --- BID for position xi is used here (approximation) if(lrevcav .AND. xold.GE.xLb) then syibid=SQRT(asigy(n-1)*asigy(n-1)+bidsq) call WAKE_CAV0(syibid,szcav0,sycav0r) AFTER change: (f) -------------- c --- Check to see if cavity data should be revised based on c --- wake sigma from previous step, with BID if(lrevcav .AND. xold.GE.xLb) then call WAKE_CAV0(asigy(n-1),szcav0,sycav0r) BEFORE change: (g) -------------- c --- Check for larger ambient sigma xzkm=(dist(n)+xzvwak)*fm2km xykm=(dist(n)+xyvwak)*fm2km tz=xzkm/ubkps ty=xykm/ubkps call WAKE_SIGA('z',asigz(n-1),dxkm,dts,Hb,xzkm,tz, & asigza,zdxdum,zdtdum) call WAKE_SIGA('y',asigy(n-1),dxkm,dts,zero,xykm,ty, & asigya,ydxdum,ydtdum) asigz(n)=AMAX1(asigz(n),asigza) asigy(n)=AMAX1(asigy(n),asigya) AFTER change: (g) -------------- c --- Check for larger ambient sigma c --- ISC-PRIME does not do this test, so restrict it with c --- a control variable if(.NOT.L_ISCPRIME) then xzkm=(dist(n)+xzvwak)*fm2km xykm=(dist(n)+xyvwak)*fm2km tz=xzkm/ubkps ty=xykm/ubkps call WAKE_SIGA('z',asigz(n-1),dxkm,dts,Hb, & xzkm,tz,asigza,zdxdum,zdtdum) call WAKE_SIGA('y',asigy(n-1),dxkm,dts,zero, & xykm,ty,asigya,ydxdum,ydtdum) asigz(n)=AMAX1(asigz(n),asigza) asigy(n)=AMAX1(asigy(n),asigya) endif BEFORE change: (h) -------------- c --- Check for larger ambient growth rate distcav=x-xbc xzkm=(distcav+xzvcav)*fm2km xykm=(distcav+xyvcav)*fm2km tz=xzkm/ubkps ty=xykm/ubkps call WAKE_SIGA('z',csigz(n-1),dxkm,dts,Hb,xzkm,tz, & csigza,zdxdum,zdtdum) call WAKE_SIGA('y',csigy(n-1),dxkm,dts,zero,xykm,ty, & csigya,ydxdum,ydtdum) csigz(n)=AMAX1(csigz(n),csigza) csigy(n)=AMAX1(csigy(n),csigya) AFTER change: (h) -------------- c --- Check for larger ambient growth rate c --- ISC-PRIME does not do this test, so restrict it with c --- a control variable if(.NOT.L_ISCPRIME) then distcav=x-xbc xzkm=(distcav+xzvcav)*fm2km xykm=(distcav+xyvcav)*fm2km tz=xzkm/ubkps ty=xykm/ubkps call WAKE_SIGA('z',csigz(n-1),dxkm,dts,Hb, & xzkm,tz,csigza,zdxdum,zdtdum) call WAKE_SIGA('y',csigy(n-1),dxkm,dts,zero, & xykm,ty,csigya,ydxdum,ydtdum) csigz(n)=AMAX1(csigz(n),csigza) csigy(n)=AMAX1(csigy(n),csigya) endif BEFORE change: (i) -------------- c --- Primary source --- if(lwak .AND. (xi.LE.xaz)) then c --- Compute sigma at xaz call WAKE_SIG(xaz,xd,xold,wakiz,wakiy,asigz(n-1), & asigy(n-1),Hb,Wb,Rb,zk,ykdum, & sigzxa,sydum,dzrate) c --- Check for larger ambient sigma xasrc=xaz+xbadj xzkm=(xasrc+xzvwak)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',asigz(n-1),dxakm,dtas,Hb, & xzkm,tz,asigzxa,dxdum,dtdum) c --- Test wake result if(asigzxa.GE.sigzxa) then sigzxa=asigzxa xzvwak=dxdum tzvwak=dtdum else call SIGTZ(sigzxa,zero,zero,Hb,szdum, & tzvwak,xzvwak) endif c --- Set virtual distance increment in m xzvwak=xzvwak*fkm2m-xasrc AFTER change: (i) -------------- c --- Primary source --- if(lwak .AND. (xi.LE.xaz)) then c --- Compute sigma at xaz xasrc=xaz+xbadj call WAKE_SIG(xaz,xd,xold,wakiz,wakiy,asigz(n-1), & asigy(n-1),Hb,Wb,Rb,zk,ykdum, & sigzxa,sydum,dzrate) call SIGTZ(sigzxa,zero,zero,Hb,szdum, & tzvwak,xzvwakkm) c --- Check for larger ambient sigma c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then xzkm=(xasrc+xzvwak)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',asigz(n-1),dxakm,dtas,Hb, & xzkm,tz,asigzxa,dxdum,dtdum) c --- Test wake result if(asigzxa.GE.sigzxa) then sigzxa=asigzxa xzvwakkm=dxdum tzvwak=dtdum endif endif c --- Set virtual distance increment in m xzvwak=xzvwakkm*fkm2m-xasrc BEFORE change: (j) -------------- c --- Cavity source --- if(lcav .AND. (xbc.LE.xaz)) then c --- Compute sigma at xaz call WAKE_SIG(xaz,xdc,xold,wakiz,wakiy,csigz(n-1), & csigy(n-1),Hb,Wb,Rb,zkc,ykdum, & sigzxa,sydum,dzrate) c --- Check for larger ambient growth rate xasrc=xaz-xbc xzkm=(xasrc+xzvcav)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',csigz(n-1),dxakm,dtas,Hb, & xzkm,tz,csigzxa,dxdum,dtdum) if(csigzxa.GE.sigzxa) then sigzxa=csigzxa xzvcav=dxdum tzvcav=dtdum else call SIGTZ(sigzxa,zero,zero,Hb,szdum, & tzvcav,xzvcav) endif c --- Set virtual distance increment in m xzvcav=xzvcav*fkm2m-xasrc AFTER change: (j) -------------- c --- Cavity source --- if(lcav .AND. (xbc.LE.xaz)) then c --- Compute sigma at xaz xasrc=xaz-xbc call WAKE_SIG(xaz,xdc,xold,wakiz,wakiy,csigz(n-1), & csigy(n-1),Hb,Wb,Rb,zkc,ykdum, & sigzxa,sydum,dzrate) call SIGTZ(sigzxa,zero,zero,Hb,szdum, & tzvcav,xzvcavkm) c --- Check for larger ambient growth rate c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then xzkm=(xasrc+xzvcav)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',csigz(n-1),dxakm,dtas,Hb, & xzkm,tz,csigzxa,dxdum,dtdum) if(csigzxa.GE.sigzxa) then sigzxa=csigzxa xzvcavkm=dxdum tzvcav=dtdum endif endif c --- Set virtual distance increment in m xzvcav=xzvcavkm*fkm2m-xasrc BEFORE change: (k) -------------- c --- Primary source --- if(lwak .AND. (xi.LE.x)) then c --- Compute sigmaz call WAKE_SIG(x,xd,xold,wakiz,wakiy,asigz(n-1), & asigy(n-1),Hb,Wb,Rb,zk,ykdum, & asigz(n),sydum,dsz(n)) c --- Check for larger ambient sigma xzkm=(dist(n)+xzvwak)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',asigz(n-1),dxkm,dts,Hb,xzkm,tz, & asigza,dxdum,dtdum) asigz(n)=AMAX1(asigz(n),asigza) AFTER change: (k) -------------- c --- Primary source --- if(lwak .AND. (xi.LE.x)) then c --- Compute sigmaz call WAKE_SIG(x,xd,xold,wakiz,wakiy,asigz(n-1), & asigy(n-1),Hb,Wb,Rb,zk,ykdum, & asigz(n),sydum,dsz(n)) c --- Check for larger ambient sigma c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then xzkm=(dist(n)+xzvwak)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',asigz(n-1),dxkm,dts,Hb, & xzkm,tz,asigza,dxdum,dtdum) asigz(n)=AMAX1(asigz(n),asigza) endif BEFORE change: (l) -------------- c --- Cavity source --- if(lcav .AND. (xbc.LE.x)) then call WAKE_SIG(x,xdc,xold,wakiz,wakiy,csigz(n-1), & csigy(n-1),Hb,Wb,Rb,zkc,ykdum, & csigz(n),sydum,dzrate) c --- Check for larger ambient growth rate distcav=x-xbc xzkm=(distcav+xzvcav)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',csigz(n-1),dxkm,dts,Hb,xzkm,tz, & csigza,zdxdum,zdtdum) csigz(n)=AMAX1(csigz(n),csigza) AFTER change: (l) -------------- c --- Cavity source --- if(lcav .AND. (xbc.LE.x)) then call WAKE_SIG(x,xdc,xold,wakiz,wakiy,csigz(n-1), & csigy(n-1),Hb,Wb,Rb,zkc,ykdum, & csigz(n),sydum,dzrate) c --- Check for larger ambient growth rate c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then distcav=x-xbc xzkm=(distcav+xzvcav)*fm2km tz=xzkm/ubkps call WAKE_SIGA('z',csigz(n-1),dxkm,dts,Hb, & xzkm,tz,csigza,zdxdum,zdtdum) csigz(n)=AMAX1(csigz(n),csigza) endif BEFORE change: (m) -------------- c --- Primary source --- if(lwak .AND. (xi.LE.xay)) then c --- Compute sigma at xay call WAKE_SIG(xay,xd,xold,wakiz,wakiy,asigz(n-1), & asigy(n-1),Hb,Wb,Rb,zkdum,yk, & szdum,sigyxa,dzrate) c --- Check for larger ambient sigma xasrc=xay+xbadj xykm=(xasrc+xyvwak)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',asigy(n-1),dxakm,dtas,zero, & xykm,ty,asigyxa,dxdum,dtdum) if(asigyxa.GE.sigyxa) then sigyxa=asigyxa xyvwak=dxdum tyvwak=dtdum else call SIGTY(sigyxa,zero,zero,sydum, & tyvwak,xyvwak) endif c --- Set virtual distance increment in m xyvwak=xyvwak*fkm2m-xasrc AFTER change: (m) -------------- c --- Primary source --- if(lwak .AND. (xi.LE.xay)) then c --- Compute sigma at xay xasrc=xay+xbadj call WAKE_SIG(xay,xd,xold,wakiz,wakiy,asigz(n-1), & asigy(n-1),Hb,Wb,Rb,zkdum,yk, & szdum,sigyxa,dzrate) call SIGTY(sigyxa,zero,zero,sydum, & tyvwak,xyvwakkm) c --- Check for larger ambient sigma c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then xykm=(xasrc+xyvwak)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',asigy(n-1),dxakm,dtas,zero, & xykm,ty,asigyxa,dxdum,dtdum) if(asigyxa.GE.sigyxa) then sigyxa=asigyxa xyvwakkm=dxdum tyvwak=dtdum endif endif c --- Set virtual distance increment in m xyvwak=xyvwakkm*fkm2m-xasrc BEFORE change: (n) -------------- c --- Cavity source --- if(lcav .AND. (xbc.LE.xay)) then c --- Compute sigma at xay call WAKE_SIG(xay,xdc,xold,wakiz,wakiy,csigz(n-1), & csigy(n-1),Hb,Wb,Rb,zkdum,ykc, & szdum,sigyxa,dzrate) c --- Check for larger ambient growth rate xasrc=xay-xbc xykm=(xasrc+xyvcav)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',csigy(n-1),dxakm,dtas,zero, & xykm,ty,csigyxa,dxdum,dtdum) if(csigyxa.GE.sigyxa) then sigyxa=csigyxa xyvcav=dxdum tyvcav=dtdum else call SIGTY(sigyxa,zero,zero,sydum, & tyvcav,xyvcav) endif c --- Set virtual distance increment in m xyvcav=xyvcav*fkm2m-xasrc AFTER change: (n) -------------- c --- Cavity source --- if(lcav .AND. (xbc.LE.xay)) then c --- Compute sigma at xay xasrc=xay-xbc call WAKE_SIG(xay,xdc,xold,wakiz,wakiy,csigz(n-1), & csigy(n-1),Hb,Wb,Rb,zkdum,ykc, & szdum,sigyxa,dzrate) call SIGTY(sigyxa,zero,zero,sydum, & tyvcav,xyvcavkm) c --- Check for larger ambient growth rate c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then xykm=(xasrc+xyvcav)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',csigy(n-1),dxakm,dtas,zero, & xykm,ty,csigyxa,dxdum,dtdum) if(csigyxa.GE.sigyxa) then sigyxa=csigyxa xyvcavkm=dxdum tyvcav=dtdum endif endif c --- Set virtual distance increment in m xyvcav=xyvcavkm*fkm2m-xasrc BEFORE change: (p) -------------- c --- Check for larger ambient sigma xykm=(dist(n)+xyvwak)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',asigy(n-1),dxkm,dts,zero, & xykm,ty,asigya,dxdum,dtdum) asigy(n)=AMAX1(asigy(n),asigya) AFTER change: (p) -------------- c --- Check for larger ambient sigma c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then xykm=(dist(n)+xyvwak)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',asigy(n-1),dxkm,dts,zero, & xykm,ty,asigya,dxdum,dtdum) asigy(n)=AMAX1(asigy(n),asigya) endif BEFORE change: (q) -------------- c --- Check for larger ambient growth rate distcav=x-xbc xykm=(distcav+xyvcav)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',csigy(n-1),dxkm,dts,zero, & xykm,ty,csigya,dxdum,dtdum) csigy(n)=AMAX1(csigy(n),csigya) AFTER change: (q) -------------- c --- Check for larger ambient growth rate c --- ISC-PRIME does not do this test, so restrict it c --- with a control variable if(.NOT.L_ISCPRIME) then distcav=x-xbc xykm=(distcav+xyvcav)*fm2km ty=xykm/ubkps call WAKE_SIGA('y',csigy(n-1),dxkm,dts,zero, & xykm,ty,csigya,dxdum,dtdum) csigy(n)=AMAX1(csigy(n),csigya) endif BEFORE change: (r) -------------- c --- Add points for portion of trajectory that lies upwind of xi, c --- using the same step interval, including the source location xwak1=xi+xbadj npup=xwak1*dxi+1 if(npup.GE.1) then nptot=npw+npup xwak(1)=zero szwak(1)=sz0 sywak(1)=sy0 drwak(1)=zero else write(io6,*)' ERROR: plume enters wake UPWIND of source' write(io6,*)' Distance from source (m) = ',xwak1 stop endif AFTER change: (r) -------------- c --- Add points for portion of trajectory that lies upwind of xi, c --- using the same step interval, including the source location xwak1=xi+xbadj npup=xwak1*dxi if(xwak1.LT.0.0) then c --- This may need more care to avoid a false halt caused by c --- precision write(io6,*)' ERROR: plume enters wake UPWIND of source' write(io6,*)' Distance from source (m) = ',xwak1 stop elseif(xwak1.GT.0.0) then c --- Keep npup=0 case only for xwak1 really = 0.0 npup=MAX(1,npup) endif if(npup.GE.1) then c --- Upwind of xi-entry so make sigmas negative to indicate c --- that BID needs to be added in WAKE_ADDBID call nptot=npw+npup xwak(1)=zero szwak(1)=-sz0 sywak(1)=-sy0 drwak(1)=zero endif BEFORE change: (s) -------------- if(nptot.LE.2*mxrise) then c --- Fill elements with nearest values deln=FLOAT(nptot)/FLOAT(nwak) do in=2,nwak-1 jn=in*deln if(jn.LE.npup) then c --- Compute ambient sigmas xtraj=jn*dx xkm=xtraj*fm2km+xz0km ts=xkm/ubkps call SIGTZ(zero,xkm,ts,Hb,szwak(in),dtdum,dxdum) xkm=xtraj*fm2km+xy0km ts=xkm/ubkps call SIGTY(zero,xkm,ts,sywak(in),dtdum,dxdum) xwak(in)=xtraj drwak(in)=zero else AFTER change: (s) -------------- if(nptot.LE.2*mxrise) then c --- Fill elements with nearest values deln=FLOAT(nptot)/FLOAT(nwak) do in=2,nwak-1 jn=(in-1)*deln if(jn.LE.npup) then c --- Compute ambient sigmas xtraj=jn*dx xkm=xtraj*fm2km+xz0km ts=xkm/ubkps call SIGTZ(zero,xkm,ts,Hb,szwak(in),dtdum,dxdum) xkm=xtraj*fm2km+xy0km ts=xkm/ubkps call SIGTY(zero,xkm,ts,sywak(in),dtdum,dxdum) xwak(in)=xtraj drwak(in)=zero c --- Upwind of xi-entry so make sigmas negative to c --- indicate that BID needs to be added in WAKE_ADDBID szwak(in)=-szwak(in) sywak(in)=-sywak(in) else BEFORE change: (t) -------------- else c --- Use sliding step-size to sample nearfield more frequently deln=2.*FLOAT(nptot-mxrise)/FLOAT(mxrise*(mxrise-1)) rn=one do in=2,nwak-1 rn=rn+one+(in-1)*deln if(INT(rn).LE.npup) then c --- Compute ambient sigmas xtraj=rn*dx xkm=xtraj*fm2km+xz0km ts=xkm/ubkps call SIGTZ(zero,xkm,ts,Hb,szwak(in),dtdum,dxdum) xkm=xtraj*fm2km+xy0km ts=xkm/ubkps call SIGTY(zero,xkm,ts,sywak(in),dtdum,dxdum) xwak(in)=xtraj drwak(in)=zero AFTER change: (t) -------------- else c --- Use sliding step-size to sample nearfield more frequently deln=2.*FLOAT(nptot-mxrise)/FLOAT(mxrise*(mxrise-1)) rn=one do in=2,nwak-1 rn=rn+one+(in-1)*deln if(INT(rn).LE.npup) then c --- Compute ambient sigmas xtraj=rn*dx xkm=xtraj*fm2km+xz0km ts=xkm/ubkps call SIGTZ(zero,xkm,ts,Hb,szwak(in),dtdum,dxdum) xkm=xtraj*fm2km+xy0km ts=xkm/ubkps call SIGTY(zero,xkm,ts,sywak(in),dtdum,dxdum) xwak(in)=xtraj drwak(in)=zero c --- Upwind of xi-entry so make sigmas negative to c --- indicate that BID needs to be added in WAKE_ADDBID szwak(in)=-szwak(in) sywak(in)=-sywak(in) BEFORE change: (u) -------------- c --- Fill only those elements used nwak=nptot do in=2,nptot if(in.LE.npup) then c --- Compute ambient sigmas xtraj=in*dx xkm=xtraj*fm2km+xz0km ts=xkm/ubkps call SIGTZ(zero,xkm,ts,Hb,szwak(in),dtdum,dxdum) xkm=xtraj*fm2km+xy0km ts=xkm/ubkps call SIGTY(zero,xkm,ts,sywak(in),dtdum,dxdum) xwak(in)=xtraj drwak(in)=zero AFTER change: (u) -------------- c --- Fill only those elements used nwak=nptot if(npup.GT.0) delta_x=xwak1/FLOAT(npup) do in=2,nptot if(in.LE.npup) then c --- Compute ambient sigmas xtraj=(in-1)*delta_x xkm=xtraj*fm2km+xz0km ts=xkm/ubkps call SIGTZ(zero,xkm,ts,Hb,szwak(in),dtdum,dxdum) xkm=xtraj*fm2km+xy0km ts=xkm/ubkps call SIGTY(zero,xkm,ts,sywak(in),dtdum,dxdum) xwak(in)=xtraj drwak(in)=zero c --- Upwind of xi-entry so make sigmas negative to c --- indicate that BID needs to be added in WAKE_ADDBID szwak(in)=-szwak(in) sywak(in)=-sywak(in) BEFORE change: (v) -------------- endif if(lcav) then c --- CAV arrays: AFTER change: (v) -------------- c --- Repeat procedure to store DRWAKISC data (as in ISC-PRIME) c --- starting at xi c --- Place initial values into first element xwakisc(1)=xi+xbadj drwakisc(1)=zero if(npw.GE.mxrise) then c --- Sample a subset of the nptot points nwakisc=mxrise xwakisc(nwakisc)=dist(np) drwakisc(nwakisc)=rtpiby2*dsz(np) if(npw.LE.2*mxrise) then c --- Fill elements with nearest values deln=FLOAT(npw)/FLOAT(nwakisc) do in=2,nwakisc-1 jn=in*deln+nws xwakisc(in)=dist(jn) drwakisc(in)=rtpiby2*dsz(jn) enddo else c --- Use sliding step-size to sample nearfield more frequently deln=2.*FLOAT(npw-mxrise)/FLOAT(mxrise*(mxrise-1)) rn=one do in=2,nwakisc-1 rn=rn+one+(in-1)*deln jn=rn+nws xwakisc(in)=dist(jn) drwakisc(in)=rtpiby2*dsz(jn) enddo endif else c --- Fill only those elements used nwakisc=npw do in=2,npw inp=in+nws xwakisc(in)=dist(inp) drwakisc(in)=rtpiby2*dsz(inp) enddo endif endif if(lcav) then c --- CAV arrays: ------------------------------------- Subroutine WAKE_DRDX (Problem Area 2) ------------------------------------- BEFORE change: (a) -------------- c Common block /WAKEDAT/ variables: c NWAK, XWAK(mxrise), DRWAK(mxrise) AFTER change: (a) -------------- c Common block /WAKEDAT/ variables: c NWAKISC, XWAKISC(mxrise), DRWAKISC(mxrise) BEFORE change: (b) -------------- if(x.gt.xwak(nwak) .OR. x.lt.xwak(1))then drdx=0.0 elseif(nwak.le.1) then AFTER change: (b) -------------- if(x.gt.xwakisc(nwakisc) .OR. x.lt.xwakisc(1))then drdx=0.0 elseif(nwakisc.le.1) then BEFORE change: (c) -------------- nwkm1=nwak-1 drdx=drwak(1) do i=nwkm1,1,-1 if(x.ge.xwak(i))then ip1=i+1 drdx=drwak(ip1)-(drwak(ip1)-drwak(i))*(xwak(ip1)-x)/ 1 (xwak(ip1)-xwak(i)) AFTER change: (c) -------------- nwkm1=nwakisc-1 drdx=drwakisc(1) do i=nwkm1,1,-1 if(x.ge.xwakisc(i))then ip1=i+1 drdx=drwakisc(ip1)-(xwakisc(ip1)-x)* 1 (drwakisc(ip1)-drwakisc(i))/(xwakisc(ip1)-xwakisc(i)) ------------------------------------ Subroutine WAKE_FIN (Problem Area 2) ------------------------------------ BEFORE change: (a) -------------- c --- Assign effective sigmas at the source which reproduce the c --- sigmas at the end of the wake region for the dispersion c --- conditions at the time of release c --- Primary source c ------------------ xzvkm=xzvwak*fm2km tzv=xzvkm/ubkps call SIGTZ(zero,xzvkm,tzv,Hb,szw0,tdum,xdum) xyvkm=xyvwak*fm2km tyv=xyvkm/ubkps call SIGTY(zero,xyvkm,tyv,syw0,tdum,xdum) AFTER change: (a) -------------- c --- Assign effective sigmas at the source which reproduce the c --- sigmas at the end of the wake region for the dispersion c --- conditions at the time of release (clip negative at zero) c --- Primary source c ------------------ szw0=0.0 syw0=0.0 if(xzvwak.GE.0.0) then xzvkm=xzvwak*fm2km tzv=xzvkm/ubkps call SIGTZ(zero,xzvkm,tzv,Hb,szw0,tdum,xdum) endif if(xyvwak.GE.0.0) then xyvkm=xyvwak*fm2km tyv=xyvkm/ubkps call SIGTY(zero,xyvkm,tyv,syw0,tdum,xdum) endif BEFORE change: (b) -------------- c --- Debug output AFTER change: (b) -------------- c --- Debug output c --- Take ABS of tabulated sigmas because they are negative c --- for x= HLID0) BEFORE change: (c) -------------- c --- VCOUP called by: VCBAR, SLUGAVE, SLUGSNP, SLUGINT c LSSLINT, CALCPF AFTER change: (c) -------------- c --- VCOUP called by: VCBAR, SLUGAVE, SLUGSNP, SLUGINT c LSSLINT, CALCPF, CAV_SAMP, PLMFOG BEFORE change: (d) -------------- c --- Treat well-mixed puff/slug codes first c ------------------------------------------ iicode=icode if(iicode.GT.10) iicode=iicode-10 if(iicode.EQ.2 .OR. iicode.EQ.6) then c --- Within past/present mixed layer and uniform vcoup=1./hlid0 return elseif(iicode .EQ. 4) then c --- Above mixed layer and uniform (no ground-level concs) vcoup=0. return endif AFTER change: (d) -------------- c --- Treat well-mixed puff/slug codes first c ------------------------------------------ iicode=icode if(iicode.GT.10) iicode=iicode-10 if(iicode.EQ.2 .OR. iicode.EQ.6) then c --- Within past/present mixed layer and uniform vcoup=0.0 if(zr.LE.hlid0) then c --- Distribution in surface layer vcoup=1./hlid0 elseif(zr.LE.hlidmax0) then c --- Distribution above surface layer vcoup=1./(hlidmax0-hlid0) endif return elseif(iicode .EQ. 4) then c --- Above mixed layer and uniform (no defined concs) vcoup=0.0 return endif BEFORE change: (e) -------------- c --- Distribution must be computed c --------------------------------- vcoup = 0.0 c --- Screen for sigma_z that is 'zero' if(sz.LT.szcut) return c c --- Remove mixing lid if source ht > mixing ht if(zs.GT.hlid0) then hlid=zunlim else hlid=hlid0 endif c c --- Check for significant trapping if((sz/hlid).gt.0.63) go to 15 AFTER change: (e) -------------- c --- Gaussian Distribution c ------------------------- vcoup = 0.0 c c --- Remove mixing lid if source ht > mixing ht c --- (Unlimited mixing for GP above current lid) if(zs.GT.hlid0) then hlid=zunlim hlidmax=zunlim else hlid=hlid0 hlidmax=AMAX1(hlid0,hlidmax0) endif c --- Case of receptor above current lid sampling mass aloft c --- (Either zero or mixed aloft) if(zr.GT.hlid) then if(zr.LE.hlidmax) then vcoup=1./(hlidmax-hlid) else vcoup=0.0 endif return endif c --- Screen for sigma_z that is 'zero' if(sz.LT.szcut) return c c --- Check for significant trapping if((sz/hlid).gt.0.63) go to 15 *********************************************************************** -- Problem Area 7 -- Trap case in which the shear adjustment to Briggs plume rise is applied in near-calm conditions. A distance to final rise of zero (calm) can produce an incorrect shear-modified rise that is zero. Also increase the number of digits used to write TSTAK and DTEMP to the list file (debug output). *********************************************************************** ------------------------------------ Subroutine POINTS1 (Problem Area 7) ------------------------------------ BEFORE change: (a) -------------- c --- Compute vertical wind shear effects if(mshear.EQ.1)then AFTER change: (a) -------------- c --- Compute vertical wind shear effects if(mshear.EQ.1 .AND. xfinb.GT.0.0)then BEFORE change: (b) -------------- 206 format(5x,'HTSTAK: ',f6.1,2x,'DIAM: ',f6.3,2x,'EXITW: ', 1 f6.3,2x,'TSTAK: ',f6.2,2x,'TEMPSS: ',f6.2,2x,'DTEMP: ',f6.3, 2 2x,'FLUXB : ',f8.3,2x,'FLUXM: ',f6.1/5x,'XF: ',f6.1,2x, 3 'HEFF: ',f6.1,2x,'TIPDW: ',f6.1) AFTER change: (b) -------------- 206 format(5x,'HTSTAK: ',f6.1,2x,'DIAM: ',f6.3,2x,'EXITW: ', 1 f6.3,2x,'TSTAK: ',f8.2,2x,'TEMPSS: ',f6.2,2x,'DTEMP: ',f8.3, 2 2x,'FLUXB : ',f8.3,2x,'FLUXM: ',f6.1/5x,'XF: ',f6.1,2x, 3 'HEFF: ',f6.1,2x,'TIPDW: ',f6.1) ------------------------------------ Subroutine POINTS2 (Problem Area 7) ------------------------------------ BEFORE change: (a) -------------- c --- Compute vertical wind shear effects if(mshear.eq.1)then AFTER change: (a) -------------- c --- Compute vertical wind shear effects if(mshear.EQ.1 .AND. xfinb.GT.0.0)then BEFORE change: (b) -------------- 206 format(5x,'HTSTAK: ',f6.1,2x,'DIAM: ',f6.3,2x,'EXITW: ', 1 f6.3,2x,'TSTAK: ',f6.2,2x,'TEMPSS: ',f6.2,2x,'DTEMP: ',f6.3, 2 2x,'FLUXB : ',f8.3,2x,'FLUXM: ',f6.1/5x,'XF: ',f6.1,2x, 3 'HEFF: ',f6.1,2x,'TIPDW: ',f6.1) AFTER change: (b) -------------- 206 format(5x,'HTSTAK: ',f6.1,2x,'DIAM: ',f6.3,2x,'EXITW: ', 1 f6.3,2x,'TSTAK: ',f8.2,2x,'TEMPSS: ',f6.2,2x,'DTEMP: ',f8.3, 2 2x,'FLUXB : ',f8.3,2x,'FLUXM: ',f6.1/5x,'XF: ',f6.1,2x, 3 'HEFF: ',f6.1,2x,'TIPDW: ',f6.1) *********************************************************************** -- Problem Area 8 -- Assign power-law exponent to default value for met file types that do not use power-law profiles (CALMET, CTDM PROFILE, and AERMET PROFILE) so that exponent is available for Briggs plume rise wind shear modification option (MSHEAR). *********************************************************************** -------------------------------------- Subroutine BLOCK DATA (Problem Area 8) -------------------------------------- BEFORE change: (a) -------------- c --- Initialize observed turbulence profiles to missing (-999.) data svprf/mxprfz*-999./ data swprf/mxprfz*-999./ AFTER change: (a) -------------- c --- Initialize observed turbulence profiles to missing (-999.) data svprf/mxprfz*-999./ data swprf/mxprfz*-999./ c --- Initialize observed wind shear power-law exponent to missing data plexp/-999./ ------------------------------------ Subroutine ADVECT (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c WSCALM AFTER change: (a) -------------- c Common block /COMPARM/ variables: c WSCALM c WSCALM, PLX0 BEFORE change: (b) -------------- if(metfm.LE.3) then AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif if(metfm.LE.3) then BEFORE change: (c) -------------- call POWLAW(ht,zgpt(1),wsold,plexp,ws) AFTER change: (c) -------------- call POWLAW(ht,zgpt(1),wsold,plx,ws) ------------------------------------ Subroutine AREAS1 (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM AFTER change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM, PLX0 BEFORE change: (b) -------------- c --- Swap emission rate for source 'i' into 1D array and scale AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c --- Swap emission rate for source 'i' into 1D array and scale BEFORE change: (c) -------------- plexp0(np)=plexp AFTER change: (c) -------------- plexp0(np)=plx ------------------------------------ Subroutine AREAS2 (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM AFTER change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM, PLX0 BEFORE change: (b) -------------- c c -------------------------------- c --- Numerical Plume Rise Section AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c c -------------------------------- c --- Numerical Plume Rise Section BEFORE change: (c) -------------- plexp0(np)=plexp AFTER change: (c) -------------- plexp0(np)=plx ------------------------------------ Subroutine LINES1 (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM AFTER change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM, PLX0 BEFORE change: (b) -------------- c --- Set logical to identify calm conditions AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c --- Set logical to identify calm conditions BEFORE change: (c) -------------- call setline(mshear,theta,ws,istab,sqrts,plexp,nlines,xl, & hbl,wml,fprimel,wsep,fptot,fbpt,ntr, & xle,xld,rzero,xtr,ztr) AFTER change: (c) -------------- call setline(mshear,theta,ws,istab,sqrts,plx,nlines,xl, & hbl,wml,fprimel,wsep,fptot,fbpt,ntr, & xle,xld,rzero,xtr,ztr) BEFORE change: (d) -------------- plexp0(np)=plexp AFTER change: (d) -------------- plexp0(np)=plx ------------------------------------ Subroutine LINES2 (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM AFTER change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM, PLX0 BEFORE change: (b) -------------- c --- Calculate the virtual times and sigmas at Heffter transition AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c --- Calculate the virtual times and sigmas at Heffter transition BEFORE change: (c) -------------- call setline(mshear,theta,ws,istab,sqrts,plexp,nl2(ig),xl2(ig), & hbl2(ig),wml2(ig),fprimel2(ig),wsep2(ig),fptot2(ig), & fbpt2(ig),ntr,xle,xld,rzero,xtr,ztr) AFTER change: (c) -------------- call setline(mshear,theta,ws,istab,sqrts,plx,nl2(ig),xl2(ig), & hbl2(ig),wml2(ig),fprimel2(ig),wsep2(ig),fptot2(ig), & fbpt2(ig),ntr,xle,xld,rzero,xtr,ztr) BEFORE change: (d) -------------- plexp0(np)=plexp AFTER change: (d) -------------- plexp0(np)=plx ------------------------------------ Subroutine NUMMET (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /NUMPARM/ variables: AFTER change: (a) -------------- c Common block /COMPARM/ variables: c PLX0 c Common block /NUMPARM/ variables: BEFORE change: (b) -------------- include 'ambient.puf' AFTER change: (b) -------------- include 'ambient.puf' include 'comparm.puf' BEFORE change: (c) -------------- c --- Impose minimum value if(ptgrad.LT.ptgrad0)ptgrad=ptgrad0 AFTER change: (c) -------------- c --- Impose minimum value if(ptgrad.LT.ptgrad0)ptgrad=ptgrad0 c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default kst6=MIN(kst,6) plx=plx0(kst) else c --- Use it plx=plexp endif BEFORE change: (c) -------------- call powlaw(zgpta(iz),zref,wsref,plexp,uamb(iz)) AFTER change: (c) -------------- call powlaw(zgpta(iz),zref,wsref,plx,uamb(iz)) BEFORE change: (d) -------------- write(io6,*)'PLEXP = ',plexp AFTER change: (d) -------------- write(io6,*)'PLX = ',plx ------------------------------------ Subroutine POINTS1 (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, PTG0(2), WSCALM, TBD, PLX0 AFTER change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, PTG0(2), WSCALM, TBD, PLX0 BEFORE change: (b) -------------- c --- Swap emission rate for source 'i' into 1D array and scale AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c --- Swap emission rate for source 'i' into 1D array and scale BEFORE change: (c) -------------- call prfinsh(fluxb,ws,plexp,htstak(i),xfinb,istab,sqrts, 1 zfinsh) AFTER change: (c) -------------- call prfinsh(fluxb,ws,plx,htstak(i),xfinb,istab,sqrts, 1 zfinsh) BEFORE change: (d) -------------- plexp0(np)=plexp AFTER change: (d) -------------- plexp0(np)=plx BEFORE change: (e) -------------- plexp0(np)=plexp AFTER change: (e) -------------- plexp0(np)=plx ------------------------------------ Subroutine POINTS2 (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, PTG0(2), WSCALM, TBD AFTER change: (a) -------------- c Common block /COMPARM/ variables: c XMXLEN, MXNEW, SYMIN, SZMIN, PTG0(2), WSCALM, TBD, PLX0 BEFORE change: (b) -------------- c --- Initialize source tabulation arrays to zero AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c --- Initialize source tabulation arrays to zero BEFORE change: (c) -------------- call prfinsh(fluxb,ws,plexp,htstak,xfinb,istab,sqrts, 1 zfinsh) AFTER change: (c) -------------- call prfinsh(fluxb,ws,plx,htstak,xfinb,istab,sqrts, 1 zfinsh) BEFORE change: (d) -------------- plexp0(np)=plexp AFTER change: (d) -------------- plexp0(np)=plx BEFORE change: (e) -------------- plexp0(np)=plexp AFTER change: (e) -------------- plexp0(np)=plx ------------------------------------ Subroutine POWLAW (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- include 'comparm.puf' AFTER change: (a) -------------- include 'comparm.puf' c --- Trap invalid power-law exponent if(p.LE.0.0) then write(io6,*) write(io6,*)'POWLAW: Invalid power-law exponent found' write(io6,*)' P = ',p stop 'Halted in POWLAW --- see list file' endif ------------------------------------ Subroutine PRFINSH (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c --- PRFINSH calls: none c---------------------------------------------------------------------- AFTER change: (a) -------------- c --- PRFINSH calls: none c---------------------------------------------------------------------- c --- Trap invalid power-law exponent if(p.LE.0.0) then write(*,*) write(*,*)'PRFINSH: Invalid power-law exponent found' write(*,*)' P = ',p stop endif ------------------------------------ Subroutine PRSHEAR (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c --- PRSHEAR calls: none c---------------------------------------------------------------------- AFTER change: (a) -------------- c --- PRSHEAR calls: none c---------------------------------------------------------------------- c --- Trap invalid power-law exponent if(p.LE.0.0) then write(*,*) write(*,*)'PRSHEAR: Invalid power-law exponent found' write(*,*)' P = ',p stop endif ------------------------------------ Subroutine SETLINE (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c --- Calculate sine & cosine, taking results in first quadrant AFTER change: (a) -------------- c --- Trap invalid power-law exponent if(plawx.LE.0.0) then write(io6,*) write(io6,*)'SETLINE: Invalid power-law exponent found' write(io6,*)' PLAWX = ',plawx stop 'Halted in SETLINE --- see list file' endif c --- Calculate sine & cosine, taking results in first quadrant ------------------------------------ Subroutine VOLS (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM AFTER change: (a) -------------- c XMXLEN, MXNEW, SYMIN, SZMIN, WSCALM, PLX0 BEFORE change: (b) -------------- c --- Adjust TYPE-1 emissions for source 'i' and obtain total AFTER change: (b) -------------- c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif c --- Adjust TYPE-1 emissions for source 'i' and obtain total BEFORE change: (c) -------------- plexp0(np)=plexp AFTER change: (c) -------------- plexp0(np)=plx ------------------------------------ Subroutine WINDSET (Problem Area 8) ------------------------------------ BEFORE change: (a) -------------- c Common block /COMPARM/ variables: c WSCALM AFTER change: (a) -------------- c Common block /COMPARM/ variables: c WSCALM, PLX0 BEFORE change: (b) -------------- data rtod/57.29578/,wdcalm/180.0/ AFTER change: (b) -------------- data rtod/57.29578/,wdcalm/180.0/ c --- Trap missing power-law exponent (substitute default) if(plexp.LT.-900.) then c --- Use default plx=plx0(istab) else c --- Use it plx=plexp endif BEFORE change: (c) -------------- call POWLAW(ht,zgpt(ilayer),wsold,plexp,ws) AFTER change: (c) -------------- call POWLAW(ht,zgpt(ilayer),wsold,plx,ws) *********************************************************************** -- Problem Area 9 -- Add call to SRCTABOUT to update the source number of the current puff in the DA file when line-source slugs are processed during the first step. The DA file of source-related tabulations was introduced in MCB-E. *********************************************************************** ----------------------------------- Subroutine COMP (Problem Area 9) ----------------------------------- BEFORE change: (a) -------------- c RDAQ, GETH2O2, CHEMBK, TCHIFLX, SRCTABIN AFTER change: (a) -------------- c RDAQ, GETH2O2, CHEMBK, TCHIFLX, c SRCTABIN, SRCTABOUT BEFORE change: (b) -------------- c --- This slug is a NEW release from segment 2+ of LINE c --- so skip call to CALCSL for 1st sampling step lskip=.TRUE. c --- Reset puff ID information isnum=isnum-1000 isrcnum(ii)=isnum AFTER change: (b) -------------- c --- This slug is a NEW release from segment 2+ of LINE c --- so skip call to CALCSL for 1st sampling step lskip=.TRUE. c --- Reset puff ID information isnum=isnum-1000 isrcnum(ii)=isnum c --- Replace source tabulations that apply to this puff c --- in the DA file, using the updated source number call SRCTABOUT(ii,isrctyp(ii),isrcnum(ii), & irlsnum(ii)) BEFORE change: (c) -------------- c --- This slug is a NEW release from segment 2+ of LINE c --- so skip call to CALCSL for 1st sampling step lskip=.TRUE. c --- Reset puff ID information isnum=isnum-1000 isrcnum(ii)=isnum AFTER change: (c) -------------- c --- This slug is a NEW release from segment 2+ of LINE c --- so skip call to CALCSL for 1st sampling step lskip=.TRUE. c --- Reset puff ID information isnum=isnum-1000 isrcnum(ii)=isnum c --- Replace source tabulations that apply to this puff c --- in the DA file, using the updated source number call SRCTABOUT(ii,isrctyp(ii),isrcnum(ii), & irlsnum(ii)) *********************************************************************** -- Problem Area 10 -- Use local source-table arrays in the subroutine used to move the puff index when inactive puffs are removed from DA file of source-related tabulations. The DA file of source-related tabulations was introduced in MCB-E. *********************************************************************** ----------------------------------- Subroutine COMP (Problem Area 10) ----------------------------------- BEFORE change: (a) -------------- c --- INPUTS: c I - integer - New record (puff) index c J - integer - Old record (puff) index c c Common block /SRCTAB/ variables: c NTR, NWK, NCV, c XTR(mxrise),ZTR(mxrise),RTR(mxrise),HTR(mxrise), c XWK(mxrise),SYWK(mxrise),SZWK(mxrise),DRWK(mxrise), c XCV(mxrise),SYCV(mxrise),SZCV(mxrise) AFTER change: (a) -------------- c --- INPUTS: c I - integer - New record (puff) index c J - integer - Old record (puff) index BEFORE change: (b) -------------- c --- OUTPUT: c Common block /SRCTAB/ variables: c NTR, NWK, NCV, c XTR(mxrise),ZTR(mxrise),RTR(mxrise),HTR(mxrise), c XWK(mxrise),SYWK(mxrise),SZWK(mxrise),DRWK(mxrise), c XCV(mxrise),SYCV(mxrise),SZCV(mxrise) AFTER change: (b) -------------- c --- OUTPUT: BEFORE change: (c) -------------- c --- Include common blocks include 'srctab.puf' AFTER change: (c) -------------- c --- Declare local arrays for DA file variables read/written integer NTR, NWK, NCV, ISRCTYP, ISRCNAM, IRLSNUM real XTR(mxrise),ZTR(mxrise),RTR(mxrise),HTR(mxrise) real XWK(mxrise),SYWK(mxrise),SZWK(mxrise),DRWK(mxrise) real XCV(mxrise),SYCV(mxrise),SZCV(mxrise) *********************************************************************** -- Problem Area 11 -- Do not allow the receptor-specific sigmas for a line source to be zero. Add a check for sy0 and sz0 when puffs/slugs are generated to assure that SYMIN and SZMIN are the smallest values used (avoids halt due to divide-by-zero). *********************************************************************** ----------------------------------- Subroutine LINES1 (Problem Area 11) ----------------------------------- BEFORE change: (a) -------------- sy0(np)=sy0dw sz0(np)=sz0dw AFTER change: (a) -------------- sy0(np)=AMAX1(sy0dw,symin) sz0(np)=AMAX1(sz0dw,szmin) ----------------------------------- Subroutine LINES2 (Problem Area 11) ----------------------------------- BEFORE change: (a) -------------- sy0(np)=sy0dw sz0(np)=sz0dw AFTER change: (a) -------------- sy0(np)=AMAX1(sy0dw,symin) sz0(np)=AMAX1(sz0dw,szmin) *********************************************************************** -- Problem Area 12 -- Add check for zero final rise from a buoyant line source when computing receptor-specific sigmas for the case of an emitting slug (attached to source). The zero rise can result in a divide-by-zero that halts the model. *********************************************************************** ------------------------------------- Subroutine RECSPEC0 (Problem Area 12) ------------------------------------- BEFORE change: (a) -------------- if(zfrise.LE.zero) then AFTER change: (a) -------------- if(zfrise.LE.zero .OR. zfrises.LE.zero) then *********************************************************************** -- Problem Area 13 -- A number of variables were not properly defined or initialized. Results are not affected. *********************************************************************** ------------------------------------- Subroutine COMP (Problem Area 13) (Undefined variables --- qratew, xlam) ------------------------------------- BEFORE change: (a) -------------- data iavg/1/ AFTER change: (a) -------------- data qratew/mxspec*0.0,mxspec*0.0/ data xlam/mxspec*0.0/ data iavg/1/ ------------------------------------- Subroutine FOGREC (Problem Area 13) (Undefined variables --- ldbhr changed to ldb) ------------------------------------- BEFORE change: (a) -------------- call ADVECT(ldbhr,ixsrc,iysrc,z0m,el,dpbl,istab,zbar, & zbot,ztop,uadv,vadv) AFTER change: (a) -------------- call ADVECT(ldb,ixsrc,iysrc,z0m,el,dpbl,istab,zbar, & zbot,ztop,uadv,vadv) ------------------------------------- Subroutine LINES1 (Problem Area 13) (Undefined variables --- thalf,xhalfkm,syhalf,szhalf) ------------------------------------- BEFORE change: (a) -------------- c --- Puffs are released, so use all MXNSEG segments nsegy=0 nsegz=0 nseg(i)=mxnseg AFTER change: (a) -------------- c --- Puffs are released, so use all MXNSEG segments thalf=0. xhalfkm=0. syhalf=0. szhalf=0. nsegy=0 nsegz=0 nseg(i)=mxnseg ------------------------------------- Subroutine MET2 (Problem Area 13) (Met station coordinates relative to grid origin) ------------------------------------- BEFORE change: (a) -------------- if(nssta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xssta(1)=xorig yssta(1)=yorig AFTER change: (a) -------------- if(nssta.gt.0)then c --- Coordinates relative to the grid origin xssta(1)=0.0 yssta(1)=0.0 BEFORE change: (b) -------------- if(nusta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xusta(1)=xorig yusta(1)=yorig AFTER change: (b) -------------- if(nusta.gt.0)then c --- Coordinates relative to the grid origin xusta(1)=0.0 yusta(1)=0.0 BEFORE change: (c) -------------- if(npsta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xpsta(1)=xorig ypsta(1)=yorig AFTER change: (c) -------------- if(npsta.gt.0)then c --- Coordinates relative to the grid origin xpsta(1)=0.0 ypsta(1)=0.0 ------------------------------------- Subroutine MET3 (Problem Area 13) (Met station coordinates relative to grid origin) ------------------------------------- BEFORE change: (a) -------------- if(nssta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xssta(1)=xorig yssta(1)=yorig AFTER change: (a) -------------- if(nssta.gt.0)then c --- Coordinates relative to the grid origin xssta(1)=0.0 yssta(1)=0.0 BEFORE change: (b) -------------- if(nusta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xusta(1)=xorig yusta(1)=yorig AFTER change: (b) -------------- if(nusta.gt.0)then c --- Coordinates relative to the grid origin xusta(1)=0.0 yusta(1)=0.0 BEFORE change: (c) -------------- if(npsta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xpsta(1)=xorig ypsta(1)=yorig AFTER change: (c) -------------- if(npsta.gt.0)then c --- Coordinates relative to the grid origin xpsta(1)=0.0 ypsta(1)=0.0 ------------------------------------- Subroutine MET4 (Problem Area 13) (Met station coordinates relative to grid origin) ------------------------------------- BEFORE change: (a) -------------- if(nssta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xssta(1)=xorig yssta(1)=yorig AFTER change: (a) -------------- if(nssta.gt.0)then c --- Coordinates relative to the grid origin xssta(1)=0.0 yssta(1)=0.0 BEFORE change: (b) -------------- if(nusta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xusta(1)=xorig yusta(1)=yorig AFTER change: (b) -------------- if(nusta.gt.0)then c --- Coordinates relative to the grid origin xusta(1)=0.0 yusta(1)=0.0 BEFORE change: (c) -------------- if(npsta.gt.0)then c --- Grab the UTM (m) coordinates from the grid origin xpsta(1)=xorig ypsta(1)=yorig AFTER change: (c) -------------- if(npsta.gt.0)then c --- Coordinates relative to the grid origin xpsta(1)=0.0 ypsta(1)=0.0 ------------------------------------- Subroutine POINTS1 (Problem Area 13) (Screen debug writes for downwash output) ------------------------------------- BEFORE change: (a) -------------- if(mbdw.EQ.2) then write(io6,*) 'hb,hw,index = ',dsbh,dsbw,index write(io6,*) 'len,[x,y]adj = ',dsbl,xadj,yadj else write(io6,*) 'hb,hw,index = ',hb,hw,index endif AFTER change: (a) -------------- if(idownw(i).GT.0) then if(mbdw.EQ.2) then write(io6,*) 'hb,hw,index = ',dsbh,dsbw,index write(io6,*) 'len,[x,y]adj = ',dsbl,xadj,yadj else write(io6,*) 'hb,hw,index = ',hb,hw,index endif endif ------------------------------------- Subroutine POINTS2 (Problem Area 13) (Screen debug writes for downwash output) ------------------------------------- BEFORE change: (a) -------------- if(mbdw.EQ.2) then write(io6,*) 'hb,hw,index = ',dsbh,dsbw,index write(io6,*) 'len,[x,y]adj = ',dsbl,xadj,yadj else write(io6,*) 'hb,hw,index = ',hb,hw,index endif AFTER change: (a) -------------- if(idownw.GT.0) then if(mbdw.EQ.2) then write(io6,*) 'hb,hw,index = ',dsbh,dsbw,index write(io6,*) 'len,[x,y]adj = ',dsbl,xadj,yadj else write(io6,*) 'hb,hw,index = ',hb,hw,index endif endif ------------------------------------- Subroutine QAPLOT1 (Problem Area 13) (Gridded receptor elevations are elevg(i,j)) ------------------------------------- BEFORE change: (a) -------------- write(ioqa,*) x,y,elev(i,j) AFTER change: (a) -------------- write(ioqa,*) x,y,elevg(i,j) ------------------------------------- Subroutine RDEMBC (Problem Area 13) (Use full MXCOL dimension for temp string) ------------------------------------- BEFORE change: (a) -------------- c Parameters: c MXBC, MXSPEC, IO15, IO6 AFTER change: (a) -------------- c Parameters: c MXBC, MXSPEC, IO15, IO6 c MXCOL BEFORE change: (a) -------------- include 'params.puf' AFTER change: (a) -------------- include 'params.puf' include 'params.cal' BEFORE change: (a) -------------- character*1 cstore1(12) AFTER change: (a) -------------- character*1 cstore1(mxcol) ------------------------------------- Subroutine RDHDBC (Problem Area 13) (BC file name is FNAMEBC) ------------------------------------- BEFORE change: (a) -------------- call QCHECK(mxspec,io6,nspec,nse,isplst,cspec,xmol,lemit, 1 fnambc,nspecbc,cspecbc,xmwtbc,ierr,ixrembc) AFTER change: (a) -------------- call QCHECK(mxspec,io6,nspec,nse,isplst,cspec,xmol,lemit, 1 fnamebc,nspecbc,cspecbc,xmwtbc,ierr,ixrembc) ------------------------------------- Subroutine RDHDBC2 (Problem Area 13) (Fix variable names FNAMEBC and I2DRHBC) ------------------------------------- BEFORE change: (a) -------------- i2drh=0 AFTER change: (a) -------------- i2drhBC=0 BEFORE change: (b) -------------- call QCHECK(mxspec,io6,nspec,nse,isplst,cspec,xmol,lemit, 1 fnambc,nspecbc,cspecbc,xmwtbc,ierr,ixrembc) AFTER change: (b) -------------- call QCHECK(mxspec,io6,nspec,nse,isplst,cspec,xmol,lemit, 1 fnamebc,nspecbc,cspecbc,xmwtbc,ierr,ixrembc) ------------------------------------- Subroutine RDMET4 (Problem Area 13) (Declare PTG as 1-D real) ------------------------------------- BEFORE change: (a) -------------- c --- Real 1-D fields real tempss(nssta),rhoss(nssta),qswss(nssta) AFTER change: (a) -------------- c --- Real 1-D fields real tempss(nssta),rhoss(nssta),qswss(nssta) real ptg(2) ------------------------------------- Subroutine RDMET5 (Problem Area 13) (Declare PTG as 1-D real) ------------------------------------- BEFORE change: (a) -------------- c --- Real 1-D fields real tempss(nssta),rhoss(nssta),qswss(nssta) AFTER change: (a) -------------- c --- Real 1-D fields real tempss(nssta),rhoss(nssta),qswss(nssta) real ptg(2) ------------------------------------- Subroutine SLGRECS (Problem Area 13) (Assign 1/slug-length for young end of emitting slugs) ------------------------------------- BEFORE change: (a) -------------- if(iage.EQ.1) then call SLGFRAC(x,y,xb1,yb1,xb2,yb2,rhoai,rhoci,fracsi,d12i,rd12i) else fracsi=fracsf rhoai=rhoaf rhoci=rhocf endif AFTER change: (a) -------------- if(iage.EQ.1) then call SLGFRAC(x,y,xb1,yb1,xb2,yb2,rhoai,rhoci,fracsi,d12i,rd12i) else d12i=d12f rd12i=rd12f fracsi=fracsf rhoai=rhoaf rhoci=rhocf endif BEFORE change: (b) -------------- if(iage.EQ.1) then call SLGFRAC(x,y,xe2,ye2,xb2,yb2,rhoay,rhocy,fracsy,d12y,rd12y) else fracsy = fracso rhoay = rhoao rhocy = rhoco endif AFTER change: (b) -------------- if(iage.EQ.1) then call SLGFRAC(x,y,xe2,ye2,xb2,yb2,rhoay,rhocy,fracsy,d12y,rd12y) else d12y=d12o rd12y=rd12o fracsy = fracso rhoay = rhoao rhocy = rhoco endif ------------------------------------- Subroutine SLUGINT (Problem Area 13) (Undefined variables --- ldbhr changed to ldb) ------------------------------------- BEFORE change: (a) -------------- ccoup = ASDF(ldbhr,dy,sy) * fcaus / (srt2pi * speedi * sy) AFTER change: (a) -------------- ccoup = ASDF(ldb,dy,sy) * fcaus / (srt2pi * speedi * sy) ------------------------------------- Subroutine TFERCF (Problem Area 13) (Reset control-file line-length limit to 200) ------------------------------------- BEFORE change: (a) -------------- character*132 aline, blank AFTER change: (a) -------------- character*132 aline, blank character*200 aline200, blank200 BEFORE change: (b) -------------- blank(100:132)=blank33 AFTER change: (b) -------------- blank200(1:100)=blank(1:100) blank200(101:200)=blank(1:100) BEFORE change: (c) -------------- if(mxcol.GT.132) then write(*,*)'ERROR in TFERCF!' write(*,*)'Control file records are assumed to be <= 132' AFTER change: (c) -------------- if(mxcol.GT.200) then write(*,*)'ERROR in TFERCF!' write(*,*)'Control file records are assumed to be <= 200' BEFORE change: (d) -------------- c --- Transfer control file records 10 aline=blank read(io5,'(a132)',end=999) aline write(iox,'(a132)') aline ncommout=ncommout+1 AFTER change: (d) -------------- 10 aline=blank aline200=blank200 read(io5,'(a200)',end=999) aline200 write(iox,'(a132)') aline200(1:132) ncommout=ncommout+1 klen=LEN_TRIM(aline200) if(klen.GT.132) then c --- Line is split aline(1:68)=aline200(133:200) write(iox,'(a132)') aline ncommout=ncommout+1 endif ------------------------------------- Subroutine WAKE_DBG (Problem Area 13) (Initialize wind speed shear factor outside of wake) ------------------------------------- BEFORE change: (a) -------------- dbuw=1.0 AFTER change: (a) -------------- dbuw=1.0 dbduw=0.0 BEFORE change: (b) -------------- dbuw=1.0 AFTER change: (b) -------------- dbuw=1.0 dbduw=0.0 BEFORE change: (c) -------------- dbuw=1.0 AFTER change: (c) -------------- dbuw=1.0 dbduw=0.0 ----------------------------------------------------------------------- 2. CALUTILS Changes (Version 2.57, Level 090202) ----------------------------------------------------------------------- ---------------------------------- Include File PARAMS.CAL ---------------------------------- BEFORE change: -------------- c --- Specify parameters parameter(mxvar=60,mxcol=132) AFTER change: -------------- c --- Specify parameters parameter(mxvar=60,mxcol=200) ---------------------------------- Subroutine READIN ---------------------------------- BEFORE change: -------------- c --- read a line of input 5 continue read(ioin,10)cstor1 10 format(150a1) if(lecho)write(ioout,7)cstor1 7 format(1x,150a1) AFTER change: -------------- c --- read a line of input 5 continue read(ioin,10)cstor1 10 format(200a1) if(lecho)write(ioout,7)cstor1 7 format(1x,200a1) ---------------------------------- Subroutine DATETM ---------------------------------- BEFORE change: -------------- c --- DATETM calls: DATE, TIME (Lahey F77 compiler utilities) c DATE_AND_TIME (COMPAQ, LF95 compilers) c ETIME (SUN CPU utility program) c YR4C AFTER change: -------------- c --- DATETM calls: DATE, TIME (Lahey F77 compiler utilities) c ETIME (SUN F77 CPU utility program) c DATE_AND_TIME (COMPAQ, LF95 compilers) c CPU_TIME (F95) c YR4C BEFORE change: -------------- c --- Get system date and time [COMPILER-SPECIFIC!] AFTER change: -------------- c --- Set initial base CPU time to -1. data rcpu0/-1./ SAVE rcpu0 c --- Get system date and time [COMPILER-SPECIFIC!] BEFORE change: -------------- cc --- Lahey F77L Compiler (end) AFTER change: -------------- cc --- Get CPU time from SUN system utility (or PC dummy) c call etime(rcpu1) cc --- Lahey F77L Compiler (end) BEFORE change: -------------- c --- COMPAQ DF90/95 and Lahey LF95 Compilers (end) AFTER change: -------------- c --- Get CPU time from LF95 system utility call CPU_TIME(rcpu1) c --- COMPAQ DF90/95 and Lahey LF95 Compilers (end) BEFORE change: -------------- c --- Get CPU time from SUN system utility (or PC dummy) call etime(rcpu) return end AFTER change: -------------- c --- Update base CPU time on first call if(rcpu0.LT.0.0) rcpu0=rcpu1 c --- Return CPU time difference from base rcpu=rcpu1-rcpu0 return end