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