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

 --- 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, RDOW
Problem 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,ykm
AFTER 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)')cset3d
AFTER 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) then
 read(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 from
 if(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') then
 read(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)then
 write(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,nspec
 xlam(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)**2
 syr=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 old
BEFORE 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

 else
 htgrise=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)) then
 call 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) then
 tiblku=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 properties
BEFORE 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 astring
c --- 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 = ',value
 write(*,*)

 stop

 endif

 return

 end