----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : C Dated : AUGUST 4, 2006 ----------------------------------------------------------------------- The following code modifications address problems identified in 2 modules of the CALPUFF Version 5 modeling system: CALMET and CALPUFF. The codes were updated in the VISTAS Recommended Version of the system and posted for distribution on August 4, 2006. Bulletin MCB-C (060804) outlines those changes and together with those in MCB-A (040716) and MCB-B (051216) represent all of the corrections made to Version 5 of of these programs since EPA first approved the CALPUFF system for regulatory applications. Related components used in the CALPUFF system for coordinate conversions have also been updated. Although these are packaged as FORTRAN include files: COORDLIB.FOR BLOCKDAT.CRD NIMA.CRD they cannot be simply swapped for the updated versions in the August 4, 2006 distribution because of one or more changes in subroutine arguments. The updates include: -- adding 5 local datums for Hawaii -- adding the WGS-72 and EMG-96 global ellipsoids -- removing the sphere based on NAD27 (the NAD27 ellipsoid remains valid) -- restricting the Lambert Azimuthal projection to sphere datums -- halting any UTM projection that is combined with the NWS-84 sphere datum -- adding the Albers Equal-Area map projection -- fixing problems with some UTM coordinate conversions that cross the equator -- recognize a projection parameter change (projection type and datum constant) These updates are not included in MCB-C. Changes to CALPUFF ------------------ Problem Area 1 -- A subroutine that computes the potential temperature gradient between 2 heights contained a bug that left an array index undefined (zero) if the upper height requested exceeds the middle of the top layer in the domain. This either results in a computed overflow error that stops the run, or in an erroneous potential temperature gradient. Computed gradients that are less than the minimum are replaced by the minimum, which is equivalent to a near-neutral atmosphere. The subroutine is called to define the potential temperature gradient: -- across the puff-depth -- at the top of the mixed layer for partial plume penetration -- in a 100m layer above stack-top -- at the top of the overwater mixing height -- at the top of a TIBL mixing height Of these, only the first is expected to involve an upper height that can exceed the mid-point of the uppermost layer, and the gradient in this case is used for the CTDM sigma-z growth option, and for terrain adjustment option MCTADJ=2 (simple CALPUFF adjustment). Modified: PTLAPS Changes to CALMET ------------------ Problem Area 1 -- Upwind averaging of mixing height and temperature includes influence of cells outside of the intended upwind cone. This issue was noted in MCB-B but no detailed code changes were included. Modified: AVEMIX, AVETMP Problem Area 2 -- The second method of vertical extrapolation of surface wind observations (power law- iextrp=2) relies on the grid cell elevation rather than the grid cell land use category to determine whether to use overland or overwater powerlaw exponents. Modified: DIAGNO Problem Area 3 -- The third method of vertical extrapolation of surface wind observations (user's extrapolation multipliers - iextrp=3) also extrapolates missing winds Modified: DIAGNO Problem Area 4 -- The solar angle array is used to determine night and day and the associated methods of vertical extrapolation of the prognostic temperatures down to the surface. In this particular night-day test, the solar angle array, which is defined on the CALMET grid, is used on the prognostic (MM4/MM5/3D) grid (i.e. at the wrong place), with solar angle values off by 1 hour (i.e. dusk/dawn occur one hour too early) Modified: RDHD4, RDHD5, MM4HD0.MET,RDMM4, RDMM5 Problem Area 5 -- The second MM5/3D record required for the time interpolation of prognostic fields to the calmet timestep is not read in if the time of the of the first record in the MM5/3D file coincides with the CALMET simulation beginning time. Modified: RDMM5 Problem Area 6 - In rare instances when prognostic records have duplicate lowest level altitudes the vertical extrapolation of the prognostic temperature down to the surface fails. Modified: RDMM4, RDMM5 Problem Area 7 - Some compilers do not keep the values of 2 internal variables in memory from one step to the next but the code assumes those values are memorized. Modified: RDMM5 Problem Area 8 -- In NOOBS-temperature mode (ITPROG>0), the mixing height growth is dictated by the current hour, rather than the previous hour, lapse rate at the top of the previous hour mixing height, which can lead to excessive growth of the mixing height. Modified: RDMM4, RDMM5 Problem Area 9 -- The surface wind observations are assumed to be measured at the first CALMET level, no matter what the anemometer height actually is. Modified: MET1.MET,OVRWAT.MET,SETUP,DIAGI,DIAGNO,RDOW,ELUSTR2,WATER2,WIND1,SIMILT Problem Area 10 -- Units of meters and kilometers are mixed up in the computation of the precipitation rate when prognostic rain data are used. Modified: INTERPQR Problem Area 11 -- Missing variable definition for output in debug NOOBS mode Modified: MIXHT2ST Problem Area 12 --- The vertical extrapolation of prognostic temperatures from the lowest prognostic level down to the surface assumes that the surface cools during the night as land typically does. However the same type of cooling is assumed over water. (Note: bug fix takes into account datasets available in V5.53 version i.e. MM4 and old MM5 datasets not including SST data- Fix for later datasets available as enhancement adjustment in VISTAS version 5.726) Modified: CALMET.INP, MM4HDO, READCF , RDHD52, RDMM4, RDMM5, RDHD53 Problem Area 13 --- The 2D-array of surface temperature in CALMET.DAT is not consistent with the surface layer of the 3D-temperature array. The latter has better temperature interpolation and upwind averaging. However CALPUFF uses the 2-D array rather than the first layer of the 3D array Modified: COMP Problem Area 14 --- Current GMT time should always be computed and not be conditional upon the value of ITPROG in subroutine COMP. Modified: COMP ----------------------------------------------------------------------- ----------------------------------------------------------------------- 1. CALPUFF ----------------------------------------------------------------------- Edit CALPUFF.FOR -------------------------------- (a) Update subroutine PTLAPS -------------------------------- BEFORE change: -------------- c --- Find level just above htbase+delz htpdz=htbase+delz klp1=kl+1 do 15 k=klp1,nz if(zgpt(k).lt.htpdz) go to 15 ku=k go to 16 15 continue 16 continue AFTER change: -------------- c --- Find level just above htbase+delz c --- (MCB-C) c --- Initialize to NZ ku=nz htpdz=htbase+delz klp1=kl+1 do 15 k=klp1,nz if(zgpt(k).lt.htpdz) go to 15 ku=k go to 16 15 continue 16 continue ----------------------------------------------------------------------- 2. CALMET ----------------------------------------------------------------------- Edit CALMET.FOR xxxxxxxxxxxxxxxx Problem Area 1 -- xxxxxxxxxxxxxxxx --------------------------------------- (a) Update subroutine AVEMIX --------------------------------------- c --- UPDATE c --- Add a check to ensure upwind averaging of mixing heights c --- is done over gridpoints within cone of influence and over c --- mnmdav square around each gridcell only BEFORE change: -------------- ja = ju - icwy jb = ju + icwy c c --- Now set the search window: ki from kilow to kihgh c kj from kjlow to kjhgh AFTER change: -------------- ja = ju - icwy jb = ju + icwy c c --- MCB-C:060804: Wind direction +/- halfang if (i.ne.ia) then alphaa=270.-atan2(float(j-ja),float(i-ia))/degrad alphaa=amod(alphaa,360.) else if(j.ge.ja) alphaa=180. if(j.lt.ja) alphaa=0. endif if (i.ne.ib) then alphab=270.-atan2(float(j-jb),float(i-ib))/degrad alphab=amod(alphab,360.) else if(j.ge.jb) alphab=180. if(j.lt.jb) alphab=0. endif c --- Now set the search window: ki from kilow to kihgh c kj from kjlow to kjhgh BEFORE change: -------------- c --- The local cell always gets a weight of one if(irel.eq.0 .and. jrel.eq.0) go to 45 c c DELA is the upwind measure of a cell's location, with (iu,ju) AFTER change: -------------- c --- The local cell always gets a weight of one if(irel.eq.0 .and. jrel.eq.0) go to 45 c c --- MCB-C:060804 - Check if the gridpoint is within mnmdav of (i,j) c --- (square) -If so, skip the check on the upwind cone (as might c --- be outside of the cone but should still be kept) if ((abs(irel).le.mnmdav).and.(abs(jrel).le.mnmdav)) goto 47 c --- MCB-C:060804 - Check if the gridpoint is within cone of influence c --- First check that within +/- halfang of wind direction c --- Direction of (ki,kj) relative to (i,j) beta=270.-atan2(float(jrel),float(irel))/degrad beta=amod(beta,360.) if ((beta.lt.alphaa).or.(beta.gt.alphab)) goto 50 c --- Then check than not beyond maximum upwind distance dist=irel**2+jrel**2 if (dist.gt.eta2) goto 50 47 continue c DELA is the upwind measure of a cell's location, with (iu,ju) c --------------------------------------- (b) Update subroutine AVETMP --------------------------------------- c --- UPDATE c --- Add a check to ensure upwind averaging of temperatures c --- is done over gridpoints within cone of influence and over c --- mnmdav square around each gridcell only BEFORE change: -------------- ja = ju - icwy jb = ju + icwy c c --- Now set the search window: ki from kilow to kihgh c kj from kjlow to kjhgh AFTER change: -------------- ja = ju - icwy jb = ju + icwy c c --- MCB-C:060804 Wind direction +/- halfang if (i.ne.ia) then alphaa=270.-atan2(float(j-ja),float(i-ia))/degrad alphaa=amod(alphaa,360.) else if(j.ge.ja) alphaa=180. if(j.lt.ja) alphaa=0. endif if (i.ne.ib) then alphab=270.-atan2(float(j-jb),float(i-ib))/degrad alphab=amod(alphab,360.) else if(j.ge.jb) alphab=180. if(j.lt.jb) alphab=0. endif c --- Now set the search window: ki from kilow to kihgh c kj from kjlow to kjhgh BEFORE change: -------------- c --- The local cell always gets a weight of one if(irel.eq.0 .and. jrel.eq.0) go to 45 c c DELA is the upwind measure of a cell's location, with (iu,ju) AFTER change: -------------- c --- The local cell always gets a weight of one if(irel.eq.0 .and. jrel.eq.0) go to 45 c c --- MCB-C:060804 - Check if the gridpoint is within mnmdav of (i,j) c --- (square) -If so, skip the check on the upwind cone (as might c --- be outside of the cone but should still be kept) if ((abs(irel).le.mnmdav).and.(abs(jrel).le.mnmdav)) goto 47 c --- MCB-C:060804 - Check if the gridpoint is within cone of influence c --- First check that within +/- halfang of wind direction c --- Direction of (ki,kj) relative to (i,j) beta=270.-atan2(float(jrel),float(irel))/degrad beta=amod(beta,360.) if ((beta.lt.alphaa).or.(beta.gt.alphab)) goto 50 c --- Then check than not beyond maximum upwind distance (050412) dist=irel**2+jrel**2 if (dist.gt.eta2)goto 50 47 continue c DELA is the upwind measure of a cell's location, with (iu,ju) xxxxxxxxxxxxxxxx Problem Area 2 -- xxxxxxxxxxxxxxxx --------------------------------------- (c) Update subroutine DIAGO --------------------------------------- c --- UPDATE c --- The second method of vertical extrapolation of surface wind observations c --- (power law - |IEXTRP|=2) uses different exponents overland and overwater. c --- For this option, the code identified overwater cells as grid cells with c --- zero elevations. However some land cells may have zero elevations and sea c --- level might not be at 0m. c --- The code change ensures that a water cell is defined by its land use category c --- and not by its elevation. BEFORE change: -------------- IF (HTOPO(IG,JG) .EQ. 0.) pexp2 = 0.286 AFTER change: -------------- c --- MCB-C (060804) c IF (HTOPO(IG,JG) .EQ. 0.) pexp2 = 0.286 if( ( (l.gt.nssta).and.(l.le.nsurf) ) .OR. : ( (ilandu(ig,jg).ge.iwat1).and.(ilandu(ig,jg).le.iwat2))) : pexp2=0.286 xxxxxxxxxxxxxxxx Problem Area 3 -- xxxxxxxxxxxxxxxx --------------------------------------- (d) Update subroutine DIAGO --------------------------------------- c --- UPDATE c --- Do not extrapolate missing wind observations when the third method c --- of vertical extrapolation of surface wind observations is selected c --- (user's extrapolation multipliers - |IEXTRP|=3) BEFORE change: -------------- ELSE IF (IABS(IEXTRP) .EQ. 3) THEN C C EXTRAPOLATE WITH USER'S EXTRAPOLATION MULTIPLIERS AFTER change: -------------- ELSE IF (IABS(IEXTRP) .EQ. 3) THEN C C EXTRAPOLATE WITH USER'S EXTRAPOLATION MULTIPLIERS c --- MCB-C (060804) Extrapolate only if not missing IF(US(1,L).GT.EDITL .OR. VS(1,L).GT.EDITL) GOTO 70 xxxxxxxxxxxxxxxx Problem Area 4 -- xxxxxxxxxxxxxxxx --------------------------------------- (e) Update include file MM4HDO.MET --------------------------------------- c --- UPDATE c --- Store 2-D array of nearest CALMET gridpoints to Mm4/MM5/3D gridpoints BEFORE change: -------------- 5 IGRAB(mxnx,mxny,4),JGRAB(mxnx,mxny,4), AFTER CHANGE ------------ 4 INEARG(mxnxp,mxnyp),JNEARG(mxnxp,mxnyp), 5 IGRAB(mxnx,mxny,4),JGRAB(mxnx,mxny,4), BEFORE change: -------------- c IGRAB(mxnx,mxny,4) - integer array - I index (1, 2, ... MXNXP) of four AFTER CHANGE ------------ c INEARG(mxnxp,mxnyp) - integer array - I index (1, 2, ... MXNX) of c closest CALMET grid point to c MM4/MM5/3D grid point c JNEARG(mxnxp,mxnyp) - integer array - J index (1, 2, ... MXNY) of c closest CALMET grid point to c MM4/MM5/3D grid point c IGRAB(mxnx,mxny,4) - integer array - I index (1, 2, ... MXNXP) of four --------------------------------------- (f) Update subroutine RDHD4 --------------------------------------- c --- UPDATE c --- Compute nearest CALMET gridpoint to all MM4 gridpoints BEFORE change: -------------- c Record first access to MM4 record ifirstpg=0 AFTER CHANGE ------------ c --- MCB-C (060804) c------------------------------------------------------------------- c --- Find the closest CALMET grid point to each MM4 c grid point (stored in MM4HDO.MET) c ------------------------------------------------------------------ do i = 1,nxp do j = 1,nyp nearii = 1 nearjj = 1 dnearg=9.9E19 do ii = 1,nx do jj = 1,ny xx = xcal + (ii - 1) * delg yy = ycal + (jj - 1) * delg pdist = sqrt ((xlcmm4(i,j) - xx) ** 2 + & (ylcmm4(i,j) - yy) **2) if (pdist .lt. dnearg) then dnearg = pdist nearii = ii nearjj = jj endif enddo enddo inearg(i,j) = nearii jnearg(i,j) = nearjj enddo enddo c Record first access to MM4 record (And first MM4.DAT file) ifirstpg=0 --------------------------------------- (g) Update subroutine RDHD5 --------------------------------------- c --- UPDATE c --- Compute nearest CALMET gridpoint to all MM5/3D gridpoints BEFORE change: -------------- c frr (09/01) check if prognostic precip data are available c NPSTA>0 : use observations (precip data) AFTER CHANGE ------------ c --- MCB-C (060804) c------------------------------------------------------------------- c --- Find the closest CALMET grid point to each MM5 c grid point (060213)v (stored in MM4HDO.MET) c ------------------------------------------------------------------ do i = 1,nxp do j = 1,nyp nearii = 0 nearjj = 0 dnearg=9.9E19 do ii = 1,nx do jj = 1,ny xx = xcal + (ii - 1) * delg yy = ycal + (jj - 1) * delg pdist = sqrt ((xlcmm4(i,j) - xx) ** 2 + & (ylcmm4(i,j) - yy) **2) if (pdist .lt. dnearg) then dnearg = pdist nearii = ii nearjj = jj endif enddo enddo inearg(i,j) = nearii jnearg(i,j) = nearjj enddo enddo c ------------------------------------------------------- c frr (09/01) check if prognostic precip data are available c NPSTA>0 : use observations (precip data) --------------------------------------- (h) Update subroutine RDMM4 --------------------------------------- c --- UPDATE c --- The solar angle array is defined on the CALMET grid but was used in c --- a loop on the MM4 gridpoints,i.e. at the wrong locations and possibly c --- exceeding array dimensions. The value at the CALMET gridpoint nearest c --- to the MM4 gridpoint of interest is now used. c --- Moreover the solar angle is defined for hours 0 to 23 with corresponding c --- array indexes between 1 to 24, a feature that was not taken into account BEFORE change: -------------- if (mhr.eq.0)then sina=sinalp(i,j,24) else sina=sinalp(i,j,mhr) endif AFTER CHANGE ------------ c --- MCB-C (060804) sina=sinalp(inearg(i,j),jnearg(i,j),mhr+1) --------------------------------------- (i) Update subroutine RDMM5 --------------------------------------- c --- UPDATE c --- The solar angle array is defined on the CALMET grid but was used in c --- a loop on the MM5 gridpoints,i.e. at the wrong locations and possibly c --- exceeding array dimensions. The value at the CALMET gridpoint nearest c --- to the MM5 gridpoint of interest is now used. The computation is done c --- only once per gridpoint (not for each level i.e. out of the k loop) c --- Moreover the solar angle is defined for hours 0 to 23 with corresponding c --- array indexes between 1 to 24, a feature that was not taken into account BEFORE change: -------------- C --- INTERPOLATE PROGNOSTIC SOUNDINGS VERTICALLY TO DIAGNOSTIC C --- MODEL LEVELS DO 75 K = 1,NZ AFTER CHANGE ------------ c --- MCB-C (060804) sina=sinalp(inearg(i,j),jnearg(i,j),mhr+1) C --- INTERPOLATE PROGNOSTIC SOUNDINGS VERTICALLY TO DIAGNOSTIC C --- MODEL LEVELS DO 75 K = 1,NZ BEFORE change: -------------- c At night: follow Stull (1983, Tellus 35Am p 219-230) c day/time determined by sinalp (passed on via gen.met) if (mhr.eq.0)then sina=sinalp(i,j,24) else sina=sinalp(i,j,mhr) endif c First record pot.temp before sunset for extrapoltion at night AFTER CHANGE ------------ c At night: follow Stull (1983, Tellus 35Am p 219-230) c day/time determined by sinalp (passed on via gen.met) c MCB-C (060804): sinalp computed before k loop c c First record pot.temp before sunset for extrapolation at night xxxxxxxxxxxxxxxx Problem Area 5 -- xxxxxxxxxxxxxxxx --------------------------------------- (j) Update subroutine RDMM5 --------------------------------------- c --- UPDATE c --- Read in a second MM5/3D record when the MM5/3D file starts at the CALMET c --- simulation beginning time. The second record is necessary to perform c --- the time interpolation of MM5/3D fields to the CALMET timestep BEFORE change: -------------- if (mdathr .lt. (ndathr-isteppg+1) ) iskip = 1 AFTER CHANGE ------------ c --- MCB-C (060804) if (mdathr.le. (ndathr-isteppg+1) ) iskip = 1 if (mdathr.eq. (ndathr-isteppg+1) .and. ifirstpg.le.1) iskip=0 xxxxxxxxxxxxxxxx Problem Area 6 -- xxxxxxxxxxxxxxxx --------------------------------------- (k) Update subroutine RDMM4 --------------------------------------- c --- UPDATE: c --- Check to ensure that the 2 temperatures used to extrapolate c --- temperatures down to the surface are at different altitudes BEFORE change: -------------- if (ifirstpg.eq.0) then call stull0(vptp(i,j,index) ,zp(i,j,index), : vptp(i,j,index+1),zp(i,j,index+1), : pt20(i,j),pt30(i,j)) endif AFTER CHANGE ------------ if (ifirstpg.eq.0) then c --- MCB-C (060804) if (zp(i,j,index).ne.zp(i,j,index+1))then index1=index+1 else index1=index+2 endif call stull0(vptp(i,j,index) ,zp(i,j,index), : vptp(i,j,index1),zp(i,j,index1), : pt20(i,j),pt30(i,j)) endif BEFORE change: -------------- c frr 030119 Extrapolate temperature downwards assuming dry adiabatic c during the day and Stull cooling profile at night if (sina .le. 0.) then c nighttime call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index+1),vptp(i,j,index+1), : pt30(i,j),cellzc(k),vptdat(i,j,k)) AFTER CHANGE ------------ c frr 030119 Extrapolate temperature downwards assuming dry adiabatic c during the day and Stull cooling profile at night if (sina .le. 0.) then c nighttime c --- MCB-C(060804) if (zp(i,j,index).ne.zp(i,j,index+1))then index1=index+1 else index1=index+2 endif call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index1),vptp(i,j,index1), : pt30(i,j),cellzc(k),vptdat(i,j,k)) BEFORE change: -------------- c frr 030119 Extrapolate temperature down to the surface if (sina .le. 0.) then c nighttime call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index+1),vptp(i,j,index+1), : pt30(i,j),0.,ptsurf) AFTER CHANGE ------------ c frr 030119 Extrapolate temperature down to the surface if (sina .le. 0.) then c nighttime c --- MCB-C(060804) if (zp(i,j,index).ne.zp(i,j,index+1))then index1=index+1 else index1=index+2 endif call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index1),vptp(i,j,index1), : pt30(i,j),0.,ptsurf) --------------------------------------- (l) Update subroutine RDMM5 --------------------------------------- c --- UPDATE: c --- Check to ensure that the 2 temperatures used to extrapolate c --- temperatures down to the surface are at different altitudes BEFORE change: -------------- c initialize pot.temp history for temp extrapolation if (ifirstpg.le.1 )then call stull0(pt(index),z(index),pt(index+1),z(index+1), : pt20(i,j),pt30(i,j)) endif AFTER CHANGE ------------ c initialize pot.temp history for temp extrapolation if (ifirstpg.le.1 )then c --- MCB-C (060804) if (z(index).ne.z(index+1))then index1=index+1 else index1=index+2 endif call stull0(pt(index),z(index),pt(index1),z(index1), : pt20(i,j),pt30(i,j)) endif BEFORE change: -------------- if (sina .le. 0.) then c nighttime call stull(zp(i,j,index) ,pt(index),pt20(i,j), : zp(i,j,index+1),pt(index+1), : pt30(i,j),cellzc(k),ptmm4) AFTER CHANGE ------------ if (sina .le. 0.) then c --- nighttime c --- MCB-C (060804) if (zp(i,j,index).ne.zp(i,j,index+1))then index1=index+1 else index1=index+2 endif call stull(zp(i,j,index) ,pt(index), : pt20(i,j),zp(i,j,index1), : pt(index1),pt30(i,j),cellzc(k) : ,ptmm4) BEFORE change: -------------- c frr 030106 Extrapolate temperature to the surface (Z=0) if (sina .le. 0.) then c nighttime call stull(zp(i,j,index) ,pt(index),pt20(i,j), : zp(i,j,index+1),pt(index+1), : pt30(i,j),0.,ptsurf) AFTER CHANGE ------------ c frr 030106 Extrapolate temperature to the surface (Z=0) if (sina .le. 0.) then c nighttime c --- MCB-C (060804) if (zp(i,j,index).ne.zp(i,j,index+1))then index1=index+1 else index1=index+2 endif call stull(zp(i,j,index) ,pt(index),pt20(i,j), : zp(i,j,index1),pt(index1), : pt30(i,j),0.,ptsurf) xxxxxxxxxxxxxxxx Problem Area 7 -- xxxxxxxxxxxxxxxx --------------------------------------- (m) Update subroutine RDMM5 --------------------------------------- c --- UPDATE: c --- Save the value of internal time counters to ensure their values are c --- memorized from one step to the next independently of the compiler. BEFORE CHANGE ------------ c --- Initialise variables at first access to MM5 data records AFTER CHANGE ------------ c --- MCB-C (060804) save it1,it2 c --- Initialise variables at first access to MM5 data records xxxxxxxxxxxxxxxx Problem Area 8 -- xxxxxxxxxxxxxxxx --------------------------------------- (n) Update subroutine RDMM4 --------------------------------------- c --- UPDATE: c --- Save the previous hour prognostic temperature soundings in order to c --- compute the previous hour lapse rates above the mixing height if ITPROG>0 c BEFORE change: -------------- common/mm5temp/nlevag1,zl(mxnx,mxny,mxnzp+1), : tz(mxnx,mxny,mxnzp+1) AFTER CHANGE ------------ c MCB-C (060804) common/mm5temp/nlevag0,zl0(mxnx,mxny,mxnzp+1), : tz0(mxnx,mxny,mxnzp+1) BEFORE change: -------------- & tsurf(mxnxp,mxnyp) AFTER CHANGE ------------ & tsurf(mxnxp,mxnyp), c --- MCB-C (060804) & zl(mxnx,mxny,mxnzp+1),tz(mxnx,mxny,mxnzp+1) BEFORE change: -------------- c USe above surface MM4 levels ( index -> top)- add one for surface nlevag1=nlevag+1 AFTER CHANGE ------------ c Use above surface MM4 levels ( index -> top)- add one for surface c --- MCB-C (060804) First save previous hour value nlevag0=nlevag1 nlevag1=nlevag+1 c --- MCB-C (060804) Save previous hour soundings to compute previous hour c --- lapse rate in MIXDT do 622 k = 1,nlevag0 tz0(i,j,k)=tz(i,j,k) zl0(i,j,k)=zl(i,j,k) 622 continue --------------------------------------- (o) Update subroutine RDMM5 --------------------------------------- c --- UPDATE: c --- Save the previous hour prognostic temperature soundings in order to c --- compute the previous hour lapse rates above the mixing height if ITPROG>0 c BEFORE change: -------------- common/mm5temp/nlevag1,zl1(mxnx,mxny,mxnzp+1), : tz1(mxnx,mxny,mxnzp+1) AFTER CHANGE ------------ c --- MCB-C (060804) common/mm5temp/nlevag0,zl0(mxnx,mxny,mxnzp+1), : tz0(mxnx,mxny,mxnzp+1) BEFORE change: -------------- & zl2(mxnx,mxny,mxnzp+1),tz2(mxnx,mxny,mxnzp+1), AFTER CHANGE ------------ & zl2(mxnx,mxny,mxnzp+1),tz2(mxnx,mxny,mxnzp+1), c --- MCB-C (060804) & zl1(mxnx,mxny,mxnzp+1),tz1(mxnx,mxny,mxnzp+1), BEFORE change: -------------- c --- Check if last step in calmet run - AFTER CHANGE ------------ c c --- MCB-C (060804) c --- last available temperature profile for use in MIXDT2 (Mix.H. growth) nlevag0= nlevag1 do 1221 i = 1,nx do 1221 j = 1,ny do 1222 k = 1,mxnzp+1 tz0(i,j,k) = tz1(i,j,k) zl0(i,j,k) = zl1(i,j,k) 1222 continue 1221 continue c c --- Check if last step in calmet run - BEFORE change: -------------- nlevag1=nlevag+1 AFTER CHANGE ------------ c --- MCB-C (060804) c --- keep track of previous step number of vertical layers nlevag0 = nlevag1 nlevag1=nlevag+1 BEFORE change: -------------- 228 continue goto 1 AFTER CHANGE ------------ 228 continue c --- MCB-C (060804) c --- Set up temperature profile for use in MIXDT2 (Mix.H. growth) nlevag0= nlevag1 do 1321 i = 1,nx do 1321 j = 1,ny do 1322 k = 1,mxnzp+1 tz0(i,j,k) = tz2(i,j,k) zl0(i,j,k) = zl2(i,j,k) 1322 continue 1321 continue goto 1 xxxxxxxxxxxxxxxx Problem Area 9 -- xxxxxxxxxxxxxxxx --------------------------------------- (p) Update common file MET1.MET --------------------------------------- c --- UPDATE: c --- Store the surface station vertical extrapolation scaling parameter c --- ZLOGSTA BEFORE change: -------------- 3 xstz(mxss),zanem(mxss),cunam(mxus),xusta(mxus), 4 yusta(mxus),xulat(mxus),xulon(mxus),xutz(mxus), AFTER CHANGE ------------ 3 xstz(mxss),zanem(mxss),zlogsta(mxss),cunam(mxus), 4 xusta(mxus),yusta(mxus),xulat(mxus),xulon(mxus),xutz(mxus), --------------------------------------- (q) Update subroutine DIAGI --------------------------------------- c --- UPDATE: c --- Compute adjustment coefficients at each surface station to ensure c --- that observed surface winds are scaled down from anemometer height c --- to the first CALMET level using a neutral log profile with roughness c --- length Zo at the gridpoint where the station is located - Store the c --- coefficients in common file MET1.MET (with associated changes in c --- calling list) BEFORE change: -------------- subroutine diagi(nssta,nusta,elev, 1 xssta,yssta,csnam,xusta,yusta,cunam,nowsta, 2 noobs) AFTER CHANGE: ------------- subroutine diagi(elev,z0,nowsta) BEFORE change: -------------- real elev(mxnx,mxny) real xssta(mxss),yssta(mxss),xusta(mxss),yusta(mxss) character*4 csnam(mxss),cunam(mxss),namst AFTER CHANGE: ------------- c --- MCB-C (060804) real elev(mxnx,mxny),z0(mxnx,mxny) character*4 namst BEFORE change: -------------- include 'd6.met' AFTER CHANGE: ------------- include 'd6.met' c --- MCB-C (060804) include 'met1.met' BEFORE change: -------------- c --- set missing data flags AFTER CHANGE: ------------- c --- MCB-C (060804) c --- Compute scaling factors to extrapolate surface wind observations c --- to first CALMET level (logarithmic profile) do n=1,nssta c --- z0: roughness length of gridpoint where station is located c --- if station outside of domain use z0 at gridpoint closest to station isc=ist(n) jsc=jst(n) if(isc.gt.nx)isc=nx if(isc.lt.1)isc=1 if(jsc.gt.ny)jsc=ny if(jsc.lt.1)jsc=1 xlnzo=alog(z0(isc,jsc)) c --- z1: first CALMET level xlnz1=alog(cellzc(1)) c --- z2: anemometer height xlnz2=alog(zanem(n)) c --- Logarithmic profile scaling factor c u(k=1)=zlogsta*u(zanem) zlogsta(n)=(xlnz1-xlnzo)/(xlnz2-xlnzo) end do c --- set missing data flags --------------------------------------- (r) Update subroutine SETUP --------------------------------------- c --- UPDATE: c --- Modify calling list to subroutine DIAGI to provide the necessary c --- parameters for the vertical extrapolation of surface wind observations c --- and to reflect the fact that common file MET1.MET is now included in c --- DIAGI BEFORE change: -------------- c --- Wind field module set-up call diagi(nssta,nusta,elev, 1 xssta,yssta,csnam,xusta,yusta,cunam, 2 nowsta,noobs) AFTER CHANGE: ------------- c --- Wind field module set-up c --- MCB-C (060804) call diagi(elev,z0,nowsta) --------------------------------------- (s) Update common file OVRWAT.MET --------------------------------------- c --- UPDATE: c --- Store the overwater station vertical extrapolation scaling parameter c --- (ZLOGWSTA) and roughness length (Z0ow) BEFORE change: -------------- 5 xowlon(mxows),wsow(mxows),wdow(mxows), 6 lremap,iutmznow,feastow,fnorthow, AFTER CHANGE ------------ 5 xowlon(mxows),wsow(mxows),wdow(mxows),zlogwsta(mxows), 6 z0ow(mxows),lremap,iutmznow,feastow,fnorthow, --------------------------------------- (t) Update subroutine RDOW --------------------------------------- c --- UPDATE: c --- Compute logarithmic adjustment coefficient ZLOGWSTRA to extrapolate c --- wind measurements from anemometer height to first CALMET level c --- (Zlogwsta stored in OVRWAT.MET) c BEFORE change: -------------- include 'map.met' AFTER CHANGE: ------------- include 'map.met' include 'grid.met' BEFORE change: -------------- character*4 c4hem AFTER CHANGE ------------ character*4 c4hem data z0min/2.e-6/ BEFORE change: -------------- 100 continue c return AFTER CHANGE ------------ 100 continue c c --- MCB-C (060804) c --- Compute scaling factors to extrapolate surface wind observations c --- to first CALMET level (neutral logarithmic profile) do n=1,nowsta c --- z0: overwater roughness length c --- Use Hosker formula (1974) here, knowing that not quite accurate c --- as it requires wind at 10m , not at Zowsta. if(wsow(n).gt.0.0 .and. wsow(n).lt.9998.)then z0ow(n)= 2.0e-6 * wsow(n) ** 2.5 z0ow(n) = amax1(z0min,z0ow(n)) else z0ow(n) = z0min endif xlnzo=alog(z0ow(n)) c --- z1: first CALMET level xlnz1=alog(zmid(1)) c --- z2: anemometer height (if missing, assume 10m) if (zowsta(n).gt.9998) then xlnz2=alog(10.) else xlnz2=alog(zowsta(n)) endif c --- Logarithmic profile scaling factor c u(k=1)=zlogssta*u(zowsta) zlogwsta(n)=(xlnz1-xlnzo)/(xlnz2-xlnzo) end do c --- return --------------------------------------- (u) Update subroutine COMP --------------------------------------- c --- UPDATE: c --- Pass on the overwater station scaling factors to subroutine c --- DIAGNO BEFORE change: -------------- call diagno(nhrind,gamma,beta2,um,vm,nowsta,zowsta, & ziconv,tprog,vptprog,ccgrid,iceilg,rho) AFTER CHANGE: ------------- call diagno(nhrind,gamma,beta2,um,vm,nowsta,zowsta, & zlogwsta,ziconv,tprog,vptprog,ccgrid,iceilg,rho) --------------------------------------- (v) Update subroutine DIAGNO --------------------------------------- c --- UPDATE: c --- Scale the surface observations (land and water) from anemometer height c --- to first CALMET level using user-defined extrapolation technique (iextrp>1) c --- or, if none has been selected (iextrp=1), using a neutral log profile BEFORE change: -------------- subroutine diagno(nh,gamma,beta2,um,vm,nowsta,zowsta, & ziconv,tprog,vptprog,ccgrid,iceilg,rho) AFTER CHANGE: ------------- subroutine diagno(nh,gamma,beta2,um,vm,nowsta,zowsta, & zlogwsta,ziconv,tprog,vptprog,ccgrid,iceilg,rho) BEFORE change: -------------- real zowsta(mxows) AFTER CHANGE: ------------- real zowsta(mxows),zlogwsta(mxows) BEFORE change: -------------- IF(IABS(IEXTRP).EQ.1) GO TO 91 AFTER CHANGE: ------------- c --- MCB-C (060804) c --- if iextrp=1, must extrapolate from anemometer to 1st CALMET level c IF(IABS(IEXTRP).EQ.1) GO TO 91 BEFORE change: -------------- 63 IF (IABS(IEXTRP) .EQ. 2) THEN AFTER CHANGE: ------------- c --- MCB-C (060804) c --- Anemometer height and neutral log adjustment coefficient if(L.le.NSOL)then c --- Station is an overland surface station zwsref=zanem(L) zlog=zlogsta(l) else c --- Station is an overwater station zlog=zlogwsta(nsow) zwsref=zowsta(nsow) c --- Use 10 m overwater anemometer height if not specified if(zwsref.ge.9998.)zwsref=10. endif IF (IABS(IEXTRP) .EQ. 1) THEN c --- Adjust from Anenometer height to first CALMET level only c --- using neutral log profile IF(US(1,L).GT.EDITL .OR. VS(1,L).GT.EDITL) GOTO 70 us(1,l)=us(1,l)*zlog vs(1,l)=vs(1,l)*zlog ENDIF 63 IF (IABS(IEXTRP) .EQ. 2) THEN BEFORE change: -------------- IF (K .EQ. 2) THEN if(L.le.NSOL)then c --- Station is an overland surface station zwsref=zanem(L) else c --- Station is an overwater station zwsref=zowsta(nsow) c --- Use 10 m overwater anemometer height if not specified if(zwsref.ge.9998.)zwsref=10. endif PADJ = (CELLZC(K)/zwsref)**pexp2 ELSE PADJ = ( CELLZC(K)/CELLZC(KM1) )**pexp2 END IF AFTER CHANGE: ------------- c --- MCB-C (060804) IF (K .EQ. 2) THEN PADJ = (CELLZC(K)/zwsref)**pexp2 ELSE PADJ = ( CELLZC(K)/CELLZC(KM1) )**pexp2 END IF BEFORE change: -------------- 67 CONTINUE AFTER CHANGE: ------------- 67 CONTINUE c --- MCB-C (060804) c --- Also adjust winds from anenometer height to first calmet level PADJ = (zwsref/CELLZC(1))**pexp2 US(1,L) = US(1,L)/PADJ VS(1,L) = VS(1,L)/PADJ BEFORE change: -------------- DO 69 K=2,NZ AFTER CHANGE: ------------- c --- MCB-C (060804) DO 69 K=NZ,1,-1 --------------------------------------- (w) Update subroutine STHEOR --------------------------------------- c --- UPDATE: c --- Pass on the anemometer height to WATER2 and ELUSTR2 subroutines c --- in order to extrapolate surface observations from anemometer height c --- to first CALMET level BEFORE change: -------------- c --- Treat SEA#.DAT station as water call water2(us,vs,ns,is,js,fcoriol,nsurf,zi) AFTER CHANGE: ------------- c --- Treat SEA#.DAT station as water c --- MCB-C call water2(us,vs,ns,9999.,is,js,fcoriol,nsurf,zi) BEFORE change: -------------- c --- If SURF.DAT station but cell is water, treat as water c --- MCB-C (060804) call water2(us,vs,ns,is,js,fcoriol,nsurf,zi) AFTER CHANGE: ------------- c --- If SURF.DAT station but cell is water, treat as water c --- MCB-C (060804) call water2(us,vs,ns,zanem(ns),is,js,fcoriol,nsurf,zi) BEFORE change: -------------- call elustr2(ziconv,us,vs,ns,is,js,tempk,icc,iupt, 1 iicloud,ccgrid,zi,itprog,tprog,rho) AFTER CHANGE: ------------- c --- MCB-C (060804) call elustr2(ziconv,us,vs,ns,zanem(ns),is,js,tempk,icc,iupt, 1 iicloud,ccgrid,zi,itprog,tprog,rho) --------------------------------------- (x) Update subroutine WATER2 --------------------------------------- c --- UPDATE: c --- Take into account that the observed surface winds are c --- at anemometer height (not necessarily the same as first CALMET level) BEFORE change: -------------- subroutine water2(us,vs,ns,is,js,fcoriol,nsurf,zi) AFTER CHANGE: ------------- subroutine water2(us,vs,ns,zzanem,is,js,fcoriol,nsurf,zi) BEFORE change: -------------- c --- INPUT: c US(MXNZ,MXWND) - real array - U component of observed winds c VS(MXNZ,MXWND) - real array - V component of observed winds c ZI(MXNX,MXNY) - real array - previous hour mixing heights c NS - integer - Location of station in surface c array AFTER CHANGE: ------------- c --- INPUT: c US(MXNZ,MXWND) - real array - U component of observed winds (for k=1: c observation is at anemometer height, not c at first CALMET level c VS(MXNZ,MXWND) - real array - V component of observed windsfor k=1: c observation is at anemometer height, not c at first CALMET level c ZI(MXNX,MXNY) - real array - previous hour mixing heights c NS - integer - Location of station in surface array c ZZANEM - real - Surface (land) station anemometer height c For overwater station, actual value is c passed on via common OVRWAT.MET BEFORE change: -------------- c --- Do not extrapolate calm winds ws10 = sqrt(us(1,ns) **2 + vs(1,ns) **2) if (ws10 .lt. 0.0001) return AFTER CHANGE: ------------- c --- Do not extrapolate calm winds c --- MCB-C (060804):us(1,n),vs(1,n) are winds measured at anemometer height c --- not necessarily at 10m wsz = sqrt(us(1,ns) **2 + vs(1,ns) **2) if (wsz .lt. 0.0001) return BEFORE change: -------------- c --- If station is not a SEA#.DAT station, skip following section AFTER CHANGE: ------------- c --- MCB-C (060804) c --- Anemometer height for actual offshore stations (in SEA.DAT) if (nwat.gt.0) then zzanem = zowsta(nwat) c --- Assume 10. m anemometer height if missing if(zzanem.ge.9998.)zzanem=10. end if c --- If station is not a SEA#.DAT station, skip following section BEFORE change: -------------- c --- compute overwater surface roughness (Hosker, 1974) if(ws10.gt.0.0)then zz0 = 2.0e-6 * ws10 ** 2.5 zz0 = amax1(z0min,zz0) else zz0 = z0min endif AFTER CHANGE: ------------- c --- compute overwater surface roughness (Hosker, 1974) c --- MCB-C (060804) if(wsz.gt.0.0)then zz0 = 2.0e-6 * wsz ** 2.5 zz0 = amax1(z0min,zz0) c --- Adjust wind speed to 10m and recompute z0 if (int(zzanem).ne.10) then ws10=wsz*alog(10./zz0)/alog(zzanem/zz0) zz0 = 2.0e-6 * ws10 ** 2.5 endif else zz0 = z0min endif BEFORE change: -------------- xlnzz0 = alog(zowsta(nwat) / zz0) AFTER CHANGE: ------------- c --- MCB-C (060804) xlnzz0=alog(zzanem/zz0) BEFORE change: -------------- ztl = zowsta(nwat) / el2 AFTER CHANGE: ------------- c --- MCB-C (060804) ztl=zzanem/el2 --------------------------------------- (y) Update subroutine ELUSTR2 --------------------------------------- c --- UPDATE: c --- Take into account that the observed surface winds are c --- at anemometer height (not necessarily the same as first CALMET level) BEFORE change: -------------- subroutine elustr2(ziconv,us,vs,ns,is,js, & tempk,icc,iupt,icloud,ccgrid,zi,itprog,tprog,rho) AFTER CHANGE: ------------- subroutine elustr2(ziconv,us,vs,ns,zzanem,is,js, & tempk,icc,iupt,icloud,ccgrid,zi,itprog,tprog,rho) BEFORE change: -------------- c NS - integer - Location of station in surface c array AFTER CHANGE: ------------- c NS - integer - Location of station in surface c array c ZZANEM - real - Anemometer height at surface station BEFORE change: -------------- xa = (1.0 - 15.0 * zmid(1) / xl) ** 0.25 AFTER CHANGE: ------------- c --- MCB-C (060804) xa = (1.0 - 15.0 * zzanem / xl) ** 0.25 BEFORE change: -------------- zzanem = 10. call similt(zzanem,el,z0(is,js),zii,ns,us,vs,zimin) AFTER CHANGE: ------------- c --- MCB-C (060804) Use actual anemometer height c zzanem = 10. call similt(zzanem,el,z0(is,js),zii,ns,us,vs,zimin) --------------------------------------- (z) Update subroutine SIMILT --------------------------------------- c --- UPDATE: c --- Use anenometer height rather than first CALMET level BEFORE change: -------------- c --- Calculate ln(z1/z0) and stability function psi-m(z1/L) at level 1 xlnz1z0 = alog(zmid(1) / z0) AFTER CHANGE: ------------- c --- Calculate ln(z1/z0) and stability function psi-m(z1/L) at c anenometer height (MCB-C 060804) xlnz1z0 = alog(zzanem / z0) BEFORE change: -------------- c --- Stable psim1 = -17. * (1. - exp(-0.29 * zmid(1) / eladj)) else c --- Unstable c ***** psim1 = (1.0 - 16.0 * zmid(1) / eladj) ** 0.25 - 1.0 x = (1. - 16. * zmid(1) / eladj) ** .25 AFTER CHANGE: ------------- c --- Stable psim1 = -17. * (1. - exp(-0.29 * zzanem / eladj)) else c --- Unstable c ***** psim1 = (1.0 - 16.0 * zmid(1) / eladj) ** 0.25 - 1.0 x = (1. - 16. * zzanem / eladj) ** .25 BEFORE change: -------------- c --- Loop over all levels above level 1 do k = 3,nzp1 AFTER CHANGE: ------------- c --- Loop over all levels above surface (including level 1 - MCB-C 060804) do k = 2,nzp1 --------------------------------------- (aa) Update subroutine WIND1 --------------------------------------- c --- UPDATE: c --- Scale the surface observations (land and water) from anemometer height c --- to first CALMET level using user-defined extrapolation technique (iextrp>1) c --- or, if none has been selected (iextrp=1), using a neutral log profile BEFORE change: -------------- real ustmp(mxnz,mxwnd),vstmp(mxnz,mxwnd),zzanem(mxwnd) AFTER CHANGE: ------------- real ustmp(mxnz,mxwnd),vstmp(mxnz,mxwnd),zzanem(mxwnd) real zlog(mxwnd) BEFORE change: -------------- c --- Define combined surface+overwater station anemometer height array if(iabs(iextrp).eq.2.or.iabs(iextrp).eq.4)then c --- Loop over surface stations first do LL=1,nssta zzanem(LL)=zanem(LL) enddo c --- Add overwater station anem. hts. to end of array indx=nssta do LL=1,nowsta indx=indx+1 zzanem(indx)=zowsta(LL) c --- If overwater anem. height missing, use 10 m if(zzanem(indx).ge.9998.)zzanem(indx)=10. enddo endif AFTER CHANGE: ------------- c --- Define combined surface+overwater station anemometer height array c --- and neutral log profile coefficient array for all extrapolation options c --- MCB-C (060804) c if(iabs(iextrp).eq.2.or.iabs(iextrp).eq.4)then c --- Loop over surface stations first do LL=1,nssta zzanem(LL)=zanem(LL) zlog(ll)=zlogsta(ll) enddo c --- Add overwater station anem. hts. to end of array indx=nssta do LL=1,nowsta indx=indx+1 zzanem(indx)=zowsta(LL) c --- If overwater anem. height missing, use 10 m if(zzanem(indx).ge.9998.)zzanem(indx)=10. zlog(indx)=zlogwsta(ll) enddo c endif BEFORE change: -------------- c --power law AFTER CHANGE: ------------- c --- MCB-C (060804) c --- if no extrapolation, still adjust from anemometer height to first c --- calmet level using neutral log profile if(iabs(iextrp).eq.1) then ustmp(1,l)=ustmp(1,l)*zlog(l) vstmp(1,l)=vstmp(1,l)*zlog(l) endif c --- power law BEFORE change: -------------- vstmp(k,l)=vstmp(km1,l)*padj end do endif AFTER CHANGE: ------------- vstmp(k,l)=vstmp(km1,l)*padj end do c --- MCB-C (060804) c --- Adjust from anemometer height to first level padj=(zzanem(L)/cellzc(1))**pexp2 ustmp(1,l)=ustmp(1,l)/padj vstmp(1,l)=vstmp(1,l)/padj endif BEFORE change: -------------- if(iabs(iextrp).eq.3) then do k=2,nz AFTER CHANGE: ------------- if(iabs(iextrp).eq.3) then c --- fextrp: nz values for scaling observations c --- => use the first one for first calmet level (MCB-C 060804) do k=nz,1,-1 BEFORE change: -------------- zml=10000. c --roughness length at gridpoint closest to station c --closest gridpoint isc=ist(l) jsc=jst(l) if(isc.gt.nx)isc=nx if(isc.lt.1)isc=1 if(jsc.gt.ny)jsc=ny if(jsc.lt.1)jsc=1 c --roughness length zr=z0(isc,jsc) AFTER CHANGE: ------------- zml=10000. c --- roughness length if (l.le.nssta) then c --- if land station use z0 at gridpoint closest to station isc=ist(l) jsc=jst(l) if(isc.gt.nx)isc=nx if(isc.lt.1)isc=1 if(jsc.gt.ny)jsc=ny if(jsc.lt.1)jsc=1 zr=z0(isc,jsc) else c --- overwater station: use current z0 (function of wind speed) c --- rather than the fixed constant from GEO.DAT (MCB-C) zr=z0ow(l-nssta) endif xxxxxxxxxxxxxxxx Problem Area 10 -- xxxxxxxxxxxxxxxx --------------------------------------- (ab) Update subroutine INTERPQR --------------------------------------- c --- UPDATE: c --- Use units of meters for the relative locations of the MM5 gridpoints BEFORE CHANGE: -------------- c 040716 (MCB-A) x(np)=xlcmm4(ip,jp)-xmap0 y(np)=ylcmm4(ip,jp)-ymap0 AFTER CHANGE: ------------- c 060804 (MCB-C) x(np)=(xlcmm4(ip,jp)-xmap0)*1000. y(np)=(ylcmm4(ip,jp)-ymap0)*1000. xxxxxxxxxxxxxxxx Problem Area 11 -- xxxxxxxxxxxxxxxx --------------------------------------- (ac) Update subroutine MIXHT2ST --------------------------------------- c --- UPDATE: c --- Include grid.met to define variables needed for output in debug NOOBS mode BEFORE CHANGE: -------------- include 'params.met' AFTER CHANGE: ------------- include 'params.met' include 'grid.met' xxxxxxxxxxxxxxxx Problem Area 12 -- xxxxxxxxxxxxxxxx --------------------------------------- (ad) Update input file CALMET.INP --------------------------------------- c --- UPDATE: c --- Insert new input parameter ILUOC3D describing the OCEAN land use c --- category in the MM4/MM5/3D dataset (Input Group 6) BEFORE change: -------------- OTHER MIXING HEIGHT VARIABLES AFTER CHANGE: ------------- Land Use category ocean in 3D.DAT datasets (ILUOC3D) Default: 16 ! ILUOC3D = 16 ! Note: if 3D.DAT from MM5 version 3.0, iluoc3d = 16 if MM4.DAT, typically iluoc3d = 7 OTHER MIXING HEIGHT VARIABLES --------------------------------------- (ae) Update subroutine READCF --------------------------------------- c --- UPDATE: c --- REad in the new input parameter ILUCO3D and store it in MM4HD0 BEFORE change: -------------- 6 'FCORIOL','CONSTW','ITPROG','IRAD','IAVET','TGDEFB','TGDEFA', 6 'JWAT1','JWAT2','TRADKM','NUMTS','NFLAGP','SIGMAP','CUTP', 6 'HA1','HA2','HB1','HB2','HC1','HC2','HC3',26*' ', AFTER CHANGE: ------------- 6 'FCORIOL','CONSTW','ITPROG','ILUOC3D', 6 'IRAD','IAVET','TGDEFB','TGDEFA', 6 'JWAT1','JWAT2','TRADKM','NUMTS','NFLAGP','SIGMAP','CUTP', 6 'HA1','HA2','HB1','HB2','HC1','HC2','HC3',25*' ', BEFORE change: -------------- 6 20*1, 2*mxwb, 12*1, 26*0, AFTER CHANGE: ------------- 6 21*1, 2*mxwb, 12*1, 25*0, BEFORE change: -------------- 6 9*1, 2*2, 1, 2, 2*1, 3*2, 2*1, 2*2, 1, 2*2, 9*1, 26*0, AFTER CHANGE: ------------- 6 9*1, 2*2, 1, 2, 2*1, 4*2, 2*1, 2*2, 1, 2*2, 9*1, 25*0, BEFORE change: -------------- 2 IAVEZI,MNMDAV,HAFANG,ILEVZI,FCORIOL,CONSTW,ITPROG, 2 IRAD,IAVET,TGDEFB, 3 TGDEFA,JWAT1,JWAT2,TRADKM,NUMTS,NFLAGP,SIGMAP,CUTP,HA1,HA2, 4 HB1,HB2,HC1,HC2,HC3,idum,idum,idum,idum,idum,idum,idum, 5 idum,idum,idum,idum,idum,idum,idum,idum,idum,idum,idum,idum, 6 idum,idum,idum,idum,idum,idum,idum) AFTER CHANGE: ------------- 2 IAVEZI,MNMDAV,HAFANG,ILEVZI,FCORIOL,CONSTW,ITPROG, 2 ILUOC3D,IRAD,IAVET,TGDEFB, 3 TGDEFA,JWAT1,JWAT2,TRADKM,NUMTS,NFLAGP,SIGMAP,CUTP,HA1,HA2, 4 HB1,HB2,HC1,HC2,HC3,idum,idum,idum,idum,idum,idum,idum, 5 idum,idum,idum,idum,idum,idum,idum,idum,idum,idum,idum,idum, 6 idum,idum,idum,idum,idum,idum) --------------------------------------- (af) Update common file MM4HDO.MET --------------------------------------- c --- UPDATE: c --- Store new input parameter iluoc3d BEFORE change: -------------- 6 IOUTMM5,IMM53D,ISTEPPG,DATUM3D, CVER3D AFTER CHANGE: ------------- 6 IOUTMM5,IMM53D,ISTEPPG,DATUM3D,ILUOC3D,CVER3D --------------------------------------- (ag) Update subroutine RDHD52 --------------------------------------- c --- UPDATE: c --- If MM5 records do not include land use, assign a water category c --- to gridpoints at zero elevations BEFORE change: -------------- 99 format(2i4,f9.4,f10.4,i5) AFTER CHANGE: ------------- 99 format(2i4,f9.4,f10.4,i5) c --- MCB-C (060804) c --- assign ocean LU to zero elevation points if (ielev4(i,j).eq.0) then ilu4(i,j)=iluoc3d else c --- make sure that not all points are flagged as ocean c --- if by chance user input iluoc3d=0 ilu4(i,j)=iluoc3d+1 endif --------------------------------------- (ah) Update subroutine RDMM4 --------------------------------------- c --- UPDATE: c --- For the vertical extrapolation of temperature from the lowest c --- MM4 level to the surface, Only use land cooling overland and c --- use linearly extrapolated profile overwater BEFORE change: -------------- call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index+1),vptp(i,j,index+1), : pt30(i,j),cellzc(k),vptdat(i,j,k)) AFTER CHANGE: ------------- c --- MCB-C (060804) if (ilu4(i,j).eq.iluoc3d) then c --- ocean zx=(cellzc(k)-zp(i,j,index))/ : (zp(i,j,index1)-zp(i,j,index)) vptdat(i,j,k)=vptp(i,j,index) : +zx*(vptp(i,j,index1)-vptp(i,j,index)) else c --- land call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index1),vptp(i,j,index1), : pt30(i,j),cellzc(k),vptdat(i,j,k)) endif BEFORE change: -------------- call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index+1),vptp(i,j,index+1), : pt30(i,j),0.,ptsurf) AFTER CHANGE: ------------- c --- MCB-C (060804) if (ilu4(i,j).eq.iluoc3d) then c --- ocean zx=-zp(i,j,index)/ : (zp(i,j,index1)-zp(i,j,index)) ptsurf=vptp(i,j,index) : +zx*(vptp(i,j,index1)-vptp(i,j,index)) else c --- land call stull(zp(i,j,index) ,vptp(i,j,index),pt20(i,j), : zp(i,j,index1),vptp(i,j,index1), : pt30(i,j),0.,ptsurf) endif --------------------------------------- (ai) Update subroutine RDMM5 --------------------------------------- c --- UPDATE: c --- For the vertical extrapolation of temperature from the lowest c --- MM5 level to the surface, Only use land cooling overland and c --- use linearly extrapolated profile overwater BEFORE change: -------------- call stull(zp(i,j,index) ,pt(index),pt20(i,j), : zp(i,j,index+1),pt(index+1), : pt30(i,j),cellzc(k),ptmm4) AFTER CHANGE: ------------- c --- MCB-C (060804) c --- Differentiate between land and water if (ilu4(i,j).eq.iluoc3d) then c --- ocean - zx=(cellzc(k)-zp(i,j,index))/ : (zp(i,j,index1)-zp(i,j,index)) ptmm4=pt(index) : +zx*(pt(index1)-pt(index)) else c --- land - surface cooling parameterized by Stull call stull(zp(i,j,index) ,pt(index), : pt20(i,j),zp(i,j,index1), : pt(index1),pt30(i,j),cellzc(k) : ,ptmm4) endif BEFORE change: -------------- call stull(zp(i,j,index) ,pt(index),pt20(i,j), : zp(i,j,index+1),pt(index+1), : pt30(i,j),0.,ptsurf) AFTER CHANGE: ------------- c --- MCB-C (060804) c --- Differentiate between land and water (060315-FRR) if (ilu4(i,j).eq.iluoc3d) then c --- ocean - c --- SSTP not available - extrapolate 1st-2nd levels zx=-zp(i,j,index)/ : (zp(i,j,index1)-zp(i,j,index)) ptsurf=pt(index) : +zx*(pt(index1)-pt(index)) else c --- land - surface cooling parameterized by Stull call stull(zp(i,j,index) ,pt(index),pt20(i,j), : zp(i,j,index1),pt(index1), : pt30(i,j),0.,ptsurf) endif ------------- --------------------------------------- (aj) Update subroutine RDMM5 --------------------------------------- c --- UPDATE: c --- Update the 2D surface temperature array (TEMP2D) to reflect the computations c --- performed on the 3D temperature array ZTEMP BEFORE CHANGE: ------------- call outhr(ndathr,nx,ny,nz,npsta,irtype,rho,qsw, 1 iomet,lcalgrd,temp2d,irh2d,ipcode2d) AFTER CHANGE -------------- c --- MCB-C (060804) Ensure Temp2D is the same as ZTEMP(k=1) c --- when ZTEMP is computed (lcalgrid=T or itprog>=1) if (LCALGRD .or. itprog.eq.1) then do j=1,ny do i=1,nx temp2d(i,j)=ztemp(i,j,1) end do end do endif call outhr(ndathr,nx,ny,nz,npsta,irtype,rho,qsw, 1 iomet,lcalgrd,temp2d,irh2d,ipcode2d) --------------------------------------- (ak) Update subroutine RDHD53 --------------------------------------- c --- UPDATE: c --- Read in prognostic land use BEFORE CHANGE: ------------- read(io20,99)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j) 99 format(2i4,f9.4,f10.4,i5) AFTER CHANGE -------------- read(io20,99)iindex,jindex,xlat4(i,j),xlong4(i,j), & ielev4(i,j),ilu4(i,j) 99 format(2i4,f9.4,f10.4,i5,i3) --------------------------------------- (al) Update BLOCK DATA --------------------------------------- BEFORE change: -------------- data datum3d/'NWS-84 '/ AFTER CHANGE: ------------- data datum3d/'NWS-84 '/ c --- MCB-C:060804: Default iluoc3d data iluoc3d/16/ xxxxxxxxxxxxxxxx Problem Area 13 -- xxxxxxxxxxxxxxxx --------------------------------------- (am) Update subroutine COMP --------------------------------------- c --- UPDATE: c --- Assign valid layer-1 temperatures to the 2D-array of surface temperature c --- array before output BEFORE CHANGE: ------------- c --- output results to disk 799 continue if(lsave)then if(iformo.eq.1)then c c FRR (01/09) !!! MUST CHANGE OUTPUT< PRTMET< CALPUFF!!! c rho, qsw,temp,rh,ipcode,and are computed at all gridpoints call surfvar(tprog,temp2d,irh2d,ipcode2d) AFTER CHANGE -------------- c --- output results to disk 799 continue if(lsave)then if(iformo.eq.1)then c c FRR (01/09) !!! MUST CHANGE OUTPUT< PRTMET< CALPUFF!!! c rho, qsw,temp,rh,ipcode,and are computed at all gridpoints call surfvar(tprog,temp2d,irh2d,ipcode2d) c c MCB-C (060804) Ensure Temp2D is the same as ZTEMP(k=1) c when LCALGRID=T or itprog.ge.1 if( LCALGRD .or. itprog.ge.1 )then do j=1,ny do i=1,nx temp2d(i,j)=ztemp(i,j,1) end do end do endif xxxxxxxxxxxxxxxx Problem Area 14 -- xxxxxxxxxxxxxxxx --------------------------------------- (an) Update subroutine COMP --------------------------------------- c --- UPDATE: c --- Assign ihrgmt before if-block BEFORE CHANGE: ------------- if (itprog .eq. 0) then ihrgmt=nhr+ibtz call mixht(el,ustar,qh,nx,ny,rho,ihrgmt,nears,nearu,iupt, 1 ilandu,iwat1,iwat2,ldbhr,zi,ziconv) AFTER CHANGE -------------- c --- MCB-C:060804: Assign ihrgmt before if-block ihrgmt=nhr+ibtz if (itprog .eq. 0) then c ihrgmt=nhr+ibtz call mixht(el,ustar,qh,nx,ny,rho,ihrgmt,nears,nearu,iupt, 1 ilandu,iwat1,iwat2,ldbhr,zi,ziconv)