c file name: heat.for subroutine heat (layrsn, bszlyd, theta, bsfcla, & bsdblk, bwtdmn, bwtdmx, bhtsmx, bhtsmn, bszlyt) c + + + PURPOSE + + + c This program simulates soil temperature based on the algorithm c described by campbell g. s. (1985) chapter 4: soil temperature c and heat flux. p. 26-39. in soil physics with basic. c The program estimates daily minimum, maximum, and average soil c temperature at the center of each simulation layer. c The inputs needed to run the program are maximum daily air c temperature, and minimum daily air temperature. Furthermore, c soil bulk density, volumetric water content, and clay c fraction are used to calculate soil thermal properties. c DATE: 10/05/93 c MODIFIED: 12/10/93 c + + + KEY WORDS + + + c soil temperature c + + + COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'p1const.inc' include 'p1unconv.inc' include 'wpath.inc' C *** include 'm1sim.inc' include 'm1flag.inc' include 'm1subr.inc' include 'w1clig.inc' include 'h1temp.inc' c + + + LOCAL COMMON BLOCKS + + + C *** include 'hydro/hlayrs.inc' *$reference c + + + ARGUMENT DECLARATIONS + + + integer layrsn real bsdblk(*) real bsfcla(*) real bszlyd(*) real theta(*) real bwtdmx real bwtdmn real bhtsmn(*) real bhtsmx(*) real bszlyt(*) c + + + ARGUMENT DEFINITIONS + + + c bsdblk - Soil bulk density, Mg/m^3 c clay - Percent clay c bszlyd - Depth to the bottom of soil layer from the soil surface (mm) c layrsn - Number of soil layers used in the simulation c theta - Volumetric soil water content for layer, m^3/m^3 c bwtdmx - Daily maximum air temperature c bhtsmn - Daily minimum soil temperature, c c bhtsmx - Daily maximum soil temperature, c c bszlyt - Layer thickness, mm c + + + PARAMETERS + + + real cm, cmjtoj, cw, dt, period parameter (cm = 2.26, cmjtoj = 1000000.0, cw = 4.18, & dt = 3600.0, period = 1.0) c cm - volumetric specific heat for soil mineral fraction, mj/m^3/c c cmjtoj - conversion factor from mj to joules c cw - volumetric specific heat for soil water fraction, mj/m^3/c c dt - time step of simulation, seconds c period - length of simulation period in days c + + + LOCAL VARIABLES + + + integer day, mo, yr integer k, l, j real a, b, c, d real claym(mnsz) real dlbc real freq real pbd real pclaym real phi real ptheta real sbd real sclaym real stheta real tamp real tdif real thermk real time(hrday) real tlbc real tsavg(mnsz) real tsoil(mnsz,hrday) real tubc real vsheat real zdamp real zdampy real dmlayr(mnsz) c + + + LOCAL DEFINITIONS + + + c claym - Clay mass fraction c day... - Day, month and year of simulation c dlbc - Depth of temperature at the lower boundry condition. c The larger of either the annual damping depth or the c depth at the lowermost simulation layer. c dmlayr - Depth to middle of layer, m c freq - Angular frequency c pbd - Weighted average bulk density of the soil profile, mg/m^3 c pclaym - Weighted average clay mass fractions of the soil profile c phi - Phase constant c ptheta - Weighted average water content of the soil profile, c m^3/m^3 c sbd - Sum of bulk densities fo the simulation layers, mg/m^3 c sclaym - Sum of clay mass fractions of the simulation layers c stheta - Sum of water contents of the simulation layers, m^3/m^3 c tamp - Amplitude of the temperature wave at the soil surface, c c tdif - Difference in temperature between the upper and lower c boundary conditions, c c thermk - Soil thermal conductivity, w/m/c c time - Time of simulation, s c tlbc - Temperature at the lower boundary condition, c c tsavg - Average daily soil temperature, calculated by c interpolation between the upper and lower boundary c conditions, c c tsoil - Hourly soil temperature, c c tubc - Temperature at the upper boundary condition, c c vsheat - Volumetric specific heat, j/m^3/c c zdamp - Diurnal damping depth, m c zdampy - Annual damping depth, m c + + + FUNCTIONS CALLED + + + c integer lentrm c + + + SUBROUTINES CALLED + + + c stat c + + + DATA INITIALIZATIONS + + + phi = -(7. * pi) / 12. c initialize various soil layer references dmlayr(1) = 0.5 * bszlyd(1) / mtomm do 100 k=2, layrsn dmlayr(k) = 0.5 * (bszlyd(k-1) + bszlyd(k)) / mtomm 100 continue c + + + OUTPUT FORMATS + + + 2000 format (' ') 2009 format (' Soil temperature output ') 2010 format (/,3x,'simulated daily minimum, and ', & 'maximum temperatures (celsius) of soil layers') 2020 format (1x,79('-')/,17x,28('-'),' layers ',27('-'),/,1x, & 'sr date',11x,'1',12x,'2',12x,'3',12x,'4',12x,'5') 2025 format (1x,22x,'6',12x,'7',12x,'8',12x,'9',11x,'10') 2026 format (1x,79('-')) 2030 format (1x,14x,5(4x,'min max'),/,1x,14x,5(2x,11('-'))) 2040 format (1x,i2,2x,2(i2,'/'),i4,2x,5(f5.1,1x,f5.1,2x), & /,17x,5(f5.1,1x,f5.1,2x)) c + + + END SPECIFICATIONS + + + if ((am0ifl .eqv. .true.).and.(am0csr .eq. 1).and.((am0hfl .eq. 4) & .or.(am0hfl .eq. 5).or.(am0hfl .eq. 6).or.(am0hfl .eq. 7))) then open(13,file= rootp(1:lentrm(rootp)) // 'temp.out', & access='sequential', status='unknown') c open(13,file='temp.out',access='sequential',status='unknown', c & mode = 'write') write(13,2009) write(13,2010) write(13,2020) if (layrsn .gt. 5) then write(13,2025) end if write(13,2026) write(13,2030) end if c calculate weighted means for soil properties of the soil profile sbd=0.0 sclaym=0.0 stheta=0.0 do 210 l=1,layrsn claym(l) = bsfcla(l) sbd=sbd+(bsdblk(l)*bszlyt(l)/mtomm) sclaym=sclaym+(claym(l)*bszlyt(l)/mtomm) stheta=stheta+(theta(l)*bszlyt(l)/mtomm) 210 continue pbd=sbd/bszlyd(layrsn) pclaym=sclaym/bszlyd(layrsn) ptheta=stheta/bszlyd(layrsn) c calculate volumetric specific heat, j/m^3/ c vsheat= cmjtoj*(((cm*pbd)/2.65)+(cw*ptheta)) c calculate soil thermal conductvity c reference: campbell (1985) p. 32 a= 0.65 - 0.78*pbd + 0.6*pbd**2 b= 1.06 * pbd c= 1.0 + (2.6/sqrt(pclaym)) d= 0.03 + (0.1*pbd**2) thermk= a+(b*ptheta)-(a-d)*exp(-(c*ptheta)**4) c calculate angular frequency of the temperature oscillation freq= (2*pi)/(86400.0*period) c calculate damping depth c diurnal damping depth zdamp= sqrt((2*(thermk/vsheat))/freq) c annual damping depth zdampy= zdamp*sqrt(365.0) c begin simulating heat flow for each day of the simulation period c calculate mean temperature and amplitude temperature at the c soil-atmosphere interface ( upper boundary condition ) c set lower boundary condition to an average air temperature tubc = (bwtdmx+bwtdmn) / 2 tamp = bwtdmx - tubc tlbc = awtyav tdif = tlbc - tubc dlbc = amax1(dmlayr(layrsn),zdampy) do 350 k=1,layrsn tsavg(k) = tubc + ( tdif * (dmlayr(k)/dlbc) ) ahtsav(k,am0csr) = tsavg(k) do 400 j = 1,24 time(j) = j * dt tsoil(k,j)= tsavg(k)+tamp*exp(-dmlayr(k)/zdamp)* & sin((freq*time(j))-(dmlayr(k)/zdamp)+phi) 400 continue c use subroutine stat to calculate the daily maximum, c and minimum temperatures at the center of each simulation c layer call stat(tsoil,bhtsmx(k),bhtsmn(k),k, mnsz, hrday) 350 continue if ((am0hfl .eq. 4) .or. (am0hfl .eq. 5) .or. (am0hfl .eq. 6) & .or. (am0hfl .eq. 7)) then call caldat (-1,day,mo,yr) if ((am0csr .eq.1) .and. (nsubr .gt. 1)) write(13,2000) write(13,2040) am0csr,day,mo,yr,(bhtsmn(k),bhtsmx(k), * k=1,layrsn) end if c close (unit = 13) return end subroutine stat (tsoil, bhtsmx, bhtsmn, k, mnsz, hrday) c This subroutine uses the simulated hourly soil temperature c values to calculate the daily maximum, and minimum soil c temperatures at the center of each simulation layer c + + + COMMON BLOCKS + + + c + + + ARGUMENT DECLARATIONS + + + integer k, j, mnsz, hrday real tsoil(mnsz, hrday) real bhtsmn real bhtsmx c c + + + ARGUMENT DEFINITIONS + + + c k - Simulation layer counter c tsoil - Hourly soil temperature, c c bhtsmn - Daily minimum soil temperature, c c bhtsmx - Daily maximum soil temperature, c c + + + END SPECIFICATIONS + + + bhtsmx= tsoil(k,1) bhtsmn= tsoil(k,1) do 100 j=1,24 if (tsoil(k,j) .gt. bhtsmx) bhtsmx= tsoil(k,j) if (tsoil(k,j) .lt. bhtsmn) bhtsmn= tsoil(k,j) 100 continue return end