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