----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : A Dated : JULY 16, 2004 ----------------------------------------------------------------------- The following code modifications address problems identified in 4 modules of the CALPUFF modeling system: CALMET, CALPUFF, READ62, and SMERGE. 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 "a" to the end of the version to indicate that the program contains the modifications identified in Model Change Bulletin A, and replace the level with the date of this bulletin. ----------------------------------------------------------------------- ----------------------------------------------------------------------- 1. CALMET (subroutines RDMM4,STULL,READCF,INTERPQR) ----------------------------------------------------------------------- Edit PARAMS.MET, PARAMSS.MET, PARAMSL.MET, PARAMSX.MET Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- c --- Specify model version character*8 mver,mlevel parameter(mver='5.53',mlevel='030709') AFTER change: -------------- c --- Specify model version character*8 mver,mlevel parameter(mver='5.53a',mlevel='040716') Edit CALMET.FOR ------------------------------------- (a) Update the MAIN program documentation ------------------------------------- Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- c c --- set version and level number of program ver='5.53' level='030709' c AFTER change: -------------- c c --- set version and level number of program ver='5.53a' level='040716' c --------------------------------------- (b) Update subroutine RDMM4 --------------------------------------- When cloud cover is estimated from humidity data in the MM4.DAT file, the relative humidity near 850mb is selected. The code had been written for pressure in tenths of millibars, but the pressure has been converted to millibars, so level 1 of the RH profile was selected instead. Reset the pressure check from 8500 to 850 (millibars). BEFORE change: -------------- c --- Store RH at (near) 850 mb for cloud cover computation c Allow for surface pressure < 850mb (happens e.g. SW Wyoming) if (n.eq.1 .or. press.gt.8500.)rh850(i,j)=rh1 AFTER change: -------------- c 040716 (MCB-A) c --- Store RH at (near) 850 mb for cloud cover computation c Allow for surface pressure < 850mb (happens e.g. SW Wyoming) if (n.eq.1 .or. press.gt.850.)rh850(i,j)=rh1 ----------------------- (c) Update subroutine STULL ----------------------- A possible exponential overflow can stop the CALMET run with an error message. Introduce an upper bound on the variable H. BEFORE change: -------------- c if Z2 and Z3 straddle the inversion, make assumption on inversion height: c assume H=z3/5 (empirical, based on inversion height~5H and z3~h) c Make this assumption also if dt2 or dt3 not "well-behaved" else H=z3/5. endif c surface temperature drop since sunset dts=min(0.,dt2*exp(alpha*z2/H)) AFTER change: -------------- c if Z2 and Z3 straddle the inversion, make assumption on inversion height: c assume H=z3/5 (empirical, based on inversion height~5H and z3~h) c Make this assumption also if dt2 or dt3 not "well-behaved" else H=z3/5. endif c 040716 (MCB-A) c Bound H to avoid exponential overflow H=max(H,z2/10.) c surface temperature drop since sunset dts=min(0.,dt2*exp(alpha*z2/H)) ------------------------ (d) Update subroutine READCF ------------------------ Incorrect variable type definitions of IAVET and TGDEFA cause the upwind temperature averaging option to be OFF, and the default temperature gradient above the mixing height over water to be ZERO. Two changes are needed in continuation line number 6 for the IVTYPE array. BEFORE change: -------------- data ivtype/7*2,3,2,51*0, 2 2*4,2*1,2,5*4,2*2,3*1,2,1, 43*0, 3 2*3, 12*2, 3, 14*2, 31*0, 4 7*2, 53*0, 5 5*2, 2*1, 2*2, 3, 9*1, 3*2, 2*1, 2, 4*1, 7*2, 1, 2, 1, 3, 2, 5 8*1, 2*2, 1, 2*2, 6*0, 6 9*1, 2*2, 1, 2, 2*1, 2*2, 2*1, 3*2, 1, 2*2, 9*1, 26*0, 7 60*0, 60*0, 60*0/ AFTER change: -------------- c 040716 (MCB-A) data ivtype/7*2,3,2,51*0, 2 2*4,2*1,2,5*4,2*2,3*1,2,1, 43*0, 3 2*3, 12*2, 3, 14*2, 31*0, 4 7*2, 53*0, 5 5*2, 2*1, 2*2, 3, 9*1, 3*2, 2*1, 2, 4*1, 7*2, 1, 2, 1, 3, 2, 5 8*1, 2*2, 1, 2*2, 6*0, c 6 9*1, 2*2, 1, 2, 2*1, 2*2, 2*1, 3*2, 1, 2*2, 9*1, 26*0, 6 9*1, 2*2, 1, 2, 2*1, 3*2, 2*1, 2*2, 1, 2*2, 9*1, 26*0, 7 60*0, 60*0, 60*0/ -------------------------- (e) Update subroutine INTERPQR -------------------------- Use consistent coordinate systems for MM5 and CALMET grid points. Precipitation assignment from the MM5 grid to CALMET grid points did not use relative CALMET grid coordinates. Subtract the CALMET grid origin from the MM5 coordinates. BEFORE change: -------------- np=np+1 x(np)=xlcmm4(ip,jp) y(np)=ylcmm4(ip,jp) p(np)=qrprog(ip,jp) AFTER change: -------------- np=np+1 c 040716 (MCB-A) x(np)=xlcmm4(ip,jp)-xmap0 y(np)=ylcmm4(ip,jp)-ymap0 p(np)=qrprog(ip,jp) ----------------------------------------------------------------------- 2. CALPUFF (subroutines XERFDIF,CALCSL) ----------------------------------------------------------------------- Edit PARAMS.PUF, PARAMSS.PUF, PARAMSL.PUF, PARAMSX.PUF Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- c --- Specify model version character*12 mver, mlevel, mmodel parameter(mver='5.711',mlevel='030625') AFTER change: -------------- c --- Specify model version character*12 mver, mlevel, mmodel parameter(mver='5.711a',mlevel='040716') Edit CALPUFF.FOR ------------------------------------- (a) Update the MAIN program documentation ------------------------------------- Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- c c --- set version and level number of program ver='5.711' level='030625' AFTER change: -------------- c c --- set version and level number of program ver='5.711a' level='040716' -------------------------------- (b) Update subroutine XERFDIF -------------------------------- The F1 integral in the analytic sampling done for attached slugs is not coded correctly. F1 is used only for slugs as they are being emitted (youngest end is at the source), and is multiplied by the change in effective emission rate during the sampling step (due to chemical transformation or deposition removal). BEFORE change: -------------- expxu = zero erfxu = one if(xu.lt.0.0) erfxu = -one c Note that for erfcut = 3, erf(erfcut) = 0.999978 if(abs(xu) .le. erfcut) then erfxu = erf(xu) expxu = exp(-xu * xu) endif c x0erf = xu * erfxu + expxu / srtpi x1erf = half * xu * x0erf - fourth * erfxu c c expxl = zero erfxl = one if(xl.lt.0.0) erfxl = -one c Note that for erfcut = 3, erf(erfcut) = 0.999978 if(abs(xl) .le. erfcut) then erfxl = erf(xl) expxl = exp(-xl * xl) endif c x0erf = x0erf - ( xl * erfxl + expxl / srtpi ) x1erf = x1erf - ( half * xl * x0erf - fourth * erfxl ) c AFTER change: -------------- expxu = zero erfxu = one if(xu.lt.0.0) erfxu = -one c Note that for erfcut = 3, erf(erfcut) = 0.999978 if(abs(xu) .le. erfcut) then erfxu = erf(xu) expxu = exp(-xu * xu) endif c expxl = zero erfxl = one if(xl.lt.0.0) erfxl = -one c Note that for erfcut = 3, erf(erfcut) = 0.999978 if(abs(xl) .le. erfcut) then erfxl = erf(xl) expxl = exp(-xl * xl) endif c c 040716 (MCB-A) c x0erf = (xu*erfxu + expxu/srtpi) - (xl*erfxl + expxl/srtpi) x1erf = (half*xu*(xu*erfxu + expxu/srtpi) - fourth*erfxu) - & (half*xl*(xl*erfxl + expxl/srtpi) - fourth*erfxl) c ------------------------ (c) Update subroutine CALCSL ------------------------ Fix array declaration in CALCSL: species mass array for wet deposition is declared using nspec (should be declared mxpec). This causes erroneous wet fluxes to be output from slug sampling. The total mass of each species in slugs is not affected so concentrations and dry deposition are not affected. BEFORE change: -------------- real qb(nspec),qe(nspec),qw(nspec,2) AFTER change: -------------- c 040716 (MCB-A) real qb(nspec),qe(nspec),qw(mxspec,2) ----------------------------------------------------------------------- 3. READ62 (subroutine COMP) ----------------------------------------------------------------------- Edit PARAMS.R62 Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- parameter(mver='5.52',mlevel='030709') AFTER change: -------------- parameter(mver='5.52a',mlevel='040716') Edit READ62.FOR ------------------------------------- (a) Update the MAIN program documentation ------------------------------------- Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- c --- Set version and level number of program (stored in /QA/ and c --- checked against values set in PARAMS.R62) ver='5.52' level='030709' AFTER change: -------------- c --- Set version and level number of program (stored in /QA/ and c --- checked against values set in PARAMS.R62) ver='5.52a' level='040716' ---------------------- (b) Update subroutine COMP ---------------------- De-activate lines using old limit of 79 levels from the TD-6201 section. These lines cause 1 or 2 dummy reads leading to "skipped soundings" in the output file. This does not affect FSL data. This change has been tested on the variable record length TD-6201 format, but not the fixed record length TD-6201 format. BEFORE change: -------------- (Change 1) c ALAT = LAT/100 + (LAT-(LAT/100*100))/60. R6201360 ALON = LON/100 + (LON-(LON/100*100))/60. R6201370 (Change 2) c 1050 IF(NUMLEV.LE.79) GO TO 1100 R6201410 READ(io8,6150,END=2000) JUNK R6201420 6150 FORMAT(A32,79(36X)) NUMLEV = MAX(NUMLEV-79,79) R6201430 GO TO 1050 R6201440 ELSEIF (JDAT.EQ.2) THEN AFTER change: -------------- (Change 1) c ALAT = LAT/100 + (LAT-(LAT/100*100))/60. R6201360 ALON = LON/100 + (LON-(LON/100*100))/60. R6201370 c 040716 (MCB-A) c --- Check for array dimensioning problems if(numlev.gt.mxlev)then write(io6,*)'ERROR -- too many levels in upper air data ', 1 'for array dimensions' write(io6,*)'YEAR: ',year,' MONTH: ',month,' DAY: ',day, 1 ' HOUR: ',hour write(io6,*)'NUMLEV in file = ',numlev write(io6,*)'MXLEV in array = ',mxlev stop endif (Change 2) c c 040716 (MCB-A) c1050 IF(NUMLEV.LE.79) GO TO 1100 R6201410 c READ(io8,6150,END=2000) JUNK R6201420 c6150 FORMAT(A32,79(36X)) c NUMLEV = MAX(NUMLEV-79,79) R6201430 c GO TO 1050 R6201440 ELSEIF (JDAT.EQ.2) THEN ----------------------------------------------------------------------- 4. SMERGE (subroutine RDWRITS) ----------------------------------------------------------------------- Edit PARAMS.SMG Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- parameter(mver='5.31',mlevel='030528') AFTER change: -------------- parameter(mver='5.31a',mlevel='040716') Edit SMERGE.FOR ------------------------------------- (a) Update the MAIN program documentation ------------------------------------- Change the version by adding the letter "a" Change the level to 040716 (example) BEFORE change: -------------- c --- Set version and level number of program ver='5.31' level='030528' AFTER change: -------------- c --- Set version and level number of program ver='5.31a' level='040716' ------------------------- (b) Update subroutine RDWRITS ------------------------- Correct the treatment of CD-144 "0" wind direction. A wind direction of ZERO is always invalid. BEFORE change: -------------- c Check for missing data if (cwd.eq.' ')then xwd(index) = 9999. xws(index) = 9999. ichkws = 1 else read (cwd,390)jwd read (cws,390)jws 390 format (i2) if(jws.GT.0)then xwd(index) = jwd * 10. else xwd(index) = 0. endif if (xwd(index).lt.0.0.or.xwd(index).gt.360.0) then xwd(index) = 9999. endif endif AFTER change: -------------- c Check for missing data if (cwd.eq.' ')then xwd(index) = 9999. xws(index) = 9999. ichkws = 1 else read (cwd,390)jwd read (cws,390)jws 390 format (i2) c 040716 (MCB-A) if(jws.GT.0)then if(jwd.GT.0) then xwd(index) = jwd * 10. else xwd(index) = 9999. xws(index) = 9999. ichkws = 1 endif else xwd(index) = 0. endif if (xwd(index).lt.0.0.or.xwd(index).gt.360.0) then xwd(index) = 9999. endif endif