-----------------------------------------------------------------------

 --- 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)