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

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