----------------------------------------------------------------------- --- CALPUFF Modeling System ---

----------------------------------------------------------------------- Model Change Bulletin : B Dated : DECEMBER 16, 2005

----------------------------------------------------------------------- The following code modifications address problems identified in 2 modules of the CALPUFF modeling system: CALMET and CALPUFF. Instructions for editing the FORTRAN files for each are provided, and the programs can then be compiled using FORTRAN 95 to create the corresponding executables. Edits to the programs update the version and level, and they introduce the required changes to the computations performed. The version and level identify the program in the output files that are produced. Users who update their current copy of each program should add the character "b" to the end of the version to indicate that the program contains the modifications identified in Model Change Bulletin B, and replace the level with the date of this bulletin. Those currently using the MCB-A versions of CALMET and CALPUFF should use the new versions of the executables, or compile from the new versions of the programs. The versions are:

MCB-A MCB-B CALMET Version 5.53a (040716) Version 5.53b (051216) CALPUFF Version 5.711a (040716) Version 5.711b (051216)

Changes to CALMET

------------------ Problem Area 1 -- Upwind averaging of mixing height and temperature includes influence of cells outside of the intended upwind cone.
Modified: AVEMIX, AVETMP

Problem Area 2 -- CALMET simulations ended one timestep after the end of a leap year due to incorrect time interval calculation with prognostic data.
Modified: RDMM5

Problem Area 3 -- Variable declarations do not match intended properties in two subroutines.
Modified: MIXHT, RDOWProblem Area 4 -- The full NOOBS option (noobs=2) means " no upper air station and no surface stations". Offshore data should be accepted (SEA.DAT) but this was not allowed.
Modified: READCF

Problem Area 5 -- Dataset version 2.0 3D.DAT files created with CALMM5 Version 3.4 Level 050413, posted on the Web, were not recognized by CALMET.

New: RDHD53 Modified: MM4HD0.MET (include file) RDHD5, RDMM5

Changes to CALPUFF

------------------ Problem Area 1 -- Wet deposition fluxes reported for SLUGs are not correct. These are computed for an effective rainfall rate that is either 0 (no precipitation) or 1 mm/hr. The actual precipitation rate should be used in this calculation. The total mass removed due to wet depletion is correct, so all concentrations and dry deposition fluxes reported for these slugs are correct. Wet fluxes reported for PUFFs are correct. Because SLUGs transition to PUFFs when sigma-y becomes large compared to the length of a slug, most receptors affected by this bug are "near-field" receptors.
Modified: WET

Problem Area 2 -- run stops due to taking the square root of a negative number, or dividing by zero. This occurs for various combinations involving gradual rise in the PRIME downwash module, use of the PDF option, and selected coastline geometries in the subgrid TIBL module.

New: WARN Modified: COMP, FIN, GRISE, NUMRISE, PDFPATH, PUFRECS, SIGTY, SIGTZ, TIBLGRO, TIBLON, WAKE_DFSN

Problem Area 3 -- An invalid file unit number is placed in several calls to subroutine YR4.
Modified: GETRCRD, RDHDBC2

Problem Area 4 -- Non-gridded receptor elevations are reported to the QA file for gridded receptors. Modified: QAPLOT1-----------------------------------------------------------------------

----------------------------------------------------------------------- 1. CALMET

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

Edit PARAMS.MET, PARAMSS.MET, PARAMSL.MET, PARAMSX.MET Change the version by adding the letter "b" Change the level to 051216 (example) BEFORE change:

-------------- c --- Specify model version character*8 mver,mlevel parameter(mver='5.53a',mlevel='040716') AFTER change:

-------------- c --- Specify model version character*8 mver,mlevel parameter(mver='5.53b',mlevel='051216')

Edit CALMET.FOR

------------------------------------- (a) Update the MAIN program documentation

------------------------------------- Change the version by adding the letter "b" Change the level to 051216 (example) BEFORE change:

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

c c --- set version and level number of program ver='5.53a' level='040716'
c AFTER change:

--------------c c --- set version and level number of program ver='5.53b' level='051216'
c

--------------------------------------- (b) Update subroutine RDMM5

--------------------------------------- c --- UPDATE c --- The time interval between prognostic records has to be c --- computed throughout the simulation period. The code did c --- not compute that interval correctly when the simulation c --- ended one timestep after the end of a leap year.

BEFORE change:

-------------- c --- Make sure dates are okay to compute final dtinc in case c --- of straddling days or years (frr-030526) if ( (ndathr-it1).gt.isteppg) then if ( (ndathr-100+24-it1).gt.isteppg) then c turn of the year it1=it1+100000-364*100-24 else c straddling days: it1=it1+100-24 endif endif AFTER change:

-------------- c --- Make sure dates are okay to compute final dtinc in case c --- of straddling days or years (frr-030526) if ( (ndathr-it1).gt.isteppg) then if ( (ndathr-100+24-it1).gt.isteppg) then c turn of the year c --- normal year it1=it1+100000-364*100-24 c --- MCB-B- leap year (051216) if ((ndathr-it1).gt.isteppg) it1=it1-100 else c straddling days: it1=it1+100-24 endif endif--------------------------------------- (c) Update subroutine MIXHT

--------------------------------------- c --- UPDATE c --- The surface air density (rho) is computed at all gridpoints c --- (2D-array). RHO was still declared as a 1D array in MIXHT. c --- Declare RHO as a 2D array in that subroutine

BEFORE change:

-------------- c RHO(mxss) - real array - Air density (kg/m**3)

AFTER CHANGE

------------ c --- MCB-B (051216) c RHO(mxnx,mxny) - real array - Air density (kg/m**3)

BEFORE change:

-------------- real rho(mxss)

AFTER CHANGE

------------ c --- MCB-B (051216) real rho(mxnx,mxny) BEFORE change:

-------------- c --- Convert heat flux (W/m**2) to w'theta' (m * deg. K/s) wt=qh(i,j)/(rho(issta)*cp)

AFTER CHANGE

------------ c --- Convert heat flux (W/m**2) to w'theta' (m * deg. K/s) c --- MCB-B (051216) wt=qh(i,j)/(rho(i,j)*cp)

--------------------------------------- (d) Update subroutine RDOW

--------------------------------------- c --- UPDATE c --- The overwater station coordinates are single-precision variables. c --- Comment out the declaration that made them double-precision

BEFORE change:

-------------- real*8 xkm,ykmAFTER CHANGE

------------ c --- MCB-B (051216) c real*8 xkm,ykm
c

------------------------- (e) Update subroutine READCF

------------------------- c --- UPDATE c --- The full NOOBS option (noobs=2) means " no upper air station and c --- no surface stations". It does not imply no offshore data. Therefore c --- a SEA.DAT file can be used in conjunction with the NOOBS=2 option c --- => allow nowsta>0 when noobs=2

BEFORE change:

-------------- elseif(noobs.EQ.2) then if(nusta.GT.0 .OR. nssta.GT.0 .OR. nowsta.GT.0) then write(io6,*) write(io6,*) 'READCF: Warning in Input Group 4' write(io6,*) 'No observations for NOOBS = ',noobs write(io6,*) 'NSSTA reset to 0 from ',nssta write(io6,*) 'NUSTA reset to 0 from ',nusta write(io6,*) 'NOWSTA reset to 0 from ',nowsta lwarncf=.TRUE. nssta=0 nowsta=0 nusta=0 endif endif

AFTER change

------------ c --- MCB-B (051216) - SEA.DAT can be used with noobs=2 elseif(noobs.EQ.2) then if(nusta.GT.0 .OR. nssta.GT.0) then write(io6,*) write(io6,*) 'READCF: Warning in Input Group 4' write(io6,*) 'No observations for NOOBS = ',noobs write(io6,*) 'NSSTA reset to 0 from ',nssta write(io6,*) 'NUSTA reset to 0 from ',nusta lwarncf=.TRUE. nssta=0 nusta=0 endif endif--------------------------------------- (f) Update external include file MM4HD0.MET

--------------------------------------- c --- UPDATE c --- Modified for 3D.DAT files version 2.0 and higher

BEFORE change:

-------------- character*8 datum3d AFTER change:

-------------- character*8 datum3d c --- MCB-B (051216) character*16 cver3d BEFORE change:

-------------- 6 IOUTMM5,IMM53D,ISTEPPG,DATUM3D AFTER change:

-------------- 6 IOUTMM5,IMM53D,ISTEPPG,DATUM3D,CVER3D

--------------------------------------- (g) Update subroutine RDHD5

--------------------------------------- c --- UPDATE c --- Modified for 3D.DAT files version 2.0 and higher

BEFORE change:

-------------- character*80 buff AFTER Cchange

------------ c --- MCB-B (051216) character*80 buff1, buff2

BEFORE change:

-------------- c --- Read first two records to determine whether file is c in 3D.DAT format (Data version 3.0) read(io20,'(a)')buff read(io20,'(a)')buff read(buff,'(a12)')cset3dAFTER change

------------ c --- MCB-B (051216) c --- Read first two records to determine whether file format c --- iMM53d = 0 for MM5.DAT file structure c --- iMM53d = 1 for 3D.DAT file structure prior to Version 2.0 c --- iMM53d = 2 for 3D.DAT file structure, Version 2.0 or later
c read(io20,'(a)')buff1 read(io20,'(a)')buff2
c c --- Step (1): Check for "3D.DAT" or "MM53D.DAT" on second c --- record of file (indication of old 3D.DAT file structure c --- prior to Version 2.0) read(buff2,'(a12)')cset3d

BEFORE change:

------------- else imm53d=0 endif AFTER change

------------ c --- MCB-B (051216)
else
c c --- Step (2) check for new standard 3D.DAT file structure c --- with dataset name, version and comments on first record if(buff1(1:6).eq.'3D.DAT')then imm53d=2 else
c c --- Structure does not fit 3D.DAT conventions -- assume c --- file is in MM5.DAT format imm53d=0 endif endif

BEFORE change:

------------- open(io20,file=mm4dat,status='old') c CALMM5.DAT in original CALMM5 formats (No data set name. c Record 1 includes title, version # and level #) if(imm53d.eq.0) call rdhd51 c CALMM5.DAT in 3D.DAT format (Record 1 is title, Record 2 is c dataset name, version # and level #)if(imm53d.eq.1) call rdhd52 AFTER change

------------ c --- MCB-B (051216) open(io20,file=mm4dat,status='old')
c c ----------------------- c --- Read header records
c ----------------------- if(imm53d.eq.0)then
c c --- Data is in original MM5 format (No data set name) c --- Record 1 includes a comment line call rdhd51
c else if(imm53d.eq.1)then
c c --- Data is in old 3D.DAT format prior to version 2.0 c --- (Record 1 is title, Record 2 is dataset name, c --- code version # and level #) call rdhd52 else if(imm53d.eq.2)then
c c --- Data is in new 3D.DAT format (dataset version 2.0 or c --- later) call rdhd53 else
c c --- Unexpected value of imm53d -- write error message write(io6,*)'Error in Subr. RDHD5 -- Unexpected value ', 1 'of IMM53D -- IMM53D = ',imm53d stop endif

--------------------------------------- (g) Update subroutine RDMM5

--------------------------------------- c --- UPDATE c --- Modified for 3D.DAT files version 2.0 and higher

BEFORE change:

------------- else if (ioutmm5.eq.83) then read(io20,83) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr 83 format (i4,i6,f6.1,i4,f5.1,i3,3(f5.2)) else if (ioutmm5.eq.84) thenread(io20,84) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr,qi,qs 84 format (i4,i6,f6.1,i4,f5.1,i3,5(f5.2)) else if (ioutmm5.eq.85) then read(io20,85) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr,qi,qs,qg 85 format (i4,i6,f6.1,i4,f5.1,i3,6(f5.2)) AFTER change:

------------- c --- MCB-B (051216) else if (ioutmm5.eq.83) then if (cver3d.eq.'2.1') then read(io20,830) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr else read(io20,83) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr endif 830 format (i4,i6,f6.1,i4,f5.1,i3,f5.2,2(f6.3)) 83 format (i4,i6,f6.1,i4,f5.1,i3,3(f5.2)) else if (ioutmm5.eq.84) then if (cver3d.eq.'2.1') then read(io20,840) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr,qi,qs else read(io20,84) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr,qi,qs endif 840 format (i4,i6,f6.1,i4,f5.1,i3,f5.2,4(f6.3)) 84 format (i4,i6,f6.1,i4,f5.1,i3,5(f5.2)) else if (ioutmm5.eq.85) then if (cver3d.eq.'2.1') then read(io20,850) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr,qi,qs,qg else read(io20,85) ipress,iz,t(n),iwd,ws(n), : irh,q,qc,qr,qi,qs,qg endif 850 format (i4,i6,f6.1,i4,f5.1,i3,f5.2,5(f6.3)) 85 format (i4,i6,f6.1,i4,f5.1,i3,6(f5.2))

BEFORE change:

-------------- c Compressed output else if (ioutmm5.ge.93) then read(io20,'(a)')buff c Read common parts and compression flag read(buff,93) ipress,iz,t(n),iwd,ws(n),w, : irh,q,fcomp 93 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,3(f5.2)) c Un-compression fromif(fcomp.gt.-0.0001) then if(ioutmm5.eq.93) then read(buff,93)ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr else if (ioutmm5.eq.94) then read(buff,94)ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr,qi,qs 94 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,5(f5.2)) else if (ioutmm5.eq.95) then read(buff,95) ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr,qi,qs,qg 95 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,6(f5.2)) endif c Compressed form else if(ioutmm5.eq.93) then qc=0 qr=0 else if (ioutmm5.eq.94) then qc=0 qr=0 qi=0 qs=0 else if (ioutmm5.eq.95) then qc=0 qr=0 qi=0 qs=0 qq=0 endif endif

endif AFTER change:

-------------- c Compressed output else if (ioutmm5.ge.93) then read(io20,'(a)')buff c Read common parts and compression flag if(cver3d.eq.'2.1') then read(buff,930) ipress,iz,t(n),iwd,ws(n),w, : irh,q,fcomp else read(buff,93) ipress,iz,t(n),iwd,ws(n),w, : irh,q,fcomp endif 93 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,3(f5.2)) 930 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,f5.2,2(f6.3)) c Un-compression from if(fcomp.gt.-0.0001) then if(ioutmm5.eq.93) then if (cver3d.eq.'2.1') thenread(buff,930)ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr else read(buff,93)ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr endif else if (ioutmm5.eq.94) then if (cver3d.eq.'2.1') then read(buff,940)ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr,qi,qs else read(buff,94)ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr,qi,qs endif 940 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,f5.2,4(f6.3)) 94 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,5(f5.2)) else if (ioutmm5.eq.95) then if (cver3d.eq.'2.1') then c --- new moisture format (6.3 instead of 5.2)(050504) read(buff,950) ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr,qi,qs,qg else read(buff,95) ipress,iz,t(n),iwd,ws(n),w, : irh,q,qc,qr,qi,qs,qg endif 95 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,6(f5.2)) 950 format (i4,i6,f6.1,i4,f5.1,f6.2,i3,f5.2,5(f6.3)) endif c Compressed form else if(ioutmm5.eq.93) then qc=0. qr=0. else if (ioutmm5.eq.94) then qc=0. qr=0. qi=0. qs=0. else if (ioutmm5.eq.95) then qc=0. qr=0. qi=0. qs=0. qq=0. endif endif

endif

---------------------------------------(h) Insert new subroutine RDHD53

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

c---------------------------------------------------------------------- subroutine rdhd53
c---------------------------------------------------------------------- c c --- CALMET Version: 5.53b Level: 051216 RDHD53 c J. Scire, Earth Tech c Adapted from RDHD52
c c --- PURPOSE: Read the header records from a file in 3D.DAT format c (dataset Version 2.0 or later)
c c --- INPUTS:
c c Parameters: c MXNXP, MXNYP, MXNZP, IO6, IO20, MXNX, MXNY c Common block /MAP/: c IUTMZN, UTMHEM, XLAT1, XLAT2, RELON0, RNLAT0, c FEAST, FNORTH, DATUM c Common block /GRID/: c NX,NY,DGRID,XORIGR,YORIGR c Common block /MM4HDO/ variables: c DATUM3D
c c --- OUTPUT:
c c Common block /MM4HDO/ variables: c IBYRM, IBJULM, IBHRM, IEYRM, IEJULM, IEHRM, c NXMM4, NYMM4, NZP, PTOPMM4, I1, J1, NXP, NYP, c SIGMA(mxnzp), XLAT4(mxnxp,mxnyp), XLONG4(mxnxp,mxnyp), c IELEV4(mxnxp,mxnyp),ILU4(mxnxp,mxnyp),XLCMM4(mxnxp,mxnyp), c YLCMM4(mxnxp,mxnyp),IGRAB(mxnx,mxny,4),JGRAB(mxnx,mxny,4), c IOUTMM5, CVER3D, c DATUM3D
c c --- RDHD53 called by: RDHD5 c --- RDHD53 calls: JULDAY, INCR, YR4, GLOBE1, GLOBE
c---------------------------------------------------------------------- c c --- Include parameters include 'params.met'
c character*128 ctemp3d character*16 cname3d character*64 ctitle3d character*132 comm3d logical lprt
c c --- Common blocks include 'map.met' include 'grid.met' include 'mm4hdo.met'COMMON /D6/ IRD,IWR,IFILE,IRDP
c --- For coordinate transformations character*8 cmapi,cmapo character*12 caction character*4 c4hem real*8 vecti(9),vecto(9)
c data lprt/.true./ c --- Scale factor for Tangential TM projection tmsone=1.00000 c --- Set translation vectors going from N.lat/E.lon c --- to projection(x,y)km iutmo=iutmzn if(utmhem.EQ.'S ' .AND. iutmzn.LT.900) iutmo=-iutmo cmapo=pmap if(cmapo.EQ.'TTM ') cmapo='TM ' cmapi='LL ' idum=0 rdum=0.0 call GLOBE1(cmapi,idum,rdum,rdum,rdum,rdum,rdum, & rdum,rdum, & cmapo,iutmo,tmsone,xlat1,xlat2,rnlat0,relon0, & feast,fnorth, & caction,vecti,vecto)
c c ------------------------------------------------------------------ c --- Read header record #1 (Dataset name, version and title
c ------------------------------------------------------------------ read(io20,10)cname3d,cver3d,ctitle3d 10 format(2a16,a64)

c ------------------------------------------------------------------ c --- Read header record #2 (Number of comment lines)
c ------------------------------------------------------------------ read(io20,11)ncomm3d 11 format(i4)

c ------------------------------------------------------------------ c --- Next "NCOMM3d" lines (comment lines)
c ------------------------------------------------------------------ if(ncomm3d.gt.0)then do i=1,ncomm3d read(io20,12)ctemp3d 12 format(a132) c --- Save first line of text for later printing if(i.eq.1)comm3d=ctemp3d enddo endif

c ------------------------------------------- c --- Next header record (MM5 output options)c ------------------------------------------- read(io20,43)ioutw,ioutq,ioutc,iouti,ioutg ioutmm5=81+10*ioutw+ioutq+ioutc+iouti+ioutg 43 format(5(i3))

c -------------------------------------------- c --- Skip next header record (map projection)
c -------------------------------------------- read(io20,*)

c ------------------------------------------------ c --- Next header record (MM5 output options) c (Add 13 output options for surface variables)
c ------------------------------------------------ read(io20,44) inhyd,imphys,icupa,ibltyp,ifrad,isoil : ,ifddaan,ifddaob : ,igrdt,ipbl,ishf,ilhf,iustr,iswdn : ,ilwdn,ist1,ist2,ist3,ist4,ist5,ist6 44 format(30(i3))
c -------------------------------------- c --- Next header record (MM5 grid data)
c -------------------------------------- read(io20,20)ibyrm,ibmom,ibdym,ibhrm,nhrsmm5, 1 nxp,nyp,nzp 20 format(i4,3i2,i5,3i4) call YR4(io6,ibyrm,ierrb) if(ierrb.NE.0)then write(io6,*)'Error encountered in Subr. YR4 ' write(io6,*)'Execution stopping in Subr. RDHD53 ', 1 '-- IERRB = ',ierrb stop endif
c c --- Calculate Julian day call julday(io6,ibyrm,ibmom,ibdym,ibjulm)
c c --- Compute ending date/time (comment out if using other format) ieyrm=ibyrm iejulm=ibjulm iehrm=ibhrm call incr(io6,ieyrm,iejulm,iehrm,nhrsmm5)
c c -------------------------------------------------- c --- Next header record (extraction subdomain data)
c -------------------------------------------------- read(io20,30)nx1,ny1,nx2,ny2,nz1,nz2, & rxmin,rxmax,rymin,rymax i1=nx1 j1=ny1 30 format(6i4,2f10.4,2f9.4)
c c --- Check that array dimensions are not exceeded if(nxp.gt.mxnxp.or.nyp.gt.mxnyp.or.nzp.gt.mxnzp)thenwrite(io6,*)'ERROR in subr. RDHD53 -- Array dimensions ', 1 'are too small for data being read' write(io6,*)'Grid being read (NXP, NYP, NZP) = ', 1 nxp,nyp,nzp write(io6,*)'Array dimensions (MXNXP, MXNYP, MXNZP) = ', 1 mxnxp,mxnyp,mxnzp stop endif c --- Check consistency between nz1,nz2, and nzp if(nzp.ne.nz2-nz1+1) then write(io6,*)'Error in RDHD53: NZ1,NZ2 and NZP not consistent' write(io6,*)'nz1,nz2,nzp:',nz1,nz2,nzp stop endif
c c --------------------------------------------- c --- Next NZP records -- MM5 half-sigma levels
c --------------------------------------------- do 40 n=1,nzp read(io20,38)sigma(n) 38 format(f6.3) 40 continue
c c -------------------------------------------------------- c --- Print the 3D.DAT header information to the list file c --- (except for the gridded fields)
c -------------------------------------------------------- if(lprt)then write(io6,101)cname3d,cver3d,ctitle3d,comm3d 101 format(//1x,'Information read from 3D.DAT file'/ 1 5x,'Dataset Name: ',a12/ 2 5x,'Dataset Version: ',a12/ 3 5x,'Dataset Title: ',a64/ 4 5x,'First line of comments: ', 5 8x,a132)
c write(io6,102) inhyd,imphys,icupa,ibltyp,ifrad,isoil, : ifddaan,ifddaob 102 format(/5x,'MM5 physics options: '/ 1 5x,' Hydrostatic: ',i2/ 1 5x,' Moisture scheme: ',i2/ 1 5x,' Convection scheme: ',i2/ 1 5x,' Boundary layer scheme: ',i2/ 1 5x,' Radiation scheme ',i2/ 1 5x,' Soil scheme: ',i2/ 1 5x,' Analysis FDDA: ',i2/ 1 5x,' Observation FDDA: ',i2)

write(io6,1021)igrdt,ipbl,ishf,ilhf,iustr,iswdn : ,ilwdn,ist1,ist2,ist3,ist4,ist5,ist6 1021 format(/5x,'MM5 surface variable options: '/ 1 5x,' Ground temperature: ',i2/ 1 5x,' PBL: ',i2/1 5x,' Sensible heat flux: ',i2/ 1 5x,' Latent heat flux: ',i2/ 1 5x,' Frictional velocity: ',i2/ 1 5x,' Downward SW radiation: ',i2/ 1 5x,' Downward LW radiation: ',i2/ 1 5x,' Soil temp at layer 1: ',i2/ 1 5x,' Soil temp at layer 2: ',i2/ 1 5x,' Soil temp at layer 3: ',i2/ 1 5x,' Soil temp at layer 4: ',i2/ 1 5x,' Soil temp at layer 5: ',i2/ 1 5x,' Soil temp at layer 6: ',i2) write(io6,103)1,ioutw,ioutq,ioutc,iouti,ioutg 103 format(/5x,'MM5 output fields (1 = YES; 0 = NO): '/ 1 5x,' Pressure, height,T,Wind speed and direction: ',i2/ 1 5x,' Vertical velocity: ',i2/ 1 5x,' RH and vapor mixing ratio: ',i2/ 1 5x,' Cloud and rain mixing ratios: ',i2/ 1 5x,' Ice and Snow mixing ratios: ',i2/ 1 5x,' Graupel: ',i2) write(io6,104)ibyrm,ibmom,ibdym,ibhrm,nhrsmm5, 1 nxp,nyp,nzp 104 format(/5x,'Date/time (YYYYMMDDHH) of MM5 data: '/ 1 5x,' Start = ',i4,3i2,' (GMT)'/ 2 5x,' No. hours = ',i4,/ 3 5x,'Extraction Subdomain in MM5 file: '/ 4 5x,' No. X cells = ',i4/ 5 5x,' No. Y cells = ',i4/ 6 5x,' No. layers = ',i4 )
c write(io6,106)nx1,ny1,nz1,nx2,ny2,nz2 106 format(/5x,' Beginning X = ',i4/ 2 5x,' Beginning Y = ',i4/ 3 5x,' Beginning Z = ',i4/ 4 5x,' Ending X = ',i4/ 5 5x,' Ending Y = ',i4/ 6 5x,' Ending Z = ',i4)
c write(io6,107)rymin,rymax,rxmin,rxmax 107 format(/5x,' Latitude range : ',f9.4,' - ',f9.4/ : 5x,' Longitude range: ',f10.4, ' - ' ,f10.4) write(io6,108) 108 format(/5x,'MM5 half-sigma levels'/5x,'Level',5x,'Sigma'/) do 110 i=1,nzp write(io6,109)i,sigma(i) 110 continue 109 format(4x,i4,6x,f6.4) endif
c c ---------------------------------------------------- c --- Next NXP * NYP records -- lat., long., elevation
c ----------------------------------------------------do 50 j=1,nyp do 50 i=1,nxp read(io20,99)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j) 99 format(2i4,f9.4,f10.4,i5) c --- Compute grid point locations from N.Lat and E.Lon call GLOBE(io6,caction,datum3d,vecti,datum,vecto, & xlong4(i,j),xlat4(i,j),xlcmm4(i,j),ylcmm4(i,j), & idum,c4hem)
c c --- QA check that I,J read match expected values
c icheck=iindex-i1+1 jcheck=jindex-j1+1 if(icheck.ne.i.or.jcheck.ne.j)then write(io6,*)'ERROR in subr. RDHD53 -- I,J do not match ', 1 'values read on header record' write(io6,*)'I, J = ',i,j write(io6,*)'ICHECK, JCHECK = ',icheck,jcheck stop endif 50 continue

return end

----------------------------------------------------------------------- 2. CALPUFF

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

Edit PARAMS.PUF, PARAMSS.PUF, PARAMSL.PUF, PARAMSX.PUF Change the version by adding the letter "b" Change the level to 051216 (example) BEFORE change:

-------------- c --- Specify model version character*12 mver, mlevel, mmodel parameter(mver='5.711a',mlevel='040716') AFTER change:-------------- c --- Specify model version character*12 mver, mlevel, mmodel parameter(mver='5.711b',mlevel='051216')

Edit CALPUFF.FOR

------------------------------------- (a) Update the MAIN program documentation

------------------------------------- Change the version by adding the letter "b" Change the level to 051216 (example) BEFORE change:

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

c c --- set version and level number of program ver='5.711a' level='040716' AFTER change:

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

c c --- set version and level number of program ver='5.711b' level='051216'

-------------------------------- (b) Update subroutine WET

-------------------------------- c --- UPDATE c --- V5.7-5.711b 051216 (MCB-B): Multiply XLAM by the precipitation c rate to produce the scavenging ratio

BEFORE change:

-------------- c --- Determine the amount of pollutant mass remaining after wet c --- removal temp1=-rmm(ixs,iys)*tsamp do i=1,nspecxlam(i)=wa(ilq,i) temp=exp(xlam(i)*temp1) fracwet(i)=temp qu(i,jp)=temp*qu(i,jp) qm(i,jp)=temp*qm(i,jp) enddo
c 101 continue

AFTER change:

-------------- c 051216 (MCB-B) c --- Determine the amount of pollutant mass remaining after wet c --- removal do i=1,nspec xlam(i)=wa(ilq,i)*rmm(ixs,iys) fracwet(i)=exp(-xlam(i)*tsamp) qu(i,jp)=fracwet(i)*qu(i,jp) qm(i,jp)=fracwet(i)*qm(i,jp) enddo
c 101 continue

-------------------------------- (c) Update subroutine GETRCRD

-------------------------------- c --- UPDATE c --- V5.711(030528) to V5.711b(051216) (MCB-B) c File unit passed to subroutine YR4 for listfile output was c not defined so that error reports attempted to write to c unit 0. Changed io1 to io6 in YR4 call argument.

BEFORE change:

-------------- c --- Enforce YYYY format for year call YR4(io1,myre,ierr) if(ierr.NE.0) stop 'Halted in GETRCRD'

AFTER change:

-------------- c --- Enforce YYYY format for year c 051216 (MCB-B) call YR4(io6,myre,ierr)if(ierr.NE.0) stop 'Halted in GETRCRD'

-------------------------------- (d) Update subroutine RDHDBC2

-------------------------------- c --- UPDATE c --- V5.711(030528) to V5.711b(051216) (MCB-B) c File unit passed to subroutine YR4 for listfile output was c not defined so that error reports attempted to write to c unit 0. Changed io1 to io6 in YR4 call argument.

BEFORE change:

-------------- c --- Set flag indicating CALPUFF data version is 'MOD5' or older ipver=1 call YR4(io1,msyrBC,ierry) if(ierry.NE.0) stop 'Halted in RDHDBC2'

AFTER change:

-------------- c --- Set flag indicating CALPUFF data version is 'MOD5' or older ipver=1 c 051216 (MCB-B) call YR4(io6,msyrBC,ierry) if(ierry.NE.0) stop 'Halted in RDHDBC2'

-------------------------------- (e) Update subroutine QAPLOT1

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (MCB-B): Fix gridded receptor file QARECG.DAT c [used discrete rec. elevations] c --- (MCB-B): Enlarge format for QATERR.DAT hts

BEFORE change:

--------------c --- Data, in rows of constant Y do j=1,ny write(ioqa,'(10000(1pe10.4,2x))') (elev(i,j),i=1,nx) enddo close(ioqa) AFTER change:

-------------- c --- Data, in rows of constant Y do j=1,ny c 051216 (MCB-B) write(ioqa,'(10000(1pe11.4,1x))') (elev(i,j),i=1,nx) enddo close(ioqa)

BEFORE change:

-------------- c --- Process gridded receptors [DAT format]
c ------------------------------------------ open(ioqa,file='qarecg.dat',status='unknown') if(LSAMP) then write(ioqa,'(a18)') '"Xkm" "Ykm" "Elev"' do j=1,nysam do i=1,nxsam x=xsamLL+float(i-1)*delsam*fm2km y=ysamLL+float(j-1)*delsam*fm2km write(ioqa,*) x,y,elevng(i) enddo enddo close(ioqa) else

AFTER change:

-------------- c --- Process gridded receptors [DAT format]
c ------------------------------------------ open(ioqa,file='qarecg.dat',status='unknown') if(LSAMP) then write(ioqa,'(a18)') '"Xkm" "Ykm" "Elev"' do j=1,nysam do i=1,nxsam x=xsamLL+float(i-1)*delsam*fm2km y=ysamLL+float(j-1)*delsam*fm2km c 051216 (MCB-B) write(ioqa,*) x,y,elevg(i,j) enddo enddo close(ioqa)-------------------------------- (f) Update subroutine PUFRECS

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (MCB-B): fix bug that did not update the rise c factor when the PRIME wake tables are c used to obtain receptor-specific c sigmas

BEFORE change:

-------------- c --- Extract sigmas if this puff was emitted with current met if(imet.EQ.1) then call WAKE_TAB(istype,isnum,idw0(ipnum),xtt1, & syr,szr,hgr,rise,xlast,idone) c --- Finish up and skip to terrain adjustments if sigmas c --- are extracted if(idone.EQ.1) then if(rise.GT.0.0) then c --- Add BID contribution to sigmas sigbidsq=(rise/3.5)**2 syr=sqrt(syr**2+sigbidsq) szr=sqrt(szr**2+sigbidsq) endif go to 500 endif endif

AFTER change:

-------------- c --- Extract sigmas if this puff was emitted with current met if(imet.EQ.1) then call WAKE_TAB(istype,isnum,idw0(ipnum),xtt1, & syr,szr,hgr,rise,xlast,idone) c --- Finish up and skip to terrain adjustments if sigmas c --- are extracted if(idone.EQ.1) then c 051216 (MCB-B) rfacsq=0.0 if(rise.GT.0.0) then c --- Add BID contribution to sigmas sigbidsq=(rise/3.5)**2syr=sqrt(syr**2+sigbidsq) szr=sqrt(szr**2+sigbidsq) c 051216 (MCB-B) c --- Compute the corresponding BIDSQ adjustment factor rfacsq=sigbidsq/bidsq endif go to 500 endif endif

-------------------------------- (g) Update subroutine FIN

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (MCB-B): Add warning function (call WARN)

BEFORE change:

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

c write(iomesg,*)'TERMINATION PHASE'

AFTER change:

-------------- c 051216 (MCB-B) c --- Process any warnings collected during the run call WARN('FIN ',0.)
c write(iomesg,*)'TERMINATION PHASE'

-------------------------------- (h) Update subroutine NUMRISE

-------------------------------- c --- UPDATE c --- V5.71-V5.711b 051216 (MCB-B): initialize first element of MXNW c arrays to stack conditions to avoid c problems with the NN pointer being c oldBEFORE change:

-------------- c --- PARAMETERS NP=NSTEP XNP=FLOAT(NP) nnp=1 AFTER change:

-------------- c --- PARAMETERS NP=NSTEP XNP=FLOAT(NP) nnp=1 c 051216 (MCB-B) c --- Initialize first location in rise storage arrays to one step (ds) c --- and no rise; this should be overwritten with actual rise results. NN=1 XT(NN)=ds ZT(NN)=h Z0T(NN)=h RT(NN)=reff

BEFORE change:

-------------- c --- Update path arrays if(float(NNP/NP).eq.float(NNP)/XNP) THEN AFTER change:

-------------- c 051216 (MCB-B) c --- Store any NNP<NP point in NN=1 array location in case rise ends c --- too quickly to reach the first scheduled point if(nnp.LT.np) then NN=1 XT(NN)=X ZT(NN)=zc Z0T(NN)=z RT(NN)=R endif c --- Update path arrays if(float(NNP/NP).eq.float(NNP)/XNP) THEN

BEFORE change:

--------------if(x15r.GT.0.0) then c --- Set stepsize dsfin=nstep*ds dx15r=x15r/(mxnw-nn) dx15r=AMAX1(dsfin,dx15r) c --- Set range for additional steps n15r=MIN(NINT(x15r/dx15r),mxnw-nn) nbeg=nn+1 nend=nn+n15r do in=nbeg,nend AFTER change:

-------------- c 051216 (MCB-B) nnleft=mxnw-nn if(x15r.GT.0.0) then c --- Check for room in arrays if(nnleft.LE.0)then write(io6,*)'Error in Subr. NUMRISE -- NN too large ', & '-- NN = ',nn,' MXNW = ',mxnw write(io6,*)'Cannot complete trajectory out to 15R' write(*,*) stop 'Halted in NUMRISE -- see list file.' endif c --- Set stepsize dsfin=nstep*ds dx15r=x15r/FLOAT(nnleft) dx15r=AMAX1(dsfin,dx15r) c --- Set range for additional steps n15r=MIN(NINT(x15r/dx15r),nnleft) nbeg=nn+1 nend=nn+n15r do in=nbeg,nend

-------------------------------- (i) Update subroutine GRISE

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (MCB-B): fix gradual rise factor for the PRIME c module. The full rise (without c streamline depression in a building c wake) was divided by the final rise c (with rather than without streamline c depression).

BEFORE change:-------------- & (htrpt1(i2,isnum)-htrpt1(i1,isnum)) if(zfrise.GT.0.0) risefac=rise/zfrise return else htgrise=zfinal(ipnum) return endif AFTER change:

-------------- & (htrpt1(i2,isnum)-htrpt1(i1,isnum)) c 051216 (MCB-B) c --- Look at full final rise without streamline depression c !!! if(zfrise.GT.0.0) risefac=rise/zfrise riselast=htrpt1(ntrpt1(isnum),isnum) if(riselast.GT.0.0) risefac=rise/riselast goto 999 else htgrise=zfinal(ipnum) goto 999 endif

BEFORE change:

-------------- & (htrpt2(i2,isnum)-htrpt2(i1,isnum)) if(zfrise.GT.0.0) risefac=rise/zfrise return else htgrise=zfinal(ipnum) return endif elseif(istype.LE.2 .AND. idw0(ipnum).EQ.4) then c --- POINT source with PRIME downwash (cavity source) c --- No rise! htgrise=zfinal(ipnum) return AFTER change:

-------------- & (htrpt2(i2,isnum)-htrpt2(i1,isnum)) c 051216 (MCB-B) c --- Look at full final rise without streamline depression c !!! if(zfrise.GT.0.0) risefac=rise/zfrise riselast=htrpt2(ntrpt2(isnum),isnum) if(riselast.GT.0.0) risefac=rise/riselast goto 999 elsehtgrise=zfinal(ipnum) goto 999 endif elseif(istype.LE.2 .AND. idw0(ipnum).EQ.4) then c --- POINT source with PRIME downwash (cavity source) c --- No rise! htgrise=zfinal(ipnum) goto 999

BEFORE change:

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

return end AFTER change:

-------------- c 051216 (MCB-B) c --- Check for valid risefac 999 if(risefac.LT.-0.0001) then write(io6,*) write(io6,*)'Fatal Error in GRISE' write(io6,*)'Computed risefac is less than 0.0' write(io6,*)' x,htgrise = ',x,htgrise write(io6,*)' zfrise,risefac = ',zfrise,risefac stop 'Halted in GRISE -- see list file' else risefac=AMAX1(risefac,0.0) endif

return end

-------------------------------- (j) Update subroutine SIGTY

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (DGS): add WARN to collect information on c the largest (most) negative c virtual time/distance; c reset very small negative virtual c time/distance to zero

BEFORE change:

--------------include 'params.puf' include 'csigma.puf'
c real c(6),d(6)
c AFTER change:

-------------- include 'params.puf' include 'csigma.puf'
c real c(6),d(6)
c c 051216 (MCB-B) c --- Declare variables to test negative virtuals c --- (threshold for warning, threshold for HALT) real zerosec(2), zerokm(2)

BEFORE change:

-------------- c --- Make sure travel time is positive. if(dt.lt.zero) then write(io6,*)'ERROR in SUBR. SIGTY -- Incremental travel ', 1 'time < 0 -- DT = ',dt write(*,*) stop 'Halted in SIGTY -- see list file.' endif
c c --- Make sure downwind distance is positive. if(dxkm.lt.zero) then write(io6,*)'ERROR in SUBR. SIGTY -- Incremental travel ', 1 'distance < 0 -- DXKM = ',dxkm,' IDOPTY = ',idoptys, 2 ' IRU = ',irus write(*,*) stop 'Halted in SIGTY -- see list file.' endif AFTER change:

-------------- c 051216 (MCB-B) c --- Set values to test for negative virtuals that are treated as zero c --- (threshold for warning, threshold for HALT) data zerosec/-.01,-1.0/, zerokm/-0.00001,-.001/
c c --- Make sure travel time is positive. if(dt.lt.zero) then if(dt.ge.zerosec(1)) then dt=zero elseif(dt.ge.zerosec(2)) thencall WARN('SIGTY-t ',dt) dt=zero else write(io6,*)'ERROR in SUBR. SIGTY -- Incremental travel ', 1 'time < 0 -- DT = ',dt write(*,*) stop 'Halted in SIGTY -- see list file.' endif endif
c c --- Make sure downwind distance is positive. if(dxkm.lt.zero) then if(dxkm.ge.zerokm(1)) then dxkm=zero elseif(dxkm.ge.zerokm(2)) then call WARN('SIGTY-x ',dxkm) dxkm=zero else write(io6,*)'ERROR in SUBR. SIGTY -- Incremental travel ', 1 'distance < 0 -- DXKM = ',dxkm,' IDOPTY = ',idoptys, 2 ' IRU = ',irus write(*,*) stop 'Halted in SIGTY -- see list file.' endif endif

-------------------------------- (k) Update subroutine SIGTZ

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (MCB-B): reset very small negative virtual c time/distance to zero

BEFORE change:

-------------- c --- Make sure travel time is positive. if(dt.lt.zero) then write(io6,*)'SR. SIGTZ: FATAL ERROR: dt = ',dt write(*,*) stop 'Halted in SIGTZ -- see list file.' endif
c c --- Make sure downwind distance is positive. if(dxkm.lt.zero) then write(io6,*)'SR. SIGTZ: FATAL ERROR: dxkm = ',dxkm write(*,*) stop 'Halted in SIGTZ -- see list file.'endif AFTER change:

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

c c 051216 (MCB-B) c --- Set values to test for negative virtuals that are treated as zero data zerosec/-.01/, zerokm/-0.00001/
c c --- Make sure travel time is positive. if(dt.lt.zero) then if(dt.ge.zerosec) then dt=zero else write(io6,*)'SR. SIGTZ: FATAL ERROR: dt = ',dt write(*,*) stop 'Halted in SIGTZ -- see list file.' endif endif
c c --- Make sure downwind distance is positive. if(dxkm.lt.zero) then if(dxkm.ge.zerokm) then dxkm=zero else write(io6,*)'SR. SIGTZ: FATAL ERROR: dxkm = ',dxkm write(*,*) stop 'Halted in SIGTZ -- see list file.' endif endif

-------------------------------- (l) Update subroutine PDFPATH

-------------------------------- c --- UPDATE c --- V5.0-V5.711b 051216 (MCB-B): Add test for a zero lofting rise c adjustment to avoid a potential c attempt to divide by zero.

BEFORE change:

-------------- c --- Set lofting rise adjustment here delh=t*dhfac/SQRT(risq+ryz*t**2)AFTER change:

-------------- c --- Set lofting rise adjustment here c 051216 (MCB-B) c --- delh=t*dhfac/SQRT(risq+ryz*t**2) xdum1=t*dhfac xdum2=risq+ryz*t**2 if(xdum1.LE.0.0) then delh=0.0 elseif(xdum2.LE.0.0) then write(*,*)'PDFPATH ERROR: Illegal Calculation' write(*,*)'risq,ryz,t = ',risq,ryz,t write(*,*)'risq+ryz*t**2 = ',xdum2 stop else delh=xdum1/SQRT(xdum2) endif

-------------------------------- (m) Update subroutine TIBLON

-------------------------------- c --- UPDATE c --- V5.7-V5.711b 051216 (MCB-B): Augment test for onshore trajectory c to exclude parallel lines c 051216 (MCB-B): Place water cell mixing ht (not ht0) c into TIBL array locations upwind of c the coast c 051216 (MCB-B): Use land cell pointer in call c to TIBLGRO so that shore properties c are consistent with coastline as puff c crosses coast

BEFORE change:

-------------- c --- Coast Vector [CROSS] Puff Vector: cxp=dxc*dyp-dxp*dyc if(cxp.GT.0.0) then c --- Trajectory is onshore, find intercept AFTER change:

--------------c --- Coast Vector [CROSS] Puff Vector: cxp=dxc*dyp-dxp*dyc c 051216 (MCB-B) c --- Screen out equal slopes (cxp=0 c --- if slopes are equal, but precision may c --- let this slip through) if(cxp.GT.0.0 .AND. slope.NE.ymxpb(1,is)) then c --- Trajectory is onshore, find intercept

BEFORE change:

-------------- c --- Compute TIBL ht during step tiblku=tiblk*u0 call TIBLGRO(ldb,xs,ys,xe,ye,ixs,iys,ixe,iye,ht1, & -.001,tiblku) c --- Case 2: Coast is within the step else c --- First entry is interval to coast ntibl=1 tstibl(ntibl)=AMIN1(frac,1.0) htibl(ntibl)=ht0 ixtibl(ntibl)=ixw iytibl(ntibl)=iyw c --- Remaining entries filled if needed if(frac.LT.1.0) then tiblku=tiblk*u0 call TIBLGRO(ldb,xs,ys,xe,ye,ixs,iys,ixe,iye,ht0, & frac,tiblku) endif AFTER change:

-------------- c --- Compute TIBL ht during step tiblku=tiblk*u0 c 051216 (MCB-B) call TIBLGRO(ldb,xs,ys,xe,ye,ixl,iyl,ht1, & -.001,tiblku) c --- Case 2: Coast is within the step else c --- First entry is interval to coast ntibl=1 tstibl(ntibl)=AMIN1(frac,1.0) c 051216 (MCB-B) htibl(ntibl)=ziwater ixtibl(ntibl)=ixw iytibl(ntibl)=iyw c --- Remaining entries filled if needed (use land cell here) if(frac.LT.1.0) thentiblku=tiblk*u0 c 051216 (MCB-B) call TIBLGRO(ldb,xs,ys,xe,ye,ixl,iyl,ht0, & frac,tiblku) endif

-------------------------------- (n) Update subroutine TIBLGRO

-------------------------------- c --- UPDATES c --- V5.711-V5.711b 051216 (MCB-B): Change puff start/end cell indices c to a cell override index ixr,iyr c (-1,-1 if no override) c (MCB-B): Test for small/negative heat flux c during step

BEFORE change:

-------------- subroutine tiblgro(ldb0,xs,ys,xe,ye,ixs,iys,ixe,iye,ht0, & frac0,tiblku0) AFTER change:

-------------- subroutine tiblgro(ldb0,xs,ys,xe,ye,ixr,iyr,ht0, & frac0,tiblku0)

BEFORE change:

-------------- ixo=xold+1 ixo=MAX(1,ixo) ixo=MIN(nx,ixo) iyo=yold+1 iyo=MAX(1,iyo) iyo=MIN(ny,iyo) ixn=xnew+1 ixn=MAX(1,ixn) ixn=MIN(nx,ixn) iyn=ynew+1 iyn=MAX(1,iyn) iyn=MIN(ny,iyn) c --- Set cell for surface properties (choose overland cell if c --- choice is needed)ix=ixo iy=iyo if(ixo.NE.ixn .OR. iyo.NE.iyn) then if(ilandu(ixo,iyo).LT.iwat1 .OR. & ilandu(ixo,iyo).GT.iwat2) then ix=ixo iy=iyo elseif(ilandu(ixn,iyn).LT.iwat1 .OR. & ilandu(ixn,iyn).GT.iwat2) then ix=ixn iy=iyn endif endif c --- Obtain cell properties AFTER change:

-------------- c 051216 (MCB-B) c --- Set cell for surface properties (choose overland cell if c --- choice is needed) if(ixr.GT.0 .AND. iyr.GT.0) then ix=ixr iy=iyr else ixo=xold+1 ixo=MAX(1,ixo) ixo=MIN(nx,ixo) iyo=yold+1 iyo=MAX(1,iyo) iyo=MIN(ny,iyo) ixn=xnew+1 ixn=MAX(1,ixn) ixn=MIN(nx,ixn) iyn=ynew+1 iyn=MAX(1,iyn) iyn=MIN(ny,iyn) ix=ixo iy=iyo if(ixo.NE.ixn .OR. iyo.NE.iyn) then if(ilandu(ixo,iyo).LT.iwat1 .OR. & ilandu(ixo,iyo).GT.iwat2) then ix=ixo iy=iyo elseif(ilandu(ixn,iyn).LT.iwat1 .OR. & ilandu(ixn,iyn).GT.iwat2) then ix=ixn iy=iyn endif endif endif c --- Obtain cell propertiesBEFORE change:

-------------- c --- Compute the effective distance from coast to reach current ht c --- using this K-factor s0=0.5*hold**2/tiblk0 AFTER change:

-------------- c --- Compute the effective distance from coast to reach current ht c --- using this K-factor c 051216 (MCB-B) c --- Note: Heat flux may be zero or negative here (set s0 large) c --- Use sensible heat flux < 5 W/m**2 as cut-off if(hbyrcp.GT.0.004) then s0=0.5*hold**2/tiblk0 else s0=x2zi1 endif

BEFORE change:

-------------- write(io6,*)'TIBLGRO --------------' write(io6,*)'Puff moves between cells:' write(io6,*)'ixs,iys,ixe,iye: ',ixs,iys,ixe,iye AFTER change:

-------------- write(io6,*)'TIBLGRO --------------' c 051216 (MCB-B) write(io6,*)'Replacement cell for Met properties ---' write(io6,*)'ixr,iyr : ',ixr,iyr

-------------------------------- (o) Update subroutine WAKE_DFSN

-------------------------------- c --- UPDATE c --- V5.711-V5.711b 051216 (MCB-B): Constrain virtual distance increment c XVWAK to be positive

BEFORE change:

--------------c --- Set virtual distance as increment xzvwak=xzvwak-xasrc xyvwak=xyvwak-xasrc AFTER change:

-------------- c --- Set virtual distance as increment xzvwak=xzvwak-xasrc xyvwak=xyvwak-xasrc c 051216 (MCB-B) c --- Constrain virtual distance to be positive c --- (Error if more negative than -1m) if(xzvwak.LE.-1. .OR. xyvwak.LE.-1.) then write(io6,*) 'WAKE_DFSN ERROR: Negative virtual distance' write(io6,*) 'xzvwak = ',xzvwak write(io6,*) 'xyvwak = ',xyvwak stop 'ERROR in WAKE_DFSN --- see list file' else xzvwak=AMAX1(zero,xzvwak) xyvwak=AMAX1(zero,xyvwak) endif

-------------------------------- (p) Update subroutine COMP

-------------------------------- c --- UPDATE c --- V5.71-V5.711b 051216 (MCB-B): TIBLGRO arguments updated and c heat flux at start of step is used c to screen TIBLGRO call c 051216 (MCB-B): Change dummy volume source array from c emvol(mxnx,mxny,mxspec) to c em4dat(mxp6,mxvol) --- array space c had been incorrect, but sufficient

BEFORE change:

-------------- real em2dat(mxp4,mxpt2),emvol(mxnx,mxny,mxspec) AFTER change:

--------------c 051216 (MCB-B) real em2dat(mxp4,mxpt2),em4dat(mxp6,mxvol)

BEFORE change:

-------------- call INITPUF(ndathr,nspec,dthr,ldb2,em2dat,emvol,em3dat, & em5dat,em5grp,metfm,newpuf) AFTER change:

-------------- c 051216 (MCB-B) call INITPUF(ndathr,nspec,dthr,ldb2,em2dat,em4dat,em3dat, & em5dat,em5grp,metfm,newpuf)

BEFORE change:

-------------- c --- Puff currently interacting with TIBL (ZITIBL>0) AFTER change:

-------------- c --- Puff currently interacting with TIBL (ZITIBL>0) c 051216 (MCB-B) c --- Do not override current land use cell index (TIBLGRO) ixr=-1 iyr=-1

BEFORE change:

-------------- hbyrcp=(wstr*wstr*wstr)*temps/ & (9.8*htmix(ixs,iys)) if(hbyrcp.LE.0.004) then AFTER change:

-------------- c 051216 (MCB-B) c --- Test for heat flux using met data at start of step, c --- and form based on u-star and M-O length c --- to conform to what is done within TIBLGRO ustars=AMAX1(0.0,ustar(ixs,iys)) hbyrcp=-ustars**3*temps/(9.8*0.4*xmonin(ixs,iys)) if(hbyrcp.LE.0.004) then

BEFORE change:

--------------c --- Compute TIBL arrays for this step ltibl=.TRUE. call TIBLGRO(ldbhr,xold,yold,xnew,ynew,ixs,iys, & ixend,iyend,zitibl(ii),-.001,-1.) c --- Store TIBL ht for next sampling step for this puff AFTER change:

-------------- c --- Compute TIBL arrays for this step ltibl=.TRUE. c 051216 (MCB-B) call TIBLGRO(ldbhr,xold,yold,xnew,ynew,ixr,iyr, & zitibl(ii),-.001,-1.) c --- Store TIBL ht for next sampling step for this puff

--------------------------------------- (q) Insert new subroutine WARN

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

c----------------------------------------------------------------------- subroutine warn(astring,value)
c----------------------------------------------------------------------- c c --- CALPUFF Version: 5.711b Level: 051216 WARN c D. Strimaitis Earth Tech, Inc.
c c --- PURPOSE: Screens data associated with potential problems and c reports outcome to list file (diagnostic)
c c --- INPUTS: c ASTRING - char*12 - subroutine/action string c VALUE - real - value to process
c c --- Parameters: IO6
c c --- OUTPUT: c (List File)
c c --- WARN called by: SIGTY, FIN c --- WARN calls: none
c---------------------------------------------------------------------- c c --- Include parameters include 'params.puf' logical LWARN character*12 astringc --- Initialize via data statements data lwarn/.FALSE./ data sigtysec/0.0/, ntsigty/0/ data sigtykm/0.0/, nxsigty/0/ if(astring.EQ.'SIGTY-t ') then c --- Process negative delta-t from SIGTY if(value.LT.sigtysec) sigtysec=value ntsigty=ntsigty+1 ntsigty=MIN(ntsigty,1000001) lwarn=.TRUE. elseif(astring.EQ.'SIGTY-x ') then c --- Process negative delta-x from SIGTY if(value.LT.sigtykm) sigtykm=value nxsigty=nxsigty+1 nxsigty=MIN(nxsigty,1000001) lwarn=.TRUE. elseif(astring.EQ.'FIN ') then c --- Report diagnostic if(LWARN) then write(io6,*) write(io6,*) write(io6,*) write(io6,*)'Diagnostic for SIGTY:' write(io6,*)' Minimum dt(sec) = ',sigtysec if(ntsigty.GT.100000) then write(io6,*)' Number less than -.01s > ',ntsigty else write(io6,*)' Number less than -.01s = ',ntsigty endif write(io6,*)' Minimum dx(km) = ',sigtykm if(nxsigty.GT.100000) then write(io6,*)' Number less than -.01m > ',nxsigty else write(io6,*)' Number less than -.01m = ',nxsigty endif write(io6,*) write(io6,*) endif

else

c --- Problem here write(*,*) write(*,*)'ERROR in subroutine WARN!' write(*,*)'Unexpected arguments:' write(*,*)' arg1 = ',astring write(*,*)' arg2 = ',valuewrite(*,*) stop
endif

return end