----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : H Dated : December 14, 2015 ----------------------------------------------------------------------- CALPUFF (From Version 5.8.4 Level 130731 to Version 5.8.5 Level 151214) ----------------------------------------------------------------------- -- Problem Area 1 -- Debug output of local species arrays (e.g., conc()) should use explicit array index since entries after NSPEC are not defined, which halts run if extended compiler checking is on. No effect on results. MODIFIED: CALCBC -- Problem Area 2 -- Initialize error flag IERR in RDHDBC to zero (PC default) No effect on results since run halts if IERR=1 is returned from call to QCHECK, and otherwise proceeds MODIFIED: RDHDBC -- Problem Area 3 -- Fix bug in boundary-source puff sampling that added the contribution of boundary-sources only when the last puff sampled in a time-step was from a boundary-source MODIFIED: COMP -- Problem Area 4 -- Fix bug in computing PRIME downwash wake-modified primary plume sigmas beyond tabulated range for receptors within the distance range of the cavity. There was an inconsistent mix of travel time and virtual-source increment for this. Sigmas for locations beyond the cavity are not affected. NEW : WAKE_CSIG MODIFIED: WAKEDAT.PUF POINTS1, POINTS2, WAKE_XSIG -- Problem Area 5 -- Fix bug in setting the high and low limits on the range of receptors that are within the range of influence of a SLUG during a sampling step. A typo in screening for the impact zone about the youngest end of the slug at the end of the step used its position at the start of the step with its spread at the end of the step. MODIFIED: SLGXLIM -- Problem Area 6 -- Rename loop index in debug output that conflicted with calling argument. Could cause model to stop if debug output is selected with more than 1 variable-emissions file for point sources. Contributed by Peter Rye. MODIFIED: RDTIEM2 -- Problem Area 7 -- Add BID to sigmas for the PRIME case where plume does not enter the wake (no wake table). Plume growth is calculated from the source rather than from the sigmas at the end of the wake (read from the wake table, including BID), so the BID due to rise must be added. Before this fix the BID sigmas were set to zero like the plumes that enter the wake, and this can over-estimate concentrations. Earlier CALPUFF versions did not include BID in the wake tables, so BID was always added in subsequent calculations including this case without a wake table. MODIFIED: POINTS1, POINTS2 -- Problem Area 8 -- Add check for Gaussian vertical distribution before calling RISEWIND to refine advection wind (may cause halt in ADVECT subroutine) MODIFIED: COMP ----------------------------------------------------------------------- CALMET (From Version 5.8.4 Level 130731 to Version 5.8.5 Level 151214) ----------------------------------------------------------------------- --- Problem Area 1 --- Fix lines that exceed FORTRAN fixed-format limit (72). These lines are active only when overwater convective mixing heights are simulated with the Batchvarova-Gryning option. Modified: WATER, WATERP --- Problem Area 2 --- Fix bug in setting the S./N. hemisphere when extrapolating the wind direction in the vertical above the surface using the similarity extrapolation option (IEXTRP=4 or -4) for surface-station data. This may produce a wind rotation error when using UTM coordinates if the latitude variable XLAT1 is not consistent with the station location. Modified: SIMILT ----------------------------------------------------------------------- 1. CALPUFF Changes ----------------------------------------------------------------------- *********************************************************************** --- Problem Area 1 --- Debug output of local species arrays (e.g., conc()) should use explicit array index since entries after NSPEC are not defined. *********************************************************************** --------------------------------- Subroutine CALCBC (Problem Area 1) --------------------------------- Ensure that explicit array index are used for conc(), concup(), dflx() and wflx() BEFORE change: -------------- if(LDBHR) then write(io6,*) 'CALCBC: conc = ',conc write(io6,*) ' concup = ',concup if(mdry.EQ.1) write(io6,*) ' dflx = ',dflx if(mwet.EQ.1) write(io6,*) ' wflx = ',wflx endif AFTER change: ------------- if(LDBHR) then write(io6,*) 'CALCBC: conc = ',(conc(i),i=1,nspec) write(io6,*) ' concup = ',(concup(i),i=1,nspec) if(mdry.EQ.1) write(io6,*) ' dflx = ',(dflx(i), & i=1,nspec) if(mwet.EQ.1) write(io6,*) ' wflx = ',(wflx(i), & i=1,nspec) endif *********************************************************************** --- Problem Area 2 --- Initialize error flag IERR in RDHDBC to zero. *********************************************************************** --------------------------------- Subroutine RDHDBC (Problem Area 2) --------------------------------- BEFORE change: -------------- c --- QA checks c ------------- AFTER change: ------------- c --- QA checks c ------------- ierr=0 *********************************************************************** --- Problem Area 3 --- Fix bug in boundary-source puff sampling that added the contribution of boundary-sources only when the last puff sampled in a time-step was from a boundary-source *********************************************************************** --------------------------------- Subroutine COMP (Problem Area 3) --------------------------------- BEFORE change: -------------- c --- Sum contribution from the boundary condition puff that passes c --- nearest each receptor this sampling step if(LBCPUF) call SUMBC AFTER change: ------------- c --- Sum contribution from the boundary condition puff that passes c --- nearest each receptor this sampling step if(nbc.GT.0) call SUMBC *********************************************************************** --- Problem Area 4 --- Fix bug in computing PRIME downwash wake-modified primary plume sigmas beyond tabulated range for receptors within the distance range of the cavity. There was an inconsistent mix of travel time and virtual-source increment for this. MODIFIED: WAKEDAT.PUF POINTS1, POINTS2, WAKE_XSIG *********************************************************************** --------------------------------- New Subroutine WAKE_CSIG (Problem Area 4) --------------------------------- BEFORE change: -------------- return end AFTER change: ------------- return end c---------------------------------------------------------------------- subroutine wake_csig(ity,itz,irur,utrans,kstab,elmo,bvfreq, & sv,sw,symin,szmin,z,zi,amode) c---------------------------------------------------------------------- c c --- CALPUFF Version: 5.8.5 Level: 151214 WAKE_CSIG c c --- PURPOSE: Save all arguments needed in SETCSIG call for release c height dispersion or for building-top dispersion c c --- INPUTS: c ITY - integer - Sigma-y method c ITZ - integer - Sigma-z method c IRUR - integer - Rural flag (rural=0 ; urban=1) c UTRANS - real - Transport wind speed (m/s) c KSTAB - integer - PGT stability class c ELMO - real - Monin-Obukhov length (m) c BVFREQ - real - Brunt-Vaisala freq (1/s) c SV - real - Sigma-v (m/s) c SW - real - Sigma-w (m/s) c SYMIN - real - Sigma-y floor (m) c SZMIN - real - Sigma-z floor (m) c Z - real - Height (m AGL) c ZI - real - Mixing height (m) c AMODE - char*8 - 'RELEASE ' for release-height data c 'BUILDING' for building-top data c c --- OUTPUT: c common/WAKECSIG/idopty_bt,idoptz_bt,iru_bt,uavg_bt,kst_bt, c & el_bt,bvf_bt,tsigv_bt,tsigw_bt,symin_bt, c & szmin_bt,zht_bt,zmix_bt, c & idopty_rh,idoptz_rh,iru_rh,uavg_rh,kst_rh, c & el_rh,bvf_rh,tsigv_rh,tsigw_rh,symin_rh, c & szmin_rh,zht_rh,zmix_rh c c --- WAKE_CSIG called by: POINTS1, POINTS2 c --- WAKE_CSIG calls: (none) c---------------------------------------------------------------------- c c --- Include common blocks include 'params.puf' include 'wakedat.puf' character*8 amode if(amode.EQ.'RELEASE ') then c --- Assign input variables to release-ht names in /WAKECSIG/ idopty_rh=ity idoptz_rh=itz iru_rh =irur uavg_rh =utrans kst_rh =kstab el_rh =elmo bvf_rh =bvfreq tsigv_rh =sv tsigw_rh =sw symin_rh =symin szmin_rh =szmin zht_rh =z zmix_rh =zi elseif(amode.EQ.'BUILDING') then c --- Assign input variables to building-top names in /WAKECSIG/ idopty_bt=ity idoptz_bt=itz iru_bt =irur uavg_bt =utrans kst_bt =kstab el_bt =elmo bvf_bt =bvfreq tsigv_bt =sv tsigw_bt =sw symin_bt =symin szmin_bt =szmin zht_bt =z zmix_bt =zi else write(*,*)'ERROR in WAKE_CSIG: invalid mode' write(*,*)' Expected RELEASE or BUILDING' write(*,*)' Found ',amode stop endif return end --------------------------------- Subroutine POINTS1 (Problem Area 4) --------------------------------- BEFORE change: -------------- c -------------------------------- c --- PRIME Downwash Section c -------------------------------- mprime=1 AFTER change: ------------- c -------------------------------- c --- PRIME Downwash Section c -------------------------------- mprime=1 c --- Store dispersion data for release height for later use zht=heff call WAKE_CSIG(idopty,idoptz,iru,ws,istab,el, & sqrts,tsigv,tsigw,symin,szmin, & zht,dpbl,'RELEASE ') BEFORE change: -------------- call SETCSIG(idoptyb,idoptzb,iru,wsb,istabb,elb,sqrtsb, & tsigvb,tsigwb,symin,szmin,zht,dpbl) AFTER change: ------------- call SETCSIG(idoptyb,idoptzb,iru,wsb,istabb,elb,sqrtsb, & tsigvb,tsigwb,symin,szmin,zht,dpbl) c --- Store dispersion data for top of building for later use call WAKE_CSIG(idoptyb,idoptzb,iru,wsb,istabb,elb, & sqrtsb,tsigvb,tsigwb,symin,szmin, & zht,dpbl,'BUILDING') --------------------------------- Subroutine POINTS2 (Problem Area 4) --------------------------------- BEFORE change: -------------- c -------------------------------- c --- PRIME Downwash Section c -------------------------------- mprime=1 AFTER change: ------------- c -------------------------------- c --- PRIME Downwash Section c -------------------------------- mprime=1 c --- Store dispersion data for release height for later use zht=heff call WAKE_CSIG(idopty,idoptz,iru,ws,istab,el, & sqrts,tsigv,tsigw,symin,szmin, & zht,dpbl,'RELEASE ') BEFORE change: -------------- call SETCSIG(idoptyb,idoptzb,iru,wsb,istabb,elb,sqrtsb, & tsigvb,tsigwb,symin,szmin,zht,dpbl) AFTER change: ------------- call SETCSIG(idoptyb,idoptzb,iru,wsb,istabb,elb,sqrtsb, & tsigvb,tsigwb,symin,szmin,zht,dpbl) c --- Store dispersion data for top of building for later use call WAKE_CSIG(idoptyb,idoptzb,iru,wsb,istabb,elb, & sqrtsb,tsigvb,tsigwb,symin,szmin, & zht,dpbl,'BUILDING') --------------------------------- Subroutine WAKE_XSIG (Problem Area 4) --------------------------------- BEFORE change: -------------- logical NOBID AFTER change: ------------- logical NOBID logical LDB LDB=.FALSE. BEFORE change: -------------- c --- Set total effective distance(km), time(s) xzkm=xkm+xzvkm xykm=xkm+xyvkm tz=t+tzv ty=t+tyv AFTER change: ------------- c --- Set total effective distance(km), time(s) xzkm=xkm+xzvkm xykm=xkm+xyvkm tz=t+tzv ty=t+tyv c --- These will be replaced below with virtuals from c --- last sigmas and release-ht data if needed --- c --- Downwind distance from last tabulated value delx=x-xwak(nwak) delxkm=x*fm2km c --- Time Increment delt=delx/uavg_rh BEFORE change: -------------- elseif(nwak.LE.1) then c --- Plume never altered by wake turbulence: c --- use HOST curves with BID adjustment call SIGTZ(zero,xzkm,tz,Hb,sz,tzvdum,xzvdum) call SIGTY(zero,xykm,ty,sy,tyvdum,xyvdum) if(.not.NOBID) then sz=SQRT(sz**2+bidsq) sy=SQRT(sy**2+bidsq) endif elseif(x.gt.xwak(nwak)) then c --- Point lies downwind of transition to ambient growth: c --- use HOST curves with the virtual offset (includes BID) call SIGTZ(zero,xzkm,tz,Hb,sz,tzvdum,xzvdum) call SIGTY(zero,xykm,ty,sy,tyvdum,xyvdum) AFTER change: ------------- elseif(nwak.LE.1) then c --- Plume never altered by wake turbulence: c --- use HOST curves with BID adjustment c --- Set selected data in /CSIGMA/ for RELEASE-HT sigma calls if(LDB) then write(io6,*)'WAKE_XSIG: nwak.LE.1, dx,dt = ',delx,delt endif call SETCSIG(idopty_rh,idoptz_rh,iru_rh,uavg_rh,kst_rh, & el_rh,bvf_rh,tsigv_rh,tsigw_rh,symin_rh, & szmin_rh,zht_rh,zmix_rh) c --- Grow sigmas from end of table with release-ht data call SIGTZ(szwak(nwak),delxkm,delt,zht_rh,sz,tzvdum,xzvdum) call SIGTY(sywak(nwak),delxkm,delt,sy,tyvdum,xyvdum) c --- Restore selected data in /CSIGMA/ for BUILDING-TOP sigma calls call SETCSIG(idopty_bt,idoptz_bt,iru_bt,uavg_bt,kst_bt, & el_bt,bvf_bt,tsigv_bt,tsigw_bt,symin_bt, & szmin_bt,zht_bt,zmix_bt) if(.not.NOBID) then sz=SQRT(sz**2+bidsq) sy=SQRT(sy**2+bidsq) endif elseif(x.gt.xwak(nwak)) then c --- Point lies downwind of transition to ambient growth: c --- Set selected data in /CSIGMA/ for RELEASE-HT sigma calls if(LDB) then write(io6,*)'WAKE_XSIG: x>xwak(nwak), dx,dt = ',delx,delt endif call SETCSIG(idopty_rh,idoptz_rh,iru_rh,uavg_rh,kst_rh, & el_rh,bvf_rh,tsigv_rh,tsigw_rh,symin_rh, & szmin_rh,zht_rh,zmix_rh) c --- Grow sigmas from end of table with release-ht data c --- Tabulated sigmas include BID call SIGTZ(szwak(nwak),delxkm,delt,zht_rh,sz,tzvdum,xzvdum) call SIGTY(sywak(nwak),delxkm,delt,sy,tyvdum,xyvdum) c --- Restore selected data in /CSIGMA/ for BUILDING-TOP sigma calls call SETCSIG(idopty_bt,idoptz_bt,iru_bt,uavg_bt,kst_bt, & el_bt,bvf_bt,tsigv_bt,tsigw_bt,symin_bt, & szmin_bt,zht_bt,zmix_bt) --------------------------------- Included File WAKEDAT.PUF (Problem Area 7) --------------------------------- BEFORE change: -------------- & xzvcav,xyvcav,fqcav,cbyqcav,istab,lrurl c AFTER change: ------------- & xzvcav,xyvcav,fqcav,cbyqcav,istab,lrurl c common/WAKECSIG/idopty_bt,idoptz_bt,iru_bt,uavg_bt,kst_bt, & el_bt,bvf_bt,tsigv_bt,tsigw_bt,symin_bt, & szmin_bt,zht_bt,zmix_bt, & idopty_rh,idoptz_rh,iru_rh,uavg_rh,kst_rh, & el_rh,bvf_rh,tsigv_rh,tsigw_rh,symin_rh, & szmin_rh,zht_rh,zmix_rh c *********************************************************************** --- Problem Area 5 --- Fix bug in setting the high and low limits on the range of receptors that are within the range of influence of a SLUG during a sampling step. A typo in screening for the impact zone about the youngest end of the slug at the end of the step used its position at the start of the step with its spread at the end of the step. MODIFIED: SLGXLIM *********************************************************************** --------------------------------- Subroutine SLGXLIM (Problem Area 5) --------------------------------- Update so the comparison is referencing XE2 and YE2 to represent the position at the end of the time step instead of XB2 and YB2 which are the position at the beginning BEFORE change: -------------- C --- NOW CONSIDER THE NEWEST END (2) OF THE SLUG AT THE END C OF THE TIME STEP. NOW SKIP EXPLICIT DEFINITION OF XTRY,YTRY C TO SAVE INSTRUCTIONS. SIGTRY = SIGCUT * SYE2 * tfaccl IF((XB2+SIGTRY) .GT. XHGH) XHGH = XE2+SIGTRY IF((XB2-SIGTRY) .LT. XLOW) XLOW = XE2-SIGTRY IF((YB2+SIGTRY) .GT. YHGH) YHGH = YE2+SIGTRY IF((YB2-SIGTRY) .LT. YLOW) YLOW = YE2-SIGTRY AFTER change: ------------- C --- NOW CONSIDER THE NEWEST END (2) OF THE SLUG AT THE END C OF THE TIME STEP. NOW SKIP EXPLICIT DEFINITION OF XTRY,YTRY C TO SAVE INSTRUCTIONS. SIGTRY = SIGCUT * SYE2 * tfaccl IF((XE2+SIGTRY) .GT. XHGH) XHGH = XE2+SIGTRY IF((XE2-SIGTRY) .LT. XLOW) XLOW = XE2-SIGTRY IF((YE2+SIGTRY) .GT. YHGH) YHGH = YE2+SIGTRY IF((YE2-SIGTRY) .LT. YLOW) YLOW = YE2-SIGTRY *********************************************************************** -- Problem Area 6 -- Rename loop index in debug output that conflicted with calling argument. MODIFIED: RDTIEM2 *********************************************************************** --------------------------------- Subroutine RDTIEM2 (Problem Area 6) --------------------------------- Change index name from 'i' to 'k' to avoid a conflict. BEFORE change: -------------- if(ldb)then write(io6,*) write(io6,*)'Time-invariant PTEMARB data' do 150 i=1,npt2 write(io6,*)'I = ',i,'CID2 = ',cid2(i) write(io6,*)'I = ',i,'TIEM2 = ',(tiem2(n,i),n=1,8) if(NINT(tiem2(6,i)).EQ.2) then write(io6,*)'I = ',i,'ZPLAT2 = ',zplat2(i) endif if(NINT(tiem2(6,i)).GE.1) then write(io6,*)'--- BHT2 = ',(bht2(n,i),n=1,36) write(io6,*)'--- BWD2 = ',(bwd2(n,i),n=1,36) if(mbdw.EQ.2) then write(io6,*)'--- BLN2 = ',(bln2(n,i),n=1,36) write(io6,*)'--- XBADJ2 = ',(xbadj2(n,i),n=1,36) write(io6,*)'--- YBADJ2 = ',(ybadj2(n,i),n=1,36) endif endif 150 continue endif AFTER change: ------------- if(ldb)then write(io6,*) write(io6,*)'Time-invariant PTEMARB data' do 150 k=1,npt2 write(io6,*)'Source = ',k,'CID2 = ',cid2(k) write(io6,*)'Source = ',k,'TIEM2 = ',(tiem2(n,k),n=1,8) if(NINT(tiem2(6,k)).EQ.2) then write(io6,*)'Source = ',k,'ZPLAT2 = ',zplat2(k) endif if(NINT(tiem2(6,k)).GE.1) then write(io6,*)'--- BHT2 = ',(bht2(n,k),n=1,36) write(io6,*)'--- BWD2 = ',(bwd2(n,k),n=1,36) if(mbdw.EQ.2) then write(io6,*)'--- BLN2 = ',(bln2(n,k),n=1,36) write(io6,*)'--- XBADJ2 = ',(xbadj2(n,k),n=1,36) write(io6,*)'--- YBADJ2 = ',(ybadj2(n,k),n=1,36) endif endif 150 continue endif *********************************************************************** -- Problem Area 7 -- Add BID to sigmas for the PRIME case where plume does not enter the wake (no wake table). *********************************************************************** --------------------------------- Subroutine POINTS1 (Problem Area 7) --------------------------------- BEFORE change: -------------- c --- Reset BID and tip downwash tipdw=0.0 c --- The BID has been moved into the wake tables bidsq=0.0 AFTER change: ------------- c --- Reset BID and tip downwash tipdw=0.0 if(nwk.LE.1 .AND. 1prm) then c Case where primary plume is not in wake (no wake table) bidsq=(rise/3.5)**2 else c --- The BID has been included in the wake tables bidsq=0.0 endif --------------------------------- Subroutine POINTS2 (Problem Area 7) --------------------------------- BEFORE change: -------------- c --- Reset BID and tip downwash tipdw=0.0 c --- The BID has been moved into the wake tables bidsq=0.0 AFTER change: ------------- c --- Reset BID and tip downwash tipdw=0.0 if(nwk.LE.1 .AND. 1prm) then c Case where primary plume is not in wake (no wake table) bidsq=(rise/3.5)**2 else c --- The BID has been included in the wake tables bidsq=0.0 endif *********************************************************************** -- Problem Area 8 -- Add check for Gaussian vertical distribution *********************************************************************** --------------------------------- Subroutine COMP (Problem Area 8) --------------------------------- BEFORE change: -------------- c --- Account for wind profile on transport during rise c --- (applies to CALMET or PROFILE winds only) if(xtotb(ii).LT.xfinal(ii)) then AFTER change: ------------- c --- Account for wind profile on transport during rise c --- (applies to CALMET or PROFILE winds only) if((MOD(icode,2).EQ.1).AND. & (xtotb(ii).LT.xfinal(ii))) then ----------------------------------------------------------------------- 2. CALMET Changes ----------------------------------------------------------------------- *********************************************************************** -- Problem Area 1 -- Fix lines that exceed FORTRAN fixed-format limit (72). *********************************************************************** --------------------------------- Subroutine WATER (Problem Area 1) --------------------------------- BEFORE change: -------------- if(imixh.eq.2) then c Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),Ilapse, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZICOLD(i,j),ZICONV, : THT,THTP) AFTER change: -------------- if(imixh.eq.2) then c Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),Ilapse, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZICOLD(i,j), : ZICONV,THT,THTP) --------------------------------- Subroutine WATERP (Problem Area 1) --------------------------------- BEFORE change: -------------- if(imixh.eq.2) then c Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),Ilapse, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZICOLD(i,j),ZICONV, : THT,THTP) AFTER change: -------------- if(imixh.eq.2) then c Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),Ilapse, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZICOLD(i,j), : ZICONV,THT,THTP) *********************************************************************** -- Problem Area 2 -- Fix bug in setting the S./N. hemisphere when extrapolating the wind direction in the vertical above the surface using the similarity extrapolation option (IEXTRP=4 or -4) for surface-station data. This may produce a wind rotation error when using UTM coordinates if the latitude variable XLAT1 is not consistent with the station location. *********************************************************************** --------------------------------- Subroutine SIMILT (Problem Area 2) --------------------------------- BEFORE change: -------------- include 'map.met' AFTER change: ------------- include 'met1.met' include 'ovrwat.met' BEFORE change: -------------- c c --- Restrict Monin-Obukhov length to be at least 5*z0 AFTER change: -------------- c --- Select latitude of the station location if(ns.LE.nssta) then idstation=ns stnlat=xslat(idstation) else idstation=ns-nssta if(idstation.GE.1 .AND. idstation.LE.nowsta) then stnlat=xowlat(idstation) else write(io6,*)'ERROR in SIMILT: invalid NS' write(io6,*)' Expected NS = 1 to NSSTA+NOWSTA' write(io6,*)' Found NS,NSSTA,NOWSTA = ',ns,nssta,nowsta stop 'ERROR in SIMILT -- see list file' endif endif c c --- Restrict Monin-Obukhov length to be at least 5*z0