----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : D Dated : JUNE 23, 2007 ----------------------------------------------------------------------- The following code modifications address problems identified in CALMET and CALPUFF. Bulletin MCB-D(070623) outlines these changes and together with those in MCB-A(040716), MCB-B(051216), and MCB-C(060804) represent all of the corrections made to these programs since EPA first approved CALPUFF for regulatory applications. All of these updates are present in CALMET version 5.8, Level 070623 and CALPUFF Version 5.8, Level 070623. Changes to CALPUFF ------------------ -- Problem Area 1 -- Pass CALM processing result (LCALM) from SETPUF/SETSLG for use later in calls to PUFRECS/SLGRECS, and turn off building downwash and transitional rise for the current puff in a calm. Previously these were turned off in the receptor-sampling routines when in a calm, and may not be triggered if receptors are not impacted in the current step. This should be done as soon as such a puff is advected into a calm region. Modified: COMP, SETPUF, SETSLG, PUFRECS, SLGRECS -- Problem Area 2 -- Initialize the following variables (results are not affected but execution may halt with some compilers) IGRDVL = 0 (READCF) HB, HW, HEFF2, ZLY, RINIT = 0. (POINTS1, POINTS2) YKDUM, SYDUM, ZKDUM, SZDUM = 0.0 (WAKE_DFSN) DBDUW = 0.0 (WAKE_DBG) XMAX, FZMAX = 0.0 (WAKE_FQC) DZDS=0.0 (NUMRISE) Modified: READCF, POINTS1, POINTS2, WAKE_DFSN, WAKE_DBG, WAKE_FQC, NUMRISE -- Problem Area 3 -- Assign PTGRAD value to default so it is always available for debug output even when not used (does not change results) Modified: NUMMET -- Problem Area 4 -- PRIME module implementation: Typo in argument IRU (0=rural, 1=urban) in call to WAKE_INI results in rural dispersion curves in PRIME calculations when making the wake-zone sigma tables (PG dispersion option) if the compiler initializes variables to zero. Without such initialization, run may halt with invalid IRU value. Modified: POINTS1, POINTS2 Local receptor height ZRC for PRIME cavity-source impact was not set in the wake downwind of the cavity. Applications that place receptors on the ground are not affected. Initialize ZRC to the actual receptor height above the ground. Modified: CAV_SAMP Cavity concentrations were not computed at receptors in cavity if they were upwind of the source. Logic updated to account for mass recirculating in cavity. Modified: CAV_CONC Puff virtual times at the time of release obtained for PRIME downwash during very strong wind speed shear in the vertical can lead to much larger or much smaller sigmas than those tabulated within the wake zone, when time-based dispersion rates are selected. Although the tabulated PRIME sigmas are used for concentrations at receptors in the wake zone, the lateral sampling screen uses the virtual-time sigmas which can lead to screening out receptors in the wake zone that should be sampled. Modified: SETPUF Changes to CALMET ------------------ -- Problem Area 1 -- Cloud cover was not updated under stable conditions for icloud=3 Modified: ELUSTR2 -- Problem Area 2 -- The non-initialization of RHOP in RDMM5 could produce wrong results with some compilers (e.g. PGI in Linux) Modified: RDMM5 -- Problem Area 3 -- The debug-print flag (ldbhr) was missing in MIXHTST and MIXHTST2, potentially causing large output files to be generated when not requested (pb on some platforms, e.g. PGI compiler in Linux) or no output to be written out when requested Additionally nx,ny were not passed on to those subroutines, potentially compounding this large debug-printing problem. Modified: COMP, ELUSTR2,MIXHTST, MIXHT2ST, OUTPT.MET -- Problem Area 4 -- For IEXTRP=1 (no vertical extrapolation), vertical extrapolation using similarity theory was performed Modified: DIAGNO -- Problem Area 5 -- Some variables were not initialized, with no consequence on the results per se, but with some compilers/platforms, un-initialized variables may stop execution Modified: PREPDI, RDMM4,RDMM5 -- Problem Area 6 -- In NOOBS-temperature mode (ITPROG>0), temperatures were adjusted to an adiabatic profile up to an undefined height (rather than up to the mixing height) Modified: TEMP3D -- Problem Area 7 -- Undefined variable ws10 may cause execution stop with some compilers/platforms if calm wind conditions are not extrapolated (icalm=0) (but this did not affect results if execution went forward) Modified: WIND1 -- Problem Area 8 -- For surface stations located on a water gridcell, the constant GEO.DAT roughness length (z0) was used rather than the Hosker z0 (function of the surface wind speed, as should be overwater) for the vertical extrapolation of surface winds with similarity theory (abs(IEXTRP)=4) Modified: WIND1 -- Problem Area 9 -- The 10m wind speed overwater was not defined when the anemometer height was exactly 10m (problem for OCD method) Modified: WATER2 -- Problem Area 10 -- Undefined index variable when skipping records might stop execution with some compilers/platforms (but not affect results) Modified: RDMM5 -- Problem Area 11 -- Prognostic Land Use, which determines how prognostic temperatures are extrapolated from the lowest M32D level to lowest CALMET level, was not read in from old MM5 format Modified: RDHD51, RDHD52 ----------------------------------------------------------------------- ----------------------------------------------------------------------- 1. CALPUFF ----------------------------------------------------------------------- -------------------------------- Update subroutine PUFRECS -------------------------------- BEFORE change: -------------- c --- If CALM period encountered: reset dispersion option for this step c --- ("calms" use time-based sigmas) and selected puff properties to c --- remove gradual rise/downwash calculations in all subsequent steps if(lcalm) then idopty=1 idoptz=1 frac=0.5 xfrise=0.0 idw0(ipnum)=0 xfinal(ipnum)=0.0 endif AFTER change: -------------- c --- If CALM period encountered: reset dispersion option for this step c --- ("calms" use time-based sigmas) if(lcalm) then idopty=1 idoptz=1 frac=0.5 xfrise=0.0 c 070623 (MCB-D) c idw0(ipnum)=0 c xfinal(ipnum)=0.0 c --- Test: Downwash should aready be off for a puff in a calm if(idw0(ipnum).GT.0) then write(io6,*)'PUFRECS: Puff with active downwash found in' write(io6,*)' a calm environment. This should have been' write(io6,*)' trapped in SETPUF!' stop 'Halted in PUFRECS -- See list file' endif endif -------------------------------- Update subroutine SLGRECS -------------------------------- BEFORE change: -------------- c --- If CALM period encountered: reset dispersion option for this step c --- ("calms" use time-based sigmas) and selected slug properties to c --- remove gradual rise/downwash calculations in all subsequent steps if(lcalm) then idopty=1 idoptz=1 xfrise=zero idw0(ipnum)=0 xfinal(ipnum)=zero endif AFTER change: -------------- c --- If CALM period encountered: reset dispersion option for this step c --- ("calms" use time-based sigmas) if(lcalm) then idopty=1 idoptz=1 xfrise=zero c 070623 (MCB-D) c idw0(ipnum)=0 c xfinal(ipnum)=zero if(idw0(ipnum).GT.0) then write(io6,*)'SLGRECS: Puff with active downwash found in' write(io6,*)' a calm environment. This should have been' write(io6,*)' trapped in SETSLG!' stop 'Halted in SLGRECS -- See list file' endif endif -------------------------------- Update subroutine COMP -------------------------------- BEFORE change: -------------- elseif(lpuff) then c c --- NOW PUFFS c --- Compute puff position and sigmas, and update stored data. call SETPUF(ii,tsamp,distm,jdstab,el,tsigv,tsigw,bvf,ws, & xnew,ynew,iru,ldbhr,uavg) c else c c --- NOW SLUGS c --- Compute position of slug-ends at start and end of step, c --- compute sigmas at these points, and update stored data. call SETSLG(ii,ifree,tsamp,iru,jdstab,el,tsigv,tsigw,bvf, & ws,dxoldm,dyoldm,dxnewm,dynewm,ldbhr,uavg) endif AFTER change: -------------- elseif(lpuff) then c c --- NOW PUFFS c --- Compute puff position and sigmas, and update stored data. c 070623 (MCB-D) c call SETPUF(ii,tsamp,distm,jdstab,el,tsigv,tsigw,bvf,ws, c & xnew,ynew,iru,ldbhr,uavg) call SETPUF(ii,tsamp,distm,jdstab,el,tsigv,tsigw,bvf,ws, & xnew,ynew,iru,ldbhr,uavg,lcalm) c else c c --- NOW SLUGS c --- Compute position of slug-ends at start and end of step, c --- compute sigmas at these points, and update stored data. c 070623 (MCB-D) c call SETSLG(ii,ifree,tsamp,iru,jdstab,el,tsigv,tsigw,bvf, c & ws,dxoldm,dyoldm,dxnewm,dynewm,ldbhr,uavg) call SETSLG(ii,ifree,tsamp,iru,jdstab,el,tsigv,tsigw,bvf, & ws,dxoldm,dyoldm,dxnewm,dynewm,ldbhr, & uavg,lcalm) endif -------------------------------- Update subroutine POINTS1 -------------------------------- BEFORE change: -------------- c --- Set minimum allowed wind speed for building downwash data wsdw0/1.0/ AFTER change: -------------- c --- Set minimum allowed wind speed for building downwash data wsdw0/1.0/ c 070623 (MCB-D) c --- Initialize local ISC building downwash variables hb=0.0 hw=0.0 heff2=0.0 zly=0.0 rinit=0.0 BEFORE change: -------------- c --- Refresh /WAKEDAT/ variables call WAKE_INI(ldbhr,istab,irub,dsbh,dsbw,dsbl,xadj,yadj, & wsb,ws,tsigvb,tsigwb,idoptyb,idoptzb) AFTER change: -------------- c --- Refresh /WAKEDAT/ variables c 070623 (MCB-D) c call WAKE_INI(ldbhr,istab,irub,dsbh,dsbw,dsbl,xadj,yadj, c & wsb,ws,tsigvb,tsigwb,idoptyb,idoptzb) call WAKE_INI(ldbhr,istab,iru,dsbh,dsbw,dsbl,xadj,yadj, & wsb,ws,tsigvb,tsigwb,idoptyb,idoptzb) BEFORE change: -------------- write(io6,208)j,np,idopty,idoptz,iru,sysrc,szsrc,dt,istab 208 format(5x,'J: ',i5,2x,'NP: ',i5,2x,'IDOPTY: ',i2,2x, 1 'IDOPTZ: ',i2,2x,'IRU: ',i2,/2x,'SIGYB: ',f8.2,2x, 2 'SIGZB: ',f8.2,2x,'DT: ',f7.2,2x,'ISTAB: ',i1) AFTER change: -------------- c 070623 (MCB-D) write(io6,208)i,np,idopty,idoptz,iru,sysrc,szsrc,dt,istab 208 format(5x,'I: ',i5,2x,'NP: ',i5,2x,'IDOPTY: ',i2,2x, 1 'IDOPTZ: ',i2,2x,'IRU: ',i2,/2x,'SIGYB: ',f8.2,2x, 2 'SIGZB: ',f8.2,2x,'DT: ',f7.2,2x,'ISTAB: ',i1) -------------------------------- Update subroutine POINTS2 -------------------------------- BEFORE change: -------------- c --- Set minimum allowed wind speed for building downwash data wsdw0/1.0/ AFTER change: -------------- c --- Set minimum allowed wind speed for building downwash data wsdw0/1.0/ c 070623 (MCB-D) c --- Initialize local ISC building downwash variables hb=0.0 hw=0.0 heff2=0.0 zly=0.0 rinit=0.0 BEFORE change: -------------- c --- Refresh /WAKEDAT/ variables call WAKE_INI(ldbhr,istab,irub,dsbh,dsbw,dsbl,xadj,yadj, & wsb,ws,tsigvb,tsigwb,idoptyb,idoptzb) AFTER change: -------------- c --- Refresh /WAKEDAT/ variables c 070623 (MCB-D) c call WAKE_INI(ldbhr,istab,irub,dsbh,dsbw,dsbl,xadj,yadj, c & wsb,ws,tsigvb,tsigwb,idoptyb,idoptzb) call WAKE_INI(ldbhr,istab,iru,dsbh,dsbw,dsbl,xadj,yadj, & wsb,ws,tsigvb,tsigwb,idoptyb,idoptzb) BEFORE change: -------------- write(io6,208)j,np,iru,sysrc,szsrc,dt,istab 208 format(5x,'J: ',i5,2x,'NP: ',i5,2x, 1 'IRU: ',i2,2x,'SIGYB: ',f8.2,2x,'SIGZB: ',f8.2,2x/ 2 5x,'DT: ',f7.2,2x,'ISTAB: ',i1) AFTER change: -------------- c 070623 (MCB-D) write(io6,208)i,np,iru,sysrc,szsrc,dt,istab 208 format(5x,'I: ',i5,2x,'NP: ',i5,2x, 1 'IRU: ',i2,2x,'SIGYB: ',f8.2,2x,'SIGZB: ',f8.2,2x/ 2 5x,'DT: ',f7.2,2x,'ISTAB: ',i1) -------------------------------- Update subroutine NUMRISE -------------------------------- BEFORE change: -------------- c --- Obtain mean streamline slopes here (within 15R of building) dxds=0.0 dzdx=0.0 dzstrm=0.0 AFTER change: -------------- c --- Obtain mean streamline slopes here (within 15R of building) dxds=0.0 dzdx=0.0 c 070623 (MCB-D) dzds=0.0 dzstrm=0.0 -------------------------------- Update subroutine NUMMET -------------------------------- BEFORE change: -------------- data zero/0.0/, half/0.5/ AFTER change: -------------- data zero/0.0/, half/0.5/ c 070623 (MCB-D) c --- Set default PTemp gradient (may not be used) kst=ipgt(ixs,iys) if(kst.LE.4) then ptgrad=ptgrad0 elseif(kst.EQ.5) then ptgrad=ptg(1) else ptgrad=ptg(2) endif c --- Impose minimum value if(ptgrad.LT.ptgrad0)ptgrad=ptgrad0 -------------------------------- Update subroutine READCF -------------------------------- BEFORE change: -------------- c --- Set "problem" marker to FALSE data problem/.FALSE./ AFTER change: -------------- c --- Set "problem" marker to FALSE data problem/.FALSE./ c 070623 (MCB-D) c --- Initialize OLD number of VOLEM volume sources to ZERO igrdvl=0 -------------------------------- Update subroutine SETPUF -------------------------------- BEFORE change: -------------- c --- Set dispersion option ("calms" use time-based sigmas) lcalm=.FALSE. if(ws.LT.wscalm) then lcalm=.TRUE. idopty=1 idoptz=1 endif AFTER change: -------------- c --- Set dispersion option ("calms" use time-based sigmas) lcalm=.FALSE. if(ws.LT.wscalm) then lcalm=.TRUE. idopty=1 idoptz=1 c 070623 (MCB-D) c --- Turn off any downwash and plume rise calculations if CALM idw0(ii)=0 xfinal(ii)=zero endif BEFORE change: -------------- c --- Try to process sigmas at start of step call WAKE_TAB(istype,isnum,idw0(ipnum),xttb1, & syb1,szb1,hgr,rise,xlast,iprime) if(iprime.EQ.1) then c --- OK, puff starts in wake region; get virtuals call WAKE_TAB(istype,isnum,idw0(ipnum),xlast, & sylast,szlast,hgr,rise,xlast2,idone) call SIGTY(sylast,zero,zero,sydum,vtylast,vdylast) call SIGTZ(szlast,zero,zero,zpb(ii),szdum,vtzlast,vdzlast) vdy0km=(vdylast-xlast*dm2km) vdz0km=(vdzlast-xlast*dm2km) c --- Process puff at start dxkm=AMAX1(0.0,xttb1*dm2km+vdy0km) dts=dxkm/ukps call SIGTY(zero,dxkm,dts,syb1,vtyb1,vdyb1) syb1=AMAX1(syb1,symin) dxkm=AMAX1(0.0,xttb1*dm2km+vdz0km) dts=dxkm/ukps call SIGTZ(zero,dxkm,dts,zpb(ii),szb1,vtzb1,vdzb1) szb1=AMAX1(szb1,szmin) c --- Process midpoint xttm1=0.5*(xttb1+xtte1) dxkm=AMAX1(0.0,xttm1*dm2km+vdy0km) dts=dxkm/ukps call SIGTY(zero,dxkm,dts,sym1,vtym1,vdym1) sym1=AMAX1(sym1,symin) dxkm=AMAX1(0.0,xttm1*dm2km+vdz0km) dts=dxkm/ukps call SIGTZ(zero,dxkm,dts,zpb(ii),szm1,vtzm1,vdzm1) szm1=AMAX1(szm1,szmin) c --- Process endpoint dxkm=AMAX1(0.0,xtte1*dm2km+vdy0km) dts=dxkm/ukps call SIGTY(zero,dxkm,dts,sye1,vtye1,vdye1) sye1=AMAX1(sye1,symin) dxkm=AMAX1(0.0,xtte1*dm2km+vdz0km) dts=dxkm/ukps call SIGTZ(zero,dxkm,dts,zpb(ii),sze1,vtze1,vdze1) sze1=AMAX1(sze1,szmin) endif AFTER change: -------------- c --- Try to process sigmas at start of step call WAKE_TAB(istype,isnum,idw0(ipnum),xttb1, & syb1,szb1,hgr,rise,xlast,iprime) if(iprime.EQ.1) then c 070623 (MCB-D) c --- OK, puff starts in wake region; get virtuals at XLAST call WAKE_TAB(istype,isnum,idw0(ipnum),xlast, & sylast,szlast,hgr,rise,xlast2,idone) call SIGTY(sylast,zero,zero,sydum,vtylast,vdylast) call SIGTZ(szlast,zero,zero,zpb(ii),szdum,vtzlast,vdzlast) vdy0km=(vdylast-xlast*dm2km) vdz0km=(vdzlast-xlast*dm2km) 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 xttm1=0.5*(xttb1+xtte1) call WAKE_TAB(istype,isnum,idw0(ipnum),xttm1, & sym1,szm1,hgr,rise,xlast,idone) if(idone.EQ.1) then c --- Point within tabulated region call SIGTY(sym1,zero,zero,sydum,vtym1,vdym1) call SIGTZ(szm1,zero,zero,zpb(ii),szdum,vtzm1,vdzm1) else c --- Point beyond tabulated region dxkm=AMAX1(0.0,xttm1*dm2km+vdy0km) dts=dxkm/ukps call SIGTY(zero,dxkm,dts,sym1,vtym1,vdym1) dxkm=AMAX1(0.0,xttm1*dm2km+vdz0km) dts=dxkm/ukps call SIGTZ(zero,dxkm,dts,zpb(ii),szm1,vtzm1,vdzm1) endif sym1=AMAX1(sym1,symin) szm1=AMAX1(szm1,szmin) c --- Process endpoint call WAKE_TAB(istype,isnum,idw0(ipnum),xtte1, & sye1,sze1,hgr,rise,xlast,idone) if(idone.EQ.1) then c --- Point within tabulated region call SIGTY(sye1,zero,zero,sydum,vtye1,vdye1) call SIGTZ(sze1,zero,zero,zpb(ii),szdum,vtze1,vdze1) else c --- Point beyond tabulated region dxkm=AMAX1(0.0,xtte1*dm2km+vdy0km) dts=dxkm/ukps call SIGTY(zero,dxkm,dts,sye1,vtye1,vdye1) dxkm=AMAX1(0.0,xtte1*dm2km+vdz0km) dts=dxkm/ukps call SIGTZ(zero,dxkm,dts,zpb(ii),sze1,vtze1,vdze1) endif sye1=AMAX1(sye1,symin) sze1=AMAX1(sze1,szmin) endif -------------------------------- Update subroutine SETSLG -------------------------------- BEFORE change: -------------- subroutine setslg(ii,iold,tsamp,iru,istab,el,sigv,sigw,bvf,ws, & dxmold,dymold,dxmnew,dymnew,ldbhr,uavg,lcalm) AFTER change: -------------- subroutine setslg(ii,iold,tsamp,iru,istab,el,sigv,sigw,bvf,ws, & dxmold,dymold,dxmnew,dymnew,ldbhr,uavg) BEFORE change: -------------- c --- OUTPUT: c UAVG - real - Mean transport speed (m/s) AFTER change: -------------- c --- OUTPUT: c UAVG - real - Mean transport speed (m/s) c LCALM - logical - CALM logical BEFORE change: -------------- logical ldbhr AFTER change: -------------- logical ldbhr,lcalm BEFORE change: -------------- c --- Set dispersion option ("calms" use time-based sigmas) if(ws.LT.wscalm) then idopty=1 idoptz=1 endif AFTER change: -------------- c --- Set dispersion option ("calms" use time-based sigmas) if(ws.LT.wscalm) then idopty=1 idoptz=1 c 070623 (MCB-D) c --- Turn off any downwash and plume rise calculations if CALM idw0(ii)=0 xfinal(ii)=zero lcalm=.TRUE. else lcalm=.FALSE. endif -------------------------------- Update subroutine CAV_CONC -------------------------------- Note: These 2 changes are implemented for EACH receptor type BEFORE change: -------------- if(xr.GT.0 .AND. xrb.LE.x115b) then AFTER change: -------------- if(xr.GT.0 .AND. xrb.LE.x115b) then c 070623 (MCB-D) c if(xr.GT.0 .AND. xrb.LE.x115b) then xtest=AMAX1(xr,xrb) if(xtest.GT.0.0 .AND. xrb.LE.x115b) then BEFORE change: -------------- call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,dpbl,cbyq0,cbyqc) AFTER change: -------------- call CAV_SAMP(ldb,xr,yr,zr,xrb,yrb,ycb,fcav,fwak, & ntr,xtr,ztr,dpbl,cbyq0,cbyqc) c 070623 (MCB-D) c --- Screen out upwind receptor impacts (if any) if(xr.LE.0.0) cbyq0=0.0 if(xrb.LE.0.0) cbyqc=0.0 -------------------------------- Update subroutine CAV_SAMP -------------------------------- BEFORE change: -------------- c --- Cavity source c ----------------- if(ipositn.LT.4 .AND. syc.GT.zero) then c --- Receptor in projected building, cavity, or wake c --- Compute concentration due to cavity if(ipositn.EQ.3) then AFTER change: -------------- c --- Cavity source c ----------------- if(ipositn.LT.4 .AND. syc.GT.zero) then c --- Receptor in projected building, cavity, or wake c --- Compute concentration due to cavity c 070623 (MCB-D) c --- Initialize receptor ht to value passed into subroutine zrc=zr if(ipositn.EQ.3) then -------------------------------- Update subroutine WAKE_DBG -------------------------------- BEFORE change: -------------- c --- Set initial values dbsz=0.0 dbsy=0.0 dbhw=0.0 dbhc=0.0 dbrsz=0.0 dbuw=1.0 ipositn=4 AFTER change: -------------- c --- Set initial values dbsz=0.0 dbsy=0.0 dbhw=0.0 dbhc=0.0 dbrsz=0.0 dbuw=1.0 ipositn=4 c 070623 (MCB-D) dbduw=0.0 BEFORE change: -------------- c --- Set initial values dbsz=0.0 dbsy=0.0 dbhw=0.0 dbhc=0.0 dbrsz=0.0 dbuw=1.0 ipositn=4 AFTER change: -------------- c --- Set initial values dbsz=0.0 dbsy=0.0 dbhw=0.0 dbhc=0.0 dbrsz=0.0 dbuw=1.0 ipositn=4 c 070623 (MCB-D) dbduw=0.0 -------------------------------- Update subroutine WAKE_DFSN -------------------------------- BEFORE change: -------------- c --- Target stepping length (m) data dx0/2.0/ AFTER change: -------------- c --- Target stepping length (m) data dx0/2.0/ c 070623 (MCB-D) c --- Initialize all dummy vars (passed but not used) szdum=0.0 zkdum=0.0 sydum=0.0 ykdum=0.0 -------------------------------- Update subroutine WAKE_FQC -------------------------------- BEFORE change: -------------- fqcav=0.0 AFTER change: -------------- fqcav=0.0 c 070623 (MCB-D) xmax=0.0 fzmax=0.0 ----------------------------------------------------------------------- 2. CALMET ----------------------------------------------------------------------- * ***************************************************************** ******************** -- Problem Area 1 -- ******************** Cloud cover was not updated under stable conditions for icloud=3 Modified: ELUSTR2 ***************************************************************** --------------------------------------- (1) Update subroutine ELUSTR2 --------------------------------------- c properly update gridded cloud cover when icloud=3 c under stable conditions BEFORE ------- c --- stable conditions -- use Venkatram(1980), Weil & Brower(1983) 100 continue c c --- Determine cloud cover (tenths) for this grid cell if(icloud.eq.2)then AFTER ----- c --- stable conditions -- use Venkatram(1980), Weil & Brower(1983) 100 continue c c --- Determine cloud cover (tenths) for this grid cell c MCB-D (070623) if(icloud.eq.2 .or.icloud.eq.3)then ***************************************************************** ******************** -- Problem Area 2 -- ******************** The non-initialization of RHOP in RDMM5 could produce wrong results with some compilers (e.g. PGI in Linux) Modified: RDMM5 ***************************************************************** --------------------------------------- (2) Update subroutine RDMM5 --------------------------------------- c Modify computation of RHOP in RDMM5 such that it is c properly (re)initialized at each timestep and is always c taken at the surface (not always maximum rho) BEFORE ------- c --- Check whether need to skip a record iskip = 0 AFTER ----- c --- Check whether need to skip a record iskip = 0 c --- initialize surface density (MCB-D 070623) rhop(i,j)=0. irho=0 BEFORE ------- if ( z(n).ge.0.) then rhop(i,j)=max(rhop(i,j), 0.3484321*ipress/t(n)) endif AFTER ----- c --- MCB-D (070623) c --- Take surface density (not always the highest rho) if ( z(n).ge.0. .and. irho.eq.0) then rhop(i,j)= 0.3484321*ipress/t(n) irho=1 endif ***************************************************************** ******************** -- Problem Area 3 -- ******************** The debug-print flag (ldbhr) was missing in MIXHTST and MIXHTST2, potentially causing large output files to be generated when not requested (pb on some platforms, e.g. PGI compiler in Linux) or no output to be written out when requested Additionally nx,ny were not passed on to those subroutines, potentially compounding this large debug-printing problem. Modified: COMP, ELUSTR2,MIXHTST, MIXHT2ST, OUTPT.MET ***************************************************************** --------------------------------------- (3a) Update subroutine MIXHTST --------------------------------------- c Remove ldbhr from calling list and include OUTPT.MET (where it is c now defined) and GRID.MET c Rename nearu to nearus to avoid conflict with nearu array in grid.met BEFORE ------- subroutine mixhtST(el,ustar,qh,rho,ihrgmt,ist,jst,nearu, 1 iupt,ilandu,iwat1,iwat2,ldbhr,imixh,zi,ziconv) AFTER ----- subroutine mixhtST(el,ustar,qh,rho,ihrgmt,ist,jst,nearus, 1 iupt,ilandu,iwat1,iwat2,imixh,zi,ziconv) BEFORE ------- include 'params.met' AFTER ----- include 'params.met' include 'outpt.met' include 'grid.met' BEFORE ------- logical ldbhr,ldate AFTER ----- logical ldate BEFORE ------- iusta=nearu AFTER ----- iusta=nearus --------------------------------------- (3b) Update subroutine MIXHT2ST --------------------------------------- c Remove ldbhr from calling list and include OUTPT.MET (where c it is now defined) BEFORE ------- subroutine mixht2ST(el,ustar,qh,rho,ihrgmt,is,js, 1 ilandu,iwat1,iwat2,ldbhr,imixh,zi,ziconv) AFTER ----- subroutine mixht2ST(el,ustar,qh,rho,ihrgmt,is,js, 1 ilandu,iwat1,iwat2,imixh,zi,ziconv) BEFORE ------- include 'params.met' include 'grid.met' AFTER ----- include 'params.met' include 'outpt.met' include 'grid.met' BEFORE ------- logical ldbhr,ldate AFTER ----- logical ldate --------------------------------------- (3c) Update subroutine COMP --------------------------------------- c Do not declare ldbhr in COMP anylonger (done OUTPT.MET) BEFORE ------- logical ldate,lprthr,ldbhr AFTER ----- logical ldate,lprthr --------------------------------------- (3d) Update subroutine ELUSTR2 --------------------------------------- c Remove ldbhr from calling lists to MIXHTST and MIXHT2ST BEFORE ------- if (itprog.eq.0)then call mixhtST(el,ustar,qqh,rho,ihr2gmt,is,js, : nearu(is,js),iupt, : ilandu(is,js),iwat1,iwat2,ldbhr,imixh,zii,ziconvi) else call mixht2ST(el,ustar,qqh,rho,ihr2gmt,is,js, : ilandu(is,js),iwat1,iwat2,ldbhr,imixh,zii, : ziconvi) endif AFTER ----- if (itprog.eq.0)then call mixhtST(el,ustar,qqh,rho,ihr2gmt,is,js, : nearu(is,js),iupt, : ilandu(is,js),iwat1,iwat2,imixh,zii,ziconvi) else call mixht2ST(el,ustar,qqh,rho,ihr2gmt,is,js, : ilandu(is,js),iwat1,iwat2,imixh,zii, : ziconvi) endif --------------------------------------- (3e) Update include file OUTPT.MET --------------------------------------- c Include ldbhr in outpt.met BEFORE ------- logical lsave,lprint,ldb AFTER ----- logical lsave,lprint,ldb,ldbhr BEFORE ------- common/outpt/lsave,iformo,lprint,iprinf,ldb,nn1,nn2,iuvout(mxnz), 1 iwout(mxnz),itout(mxnz),imtout(8),iomet AFTER ----- common/outpt/lsave,iformo,lprint,iprinf,ldb,ldbhr,nn1,nn2, 1 iuvout(mxnz),iwout(mxnz),itout(mxnz),imtout(8),iomet ***************************************************************** ******************** -- Problem Area 4 -- ******************** For IEXTRP=1 (no vertical extrapolation), vertical extrapolation using similarity theory was performed Modified: DIAGNO ***************************************************************** --------------------------------------- (4) Update subroutine DIAGNO --------------------------------------- c restrict vertical extrapolation using similarity theory to c ABS(IEXTRP) = 4 BEFORE ------- ELSE C C EXTRAPOLATE WITH SIMILARITY THEORY AFTER ----- ELSE IF (IABS(IEXTRP) .EQ. 4) THEN C C EXTRAPOLATE WITH SIMILARITY THEORY ***************************************************************** ******************** -- Problem Area 5 -- ******************** Some variables were not initialized, with no consequence on the results per se, but with some compilers/platforms, un-initialized variables may stop execution Modified: PREPDI, RDMM4,RDMM5 ***************************************************************** --------------------------------------- (5a) Update subroutine PREPDI --------------------------------------- c Initialize um and vm BEFORE ------- htfac=1.0 AFTER ----- htfac=1.0 c --- MCB-D (070623) um=0. vm=0. --------------------------------------- (5b) Update subroutine RDMM4 --------------------------------------- c Initialize nlevag1,TZ,ZL at first call to subroutine BEFORE ------- c --- Initialize constants deg2rad=0.0174533 AFTER ----- c --- Initialize constants deg2rad=0.0174533 c --- MCB-D (070623) c --- Initialize variables if (ifirstpg.eq.0) then nlevag1=0 do k=1,mxnzp+1 do j=1,mxny do i=1,mxnx tz(i,j,k)=0. zl(i,j,k)=0. end do end do end do endif --------------------------------------- (5c) Update subroutine RDMM5 --------------------------------------- c Initialize nlevag1,TZ1,ZL1,TZ2,ZL2 at fist call to subroutine BEFORE ------- if (ifirstpg.eq.0) then it1=0 it2=0 AFTER ----- if (ifirstpg.eq.0) then it1=0 it2=0 c --- MCB-D (070623) nlevag1=0 BEFORE ------- 58 continue endif AFTER ----- 58 continue c --- MCB-D (070623) do 59 k = 1,mxnzp+1 do 59 j = 1,mxny do 59 i = 1,mxnx tz1(i,j,k)=0. zl1(i,j,k)=0. tz2(i,j,k)=0. zl2(i,j,k)=0. 59 continue endif ***************************************************************** ******************** -- Problem Area 6 -- ******************** In NOOBS-temperature mode (ITPROG>0), temperatures were adjusted to an adiabatic profile up to an undefined height (rather than up to the mixing height) Modified: TEMP3D ***************************************************************** --------------------------------------- (6) Update subroutine TEMP3D --------------------------------------- c Replace variable zic by ziconv(i,j) as zic is not defined c in NOOBS mode BEFORE ------- if (zmid(k).le.zic) & ztemp(i,j,k) = ztemp(i,j,1)-0.0098*(zmid(k)-zmid(1)) AFTER ----- c --- MCB-D (070623) if (zmid(k).le.ziconv(i,j)) & ztemp(i,j,k) = ztemp(i,j,1)-0.0098*(zmid(k)-zmid(1)) ***************************************************************** ******************** -- Problem Area 7 -- ******************** Undefined variable ws10 may cause execution stop with some compilers/platforms if calm wind conditions are not extrapolated (icalm=0) (but this did not affect results if execution went forward) Modified: WIND1 ***************************************************************** --------------------------------------- (7) Update subroutine WIND1 --------------------------------------- c Only readjust winds (under calm wind conditions)if valid data BEFORE ------- IF (ICALM.EQ.0) THEN do l=1,nssta+nowsta if (ws10(l).lt.0.0001)then do k=2,nz c --- if only partial weight of station with calm winds if (vvw(k,l).lt.0.9999) then do ll=1,nstat if (ll.ne.l) vvw(k,ll)=vvw(k,ll)/(1.-vvw(k,l)) end do c --- if full weight (i.e.only surface station with valid data c --- interpolate upper air data only (1/R**2) else do ll=nssta+nowsta+1,nstat vvw(k,ll)=rsqwt(ll)/(1.-rsqwt(l)) end do endif vvw(k,l)=0. end do endif end do ENDIF AFTER ----- IF (ICALM.EQ.0) THEN do l=1,nssta+nowsta c --- only readjust for valid data (MCB-D 070623) - if((us(1,L).lt.editl).and.(vs(1,L).lt.editl)) then if (ws10(l).lt.0.0001)then do k=2,nz c --- if only partial weight of station with calm winds if (vvw(k,l).lt.0.9999) then do ll=1,nstat if (ll.ne.l) vvw(k,ll)=vvw(k,ll)/(1.-vvw(k,l)) end do c --- if full weight (i.e.only surface station with valid data c --- interpolate upper air data only (1/R**2) else do ll=nssta+nowsta+1,nstat vvw(k,ll)=rsqwt(ll)/(1.-rsqwt(l)) end do endif vvw(k,l)=0. end do endif endif end do ENDIF ***************************************************************** ********************* -- Problem Area 8 -- ********************* For surface stations located on a water gridcell, the constant GEO.DAT roughness length (z0) was used rather than the Hosker z0 (function of the surface wind speed, as should be overwater) for the vertical extrapolation of surface winds with similarity theory (abs(IEXTRP)=4) Modified: WIND1 ***************************************************************** --------------------------------------- (8) Update subroutine WIND1 --------------------------------------- c Update roughness length and log profile coefficients for surface c stations located on a water gridcell, with Hosker z0 c Define minimum roughness length BEFORE ------- real vvw(mxnz,mxwnd) AFTER ----- real vvw(mxnz,mxwnd) c --- MCB-D (070623) data z0min/2.0e-6/ BEFORE ------- c --- Define combined surface+overwater station anemometer height array c --- and neutral log profile coefficient array AFTER ----- c --- MCB-D (070623) do L=1,nssta Isc = IST(L) Jsc = JST(L) if(isc.gt.nx)isc=nx if(isc.lt.1)isc=1 if(jsc.gt.ny)jsc=ny if(jsc.lt.1)jsc=1 if (ilandu(isc,jsc).ge.iwat1.and.ilandu(isc,jsc).le.iwat2) then if(us(1,l).lt.editl.and.vs(1,l).lt.editl)then c --- update z0 wsz = sqrt(us(1,L) **2 + vs(1,L) **2) c --- compute overwater surface roughness (Hosker, 1974) if(wsz.gt.0.0)then zz0 = 2.0e-6 * wsz ** 2.5 z0(isc,jsc) = amax1(z0min,zz0) else z0(isc,jsc) = z0min endif c --- update zlog xlnzo=alog(z0(isc,jsc)) c --- z1: first CALMET level xlnz1=alog(cellzc(1)) c --- z2: anemometer height xlnz2=alog(zanem(l)) c --- Logarithmic profile scaling factor c u(k=1)=zlogsta*u(zanem) zlogsta(l)=(xlnz1-xlnzo)/(xlnz2-xlnzo) endif endif enddo c --- Define combined surface+overwater station anemometer height array c --- and neutral log profile coefficient array ***************************************************************** ******************** -- Problem Area 9 -- ******************** The 10m wind speed overwater was not defined when the anemometer height was exactly 10m (problem for OCD method) Modified: WATER2 ***************************************************************** --------------------------------------- (9) Update subroutine WATER2 --------------------------------------- c Initialize ws10 BEFORE ------- c --- compute overwater surface roughness (Hosker, 1974) if(wsz.gt.0.0)then AFTER ----- c --- compute overwater surface roughness (Hosker, 1974) c --- MCB-D (070623) ws10=wsz if(wsz.gt.0.0)then ***************************************************************** ********************* -- Problem Area 10 -- ********************* Undefined index variable when skipping records might stop execution with some compilers/platforms (but not affect results) Modified: RDMM5 ***************************************************************** --------------------------------------- (10) Update subroutine RDMM5 --------------------------------------- c Fill in indx array only for valid records (not when skipping) BEFORE ------- else do 65 n=1,nlev c --- "Read" (skip) data read(io20,*) idum,idum 65 continue endif c frr (030106) c store first level above ground for each MM5 gridpoint indx(i,j)=index c end of loop on nxp,nyp 5 continue AFTER ----- c --- MCB-D (070623) c store first level above ground for each MM5 gridpoint indx(i,j)=index else do 65 n=1,nlev c --- "Read" (skip) data read(io20,*) idum,idum 65 continue endif c end of loop on nxp,nyp 5 continue ***************************************************************** ********************* -- Problem Area 11 -- ********************* Prognostic Land Use, which determines how prognostic temperatures are extrapolated from the lowest M32D level to lowest CALMET level, was not read in from old MM5 format Modified: RDHD51, RDHD52 ***************************************************************** --------------------------------------- (15a) Update subroutine RDHD51 --------------------------------------- c read in prognostic LU BEFORE ----- data lprt/.true./ AFTER ----- c --- Prognostic land use (MCB-D 070623) integer ilu4p(mxnxp,mxnyp) data lprt/.true./ BEFORE ----- if(cline(13:13).EQ.' ') then c --- Early MM5.DAT format for these records read(cline,47)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j) else c --- Revised MM5.DAT format for these records read(cline,48)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j) endif 47 format(2i3,f6.2,f8.2,i5) 48 format(2i3,f7.3,f8.3,i5) AFTER ----- if(cline(13:13).EQ.' ') then c --- MCB-D (070623 - LU read in) c --- Early MM5.DAT format for these records read(cline,47)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j), ilu4p(i,j) else c --- Revised MM5.DAT format for these records read(cline,48)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j),ilu4p(i,j) endif 47 format(2i3,f6.2,f8.2,i5,i3) 48 format(2i3,f7.3,f8.3,i5,i3) c --- Fill in ILU4 array stored in MM4HDO.MET ilu4(i,j)=ilu4p(i,j) --------------------------------------- (15b) Update subroutine RDHD52 --------------------------------------- c read in prognostic LU rather than inferring it from zero elevation BEFORE ----- data lprt/.true./ AFTER ----- c --- Prognostic land use (MCB-D 070623) integer ilu4p(mxnxp,mxnyp) data lprt/.true./ BEFORE ----- read(io20,99)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j) 99 format(2i4,f9.4,f10.4,i5) c --- MCB-C (070501) c --- assign ocean LU to zero elevation points if (ielev4(i,j).eq.0) then ilu4(i,j)=iluoc3d else c --- make sure that not all points are flagged as ocean c --- if by chance user input iluoc3d=0 ilu4(i,j)=iluoc3d+1 endif AFTER ----- c --- MCB-D (070623) c --- Land use is in fact recorded in MM5.DAT so read it read(io20,99)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j),ilu4p(i,j) 99 format(2i4,f9.4,f10.4,i5,i3) c --- Fill in ILU4 array stored in MM4HDO.MET ilu4(i,j)=ilu4p(i,j)