subroutine tmpadj(hour,halfdy) c c +++PURPOSE+++ c This function adjusts the surface temperature on an hourly c basis by using the daily minumum and maximum temperatures. c It also uses a similar method for computing the hourly c dewpoint based on the daily average. c c Authors(s): John Witte, UofMn WCES @ USDA-ARS-NCSRL c Date: 04/01/93 c Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051 c August 1994 cd The program is modiffied by S. Dun Aug 14, 2008 to direcly use SI units c c +++ARGUMENT DECLARATIONS+++ real halfdy integer hour c c +++ARGUMENT DEFINITIONS+++ c halfdy - The amount of time from sunrise until noon (hrs). c hour - The hour of the day that we are calculating. c c +++PARAMETERS+++ include 'pmxpln.inc' include 'pmxnsl.inc' include 'pmxhil.inc' include 'pmxtls.inc' include 'pmxtil.inc' include 'pmxcrp.inc' include 'pmxelm.inc' c c +++COMMON BLOCKS+++ include 'cclim.inc' c read: tave,vwind,hradmj,hrtemp,radmj c include 'cwint.inc' c read: snodpy(mxplane),tfrdp(mxplan),tthawd(mxplan),surtmp c frdp(mxplan),thdp(mxplan),densg(mxplan) c include 'cwater.inc' c read: salb(mxplan), c include 'ccover.inc' c read: canhgt(mxplan) c include 'ccrpout.inc' c read: rrc(mxplan) c include 'cstruc.inc' c read: iplane cd Added by S. Dun, Nov. 21, 2006 include 'cangie.inc' c read : solef (solar energy adjust fator for sloped surface) cd End adding c c +++LOCAL VARIABLES+++ c save real alb,rads,lytomj,clouds,ktemp,aemiss,displ,wsrgh,etrgh, 1 convht,vkcons,denair,hcpair,htwind,radcof,semiss,sbcons, 2 netrad,grdp,gutdp,sysdep,ksnow,kres,kftill,kfutil, 3 numer,denom,surfcn,effk,bakup, gtdp, hrtave,factor c 2 netrad,radadj,grdp,gutdp,sysdep,ksnow,kres,kftill,kfutil, c c +++LOCAL DEFINITIONS+++ c alb - Albedo of the top layer for the system (Dec %). c rads - Radiation on a sloping surface (W/m^2). c lytomj - Conversn factor of solar radiation (Langly's to MJ/m^2). c clouds - Estimated cloud cover (Dec %). c ktemp - Average temperature in degrees Kelvin units. c aemiss - Atmospheric emissivity. c displ - Zero-plane displacment, a function of canopy cover (m). c wsrgh - Roughness parameter for wind speed (m). c etrgh - Roughness parameter for energy transfer (m). c convht - Convective heat transfer coeff (J/m^3-K). c vkcons - Von Karman constant. c denair - Approximate density of air (Kg/m^3). c hcpair - Heat capacity of air (J/Kg-K). c htwind - Height at which wind velocity is measured (m). c radcof - Long_wave radiation transfer coefficient (W/m^2-K). c semiss - Surface emissivity for all layers. c sbcons - Stefan-Bolzman constant (W/m^2-K^4). c netrad - Calculated net radiation (W/m^2). c radadj - Daily radiation adjusted to hourly (MJ/m^2). c grdp - Depth from soil's surface to first 0 degree isotherm (m). c gutdp - Temperature gradient depth of untilled layer (m). c sysdep - Depth of the frozen layer system from surface to c first 0 degree isotherm (m). c ksnow - These are temporary variables representing the thermal c kres --_ conductivity of the snow, residue, frozen tilled and c kftill -- frozen untilled layers. Reason for temps is because c kfutil - if no such layer exists, the values are set to zero. c numer - Temporary numerator. c denom - Temporary denominator. c surfcn - Surface conductivity of top layers (W/m-C). c effk - Effective thermal conductivity of the entire layered c system to 0 degree isotherm (W/m-C). c bakup - Temporary variable used in simplifying equations. c hrtave - hourly air temperature to aviod messing up with daily average temperature stored c c +++DATA INITIALIZATIONS+++ c data sbcons/8.1247e-11/ c data vkcons/0.4/ c data hcpair/101.0/ c data denair/1.0/ c data htwind/2.0/ c data lytomj/0.04184/ c data kres/0.0232/ c c +++END SPECIFICATIONS+++ cd sbcons = 8.1247e-11 sbcons = 5.6697e-8 vkcons = 0.4 hcpair = 1012.0 denair = 1.2 htwind = 2.0 lytomj = 0.04184 c c XXX WHY is the value for KRES here set to 0.0232, while in c subroutine FROST, kres = 0.168 ?????? dcf 2/23/94 c cd kres = 0.0232 kres = 0.05 c kres = 0.168 kres = kres*kresf grdp = 0.0 gtdp = 0.0 gutdp = 0.0 alb = salb(iplane) if (snodpy(iplane) .gt. 0.01) alb = 0.5 c -- Convert from MJ/m^2/hr to L/min. cd rads = (hradmj / 60) * (1 / lytomj) c -- Convert from MJ/m^2/hr to W/m^2 rads = (hradmj / 3600) *1E+6 c Modified by S. Dun, Jan 12, 2007 clouds = CloudC c cd if (rpoth .lt. 0.0001) then cd clouds = 1.0 cd else cd clouds = (1.0 - radmj /rpoth) / 0.7 cd endif cd if (clouds .gt. 0.999) then cd clouds = 1.0 cd endif cd End modifying call hrtmp(hour,halfdy) hrtave = hrtemp ktemp = hrtave + 273.16 semiss = 1.0 aemiss = (1.0 - 0.84 * clouds) * (1.0 - 0.261 * 1 exp(7.77e-4 * (hrtave * hrtave))) + (0.84 * clouds) displ = 0.77 * canhgt(iplane) c reza added next line 7/27/94 c if(displ.gt.htwind)displ=0.77*htwind c reza added the above line 7/27/94 c fix from John Laflen to prevent divide by zero below if(displ.ge.htwind)displ=0.77*htwind c fix from John Laflen above - 10/2/2001 if((snodpy(iplane) .lt. 0.01).and.(canhgt(iplane) .gt. 0.0)) then wsrgh = 0.13 * canhgt(iplane) elseif (snodpy(iplane) .lt. 0.01) then wsrgh = rrc(iplane) elseif (snodpy(iplane) .gt. canhgt(iplane)) then wsrgh = 0.0002 else wsrgh = 0.13 * (canhgt(iplane) - snodpy(iplane)) endif if (wsrgh .lt. 0.001) wsrgh = 0.001 c Reza added 7/27/94 if(wsrgh.gt..26)wsrgh=0.26 c end Reza 7/27/94 etrgh = 0.2 * wsrgh convht =((vkcons * vkcons) * denair * hcpair) / 1 (log((htwind - displ + etrgh) / etrgh) * 2 log((htwind - displ + wsrgh) / wsrgh)) radcof = 4.0 * semiss * sbcons * (ktemp**3) c c XXX WHY is RADADJ set equal to RADCOF here, then never used??????? c Is the correct variable being used below (RADCOF) or should c RADADJ be used??? Can RADADJ be deleted???? dcf 2/22/94 c c radadj = radcof cd netrad = (1.0 - alb)*rads + (aemiss - semiss)*sbcons*hrtave**4 netrad = (1.0 - alb)*rads + (aemiss - semiss)*sbcons*ktemp**4 c -- Now we calculate the gradient depth for both top frost and none... c -- If average temperature is below zero. if (hrtave .lt. 0.0) then c -- If top frost is present... if (tfrdp(iplane) .gt. 0.001) then if (thdp(iplane) .gt. 0.001) then grdp = 0.0 else grdp = tfrdp(iplane) endif else c -- If no top frost is present... if (frdp(iplane) .gt. 0.001) then if (thdp(iplane) .gt. 0.001) then grdp = 0.0 else grdp = frdp(iplane) endif endif endif c -- The above assumption comes directly out of Flershinger's notes c -- and assumes frost exists all the way from the frost depth to c -- the surface. c -- Now, we consider the gradient depth is avg temperature is above 0. else c -- If top frost is present... if (tfrdp(iplane) .gt. 0.001) then grdp = thdp(iplane) else c -- If no top frost is present... if (frdp(iplane) .gt. 0.001) then if (thdp(iplane) .gt. 0.001) then grdp = thdp(iplane) else grdp = 0.0 endif endif endif endif if ((grdp .le. tilld(iplane)) .or. (tilld(iplane) .gt. 1.0)) then gtdp = grdp gutdp = 0.0 else gtdp = tilld(iplane) gutdp = grdp - tilld(iplane) endif if (gtdp .lt. 0.001) gtdp = 0.0 if (gutdp .lt. 0.001) gutdp = 0.0 if (grdp .lt. 0.001) grdp = 0.0 sysdep = snodpy(iplane) + resdep(iplane) + grdp c -- Now we make sure that each layer exists, if not we simply c -- set the "tc" value of that layer to 1.0 so that multiplying c -- by the "tc" will leave the final value unaffected by that c -- absent layer. if (snodpy(iplane) .lt. 0.0001) then ksnow = 1.0 else c Considering thermal conductivity of the snow pack c model from Sturm et al., 1997 if (densg(iplane).lt. 156) then ksnow = 0.023 + 0.234 * (densg(iplane)/1000.) else ksnow = 0.138 - 1.01*(densg(iplane)/1000.) 1 + 3.233* (densg(iplane)/1000.)** 2 endif cd ksnow = 0.046 + 0.0451 * (densg(iplane)/1000) + 0.3447 * cd 1 ((densg(iplane)/1000) * (densg(iplane)/1000)) c ksnow = ksnow * ksnowf endif if (resdep(iplane) .lt. 0.0001) kres = 1.0 c c XXX IS the following right??? Values for KFTILL when gtdp.ge.0.0001 c were not defined. I got the value of 1.75 from subroutine FROST c REZA ??? if (gtdp .lt. 0.0001) then kftill = 1.0 else kftill = 1.75 endif c c XXX IS the following right??? Values for KFUTIL when gutdp.ge.0.0001 c were not defined. I got the value of 2.1 from subroutine FROST c REZA ??? if (gutdp .lt. 0.0001) then kfutil = 1.0 else kfutil = 2.1 endif cd Added by S. Dun, May 25, 2009 c When there is no frost layers in the soil, we assume 0c isotherm is at 0.1 mm if (sysdep .lt. 0.0001) then kftill = 1.75 grdp = 0.001 gtdp = 0.001 sysdep = snodpy(iplane) + resdep(iplane) + grdp endif cd End adding c -- Now we use the harmonic mean method to calculate the surface c -- temperature for the hour. NOTE that the final calculation of c -- numerator/denominator has been broken up into two parts... the c -- "numer" and the "denom". numer = (ksnow * kres * kftill * kfutil) * 1 (snodpy(iplane) + resdep(iplane) + grdp) denom = (ksnow * kres * kftill * gutdp) + 1 (ksnow * kres * gtdp * kfutil) + 2 (ksnow * resdep(iplane) * kftill * kfutil) + 3 (snodpy(iplane) * kres * kftill * kfutil) c -- Now we must make sure we have no floating point errors... if ((denom .lt. -0.0001) .or. (denom .gt. 0.0001)) then surfcn = numer / denom else surfcn = 0.0 endif c -- If the second case is chosen, an error has occurred! c -- The following line converts from MJoules to Langly's cd effk = surfcn / 0.0697333 c Directly use the SI units effk = surfcn if (sysdep .gt. 0.0) then cd surtmp(hour) = (netrad+(radcof+convht*(vwind*100))*hrtave+0.0)/ cd 1 (radcof+convht*(vwind*100)+(effk/sysdep)) surtmp(hour) = (netrad+(radcof+convht*(vwind))*hrtave)/ 1 (radcof+convht*(vwind)+(effk/sysdep)) else cd surtmp(hour) = (netrad+(radcof+convht*(vwind*100))*hrtave+0.0)/ cd 1 (radcof+convht*(vwind*100)) surtmp(hour) = (netrad+(radcof+convht*(vwind))*hrtave)/ 1 (radcof+convht*(vwind)) endif cd Added by S. Dun, May 25, 2009 c When there is snow, we assume surface temperture is not greater than 0C if ((surtmp(hour).gt. 0.0).and.(snodpy(iplane) .gt. 0.001)) then surtmp(hour) = 0.0 endif cd End adding return end