!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOILEBM ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE soilebm(nx,ny,nz,soiltyp,vegtyp,lai,veg, & 1,5
tsfc,tdeep,wetsfc,wetdp,wetcanp,snowdpth, &
qvsfc,windsp,psfc,rhoa,precip, &
tair,qvair,cdha,cdqa,radsw,rnflx,shflx,lhflx,gflx,ct, &
evaprg,evaprtr,evaprr,qvsat,qvsata,f34)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Predict the soil surface temperature and moisture contents by solving
! the surface energy and moisture budget equtions:
!
! 1. the ground surface temperature, Ts -- tsfc
! 2. the deep ground temperature, T2 -- tdeep
! 3. the surface soil moisture, wg -- wetsfc
! 4. the deep soil moisture, w2 -- wetdp
! 5. the canopy moisture, wr -- wetcanp
!
!-----------------------------------------------------------------------
!
! The equations are listed as follows.
!
!
! d Ts 2PI
! ------ = Ct (Rn - H - LE) - ----- (Ts - T2)
! d t Tau
!
!
! d T2 1
! ------ = ----- (Ts - T2)
! d t Tau
!
!
! d Wg C1 C2
! ------ = ------ (Pg - Eg) - ----- (Wg - Wgeq)
! d t ROw d1 Tau
!
!
! d W2 1
! ------ = ------ (Pg - Eg - Etr)
! d t ROw d2
!
!
! d Wr
! ------ = Veg P - Er
! d t
!
!
! where
!
! Tau -- 1 day in seconds = 86400 seconds
! PI -- number of PI = 3.141592654
! ROw -- Density of liquid water
! d1 -- Top layer depth of soil column, 0.01 m
! d2 -- Deep layer depth of soil column, 1 m
! Veg -- Vegetation fraction
! Ct -- Thermal capacity
! Rn -- Radiation flux, rnflx
! H -- Sensible heat flux, shflx
! LE -- Latent heat flux, lhflx = latent*(Eg + Ev)
! Eg -- Evaporation from ground
! Ev -- Evapotranspiration from vegetation, Ev = Etr + Er
! Etr -- Transpiration of the remaining part of the leaves
! Er -- Evaporation directly from the foliage covered by
! intercepted water
! P -- Precipitation rates
! Pg -- Precipitation reaching the ground,
! Pg = (1 - Veg) P
! Wgeq -- Surface volumetric moisture
! C1, C2 -- Coefficients
!
! For detailed information about the surface energy budget model,
! see the articles in the reference list.
!
! The second-order Rouge-Kutta time integration scheme is used,
! which is described below.
!
! Assume a equation in the form of
!
! d X
! ----- = F(X, t)
! d t
!
! In the forward scheme, we have
!
! X(1) = X(0) + dt * F[X(0), t0]
!
! We split one time step into two halves, dt2 = dt/2, and use the
! forward scheme to calculate the first half step X(1/2).
!
! X(1/2) = X(0) + dt2 * F[X(0), t0]
!
! Then we can calculate the Right Hand Side (RHS) of the equation at
! the half step, F[X(1/2), t(1/2)]. Finally, we calculate the one
! step prediction, X(t1), by use of the average of F[X(0), t0] and
! F[X(1/2), t(1/2)].
!
! X(1) = X(0) + dt * 0.5 * { F[X(0),t0] + F[X(1/2), t(1/2)] }
! REFERENCES:
!
! Jacquemin, B. and J. Noilhan, 1990: Sensitivity Study and
! Validation of a Land Surface Parameterization Using the
! HAPEX-MOBILHP Data Set, Boundary-Layer Meteorology, 52,
! 93-134, (JN).
!
! Noilhan, J. and S. Planton, 1989: A Simple Parameterization of
! Land Surface Processes for Meteorological Model, Mon. Wea.
! Rev., 117, 536-549, (NP).
!
! Pleim, J. E. and A. Xiu, 1993: Development and Testing of a
! Land-Surface and PBL Model with Explicit Soil Moisture
! Parameterization, Preprints, Conf. Hydroclimat., AMS, 45-51,
! (PX).
!
! Bougeault, P., et al., 1991: An Experiment with an Advanced
! Surface Parameterization in a Mesobeta-Scale Model. Part I:
! Implementation, Monthly Weather Review, 119, 2358-2373.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu and Vince Wong
! 11/16/93
!
! MODIFICATION HISTORY:
!
! 10/30/94 (Y. Liu)
! Fixed a bug reported by Richard Carpenter.
!
! 11/01/1994 (Y. Liu)
! Subroutine name of soil model were changed from SFCEBM to SOILEBM
! and the arguments passed were changed to 2-D arrays instead of 3-D
! temporary arrays.
!
! 12/12/1994 (Y. Liu)
! Fixed a bug for the final calcultaion of qvsfc.
!
! 12/14/1994 (Y. Liu)
! Fixed a bug in phycst.inc for the water density value which was
! previously mistakingly set to 1 kg/m**3. The correct value should
! be 1000 kg/m**3. This bug largely influenced the integration of Wg
! and W2.
!
! 12/23/1994 (Y. Liu)
! Added the runoff calculation for Wr.
!
! 02/07/1995 (Yuhe Liu)
! Added a new 2-D permanent array, veg(nx,ny), to the soil model and
! at the same time delete the table data array veg(14).
!
! 03/27/1995 (Yuhe Liu)
! Changed the solor radiation used in the calculation of surface
! resistence factor F1 from the one at the top of atmosphere to the
! one at the surface.
!
! Changed the formula of calculating the surface resistence factor
! F3 to F3=1, instead of varying with qvsat(Tair) and qvair.
!
! 12/8/1998 (Donghai Wang and Vince Wong)
! Added a new 2-D permanent array, snowcvr(nx,ny), for snow cover.
! We just used a simple scheme to consider the snow cover process.
!
! 2000/01/10 (Gene Bassett)
! Snow cover (0 or 1) changed to snow depth (snowdpth). For simplicity
! a fractional value for snow cover is not used (simply say the grid
! point is completely covered with snow if snowdpth > snowdepth_crit,
! otherwise no snow).
!
! 2000/02/04 (Gene Bassett, Yang Kun)
! Fixed an error in tsoil integration (rhst2) causing tsoil to change
! by only 50% of what it should.
!
! 2001/12/07 (Diandong Ren, Ming Xue)
! Re-structured the code, moved the calculations of the right hand
! side terms of the soil-vegetation model into subroutine soilebm_frc.
!
! Soil seasonal temperature trend to be added.
!
! 2002/02/15 (Yunheng Wang)
! Changed wrmax to a 2-D array which was a bug found during mpi testing.
!
! 2002/06/9 (Ming Xue)
! Revoked some modifications to the soil moisture related caculations
! that Diandong Ren put in since IHOP_2 - the mods need more testing.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the z-direction (sfc/top)
!
! soiltyp Soil type at the horizontal grid points
! vegtyp Vegetation type at the horizontal grid points
! lai Leaf Area Index
! veg Vegetation fraction
!
! windsp Wind speed just above the surface (m/s)
! psfc Surface pressure (Pascal)
! rhoa Near sfc air density
! precip Precipitation flux reaching the surface (m/s)
!
! cdha Array for cdh, surface drag coefficient for heat
! cdqa Array for cdq, surface drag coefficient for moisture
!
! pres 3-dimensional pressure
! temp 3-dimensional temperature
! qv 3-dimensional specific humidity
!
! INPUT/OUTPUT:
!
! tsfc Temperature at ground surface (K)
! tdeep Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
!
! OUTPUT:
!
! qvsfc Effective S. H. at sfc.
!
! Local automatic work arrays for storing the right hand forcing terms
! of the soil model equations and for storing intermediate values of the
! soil state variables.
!
! frc_tsfc Temporary array
! frc_tdeep Temporary array
! frc_wsfc Temporary array
! frc_wdp Temporary array
! frc_wcnp Temporary array
!
! tsfcn Temporary array
! tdeepn Temporary array
! wgn Temporary array
! w2n Temporary array
! wrn Temporary array
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: windsp(nx,ny) ! Wind speed
REAL :: psfc (nx,ny) ! Surface pressure
REAL :: rhoa (nx,ny) ! Air density near the surface
REAL :: precip(nx,ny) ! Precipitation rate at the surface
INTEGER :: soiltyp(nx,ny) ! Soil type at each point
INTEGER :: vegtyp (nx,ny) ! Vegetation type at each point
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny) ! Temperature at ground surface (K)
REAL :: tdeep (nx,ny) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny) ! Surface soil moisture
REAL :: wetdp (nx,ny) ! Deep soil moisture
REAL :: wetcanp(nx,ny) ! Canopy water amount
REAL :: qvsfc (nx,ny) ! Effective humidity at surface
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: tair (nx,ny) ! Surface air temperature (K)
REAL :: qvair (nx,ny) ! Surface air specific humidity (kg/kg)
REAL :: cdha (nx,ny) ! Surface drag coeff. for heat
REAL :: cdqa (nx,ny) ! Surface drag coeff. for moisture
REAL :: radsw (nx,ny) ! Solar radiation reaching the surface
REAL :: rnflx (nx,ny) ! Radiation flux at surface
REAL :: shflx (nx,ny) ! Sensible heat flux at surface
REAL :: lhflx (nx,ny) ! Latent heat flux at surface
REAL :: gflx (nx,ny) ! Ground diffusive heat flux
REAL :: ct (nx,ny) ! Soil thermal coefficient
REAL :: evaprg (nx,ny) ! Evaporation
REAL :: evaprtr(nx,ny) ! Transpiration from leaves
REAL :: evaprr (nx,ny) ! Direct evaporation from leaves
REAL :: f34 (nx,ny) ! Resistance factor of F3*F4
REAL :: qvsata (nx,ny) ! qvsat(tair) (kg/kg)
REAL :: qvsat (nx,ny) !
REAL :: frc_tsfc(nx,ny) ! Right hand side forcing for tsfc eq.
REAL :: frc_tdeep(nx,ny) ! Right hand side forcing for tdeep eq.
REAL :: frc_wsfc(nx,ny) ! Right hand side forcing for wetsfc eq.
REAL :: frc_wdp(nx,ny) ! Right hand side forcing for wetsp eq.
REAL :: frc_wcnp(nx,ny) ! Right hand side forcing for wetcanp eq.
REAL :: tsfcn(nx,ny) ! Temporary array, tsn
REAL :: tdeepn(nx,ny) ! Temporary array, t2n
REAL :: wgn(nx,ny) ! Temporary array, wetsfcNEW
REAL :: w2n(nx,ny) ! Temporary array, w2n
REAL :: wrn(nx,ny) ! Temporary array, wrn
REAL :: relief ! Difference between seasonal average skin and deep soil
!
!-----------------------------------------------------------------------
!
! Include files: globcst.inc and phycst.inc
!
!-----------------------------------------------------------------------
!
! Parameters and variables are defined in globcst.inc:
!
! dtsfc Surface model time step
! nsfcst # of surface model time steps
!
! moist Moist flag
!
! year Reference year
! month Reference month
! day Reference day
! jday Reference Julian day
! hour Hour of reference time
! minute Minute of reference time
! second Second of reference time
!
! latitud Latitude at the domain center
! longitud Longitude at the domain center
!
! curtim Current model time
! dtbig Length of big time step
!
! bslope Slope of the retention curve
! cgsat Soil thermal coefficient for bare ground at saturation
! cgv Soil thermal coef. for totally shielded ground by veg.
! pwgeq Coefficient of Wgeq formula. NP, Tab. 2
! awgeq Coefficient of Wgeq formula. NP, Tab. 2
! c1sat Value of C1 at saturation. NP, Tab. 2
! c2ref Value of C2 for W2 = .5 * Wsat. NP, Tab. 2
! wsat Saturated volumetric moisture content. JN, Tab. 1
! wfc Field capacity moisture. JN, Tab. 1
! wwlt Wilting volumetric moisture content. JN, Tab. 1
!
! Parameters and variables are defined in phycst.inc:
!
! solarc Solar constant (W/m**2)
! emissg Emissivity of the ground
! emissa Emissivity of the atmosphere
! sbcst Stefen-Boltzmann constant
!
! rhow Liquid water reference density (kg/m**3)
! rd Gas constant for dry air (kg/(m s**2))
! cp Gas heat capacity at constant pressure
! cv Gas heat capacity at constant volume
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
REAL :: pi ! Pi
PARAMETER (pi = 3.141592654)
! REAL :: tau ! Seconds of a day = 24. * 3600.
! PARAMETER (tau = 86400.)
REAL :: dtsfc2 ! Length of half time step in SFCEBM,
! dtsfc2 = dtsfc/2.
REAL :: log100 ! Constant: alog(100)
! dependent distance from the earth to the sun
REAL :: cg ! Soil thermal coefficient for bare ground
REAL :: rhsts(nx,ny) ! Right hand side of Eq. for Ts at current time
REAL :: rhst2(nx,ny) ! Right hand side of Eq. for T2 at current time
REAL :: rhswg(nx,ny) ! Right hand side of Eq. for Wg at current time
REAL :: rhsw2(nx,ny) ! Right hand side of Eq. for W2 at current time
REAL :: rhswr(nx,ny) ! Right hand side of Eq. for Wr at current time
REAL :: wrmax(nx,ny) ! Maximum value for canopy moisture, wetcanp
REAL :: c1wg ! Coefficient in the surface moisture Eq. of Wg
REAL :: c2wg ! Coefficient in the surface moisture Eq. of Wg
REAL :: wgeq ! Surface moisture when gravity balances the capillarity
REAL :: wr2max ! Tendency to reach the maximum wrmax
REAL :: runoff ! Runoff of the interception reservoir.
REAL :: pnet ! Residual of precip. and evap.
REAL :: vegp ! Precip. intercepted by vegetation
INTEGER :: i, j, it
REAL :: tema
LOGICAL :: firstcall ! First call flag of this subroutine
! SAVE firstcall, log100, dtsfc2, tauinv
SAVE firstcall, log100, dtsfc2
DATA firstcall/.true./
INTEGER :: jday_min ! offset value from Jan 01.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (firstcall) THEN
log100 = ALOG(100.0)
dtsfc2 = dtsfc/2.0
firstcall = .false.
END IF
IF ( moist /= 0 ) THEN
! Calculate saturated specific humidity near the surface, qvsata.
CALL getqvs
(nx,ny,1, 1,nx-1,1,ny-1,1,1, psfc,tair,qvsata)
DO j = 1, ny-1
DO i = 1, nx-1
f34(i,j) = MAX( 0., 1.0 - 0.0016 * (298.0-tair(i,j))**2 )
END DO
END DO
END IF
jday_min = 61 ! may change with latitude
relief=tsoil_offset_amplitude*sin((jday-jday_min)/365.0*2.0*PI)
!
!-----------------------------------------------------------------------
!
! Start time integration loop
!
!-----------------------------------------------------------------------
!
DO it = 1, nsfcst
CALL soilebm_frc(nx,ny,nz,soiltyp,vegtyp,lai,veg, &
tsfc,tdeep,wetsfc,wetdp,wetcanp,snowdpth, &
qvsfc,windsp,psfc,rhoa,precip,tair,qvair,cdha,cdqa, &
radsw,rnflx,shflx,lhflx,gflx,ct,evaprg,evaprtr, &
evaprr,qvsat,qvsata,f34, &
frc_tsfc,frc_tdeep,frc_wsfc,frc_wdp,frc_wcnp,wrmax,relief)
DO j = 1, ny-1
DO i = 1, nx-1
tsfcn (i,j) = tsfc(i,j) + dtsfc2 * frc_tsfc(i,j)
tdeepn(i,j) = tdeep(i,j)+ dtsfc2 * frc_tdeep(i,j)
END DO
END DO
IF ( moist /= 0 ) THEN
DO j = 1, ny-1
DO i = 1, nx-1
wgn(i,j) = wetsfc(i,j) + dtsfc2 * frc_wsfc(i,j)
wgn(i,j) = MAX(wgn(i,j), 0.0 )
wgn(i,j) = MIN(wgn(i,j), wsat(soiltyp(i,j)) )
w2n(i,j) = wetdp(i,j) + dtsfc2 * frc_wdp(i,j)
w2n(i,j) = MAX( w2n(i,j), 0.0 )
w2n(i,j) = MIN(w2n(i,j), wsat(soiltyp(i,j)) )
wrn(i,j) = wetcanp(i,j) + dtsfc2 * frc_wcnp(i,j)
wrn(i,j) = MAX(wrn(i,j), 0.0 )
wrn(i,j) = MIN(wrn(i,j), wrmax(i,j) )
END DO
END DO
END IF
DO j = 1, ny-1
DO i = 1, nx-1
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 ) THEN
rhsts(i,j) = 0.0
rhst2(i,j) = 0.0
ELSE
rhsts(i,j)=frc_tsfc(i,j)
rhst2(i,j)=frc_tdeep(i,j)
END IF
END DO
END DO
IF ( moist /= 0 ) THEN
DO j = 1, ny-1
DO i = 1, nx-1
rhswg(i,j)=frc_wsfc(i,j)
rhsw2(i,j)=frc_wdp(i,j)
rhswr(i,j)=frc_wcnp(i,j)
END DO
END DO
END IF
CALL soilebm_frc(nx,ny,nz,soiltyp,vegtyp,lai,veg, &
tsfcn,tdeepn,wgn,w2n,wrn,snowdpth, &
qvsfc,windsp,psfc,rhoa,precip,tair,qvair,cdha,cdqa, &
radsw,rnflx,shflx,lhflx,gflx,ct, &
evaprg,evaprtr,evaprr,qvsat,qvsata,f34, &
frc_tsfc,frc_tdeep,frc_wsfc,frc_wdp,frc_wcnp,wrmax,relief)
! Integration for Ts and T2 at one time step.
DO j = 1, ny-1
DO i = 1, nx-1
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 ) THEN
rhsts(i,j) = 0.0
rhst2(i,j) = 0.0
ELSE
rhsts(i,j) = 0.5 * (frc_tsfc(i,j) +rhsts(i,j))
rhst2(i,j) = 0.5 * (frc_tdeep(i,j)+rhst2(i,j))
END IF
tsfc(i,j) = tsfc(i,j) + dtsfc * rhsts(i,j)
tdeep(i,j) = tdeep(i,j)+ dtsfc * rhst2(i,j)
END DO
END DO
IF ( moist /= 0 ) THEN
DO j = 1, ny-1
DO i = 1, nx-1
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 .OR. &
snowdpth(i,j) >= snowdepth_crit ) THEN
rhswg(i,j) = 0.0
rhsw2(i,j) = 0.0
rhswr(i,j) = 0.0
ELSE
rhswg(i,j) = 0.5 * (frc_wsfc(i,j)+rhswg(i,j))
rhsw2(i,j) = 0.5 * (frc_wdp (i,j)+rhsw2(i,j))
rhswr(i,j) = 0.5 * (frc_wcnp(i,j)+rhswr(i,j))
END IF
wetsfc(i,j) = wetsfc(i,j) + dtsfc * rhswg(i,j)
wetsfc(i,j) = MAX( wetsfc(i,j), 0.0 )
wetsfc(i,j) = MIN( wetsfc(i,j), wsat(soiltyp(i,j)) )
wetdp(i,j) = wetdp(i,j) + dtsfc * rhsw2(i,j)
wetdp(i,j) = MAX( wetdp(i,j), 0.0 )
wetdp(i,j) = MIN( wetdp(i,j), wsat(soiltyp(i,j)) )
wetcanp(i,j) = wetcanp(i,j) + dtsfc * rhswr(i,j)
wetcanp(i,j) = MAX( wetcanp(i,j), 0.0 )
wetcanp(i,j) = MIN( wetcanp(i,j), wrmax(i,j) )
END DO
END DO
END IF
END DO ! TIME INTEGRATION
DO j = 1, ny-1 !SOIL
DO i = 1, nx-1
tema = MAX( wetsfc(i,j), wwlt(soiltyp(i,j)) )
IF (snowdpth(i,j) >= snowdepth_crit) THEN
ct(i,j) = cg_snow
gflx(i,j) = 2.0*pi*(tsfc(i,j)-tdeep(i,j)-relief) &
*snowflxfac/(tau*ct(i,j))
! Snow cover
ELSE
cg = cgsat(soiltyp(i,j)) &
* ( wsat(soiltyp(i,j))/tema ) &
**( bslope(soiltyp(i,j))/log100 )
ct(i,j) = cg * cgv / ( (1.0-veg(i,j)) * cgv &
+ veg(i,j) * cg ) ! NP, Eq. 8
gflx(i,j) = 2.0*pi*(tsfc(i,j)-tdeep(i,j)-relief)/(tau*ct(i,j))
! Ground diffusive heat flux
END IF
shflx(i,j) =rhoa(i,j) * cp * cdha(i,j) * windsp(i,j) &
* ( tsfc(i,j) - tair(i,j) *(psfc(i,j)/1.0E5)**rddcp ) !RDDRDD
END DO
END DO
CALL getqvs(nx,ny,1, 1,nx-1,1,ny-1,1,1, psfc,tsfc,qvsat) !HYDROLOGY
CALL evapflx(nx,ny,radsw,f34,cdqa,windsp,psfc,rhoa,qvair, &
soiltyp,vegtyp,lai,veg,tsfc,wetsfc,wetdp,wetcanp, &
snowdpth,evaprg,evaprtr,evaprr,lhflx,qvsat)
DO j=1, ny-1
DO i=1, nx-1
qvsfc(i,j) = lhflx(i,j) + qvair(i,j)
evaprg (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprg(i,j)
evaprtr(i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprtr(i,j)
evaprr (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprr(i,j)
lhflx (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*lhflx(i,j) !MOISTURE
lhflx (i,j) = lhflx(i,j) * lathv !QE
END DO
END DO
RETURN
END SUBROUTINE soilebm
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOILEBM_FRC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE soilebm_frc(nx,ny,nz,soiltyp,vegtyp,lai,veg, & 2,2
tsfc,tdeep,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,windsp,psfc, &
rhoa,precip,tair,qvair, cdha,cdqa,radsw, rnflx,shflx,lhflx, &
gflx,ct,evaprg,evaprtr,evaprr,qvsat,qvsata,f34, &
frc_tsfc,frc_tdeep,frc_wsfc,frc_wdp,frc_wcnp,wrmax,relief)
!------------------------------------------------------------------
!
! AUTHOR: Diandong Ren (dd_ren@rossby.metr.ou.edu)
! Ming Xue (mxue@ou.edu)
! 12/08/2001
!
! 2002/02/15 (Yunheng Wang)
! Changed wrmax to a 2-D array which was a bug found during mpi testing.
!
!------------------------------------------------------------------
!
! PURPOSE:
! To calculate the right hand side forcing terms in the
! soil-vegetation model.
!------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: eta_fac,wmax,c1max,sigma_sqr
REAL :: windsp(nx,ny) ! Wind speed
REAL :: psfc (nx,ny) ! Surface pressure
REAL :: rhoa (nx,ny) ! Air density near the surface
REAL :: precip(nx,ny) ! Precipitation rate at the surface
INTEGER :: soiltyp(nx,ny) ! Soil type at each point
INTEGER :: vegtyp (nx,ny) ! Vegetation type at each point
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny) ! Temperature at ground surface (K)
REAL :: tdeep (nx,ny) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny) ! Surface soil moisture
REAL :: wetdp (nx,ny) ! Deep soil moisture
REAL :: wetcanp(nx,ny) ! Canopy water amount
REAL :: qvsfc (nx,ny) ! Effective humidity at surface
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: tair (nx,ny) ! Surface air temperature (K)
REAL :: qvair (nx,ny) ! Surface air specific humidity (kg/kg)
REAL :: cdha (nx,ny) ! Surface drag coeff. for heat
REAL :: cdqa (nx,ny) ! Surface drag coeff. for moisture
REAL :: radsw (nx,ny) ! Solar radiation reaching the surface
REAL :: rnflx (nx,ny) ! Radiation flux at surface
REAL :: shflx (nx,ny) ! Sensible heat flux at surface
REAL :: lhflx (nx,ny) ! Latent heat flux at surface
REAL :: gflx (nx,ny) ! Ground diffusive heat flux
REAL :: ct (nx,ny) ! Soil thermal coefficient
REAL :: evaprg (nx,ny) ! Evaporation
REAL :: evaprtr(nx,ny) ! Transpiration from leaves
REAL :: evaprr (nx,ny) ! Direct evaporation from leaves
REAL :: f34 (nx,ny) ! Resistance factor of F3*F4
REAL :: qvsata (nx,ny) ! qvsat(tair) (kg/kg)
REAL :: qvsat (nx,ny) !
REAL :: frc_tsfc(nx,ny) ! Right hand side forcing for tsfc eq.
REAL :: frc_tdeep(nx,ny) ! Right hand side forcing for tsoil eq.
REAL :: frc_wsfc(nx,ny) ! Right hand side forcing for wetsfc eq.
REAL :: frc_wdp(nx,ny) ! Right hand side forcing for wetsp eq.
REAL :: frc_wcnp(nx,ny) ! Right hand side forcing for wetcanp eq.
REAL :: wrmax(nx, ny) ! Maximum value for canopy moisture, wetcanp
REAL :: c1wg ! Coefficient in the surface moisture Eq. of Wg
REAL :: c2wg ! Coefficient in the surface moisture Eq. of Wg
REAL :: pi ! Pi
PARAMETER (pi = 3.141592654)
! REAL :: tau ! Seconds of a day = 24. * 3600.
! PARAMETER (tau = 86400.)
REAL :: dtsfc2 ! Length of half time step in SFCEBM,
! dtsfc2 = dtsfc/2.
REAL :: log100 ! Constant: alog(100)
! dependent distance from the earth to the sun
REAL :: cg ! Soil thermal coefficient for bare ground
REAL :: wgeq ! Surface moisture when gravity balances the capillarity
REAL :: wr2max ! Tendency to reach the maximum wrmax
REAL :: runoff ! Runoff of the interception reservoir.
REAL :: pnet ! Residual of precip. and evap.
REAL :: vegp ! Precip. intercepted by vegetation
INTEGER :: i, j, it
REAL :: tema
REAL :: relief ! Difference between seasonal average skin and deep soil
! temperature (t_skin - t_deeplayer).
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
log100 = ALOG(100.0)
dtsfc2 = dtsfc/2.0
DO j = 1, ny-1
DO i = 1, nx-1
tema = MAX( wetdp(i,j), wwlt(soiltyp(i,j)) )
IF (snowdpth(i,j) >= snowdepth_crit) THEN !snow cover
ct(i,j) = cg_snow
gflx(i,j) = 2.0*pi*(tsfc(i,j)-tdeep(i,j)-relief) &
*snowflxfac/(tau*ct(i,j))
ELSE
cg = cgsat(soiltyp(i,j)) &
* ( wsat(soiltyp(i,j))/tema ) &
**( bslope(soiltyp(i,j))/log100 )
ct(i,j) = cg * cgv / ( (1.0-veg(i,j)) * cgv + veg(i,j) * cg )
gflx(i,j) = 2.0*pi*(tsfc(i,j)-tdeep(i,j)-relief)/(tau*ct(i,j))
END IF
shflx(i,j) = rhoa(i,j) * cp * cdha(i,j) * windsp(i,j) &
* ( tsfc(i,j) - tair(i,j) )
END DO
END DO
IF ( moist == 0 ) THEN
DO j = 1, ny-1
DO i = 1, nx-1
evaprg(i,j) = 0.0
evaprtr(i,j) = 0.0
evaprr(i,j) = 0.0
lhflx(i,j) = 0.0
END DO
END DO
ELSE
CALL getqvs(nx,ny,1, 1,nx-1,1,ny-1,1,1, psfc,tsfc,qvsat)
CALL evapflx(nx,ny,radsw,f34,cdqa,windsp,psfc,rhoa,qvair, &
soiltyp,vegtyp,lai,veg,tsfc,wetsfc,wetdp,wetcanp, &
snowdpth,evaprg,evaprtr,evaprr,lhflx,qvsat)
END IF
DO j = 1, ny-1
DO i = 1, nx-1
evaprg (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprg(i,j)
evaprtr(i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprtr(i,j)
evaprr (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprr(i,j)
lhflx (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*lhflx(i,j)
lhflx(i,j) = lhflx(i,j) * lathv ! Latent heat flux
END DO
END DO
DO j = 1, ny-1
DO i = 1, nx-1
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 ) THEN
frc_tsfc (i,j) = 0.
frc_tdeep(i,j) = 0.
ELSE
frc_tsfc(i,j) = ct(i,j)*(rnflx(i,j)-shflx(i,j)-lhflx(i,j)-gflx(i,j))
IF ( snowdpth(i,j) >= snowdepth_crit ) THEN
frc_tdeep(i,j)= 1.03* (tsfc(i,j)-tdeep(i,j)-relief)*tauinv*snowflxfac ! Snow cover
ELSE
frc_tdeep(i,j)= 1.03*(tsfc(i,j) - tdeep(i,j)-relief) * tauinv
END IF
END IF
END DO
END DO
IF ( moist /= 0 ) THEN
DO j = 1, ny-1 !HYDROLOGY
DO i = 1, nx-1
wrmax(i,j) = 0.2 * 1.e-3 * veg(i,j) * lai(i,j) ! meter
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 .OR. &
snowdpth(i,j) >= snowdepth_crit ) THEN
frc_wsfc(i,j) = 0.
frc_wdp(i,j) = 0.
frc_wcnp(i,j) = 0.
ELSE
tema = MAX( wetsfc(i,j), wwlt(soiltyp(i,j)) )
c1wg = 0.4* c1sat(soiltyp(i,j)) &
* ( wsat(soiltyp(i,j)) / tema ) &
**( bslope(soiltyp(i,j)) / 2.0 + 1.0)
c2wg = c2ref(soiltyp(i,j)) * wetdp(i,j) &
/ ( wsat(soiltyp(i,j)) - wetdp(i,j) + wetsml )
wgeq = wetdp(i,j) - wsat(soiltyp(i,j)) &
* awgeq(soiltyp(i,j)) &
* ( wetdp(i,j) / wsat(soiltyp(i,j)) ) &
** pwgeq(soiltyp(i,j)) &
* ( 1.0 - ( wetdp(i,j) / wsat(soiltyp(i,j)) ) &
** ( 8 * pwgeq(soiltyp(i,j)) ) )
frc_wsfc(i,j) = c1wg * ( ( 1.0 - veg(i,j) ) &
* precip(i,j) - evaprg(i,j) ) &
/(rhow * d1) &
- c2wg * ( wetsfc(i,j) - wgeq ) * tauinv
frc_wdp(i,j) = ( ( 1.0 - veg(i,j) ) * precip(i,j) &
- evaprg(i,j) - evaprtr(i,j) ) &
/ ( rhow * d2 )
wr2max = ( wrmax(i,j) - wetcanp(i,j) ) / dtsfc2
vegp = veg(i,j) * precip(i,j)
pnet = vegp - evaprr(i,j)
tema = pnet - wr2max * rhow
runoff = MAX( tema, 0.0 )
vegp = vegp - runoff
frc_wcnp(i,j) = ( vegp - evaprr(i,j) ) / rhow
END IF
END DO
END DO
END IF
RETURN
END SUBROUTINE soilebm_frc
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EVAPFLX ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE evapflx(nx,ny, radsw, f34, cdqa, & 2
windsp,psfc,rhoa,qvair, &
soiltyp,vegtyp,lai,veg, &
tsfc,wetsfc,wetdp,wetcanp,snowdpth, &
evaprg,evaprtr,evaprr,qvflx,qvsat)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate:
!
! 1. Evaporation from ground surface,
!
! evaprg = rhoa * cdq * windsp * evaprg'
!
! where
!
! evaprg' = (1.-veg) * (rhgs * qvsats - qvair)
! Evaporation from the ground
! NP, Eq. 27
!
! 2. Direct evaporation from the fraction delta of the foliage covered
! by intercepted water.
!
! Er = rhoa * cdq * windsp * Er'
!
! Er' = delta * veg * (qvsats-qvair)
!
! 3. Transpiration of the remaining part (1-delta) of leaves,
!
! Etr = rhoa * cdq * windsp * Etr'
!
! where
!
! Etr' = veg * (1-delta) * Ra/(Ra+Rs) * (qvsats-qvair )
!
! and Ra is aerodynamic resistance and Rs is the surface resistance
!
! 1
! Ra = ----------------
! cdq * windsp
!
! Rsmin
! Rs = ------------------
! LAI*F1*F2*F3*F4
!
!
! f + Rsmin/Rsmax
! F1 = -------------------
! f + 1
!
! Rg 2
! f = 0.55 ----- -----
! Rgl LAI
!
! - 1, Wfc < W2
! |
! | W2 - Wwlt
! F2 = - ------------, Wwlt <= W2 <= Wfc
! | Wfc - Wwlt
! |
! - 0, W2 < Wwlt
!
!
! 1-0.06*(qvsats-qvair), qvsats-qvair <= 12.5 g/kg
! F3 = {
! 0.25, otherwise
!
!
! F4 = 1 - 0.0016 * (298-tair)**2
!
! 4. Water vapor flux, qvflx,
!
! qvflx = rhoa * cdq * windsp * qvflx'
!
! qvflx' = (Eg' + Etr' + Er') = (qvsfc - qvair)
!
! where qvsfc is the effective surface specific humidity
!
! (qvsfc - qvair) = (Eg' + Etr' + Er')
!
! This subroutine will solve Eg', Etr', Er', and qvflx'
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu and Vince Wong
! 4/20/94
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
!
! radsw Incoming solar radiation
!
! f34 f3*f4;
! f3 = Fractional conductance of atmospheric vapor pressure
! f4 = Fractional conductance of air temperature
!
! cdqa Array for surface drag coefficient for water vapor
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! veg Vegetation fraction
!
! tsfc Surface soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Vegetation moisture
!
! psfc Surface pressure
! qvair Specific humidity near the surface
! windsp Wind speed near the surface
! rhoa Air density near the surface
!
! OUTPUT:
!
! evaprp Evaporation from groud surface
! evaprtr Transpiration of the remaining part (1-delta) of leaves
! evaprr Direct evaporation from the fraction delta
! qvflx Water vapor flux
!
! WORK ARRAY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny
!
REAL :: radsw (nx,ny)
REAL :: f34 (nx,ny)
REAL :: cdqa (nx,ny)
INTEGER :: soiltyp(nx,ny)
INTEGER :: vegtyp (nx,ny)
REAL :: lai (nx,ny)
REAL :: veg (nx,ny)
REAL :: tsfc (nx,ny)
REAL :: wetsfc (nx,ny)
REAL :: wetdp (nx,ny)
REAL :: wetcanp(nx,ny)
REAL :: snowdpth(nx,ny)
!
REAL :: psfc (nx,ny)
REAL :: qvair (nx,ny)
REAL :: windsp (nx,ny)
REAL :: rhoa (nx,ny)
REAL :: evaprg (nx,ny)
REAL :: evaprtr(nx,ny)
REAL :: evaprr (nx,ny)
REAL :: qvflx (nx,ny)
REAL :: qvsat (nx,ny)
!
!-----------------------------------------------------------------------
!
! Include files: globcst.inc and phycst.inc
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
REAL :: pi
PARAMETER ( pi = 3.141592654 )
!
INTEGER :: i, j
REAL :: wrmax ! Maximum value for canopy moisture, wetcanp
REAL :: rstcoef ! Coefficient of resistance
REAL :: delta
REAL :: rhgs ! R.H. at ground surface
!
REAL :: tema
REAL :: temb
REAL :: pterm
REAL :: mterm
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j = 1, ny-1
DO i = 1, nx-1
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 ) THEN
!
!-----------------------------------------------------------------------
!
! Over water and ice, qvsfc should be saturated
!
!-----------------------------------------------------------------------
!
qvflx(i,j) = qvsat(i,j) - qvair(i,j)
evaprg (i,j) = qvflx(i,j)
evaprtr(i,j) = 0.0
evaprr (i,j) = 0.0
ELSE ! over land
wrmax = 0.2 * 1.e-3 * veg(i,j) * lai(i,j) ! in meter
!
!-----------------------------------------------------------------------
!
! In order to calculate the qv flux at the surface, we need to
! calculate some parameters like the resistance coefficient,
!
! rstcoef = rsta / (rsts + rsta)
!
! where rsta is the aerodynamic resistance and rsts is the surface
! resistances.
!
! rsta = 1 / ( cdq * Va )
!
! rsts = rsmin(vtyp) / ( lai(vtyp) * f1 * f2 * f3 * f4 )
!
! f3 * f4 is time-independent and has been calculated previously
! and stored in f34(i,j)
!
!-----------------------------------------------------------------------
!
! Calculate f1.
!
! f = 0.55*(radsw/rgl(vtyp))*(2./lai(vtyp)) ! NP, Eq. 34
! f1 = ( rsmin(vtyp)/rsmax + f ) / ( 1. + f ) ! NP, Eq. 34
!
! Note: the incoming solar radiation radsw is stored in radsw(i,j).
!
!-----------------------------------------------------------------------
!
IF ( lai(i,j) == 0. ) THEN ! No vegetation, I/W
rstcoef = 1.
ELSE
temb = 0.55 * ( radsw(i,j) / rgl(vegtyp(i,j)) ) &
* ( 2.0 / lai(i,j) )
rstcoef = ( rsmin(vegtyp(i,j)) / rsmax + temb ) &
/ ( 1.0 + temb )
END IF
!
!-----------------------------------------------------------------------
!
! Calculate f2 and f1*f2.
!
!-----------------------------------------------------------------------
!
pterm = .5 + SIGN( .5, wetdp(i,j) - wfc(soiltyp(i,j)) )
mterm = SIGN( .5, wetdp(i,j) - wwlt(soiltyp(i,j)) ) &
- SIGN( .5, wetdp(i,j) - wfc(soiltyp(i,j)) )
rstcoef = rstcoef * ( pterm + mterm &
* ( wetdp(i,j)-wwlt(soiltyp(i,j)) ) &
/ ( wfc(soiltyp(i,j)) - wwlt(soiltyp(i,j)) ) )
!
!-----------------------------------------------------------------------
!
! Calculate lai*f1*f2*f3*f4 where f3*f4 is stored in f34(i,j).
!
!-----------------------------------------------------------------------
!
rstcoef = lai(i,j)*rstcoef*f34(i,j) ! lai*f1*f2*f3*f4
!
!-----------------------------------------------------------------------
!
! Calculate the resistance coefficient, rsta/(rsts+rsta)
!
! rsts = rsmin(vtyp)/(lai(i,j)*f1*f2*f3*f4) ! Sfc. resistance
! rsta = 1./(cdh*va) ! NP, between Eq. 32 & 33
! rstcoef = rsta/(rsta+rsts)
! = 1/(1+rsts/rsta)
!
!-----------------------------------------------------------------------
!
tema = rsmin(vegtyp(i,j)) * cdqa(i,j) * windsp(i,j)
IF ( ABS(rstcoef) > 1.0E-30 ) THEN
rstcoef = 1.0 / (1.0 + tema/rstcoef)
END IF
!
!-----------------------------------------------------------------------
!
! 1. evaprg'
!
! evaprg' = (1.-veg) * (rhgs*qvsats - qvair)
! Evaporation from the ground
! NP, Eq. 27
!
! evaprg will be stored for current and future use to
! calculate the latent heat flux and soil moisture transports.
!
!-----------------------------------------------------------------------
!
pterm = .5 + SIGN( .5, wetsfc(i,j)-1.1*wfc(soiltyp(i,j)) )
rhgs = pterm &
+ (1.-pterm) * ( 0.25 * ( 1.0 - COS( wetsfc(i,j) &
* pi / (1.1*wfc(soiltyp(i,j))))) ** 2 )
! IF (snowdpth(i,j) .ge. snowdepth_crit) rhgs=1.0 !Snow cover
evaprg(i,j) = ( 1.0 - veg(i,j) ) &
* ( rhgs * qvsat(i,j) - qvair(i,j) )
!
!-----------------------------------------------------------------------
!
! 2. Transpiration of the remaining part (1-delta) of leaves, Etr',
!
! Etr' = (1-delta) * veg
! * Ra/(Ra+Rs) * ( qvsats - qvair )
!
! 3. Direct evaporation from the fraction delta, Er'
!
! Er' = delta * veg * ( qvsats - qvair )
!
! Er' and Etr' are stored for future use to calculate the latent heat
! flux and soil moisture transports.
!
!-----------------------------------------------------------------------
!
IF ( wrmax == 0.0 ) THEN
delta = 0.0
ELSE
delta = ( wetcanp(i,j) / wrmax ) ** 0.66666667
END IF
pterm = .5 + SIGN( .5, qvsat(i,j) - qvair(i,j) )
delta = pterm * delta + ( 1. - pterm )
tema = veg(i,j) * ( qvsat(i,j) - qvair(i,j) )
evaprtr(i,j) = ( 1.0 - delta ) * rstcoef * tema
evaprr (i,j) = delta * tema
!
!-----------------------------------------------------------------------
!
! 4. Water vapor flux, qvflx',
!
! qvflx' = evaprg' + evaprtr' + evaprr'
!
! qvflx will be saved for the future use to calculate the latent
! heat flux.
!
!-----------------------------------------------------------------------
!
qvflx(i,j) = evaprg(i,j) + evaprtr(i,j) + evaprr(i,j)
END IF
! NP, expl. Eq. 26-27
END DO
END DO
RETURN
END SUBROUTINE evapflx
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE OUSOIL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE ousoil(nx,ny,nz,nzsoil,zpsoil,j3soil,j3soilinv, & 1,12
soiltyp,vegtyp,lai,veg, &
tsoil,qsoil,wetcanp,snowdpth,qvsfc, &
windsp,psfc,rhoa,precip, &
tair,qvair,cdma, cdha,cdqa, &
radsw, rnflx, radswnet, radlwin, &
shflx,lhflx,gflx,ct,evaprg,evaprtr,evaprr,qvsat,qvsata,f34, &
tem1soil,tem2soil,tem3soil,tem4soil,tsdiffus,deltem,rrtem,temple)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Predict the soil surface temperature and moisture contents by solving
! the surface energy and moisture budget equtions:
!
! 1. the soil temperatures, tsoil(nx,ny,nz)
! 2. the soil moisture, qsoil(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! AUTHOR: Jerry Brotzge and Dan Weber
! 1/22/02
!
! Updated code to allow for implicit soil temp/moisture scheme
! and to allow for multiple options of soil/temp schemes.
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! where
!
! Tau -- 1 day in seconds = 86400 seconds
! PI -- number of PI = 3.141592654
! ROw -- Density of liquid water
! Veg -- Vegetation fraction
! Ct -- Thermal capacity
! Rn -- Radiation flux, rnflx
! H -- Sensible heat flux, shflx
! LE -- Latent heat flux, lhflx = latent*(Eg + Ev)
!
! Eg -- Evaporation from ground
! Ev -- Evapotranspiration from vegetation, Ev = Etr + Er
! Etr -- Transpiration of the remaining part of the leaves
! Er -- Evaporation directly from the foliage covered by
! intercepted water
!
! P -- Precipitation rates
! Pg -- Precipitation reaching the ground,
! Pg = (1 - Veg) P
!
! Wgeq -- Surface volumetric moisture
! C1, C2 -- Coefficients
!
!
! For detailed information about the surface energy budget model,
! see the articles in the reference list.
!
!-----------------------------------------------------------------------
!
! REFERENCES:
!
! Bougeault, P., et al., 1991: An Experiment with an Advanced
! Surface Parameterization in a Mesobeta-Scale Model. Part I:
! Implementation, Monthly Weather Review, 119, 2358-2373.
!
! Chen, F., and J. Dudhia, 2001: Coupling an Advanced Land Surface-
! Hydrology Model with the Penn State-NCAR MM5 Modeling System.
! Part 1: Model Implementation and Sensitivity, Mon Wea Rev.,
! 129, 569-585.
!
! Ek, M., and L. Mahrt, 1991: OSU 1-D PBL model user's guide.
! Version 1.04, 120 pp. [Available from the Dept of
! Atmospheric Sciences, Oregon State Univ., Corvallis, OR
! 97331-2209].
!
! Jacquemin, B. and J. Noilhan, 1990: Sensitivity Study and
! Validation of a Land Surface Parameterization Using the
! HAPEX-MOBILHP Data Set, Boundary-Layer Meteorology, 52,
! 93-134, (JN).
!
! Noilhan, J. and S. Planton, 1989: A Simple Parameterization of
! Land Surface Processes for Meteorological Model, Mon. Wea.
! Rev., 117, 536-549, (NP).
!
! Peters-Lidard, C.D., E. Blackburn, X. Liang, and E.F. Wood,
! 1998: The effect of soil thermal conductivity parameterization
! on surface energy fluxes and temperatures, J. Atmos. Sci.,
! 55, 1209-1224.
!
! Pleim, J. E. and A. Xiu, 1993: Development and Testing of a
! Land-Surface and PBL Model with Explicit Soil Moisture
! Parameterization, Preprints, Conf. Hydroclimat., AMS, 45-51,
! (PX).
!
! Smirnova, T. G., et al., 1997: Performance of Different Soil
! Model Configurations in Simulating Ground Surface Temperatures
! and Surface Fluxes, Mon. Wea. Rev., 125, 1870-1884.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the z-direction (sfc/top)
! nzsoil Number of grid points in the z-direction (sfc/soil substrate)
! zpsoil The physical depth defined at w-point in soil
! j3soil Coordinate transformation Jacobian, d(zpsoil)/d(zsoil)
! j3soilinv Inverse of j3soil
!
! soiltyp Soil type at the horizontal grid points
! vegtyp Vegetation type at the horizontal grid points
! lai Leaf Area Index
! veg Vegetation fraction
!
! windsp Wind speed just above the surface (m/s)
! psfc Surface pressure (Pascal)
! rhoa Near sfc air density
! precip Precipitation flux reaching the surface (m/s)
!
!
! cdma Array for cdm, surface drag coefficient for momentum
! cdha Array for cdh, surface drag coefficient for heat
! cdqa Array for cdq, surface drag coefficient for moisture
!
! pres 3-dimensional pressure
! temp 3-dimensional temperature
! qv 3-dimensional specific humidity
!
! INPUT and OUTPUT:
!
! tsfc Temperature at skin surface (K)
! qsfc Moisture at skin surface
! tsoil Soil temperatures (K)
! qsoil Soil moisture
! wetcanp Canopy wetness
!
! OUTPUT:
!
! qvsfc Effective S. H. at sfc.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: nzsoil ! Grid points in the soil
REAL :: zpsoil (nx,ny,nzsoil) ! The depth of the soil.
REAL :: j3soil (nx,ny,nzsoil) ! Coordinate transformation
REAL :: j3soilinv (nx,ny,nzsoil) ! Inverse of j3soil
REAL :: windsp(nx,ny) ! Wind speed
REAL :: psfc (nx,ny) ! Surface pressure
REAL :: rhoa (nx,ny) ! Air density near the surface
REAL :: precip(nx,ny) ! Precipitation rate at the surface
INTEGER :: soiltyp(nx,ny) ! Soil type at each point
INTEGER :: vegtyp (nx,ny) ! Vegetation type at each point
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny) ! Temperature at ground surface (K)
REAL :: qsfc (nx,ny) ! Moisture at ground surface
REAL :: qvsfc (nx,ny) ! Effective humidity at surface
REAL :: tsoil (nx,ny,nzsoil) ! Soil temperature at nzsoil levels
REAL :: qsoil (nx,ny,nzsoil) ! Soil moisture at nzsoil levels
REAL :: wetcanp(nx,ny) ! Canopy wetness
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: tair (nx,ny) ! Surface air temperature (K)
REAL :: qvair (nx,ny) ! Surface air specific humidity (kg/kg)
REAL :: cdma (nx,ny) ! Surface drag coeff. for momentum
REAL :: cdha (nx,ny) ! Surface drag coeff. for heat
REAL :: cdqa (nx,ny) ! Surface drag coeff. for moisture
REAL :: radsw (nx,ny) ! Solar radiation reaching the surface
REAL :: rnflx (nx,ny) ! Radiation flux at surface
REAL :: radswnet(nx,ny) ! Net solar radiation, SWin - SWout
REAL :: radlwin(nx,ny) ! Incoming longwave radiation to sfc
REAL :: shflx (nx,ny) ! Sensible heat flux at surface
REAL :: lhflx (nx,ny) ! Latent heat flux at surface
REAL :: gflx (nx,ny) ! Ground diffusive heat flux
REAL :: ct (nx,ny) ! Soil thermal coefficient
REAL :: evaprg (nx,ny) ! Evaporation
REAL :: evaprtr(nx,ny) ! Transpiration from leaves
REAL :: evaprr (nx,ny) ! Direct evaporation from leaves
REAL :: f34 (nx,ny) ! Resistance factor of F3*F4
REAL :: qvsata (nx,ny) ! qvsat(tair) (kg/kg)
REAL :: qvsat (nx,ny) !
REAL :: aatem ! Temporary array for est'g Pot Evapo.
REAL :: fftem (nx,ny) ! Temporary array for est'g Pot Evapo.
REAL :: rrtem (nx,ny) ! Temporary array for est'g Pot Evapo.
REAL :: deltem (nx,ny) ! Temporary array for est'g Pot Evapo.
REAL :: rnettot ! Temporary array for est'g Pot Evapo.
REAL :: evappot (nx,ny) ! Temporary array for est'g Pot Evapo.
REAL :: temple (nx,ny) ! Temporary array for est'g Pot Evapo.
REAL :: rsttemp (nx,ny) ! Canopy resistance used for LE,tr
REAL :: temwra (nx,ny) ! Temporary array for est'd wr
REAL :: temwrb (nx,ny) ! Temporary array for est'd wr
REAL :: tempbc ! Temporary array for est'd LE
!
!-----------------------------------------------------------------------
!
! Include files: globcst.inc and phycst.inc
!
!-----------------------------------------------------------------------
!
! Parameters and variables are defined in globcst.inc:
!
! dtsfc Surface model time step
! nsfcst # of surface model time steps
!
! moist Moist flag
!
! year Reference year
! month Reference month
! day Reference day
! jday Reference Julian day
! hour Hour of reference time
! minute Minute of reference time
! second Second of reference time
!
! latitud Latitude at the domain center
! longitud Longitude at the domain center
!
! curtim Current model time
! dtbig Length of big time step
!
! bslope Slope of the retention curve
! cgsat Soil thermal coefficient for bare ground at saturation
! cgv Soil thermal coef. for totally shielded ground by veg.
! wsat Saturated volumetric moisture content. JN, Tab. 1
! wfc Field capacity moisture. JN, Tab. 1
! wwlt Wilting volumetric moisture content. JN, Tab. 1
!
! Parameters and variables are defined in phycst.inc:
!
! solarc Solar constant (W/m**2)
! emissg Emissivity of the ground
! emissa Emissivity of the atmosphere
! sbcst Stefen-Boltzmann constant
!
! rhow Liquid water reference density (kg/m**3)
! rd Gas constant for dry air (kg/(m s**2))
! cp Gas heat capacity at constant pressure
! cv Gas heat capacity at constant volume
! lathv LH of vaporization at 0 deg C (Lv = 2.50, [m^2/s^2])
!
! tsoil0 Initial skin temperature
! qsoil0 Initial skin wetness
!
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
REAL :: pi ! Pi
PARAMETER (pi = 3.141592654)
REAL :: dtsfc2 ! Length of half time step in SFCEBM,
! dtsfc2 = dtsfc/2.
REAL :: log100 ! Constant: alog(100)
! dependent distance from the earth to the sun
REAL :: cg ! Soil thermal coefficient for bare ground
REAL :: wrmax ! Maximum value for canopy moisture, wetcanp
REAL :: wr2max ! Tendency to reach the maximum wrmax
REAL :: runoff ! Runoff of the interception reservoir.
REAL :: pnet ! Residual of precip. and evap.
REAL :: vegp ! Precip. intercepted by vegetation
INTEGER :: i, j, k, it
INTEGER :: temcount
REAL :: temsum
REAL :: psi ! Matric water potential
REAL :: psitemp ! Temporary array
REAL :: pf ! log(psi)
REAL :: cisoil(nx,ny) ! Volumetric heat capacity
REAL :: ktsoil ! Thermal conductivity
REAL :: ktsoiltop(nx,ny) ! Thermal conductivity, top soil level
REAL :: ghterm2 ! Ground storage term
REAL :: dtsoildt(nx,ny) ! Storage heat change w/time
REAL :: t2old(nx,ny) ! 2nd soil temp, stored
REAL :: totaldepth ! Thickness of total soil column
REAL :: tem1(nx,ny,nzsoil) ! Used for stretching term
REAL :: tsdiffus(nx,ny,nzsoil) ! Thermal diffusivity
REAL :: tem1soil (nx,ny,nzsoil) ! Tridiagonal parameters
REAL :: tem2soil (nx,ny,nzsoil) ! Tridiagonal parameters
REAL :: tem3soil (nx,ny,nzsoil) ! Tridiagonal parameters
REAL :: tem4soil (nx,ny,nzsoil) ! Tridiagonal parameters
REAL :: dampdepth(nx,ny,nzsoil) ! Damping depth (Garratt, pg 118)
REAL :: tempcoefa ! Temporary arrays for tridiag.
REAL :: tempcoefb ! Temporary arrays for tridiag.
REAL :: tempcoefc ! Temporary arrays for tridiag.
REAL :: tempcoefx1 ! Temporary arrays for tridiag.
REAL :: tempcoefx2 ! Temporary arrays for tridiag.
REAL :: tempcoefz1 ! Temporary arrays for tridiag.
REAL :: tempcoefz2 ! Temporary arrays for tridiag.
REAL :: tempbeta ! Used to estimate skin temperature
REAL :: tempbeta2 ! Used to estimate skin temperature
REAL :: potair (nx,ny) ! Potential air temperature
REAL :: tempcg ! Used to compute cg
REAL :: tskintema
REAL :: tskintemb
REAL :: tema, temb, temc
REAL :: rhswr
REAL :: desdt
REAL :: dqsdt
REAL :: dew ! Dew, forms when pot E < 0
LOGICAL :: firstcall ! First call flag of this subroutine
SAVE firstcall, log100, dtsfc2
DATA firstcall/.true./
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (firstcall) THEN
log100 = ALOG(100.0)
dtsfc2 = dtsfc/2.0
firstcall = .false.
END IF
!
!-----------------------------------------------------------------------
!
! Calculate saturated specific humidity near the surface, qvsata.
!
!-----------------------------------------------------------------------
!
IF ( moist /= 0 ) THEN
CALL getqvs(nx,ny,1, 1,nx-1,1,ny-1,1,1, psfc,tair,qvsata)
END IF
!-----------------------------------------------------------------------
!
! Converts standard soil temperatures to potential temperatures.
!
!-----------------------------------------------------------------------
DO k=1,nzsoil
DO j=1,ny-1
DO i=1,nx-1
tsoil(i,j,k) = tsoil(i,j,k) * (100000.0/psfc(i,j))**0.286
END DO
END DO
END DO
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of Implicit Scheme (JAB/DBW)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
CALL getqvs(nx,ny,1, 1,nx-1,1,ny-1,1,1, psfc,tsoil(1,1,1),qvsat)
!-----------------------------------------------------------------------
! Initialize heat capacity (C) and thermal conductivity (K) for each
! level
!-----------------------------------------------------------------------
! Compute soil thermal conductivity and heat capacity at each level k
CALL getcon(nx, ny, nzsoil, soiltyp, vegtyp, tsoil(1,1,1), &
qsoil(1,1,1), ktsoiltop, cisoil, tsdiffus)
!--------------------------------------------------------------------------
! Compute boundary conditions using the OSU Scheme
!--------------------------------------------------------------------------
!
!
! Compute the potential evaporation
!
! Ep = [ (Rn/rho/Cp/Ch) + (Ta-Ts)] * delta + (r + 1)*A
! ----------------------------------------------- * (rho*Cp*Ch/Lv)
! delta + r + 1.0
!
! delta = (qsfcsat - qairsat) / (Ts - Tair) * Lv / Cp
!
!
! r = 4.0 * sigma * Tair**4.0 * Rd
! -----------------------------
! Psfc * Cp * Ch
!
!-------------------------------------------------------------------------
temb = 4.0 * sbcst * rd/cp
temc = lathv * cpinv
DO j=1,ny-1
DO i=1,nx-1
tsfc(i,j) = tsoil(i,j,1)
cdha(i,j) = cdha(i,j) * MAX(windsp(i,j),0.1)
!-------------------------------------------
! Estimate GH flux term
!-------------------------------------------
gflx(i,j)=ktsoiltop(i,j)*dzsoilinv*j3soilinv(i,j,2)* &
(tsoil(i,j,1)-tsoil(i,j,2))
fftem(i,j) = radswnet(i,j) + radlwin(i,j)
rnettot = fftem(i,j) - (sbcst * tair(i,j)**4.0) - gflx(i,j)
rrtem(i,j) = temb * (tair(i,j)**4.0) / (psfc(i,j)* cdha(i,j))
aatem = temc * (qvsata(i,j) - qvair(i,j))
potair(i,j) = tair(i,j) * ((100000.0 / psfc(i,j)) ** 0.286)
!---------------------------------------------------------
! Compute Clausius-Clapyron Eqtn.
!---------------------------------------------------------
dqsdt = 0.0
desdt = 0.0
Do k = 7,2,-1
desdt = alpha(k)*(k-1) + (tair(i,j)-273.16)*desdt
END DO
dqsdt = 0.622 * desdt / psfc(i,j)
deltem(i,j) = temc * dqsdt
! potair(i,j) = tair(i,j) * ((100000.0 / psfc(i,j)) ** 0.286)
! evappot(i,j)=((rnettot/(rhoa(i,j)*cp*cdha(i,j))+(potair(i,j)-tair(i,j))) &
! * deltem(i,j) + ((rrtem(i,j) + 1.0) * aatem)) / ((deltem(i,j) + &
! rrtem(i,j) + 1.0)*lathv) * (rhoa(i,j) * cp * cdha(i,j))
evappot(i,j)=( ((rnettot/(rhoa(i,j)*cp*cdha(i,j)))+(potair(i,j)-tair(i,j))) &
* deltem(i,j) + ((rrtem(i,j) + 1.0) * aatem) + (aatem - deltem(i,j)* &
tair(i,j)) * (((100000.0/psfc(i,j))**0.286)-1.0) ) &
/ ( (deltem(i,j) + rrtem(i,j) + 1.0) &
+ (((100000.0/psfc(i,j))**0.286)-1.0) ) &
* (rhoa(i,j) * cp *cdha(i,j)/lathv)
!------------ Set upper limit on potential evaporation.
! tempcoefc = evappot(i,j) * 2500000.0
! IF (tempcoefc > rnettot) evappot(i,j) = rnettot/2500000.0
!------------ Computes dew; adds to precip total ----------------
IF (evappot(i,j) < 0) THEN
dew = -evappot(i,j)/rhow
precip(i,j) = precip(i,j) + dew
END IF
IF (soiltyp(i,j) /= 13) THEN
tempbeta = (qsoil(i,j,1) - wwlt(soiltyp(i,j))) / &
(wfc(soiltyp(i,j)) - wwlt(soiltyp(i,j)) )
ELSE IF (soiltyp(i,j) == 13) THEN
tempbeta = 1.0
END IF
!---------------Set upper limit on soil wetness
IF (tempbeta > 1.0) tempbeta = 1.0
!---------------Set lower limit on soil wetness
IF (tempbeta < 0.0) tempbeta = 0.0
temple(i,j) = (1.0 - veg(i,j)) * tempbeta * evappot(i,j)
END DO
END DO
!-------------------------------------------------------------
! Computes initial LE/H flux estimate
!-------------------------------------------------------------
CALL canres(nx,ny, radsw, f34, cdqa, &
windsp,psfc,rhoa,qvair, &
soiltyp,vegtyp,lai,veg,tair,tsoil(1,1,1), &
qsoil,wetcanp,nzsoil,zpsoil,snowdpth, &
qvsata,rsttemp)
CALL leflx(nx,ny,cdha,cdqa,windsp, rhoa, qvair, veg, wetcanp, &
deltem,rrtem,temple, lhflx, evappot,soiltyp,evaprg, &
evaprtr,evaprr,qvsfc,qvsat,rsttemp)
!--------------------------------------------------------------
!-------------------------------------------------------------
! Computes skin temperature.
!-------------------------------------------------------------
CALL tskin(nx,ny, nzsoil, cdha, windsp, rhoa, tair, potair, &
rrtem, fftem, ktsoiltop, psfc, evappot, lhflx, tsoil)
!----------------------------------------------------------------------
! Estimate Soil Moisture
!-----------------------------------------------------------------------
!
CALL qdiff(nx,ny,nzsoil, j3soilinv, soiltyp, lai, veg, &
tsoil, qsoil, wetcanp, windsp, rhoa, precip, &
qvair, qvsata, evaprr, evaprg, tem1soil, tem2soil, tem3soil, &
tem4soil, dtsfc2)
!
!-----------------------------------------------------------------------
!
! Call the tridiagonal solver to solve the diffusion of soil moisture
! from k= 2, nzsoil-1.
!
!-----------------------------------------------------------------------
!
CALL tridiag2(nx,ny,nzsoil,1,nx-1,1,ny-1,2,nzsoil-1,tem1soil,tem2soil, &
tem3soil,tem4soil)
DO k=2,nzsoil-1 ! load the final result from tridiag2 into qsoil
DO j=1,ny-1
DO i=1,nx-1
qsoil(i,j,k) = tem4soil(i,j,k)
END DO
END DO
END DO ! qsoil is updated.........
!--------------------------------------------------------------------
! Compute bottom boundary condition
!---------------------------------------------------------------------
DO j=1,ny-1
DO i=1,nx-1
qsoil(i,j,nzsoil) = qsoil(i,j,nzsoil-1)
END DO
END DO
!------------------------------------------------------------------------
! Estimate Soil Temperature
!-----------------------------------------------------------------------
!
! C * dT/dt = d/dz (K * dT/dz)
!
!-----------------------------------------------------------------------
!
! Calculate coefficients of the tridigonal equation for soil temperature
!
!-----------------------------------------------------------------------
!
DO k=2,nzsoil-1
DO j=1,ny-1
DO i=1,nx-1
tempcoefa = cnbeta * dtsfc * 0.5 * dzsoilinv2 * j3soilinv(i,j,k)
tempcoefb = (1.0 - cnbeta) * dtsfc * 0.5 * dzsoilinv2 * j3soilinv(i,j,k)
tempcoefx1 = (tsdiffus(i,j,k-1) * j3soilinv(i,j,k-1)) + &
(tsdiffus(i,j,k ) * j3soilinv(i,j,k ))
tempcoefx2 = (tsdiffus(i,j,k ) * j3soilinv(i,j,k )) + &
(tsdiffus(i,j,k+1) * j3soilinv(i,j,k+1))
tem1soil(i,j,k) = -tempcoefa * tempcoefx1
tem2soil(i,j,k) = 1.0 + tempcoefa * (tempcoefx1 + tempcoefx2)
tem3soil(i,j,k) = -tempcoefa * tempcoefx2
tem4soil(i,j,k) = tempcoefb * tempcoefx1 * tsoil(i,j,k-1) + &
(1.0-tempcoefb*(tempcoefx1 + tempcoefx2))*tsoil(i,j,k) &
+ tempcoefb * tempcoefx2 *tsoil(i,j,k+1)
END DO
END DO
END DO
! Set the lower temperature boundary conditions
! note lower boundary is zero gradient
! (tsoil(i,j,nzsoil) = tsoil(i,j,nz-1)
DO j=1,ny-1 ! must set the boundry condition for tem2soil(i,j,nzsoil-1)
DO i=1,nx-1
tem4soil(i,j,2) = tem4soil(i,j,2) - tem1soil(i,j,2)*tsoil(i,j,1) !Top BC
tem2soil(i,j,nzsoil-1) = tem2soil(i,j,nzsoil-1) + tem3soil(i,j,nzsoil-1)
!Bottom, zero gradient
END DO
END DO ! coefficients are ready for the crank-nicholson
! solver.
!
!-----------------------------------------------------------------------
!
! Call the tridiagonal solver.
!
!-----------------------------------------------------------------------
!
CALL tridiag2(nx,ny,nzsoil,1,nx-1,1,ny-1,2,nzsoil-1,tem1soil, &
tem2soil,tem3soil,tem4soil)
DO k=2,nzsoil-1
DO j=1,ny-1
DO i=1,nx-1
tsoil(i,j,k) = tem4soil(i,j,k)
END DO
END DO
END DO ! soil temperature is up to date....
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! End of qsoil and tsoil implicit solving technique
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!---------------------------------------------------------------------
!
! Calculate the canopy wetness, wr
!
!----------------------------------------------------------------------
IF ( moist /= 0 ) THEN
DO j = 1, ny-1
DO i = 1, nx-1
wrmax = 0.2 * 1.e-3 * veg(i,j) * lai(i,j) ! meter
wr2max = (wrmax - wetcanp(i,j) )/ dtsfc
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 .OR. &
snowdpth(i,j) >= snowdepth_crit ) THEN
rhswr = 0.0
ELSE
vegp = veg(i,j) * precip(i,j)
pnet = vegp - evaprr(i,j)
tema = pnet - wr2max * rhow
runoff = MAX( tema, 0.0 )
vegp = vegp - runoff
rhswr = (vegp - evaprr(i,j) ) / rhow
END IF
wetcanp(i,j) = wetcanp(i,j) + dtsfc * rhswr
wetcanp(i,j) = MAX( wetcanp(i,j), 0.0 )
wetcanp(i,j) = MIN( wetcanp(i,j), wrmax )
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Calculate the heat capacity cg and ct, sensible heat flux shflx,
! and ground heat flux gflx
!
!-----------------------------------------------------------------------
!
DO j = 1, ny-1
DO i = 1, nx-1
tempcg = MAX( qsoil(i,j,1), wwlt(soiltyp(i,j)) ) !(JAB)
IF (snowdpth(i,j) >= snowdepth_crit) THEN !snow cover
ct(i,j) = cg_snow
gflx(i,j) = 2.0*pi*(tsoil(i,j,1)-tsoil(i,j,2)) &
*snowflxfac/(tau*ct(i,j))
! Snow cover
ELSE
!-------------------------------------------
! Estimate 2nd G flux for plotting
!-------------------------------------------
totaldepth = zrefsfc - zpsoil(i,j,nzsoil)
IF (totaldepth > 0.10) THEN
temsum = 0.0
temcount = 0
DO k=2,nzsoil
zrefsoil = zrefsfc - zpsoil(i,j,k)
IF (zrefsoil <= 0.10) THEN
temsum = temsum + (tsoil(i,j,k-1) - tsoil(i,j,k))
temcount = temcount + 1
END IF
END DO
IF (temcount > 0) THEN
gflx(i,j) = ktsoiltop(i,j)*j3soilinv(i,j,2)*dzsoilinv* &
(1.0/temcount) * temsum
ELSE IF (temcount == 0) THEN
gflx(i,j) = ktsoiltop(i,j)*j3soilinv(i,j,2)*dzsoilinv* &
0.5*(tsoil(i,j,1)-tsoil(i,j,2))
END IF
ELSE IF (totaldepth <= 0.10) THEN
gflx(i,j) = ktsoiltop(i,j)*j3soilinv(i,j,2)*dzsoilinv* &
0.5*(tsoil(i,j,1)-tsoil(i,j,2))
END IF
END IF
shflx(i,j) = rhoa(i,j)*cp*cdha(i,j)* &
( tsoil(i,j,1) - potair(i,j) ) ! Sensible heat flux
! NP, Eq. 26
tsfc(i,j) = tsoil(i,j,1)
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate the water vapor flux, qvflx, and hence the latent heat
! flux, lhflx. They share the same array lhflx(i,j).
!
! If moisture flag is off, all moisture fields are set to zero.
!
!-----------------------------------------------------------------------
!
! Calculate saturated specific humidity at the ground surface.
!
!-----------------------------------------------------------------------
CALL getqvs(nx,ny,1, 1,nx-1,1,ny-1,1,1, psfc,tsfc,qvsat)
!
!-----------------------------------------------------------------------
! Calculates new estimate of LE
!-----------------------------------------------------------------------
!
CALL canres(nx,ny, radsw, f34, cdqa, &
windsp,psfc,rhoa,qvair, &
soiltyp,vegtyp,lai,veg,tair,tsoil(1,1,1), &
qsoil,wetcanp,nzsoil,zpsoil,snowdpth, &
qvsata,rsttemp)
CALL leflx(nx,ny,cdha,cdqa, windsp, rhoa, qvair, veg, wetcanp, &
deltem,rrtem,temple, lhflx, evappot, soiltyp, &
evaprg,evaprtr,evaprr,qvsfc,qvsat,rsttemp)
DO j = 1, ny-1
DO i = 1, nx-1
cdha(i,j) = cdha(i,j) / MAX(windsp(i,j),0.1)
IF (lhflx(i,j) > (evappot(i,j)*lathv)) lhflx(i,j)=evappot(i,j)*lathv
qvsfc(i,j) = lhflx(i,j)/(lathv*rhoa(i,j)*cdqa(i,j)*windsp(i,j)) &
+ qvair(i,j)
IF (qvsfc(i,j) > qvsat(i,j)) qvsfc(i,j) = qvsat(i,j)
!-----------------------------------------------------------------------
!
! Converts potential soil temperatures to standard soil temperatures.
!
!-----------------------------------------------------------------------
DO k=1,nzsoil
tsoil(i,j,k) = tsoil(i,j,k) * ((100000.0 / psfc(i,j)) **(-0.286))
tsfc(i,j) = tsoil(i,j,1)
END DO
END DO
END DO
RETURN
END SUBROUTINE ousoil
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETCON ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE getcon(nx,ny,nzsoil, soiltyp, vegtyp, tsoil, qsoil, & 1
ktsoiltop, cisoil, tsdiffus)
!
!-----------------------------------------------------------------------
!
! PURPOSE: To calculate the thermal conductivity and diffusivity
!-----------------------------------------------------------------------
!
! AUTHOR: Jerry Brotzge
! 7/18/02
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nzsoil Number of grid points in the z-direction (soil)
! soiltyp Soil type
! vegtyp Vegetation type
! tsoil Soil temperatures (K)
! qsoil Soil moisture
!
! OUTPUT:
!
! cisoil ! Volumetric heat capacity
! ktsoiltop ! Thermal conductivity, top soil level
! tsdiffus ! Thermal diffusivity
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nzsoil
INTEGER :: soiltyp(nx,ny)
INTEGER :: vegtyp(nx,ny)
REAL :: qsoil (nx,ny,nzsoil)
REAL :: tsoil (nx,ny,nzsoil)
REAL :: tsdiffus(nx,ny,nzsoil) ! Thermal diffusivity
REAL :: cisoil(nx,ny) ! Volumetric heat capacity
REAL :: ktsoiltop(nx,ny) ! Thermal conductivity, top soil level
!
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
!
INTEGER :: i, j, k
REAL :: psi ! Matric water potential
REAL :: psitemp ! Temporary array
REAL :: pf ! log(psi)
REAL :: ktsoil ! Thermal conductivity
REAL :: totaldepth ! Thickness of total soil column
REAL :: satker ! Saturation of soil, qsoil/porosity
REAL :: kersten ! Kersten Number
REAL :: dryden ! Dry soil density [kg/m3]
REAL :: liqfrc ! Nonfrozen liquid fraction
REAL :: tcs ! Solids thermal conductivity
REAL :: tcko ! Minerals thermal conductivity
REAL :: tcdry ! Dry thermal conductivity
REAL :: tcsat ! Saturated thermal conductivity
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!----------------------------------------------------------------------
DO k=1,nzsoil ! set initial ground heat flux
DO j=1,ny-1
DO i=1,nx-1
!----------------------------------------------------------------------
! Option #1 for computing soil thermal conductivity
! Compute soil thermal conductivity using McCumber and Pielke (1981)
!
! Kt = 420 * EXP[ -(2.7 + Pf)] Pf <= 5.1
! = 0.1744 Pf > 5.1, Pf < 0
!
! Pf = log[ psi,sat * (Qsat / Q)**b ]
!
!----------------------------------------------------------------------
! psitemp = 0.0
!
! IF (qsoil(i,j,k) /= 0.0) THEN
! psitemp = wsat(soiltyp(i,j)) / qsoil(i,j,k)
! END IF
!
! psi = psisat(soiltyp(i,j))*(psitemp**bslope(soiltyp(i,j)) )
!
! pf = ALOG10(psi)
!
! IF (pf <= 5.1) ktsoil = 418.46*EXP(-(pf + 2.7))
!
! IF (pf > 5.1) ktsoil = 0.172
! IF (pf < 0.0) ktsoil = 0.172
!
!
! IF (ktsoil >= 1.9) ktsoil = 1.90
!
! IF (k == 1) ktsoiltop(i,j) = ktsoil
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Option #2 for computing soil thermal conductivity
! Compute soil thermal cond. using Johansen (1975); Peters-Lidard (1998)
!
! Kt = KE*(Ksat - Kdry) + Kdry
! Kdry = (0.135*gammadry + 64.7)/(2700 - 0.947*gammadry)
! gammadry = (1 - porosity) * 2700
!
! Ksat = Ks**(1-porosity) * Kice**(porosity-xu) * Kh2o**(xu)
! Ks = Kq**quartz * Ko**(1-quartz)
!
! KE = 0.7 * log(SR) + 1.0 SR > 0.05 (coarse)
! = log(SR) + 1.0 SR > 0.10 (fine)
!-------------------------------------------------------------------------
IF (soiltyp(i,j) < 12) THEN
satker = qsoil(i,j,k) / porosity(soiltyp(i,j))
ELSE
satker = 1.0
END IF
IF (satker > 1.0) satker = 1.0
IF (satker <= 0.1) satker = 0.11
kersten = ALOG10(satker) + 1.0 !Kersten Number, fine soils
IF (soiltyp(i,j) == 12) kersten = satker
dryden = (1.0 - porosity(soiltyp(i,j))) * 2700.0 ! [kg/m3]
tcdry = (0.135*dryden + 64.7) / (2700.0 - 0.947*dryden)
tcko = 3.0 ! [ W/m/K]
IF (quartz(soiltyp(i,j)) > 0.2) tcko = 2.0 ! [ W/m/K]
tcs = (7.7**quartz(soiltyp(i,j))) * (tcko**(1.0-quartz(soiltyp(i,j)) ))
liqfrc = 1.0 !Fraction of liquid water
IF (tsoil(i,j,k) < 273.15) liqfrc = 0.0
tcsat = (tcs ** (1.0 - porosity(soiltyp(i,j)) )) * &
(2.2 ** (porosity(soiltyp(i,j)) - liqfrc) ) * &
(0.57** (liqfrc))
ktsoil = kersten * (tcsat-tcdry) + tcdry
IF (ktsoil > 1.9) ktsoil = 1.90
IF (ktsoil < 0.2) ktsoil = 0.20
IF (k == 1) ktsoiltop(i,j) = ktsoil
!---------------------------------------------------------------------------
! Compute heat capacity as a function of soil moisture
!
! C = Q * C,h20 + (1 - Qsat)*C,soil + (Qsat - Q)*C,air
!
!-----------------------------------------------------------------------
IF (qsoil(i,j,k) > wfc(soiltyp(i,j))) qsoil(i,j,k) = wfc(soiltyp(i,j))
IF (qsoil(i,j,k) > 1.0) qsoil(i,j,k) = 1.0
cisoil(i,j) = (qsoil(i,j,k)*cwater) + &
(1.0 - wsat(soiltyp(i,j))) * csoil + &
(wsat(soiltyp(i,j)) - qsoil(i,j,k)) * cair
tsdiffus(i,j,k) = ktsoil / cisoil(i,j)
END DO
END DO
END DO
RETURN
END SUBROUTINE getcon
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TSKIN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE tskin(nx,ny, nzsoil, cdha, windsp, rhoa, tair, potair, & 1
rrtem, fftem, ktsoiltop, psfc, evappot, lhflx, tsoil)
!
!-----------------------------------------------------------------------
!
! PURPOSE: To estimate the top surface (skin) temperature, (Ek and Mahrt, 1991)
!
! Ts = (YY + ZZ * Tsoil(x,y,2)) / (ZZ + 1.0)
!
! FF = (1 - albedo) * SW,in + LW,in
!
! RCH = rho * Cp * Ch
!
! YY = Tair + [(FF-sigma*Tair**4)/RCH + (Tp-Tair)- Beta*(Lv*Ep/RCH)
! ----------------------------------------------------
! r + 1
!
! ZZ = Kt / (-0.5*zsoil(1) * RCH * (r + 1) )
!
!-------------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! AUTHOR: Jerry Brotzge
! 7/19/02
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! tsoil Soil temperatures (K)
! cdha Drag coefficient for heat
! windsp Wind speed
! ktsoiltop Soil conductivity
! evappot Potential evaporation
! lhflx Initial estimates of LE flux
! fftem Radiative forcing (SW,net + LW,in)
!
! OUTPUT:
!
! tsoil(,,1) ! Potential skin temperature (K)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nzsoil
REAL :: tsoil (nx,ny,nzsoil)
REAL :: ktsoiltop(nx,ny) ! Thermal conductivity, top soil level
REAL :: cdha(nx,ny) ! Drag coefficient
REAL :: windsp(nx,ny) ! Wind speed
REAL :: rhoa(nx,ny) ! Density of air
REAL :: tair(nx,ny) ! Air temperature (K)
REAL :: potair(nx,ny) ! Potential temp of air (K)
REAL :: rrtem(nx,ny) !
REAL :: fftem(nx,ny) ! Radiative forcing, SW net + LW in
REAL :: psfc(nx,ny) ! Atmospheric pressure (pa)
REAL :: evappot(nx,ny) ! Potential evaporation
REAL :: lhflx(nx,ny) ! Latent heat flux (W/m2)
!
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
!
INTEGER :: i, j, k
REAL :: aatemp ! Temporary array
REAL :: bbtemp ! Temporary array
REAL :: cctemp ! Temporary array
REAL :: ddtemp ! Temporary array
REAL :: eetemp ! Temporary array
REAL :: ggtemp ! Temporary array
REAL :: hhtemp ! Temporary array
REAL :: jjtemp ! Temporary array
REAL :: rrtemp ! Temporary array
REAL :: sstemp ! Temporary array
REAL :: zztemp ! Temporary array
! REAL :: dzsoilinv
REAL :: tempbeta2
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j=1,ny-1
DO i=1,nx-1
IF (lhflx(i,j) > (evappot(i,j)*lathv)) lhflx(i,j)=evappot(i,j)*lathv
tempbeta2 = 0.0
IF ((lhflx(i,j) /= 0.0).and.(evappot(i,j) /= 0.0)) THEN
tempbeta2 = lhflx(i,j) / (evappot(i,j)*lathv)
END IF
aatemp = fftem(i,j)
bbtemp = sbcst*(tair(i,j)**4.0)
cctemp = tempbeta2*lathv*evappot(i,j)
ddtemp = ktsoiltop(i,j) * tsoil(i,j,2) * dzsoilinv
eetemp = potair(i,j) - tair(i,j)
! Note that rrtem(i,j) is r
sstemp = 1.0 / (rhoa(i,j) * cp * cdha(i,j) )
rrtemp = 1.0 / (rrtem(i,j) + 1.0)
ggtemp = rrtemp*(eetemp + sstemp*(aatemp-bbtemp-cctemp+ddtemp))
zztemp = (rrtem(i,j) - 1.0 + ((psfc(i,j)/1.0E5)**(-0.286)) ) / rrtem(i,j)
! jjtemp = rrtemp *sstemp * (ktsoiltop(i,j) * dzsoilinv ) + zztemp
jjtemp = 1.0 + rrtemp *sstemp * (ktsoiltop(i,j) * dzsoilinv )
tsoil(i,j,1) = (tair(i,j) + ggtemp) / jjtemp !Original
! tsfc(i,j) = tsoil(i,j,1)
END DO
END DO
RETURN
END SUBROUTINE tskin
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE QDIFF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE qdiff(nx,ny,nzsoil, j3soilinv, soiltyp, lai, veg, & 1
tsoil, qsoil, wetcanp, windsp, rhoa, precip, &
qvair, qvsata, evaprr, evaprg, tem1soil, tem2soil, &
tem3soil, tem4soil, dtsfc2)
!
!-----------------------------------------------------------------------
!
! PURPOSE: To calculate the soil moisture profile and input matrix
! for tridiagonalization.
!-----------------------------------------------------------------------
!
! AUTHOR: Jerry Brotzge
! 3/06/02
!
! MODIFICATION HISTORY:
!
! (4/11/02) Dan Weber
! Cleaned up code and modified loop calculations to improve
! optimization.
!
!-----------------------------------------------------------------------
!
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nzsoil Number of grid points in the z-direction (soil)
!
! cdqa Array for surface drag coefficient for water vapor
! cdha Array for surface drag coefficient for heat
! veg Vegetation fraction
!
! qvair Specific humidity near the surface
! windsp Wind speed near the surface
! rhoa Air density near the surface
!
! OUTPUT:
!
! tem1soil Tridiagonal matrix
! tem2soil Tridiagonal matrix
! tem3soil Tridiagonal matrix
! tem4soil Tridiagonal matrix
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nzsoil
INTEGER :: soiltyp(nx,ny)
REAL :: precip (nx,ny)
REAL :: lai (nx,ny)
REAL :: veg (nx,ny)
REAL :: wetcanp(nx,ny)
REAL :: windsp (nx,ny)
REAL :: rhoa (nx,ny)
REAL :: evaprr (nx,ny)
REAL :: evaprg (nx,ny)
REAL :: qsdiffustop (nx,ny)
REAL :: qsdiffusm (nx,ny)
REAL :: qsdiffus (nx,ny)
REAL :: qsdiffusa (nx,ny)
REAL :: qsconducttop (nx,ny)
REAL :: qsconductm (nx,ny)
REAL :: qsconduct (nx,ny)
REAL :: qsconducta (nx,ny)
REAL :: qvsata (nx,ny)
REAL :: qvair (nx,ny)
REAL :: j3soilinv(nx,ny,nzsoil)
REAL :: qsoil (nx,ny,nzsoil)
REAL :: tsoil (nx,ny,nzsoil)
REAL :: tem1soil (nx,ny,nzsoil)
REAL :: tem2soil (nx,ny,nzsoil)
REAL :: tem3soil (nx,ny,nzsoil)
REAL :: tem4soil (nx,ny,nzsoil)
!
!-----------------------------------------------------------------------
! Include file: phycst.inc
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
!
INTEGER :: i, j, k
REAL :: tempcoefa
REAL :: tempcoefb
REAL :: tempcoefc
REAL :: tempcoefx1
REAL :: tempcoefx2
REAL :: tempcoefz1
REAL :: tempcoefz2
REAL :: tempbeta
REAL :: tempkq
REAL :: tempel
REAL :: tempws
REAL :: tempinf
REAL :: tema, temb, temc, temd
REAL :: qcond0
REAL :: qdiff0
REAL :: wrmax
REAL :: wr2max
REAL :: vegp
REAL :: pnet
REAL :: runoff
REAL :: dtsfc2
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!----------------------------------------------------------------------
! Estimate Soil Moisture
!-----------------------------------------------------------------------
!
! Compute moisture diffusional coefficient (Dn) and hydraulic (Kn)
! conductivity
!
! Note that the constants (Ksat, Psisat, Qsat, b) are all a f(Veg type)
!
! dQ/dt = d/dz( D * dQ/dz) + dK/dz
!
! D = - (b * Ksat * Psisat / Q) * (Q / Qsat)**(b+3)
!
!
! K = Ksat * (Q / Qsat )**(2*b+3)
!
!-----------------------------------------------------------------------
!------------------------------------------------------------------
! Compute top moisture boundary conditions
!------------------------------------------------------------------
tema = 0.2 * 1.e-3
temb = 1.0/dtsfc
temc = dtsfc * dzinv
DO j=1,ny-1
DO i=1,nx-1
qsdiffustop(i,j) = - (bslope(soiltyp(i,j)) * kns(soiltyp(i,j)) * &
wsat(soiltyp(i,j)) / qsoil(i,j,1)) * &
( (qsoil(i,j,1)/wsat(soiltyp(i,j))) ** (bslope(soiltyp(i,j)) &
+ 3.0) )
qsconducttop(i,j) = (kns(soiltyp(i,j)) / wsat(soiltyp(i,j))) * &
( (qsoil(i,j,1)/wsat(soiltyp(i,j)) ) ** (2.0 * &
bslope(soiltyp(i,j)) + 2.0) )
tempel = evaprg(i,j)
tempws = rhow * ( (qsdiffustop(i,j) * (qsoil(i,j,1) - &
qsoil(i,j,2)) * dzsoilinv) + qsconducttop(i,j) )
wrmax = tema * veg(i,j) * lai(i,j) ! m
wr2max = temb * ( wrmax - wetcanp(i,j)) ! m/s
vegp = veg(i,j) * precip(i,j) ! m/s
pnet = vegp - evaprr(i,j) ! m/s
temd = pnet - wr2max * rhow ! m/s
runoff = MAX( temd, 0.0 ) ! m/s
tempinf = (precip(i,j)-runoff) * rhow ! kg/m/m/s
qsoil(i,j,1) = qsoil(i,j,1) + temc * (tempws + &
tempel + tempinf)
IF (qsoil(i,j,1) > wfc(soiltyp(i,j))) qsoil(i,j,1) = wfc(soiltyp(i,j))
IF (qsoil(i,j,1) > 1.0) qsoil(i,j,1) = 1.0
END DO
END DO
!------------------------------------------------------------------------
!
! Compute moisture diffusivity and conductivity terms
! See Smirnova et al.(1997)
!
! Qdiff = (-b*Knsat*psisat/theta)*(theta/thetasat)**(b+3)
!
! Qcond = Knsat * (theta/thetasat)**(2b+3)
!
!------------------------------------------------------------------------
tema = dtsfc*0.5*dzsoilinv2
DO k=2,nzsoil-1
tempcoefc = dtsfc * dzsoilinv
DO j=1,ny-1
DO i=1,nx-1
tempcoefa = cnbeta*tema*j3soilinv(i,j,k)
tempcoefb = tema*j3soilinv(i,j,k)-tempcoefa
qsdiffusm(i,j) = (bslope(soiltyp(i,j)) * kns(soiltyp(i,j)) * &
psisat(soiltyp(i,j)) / qsoil(i,j,k-1)) * &
( (qsoil(i,j,k-1)/wsat(soiltyp(i,j))) ** (bslope(soiltyp(i,j)) &
+ 3.0) )
qsdiffus (i,j) = (bslope(soiltyp(i,j)) * kns(soiltyp(i,j)) * &
psisat(soiltyp(i,j)) / qsoil(i,j,k)) * &
( (qsoil(i,j,k)/wsat(soiltyp(i,j))) ** (bslope(soiltyp(i,j)) &
+ 3.0) )
qsdiffusa(i,j) = (bslope(soiltyp(i,j)) * kns(soiltyp(i,j)) * &
psisat(soiltyp(i,j)) / qsoil(i,j,k+1)) * &
( (qsoil(i,j,k+1)/wsat(soiltyp(i,j))) ** (bslope(soiltyp(i,j)) &
+ 3.0) )
qsconductm(i,j) = (kns(soiltyp(i,j)) / wsat(soiltyp(i,j))) * &
( (qsoil(i,j,k-1)/wsat(soiltyp(i,j)) ) ** (2.0 * &
bslope(soiltyp(i,j)) + 2.0) )
qsconduct (i,j) = (kns(soiltyp(i,j)) / wsat(soiltyp(i,j))) * &
( (qsoil(i,j,k)/wsat(soiltyp(i,j)) ) ** (2.0 * &
bslope(soiltyp(i,j)) + 2.0) )
qsconducta(i,j) = (kns(soiltyp(i,j)) / wsat(soiltyp(i,j))) * &
( (qsoil(i,j,k+1)/wsat(soiltyp(i,j)) ) ** (2.0 * &
bslope(soiltyp(i,j)) + 2.0) )
END DO
END DO
!------------------------------------------------------------------------
!
! Compute coefficients for tridiagonalization of moisture variables
!
!------------------------------------------------------------------------
DO j=1,ny-1
DO i=1,nx-1
tempcoefx1 = (qsdiffusm(i,j) * j3soilinv(i,j,k-1)) + &
(qsdiffus(i,j) * j3soilinv(i,j,k ))
tempcoefx2 = (qsdiffus(i,j) * j3soilinv(i,j,k )) + &
(qsdiffusa(i,j) * j3soilinv(i,j,k+1))
tempcoefz1 = tempcoefc * qsconductm(i,j) * j3soilinv(i,j,k-1)
tempcoefz2 = tempcoefc * qsconducta(i,j) * j3soilinv(i,j,k+1)
tem1soil(i,j,k) = - tempcoefa * tempcoefx1 + tempcoefz1
tem2soil(i,j,k) = 1.0 + tempcoefa * (tempcoefx1 + tempcoefx2)
tem3soil(i,j,k) = - tempcoefa * tempcoefx2 - tempcoefz2
tem4soil(i,j,k) = tempcoefb * tempcoefx1 * qsoil(i,j,k-1) + &
(1.0-tempcoefb*(tempcoefx1 + tempcoefx2)) * qsoil(i,j,k) &
+ (tempcoefb * tempcoefx2) * qsoil(i,j,k+1)
END DO
END DO
END DO ! end of k loop, the moisture coefficients are set
!-----------------------------------------------------------------
! Set upper and lower moisture boundary conditions
! note lower boundary is zero gradient
! (qsoil(i,j,nzsoil) = qsoil(i,j,nz-1)
!-----------------------------------------------------------------
DO j=1,ny-1 ! must set the boundry condition for tem2soil(i,j,nzsoil-1)
DO i=1,nx-1
tem4soil(i,j,2) = tem4soil(i,j,2)-tem1soil(i,j,2)*qsoil(i,j,1) !Top BC
tem2soil(i,j,nzsoil-1) = tem2soil(i,j,nzsoil-1) + &
tem3soil(i,j,nzsoil-1) !Bottom, zero gradient
END DO
END DO
RETURN
END SUBROUTINE qdiff
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EVAPFLX2 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE canres(nx,ny, radsw, f34, cdqa, & 2
windsp,psfc,rhoa,qvair, &
soiltyp,vegtyp,lai,veg, &
tair,tsoil,qsoil,wetcanp,nzsoil,zpsoil,snowdpth, &
qvsata,rstcoef)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate canopy resistance.
!
!
! 1
! Ra = ----------------
! cdq * windsp
!
! Rsmin
! Rs = ------------------
! LAI*F1*F2*F3*F4
!
!
! f + Rsmin/Rsmax
! F1 = -------------------
! f + 1
!
! Rg 2
! f = 0.55 ----- -----
! Rgl LAI
!
! - 1, Wfc < W2
! |
! | W2 - Wwlt
! F2 = - ------------, Wwlt <= W2 <= Wfc
! | Wfc - Wwlt
! |
! - 0, W2 < Wwlt
!
!
! 1-0.06*(qvsats-qvair), qvsats-qvair <= 12.5 g/kg
! F3 = {
! 0.25, otherwise
!
!
! F4 = 1 - 0.0016 * (298-tair)**2
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu and Vince Wong
! 4/20/94
!
! MODIFICATION HISTORY:
!
! 07/16/02 Jerry Brotzge
! Modified according to Chen and Dudhia (2001, MWR)
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nzsoil Number of grid points in the soil
!
! radsw Incoming solar radiation
!
! f34 f3*f4;
! f3 = Fractional conductance of atmospheric vapor pressure
! f4 = Fractional conductance of air temperature
!
! cdqa Array for surface drag coefficient for water vapor
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! veg Vegetation fraction
!
! tsfc Surface soil temperature (K)
! tsoil Soil temperatures
! qsoil Soil moistures
!
! psfc Surface pressure
! qvair Specific humidity near the surface
! windsp Wind speed near the surface
! rhoa Air density near the surface
!
! OUTPUT:
!
! rstcoef Canopy resistance (s/m)
!
! WORK ARRAY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz
INTEGER :: nzsoil
!
REAL :: radsw (nx,ny)
REAL :: f34 (nx,ny)
REAL :: cdqa (nx,ny)
INTEGER :: soiltyp(nx,ny)
INTEGER :: vegtyp (nx,ny)
REAL :: lai (nx,ny)
REAL :: veg (nx,ny)
REAL :: tsoil (nx,ny,nzsoil)
REAL :: qsoil (nx,ny,nzsoil)
REAL :: wetcanp(nx,ny)
REAL :: sumgdz(nx,ny)
REAL :: zpsoil (nx,ny,nzsoil)
REAL :: tsfc (nx,ny)
REAL :: tair (nx,ny)
REAL :: snowdpth(nx,ny)
REAL :: psfc (nx,ny)
REAL :: qvair (nx,ny)
REAL :: windsp (nx,ny)
REAL :: rhoa (nx,ny)
REAL :: evaprg (nx,ny)
REAL :: evaprtr(nx,ny)
REAL :: evaprr (nx,ny)
REAL :: qvflx (nx,ny)
REAL :: qvsat (nx,ny)
REAL :: qvsata (nx,ny)
REAL :: rstcoef(nx,ny)
!
!-----------------------------------------------------------------------
!
! Include files: globcst.inc and phycst.inc
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
REAL :: pi
PARAMETER ( pi = 3.141592654 )
!
INTEGER :: i, j, k
REAL :: wrmax ! Maximum value for canopy moisture, wetcanp
! REAL :: rstcoef (nx,ny) ! Coefficient of resistance
REAL :: delta
REAL :: rhgs ! R.H. at ground surface
!
REAL :: tema ! Used to compute Cg, heat capacity
REAL :: temb ! Used to compute f1*f2
REAL :: f2 ! Temporary variable for f2
REAL :: pterm
REAL :: mterm
REAL :: waf ! Available soil moisture fraction
REAL :: bw ! Half point of the available soil mstr frctn curve
REAL :: ps ! Shelter factor
REAL :: beta2 ! Used to estimate f1*f2
REAL :: delzneg
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j=1,ny-1
DO i=1,nx-1
sumgdz(i,j) = 0.0
rstcoef(i,j) = 0.0
END DO
END DO
DO j = 1, ny-1
DO i = 1, nx-1
IF ( soiltyp(i,j) /= 12 .AND. soiltyp(i,j) /= 13) THEN
wrmax = 0.2 * 1.e-3 * veg(i,j) * lai(i,j) ! in meter
!
!-----------------------------------------------------------------------
!
! In order to calculate the qv flux at the surface, we need to
! calculate some parameters like the resistance coefficient,
!
! rstcoef = rsta / (rsts + rsta)
!
! where rsta is the aerodynamic resistance and rsts is the surface
! resistances.
!
! rsta = 1 / ( cdq * Va )
!
! rsts = rsmin(vtyp) / ( lai(vtyp) * f1 * f2 * f3 * f4 )
!
! f3 * f4 is time-independent and has been calculated previously
! and stored in f34(i,j)
!
!-----------------------------------------------------------------------
!
! Calculate f1.
!
! f = 0.55*(radsw/rgl(vtyp))*(2./lai(vtyp)) ! NP, Eq. 34
! f1 = ( rsmin(vtyp)/rsmax + f ) / ( 1. + f ) ! NP, Eq. 34
!
! Note: the incoming solar radiation radsw is stored in radsw(i,j).
!
!-----------------------------------------------------------------------
!
IF ( lai(i,j) == 0. ) THEN ! No vegetation, I/W
rstcoef(i,j) = 1.
ELSE
IF (radsw(i,j) < 0.0 ) radsw(i,j) = 0.0
temb = 0.55 * ( radsw(i,j) / rgl(vegtyp(i,j)) ) &
* ( 2.0 / lai(i,j) )
rstcoef(i,j) = ( rsmin(vegtyp(i,j)) / rsmax + temb ) &
/ ( 1.0 + temb )
rstcoef(i,j) = MAX( 0.0001, rstcoef(i,j))
END IF
!---------------------------------------------------------------------------
!
! Calculate F4.
! Replaced force-restore wetdp estimate with integrated soil moisture (JAB)
! estimated as a function of root zone depth
! (See Chen and Dudhia MWR 2001)
!
!---------------------------------------------------------------------------
DO k = 2,nzsoil-1
! delzneg = (zpsoil(i,j,k)-zpsoil(i,j,k-1))/ &
! (zpsoil(i,j,nzsoil)-zpsoil(i,j,1))
IF ((zpsoil(i,j,1)-zpsoil(i,j,nzsoil)) <= rootzone(vegtyp(i,j))) &
rootzone(vegtyp(i,j)) = zpsoil(i,j,1) - zpsoil(i,j,nzsoil)
IF (zpsoil(i,j,k) >= (zpsoil(i,j,1)-rootzone(vegtyp(i,j))) ) THEN
delzneg = (zpsoil(i,j,k-1)-zpsoil(i,j,k))/ &
rootzone(vegtyp(i,j))
ELSE
delzneg = 0.0
END IF
beta2 = (qsoil(i,j,k) - wwlt(soiltyp(i,j))) / &
(wfc(soiltyp(i,j)) - wwlt(soiltyp(i,j)) )
IF (qsoil(i,j,k) > wfc(soiltyp(i,j)) ) beta2 = 1.0
IF (qsoil(i,j,k) <= wwlt(soiltyp(i,j))) beta2 = 0.0
sumgdz(i,j) = sumgdz(i,j) + (beta2 * delzneg)
END DO
sumgdz(i,j) = MAX( 0.0001, sumgdz(i,j))
rstcoef(i,j) = rstcoef(i,j) * sumgdz(i,j)
!
!-----------------------------------------------------------------------
!
! Calculate F2.
!
!-----------------------------------------------------------------------
f2 = 1.0 + hsf2(vegtyp(i,j))*(qvsata(i,j) - qvair(i,j) )
IF (qvair(i,j) > qvsata(i,j) .OR. f2 < 0) f2 = 1.0
IF (f2 > 1.0E-30) THEN
f2 = 1.0/f2
f2 = MAX( 0.01, f2)
ELSE IF (f2 <= 1.0E-30 .AND. f2 /= 1.0) THEN
f2 = 0.01
END IF
!
!-----------------------------------------------------------------------
!
! Calculate f3 * f4, stored in f34(i,j), where
!
! f3 -- Fractional conductance of atmospheric vapor pressure
! f3 = 1 - 0.06 * ( qvsata(i,j) - qvair ) * 1.e3 ! use kg/kg for qv
!
! (We use f3 = 1 in the code instead of the above formula.)
!
! f4 -- Fractional conductance of air temperature
! f4 = 1 - 0.0016 * ( 298 - tanem ) ** 2
!
! f3 * f4 will be used to calculate the resistance coefficient.
!
!-----------------------------------------------------------------------
!
! f3 = 1.0 ! f3 set to 1
! f34(i,j) = f3 * ( 1.0 - 0.0016 * ( 298. - tair(i,j) ) ** 2 )
!
!-----------------------------------------------------------------------
!
f34(i,j) = MAX( 0.0001, 1.0 - 0.0016 * (298.0-tair(i,j))**2 )
! f3 * f4, JN, Eq. 11
!
!-----------------------------------------------------------------------
!
! Calculate lai*f1*f2*f3*f4 where f3*f4 is stored in f34(i,j).
!
!-----------------------------------------------------------------------
!
rstcoef(i,j) = lai(i,j)*rstcoef(i,j)*f2*f34(i,j) ! lai*f1*f2*f3*f4
!
!-----------------------------------------------------------------------
!
! Calculate the resistance coefficient, rsta/(rsts+rsta)
!
! rsts = rsmin(vtyp)/(lai(i,j)*f1*f2*f3*f4) ! Sfc. resistance
! rsta = 1./(cdh*va) ! NP, between Eq. 32 & 33
! rstcoef = rsta/(rsta+rsts)
! = 1/(1+rsts/rsta)
!
!-----------------------------------------------------------------------
!
IF ( ABS(rstcoef(i,j)) > 1.0E-30) THEN
rstcoef(i,j) = rsmin(vegtyp(i,j)) / rstcoef(i,j)
END IF
END IF !soiltyp /= 12 and 13
END DO
END DO
RETURN
END SUBROUTINE canres
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE LEFLX ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE leflx(nx,ny, cdha, cdqa, windsp, rhoa, qvair, veg, wetcanp, & 2
deltem, rrtem, temple, lhflx, evappot, soiltyp, &
evaprg, evaprtr, evaprr, qvsfc, qvsat, rsttemp)
!
!-----------------------------------------------------------------------
!
! PURPOSE: To calculate the latent heat flux.
!
! Calculate:
!
! 1. Evaporation from ground surface,
!
! evaprg = rhoa * cdq * windsp * evaprg'
!
! where
!
! evaprg' = (1.-veg) * (rhgs * qvsats - qvair)
! Evaporation from the ground
! NP, Eq. 27
!
! 2. Direct evaporation from the fraction delta of the foliage covered
! by intercepted water.
!
! Er = rhoa * cdq * windsp * Er'
!
! Er' = delta * veg * (qvsats-qvair)
!
! 3. Transpiration of the remaining part (1-delta) of leaves,
!
! Etr = rhoa * cdq * windsp * Etr'
!
! where
!
! Etr' = veg * (1-delta) * Ra/(Ra+Rs) * (qvsats-qvair )
!
! and Ra is aerodynamic resistance and Rs is the surface resistance
!
! 1
! Ra = ----------------
! cdq * windsp
!
! Rsmin
! Rs = ------------------
! LAI*F1*F2*F3*F4
!
! 4. Water vapor flux, lhflx,
!
! lhflx = rhoa * cdq * windsp * qvflx'
!
! lhflx' = (Eg' + Etr' + Er') = (qvsfc - qvair)
!
! where qvsfc is the effective surface specific humidity
!
! (qvsfc - qvair) = (Eg' + Etr' + Er')
!
! This subroutine will solve Eg', Etr', Er', and leflx'
!
!-----------------------------------------------------------------------
!
! AUTHOR: J. Brotzge
! 3/06/02
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
!
! cdqa Array for surface drag coefficient for water vapor
! cdha Array for surface drag coefficient for heat
! veg Vegetation fraction
!
! qvair Specific humidity near the surface
! windsp Wind speed near the surface
! rhoa Air density near the surface
!
! OUTPUT:
!
! evaprp Evaporation from groud surface
! evaprtr Transpiration of the remaining part (1-delta) of leaves
! evaprr Direct evaporation from the fraction delta
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny
INTEGER :: soiltyp(nx,ny)
REAL :: cdqa (nx,ny)
REAL :: cdha (nx,ny)
REAL :: veg (nx,ny)
REAL :: wetcanp(nx,ny)
REAL :: qvair (nx,ny)
REAL :: qvsfc (nx,ny)
REAL :: qvsat (nx,ny)
REAL :: windsp (nx,ny)
REAL :: rhoa (nx,ny)
REAL :: rsttemp(nx,ny)
REAL :: evappot(nx,ny)
REAL :: evaprg (nx,ny)
REAL :: evaprtr(nx,ny)
REAL :: evaprr (nx,ny)
REAL :: lhflx (nx,ny)
REAL :: deltem (nx,ny)
REAL :: rrtem (nx,ny)
REAL :: temple (nx,ny)
!
!-----------------------------------------------------------------------
! Include file: phycst.inc
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!
!-----------------------------------------------------------------------
!
! Local variables:
!
!-----------------------------------------------------------------------
!
!
INTEGER :: i, j
REAL :: tempbc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j = 1, ny-1
DO i = 1, nx-1
!-----------------------------------------------------------------------
!
! Over water and ice, qvsfc should be saturated
!
!-----------------------------------------------------------------------
!
IF ( soiltyp(i,j) == 12 .OR. soiltyp(i,j) == 13 ) THEN
! evaprg(i,j) = qvsat(i,j) - qvair(i,j)
evaprg(i,j) = (qvsat(i,j) - qvair(i,j)) * cdqa(i,j) * &
rhoa(i,j) * windsp(i,j)
evaprtr(i,j) = 0.0
evaprr (i,j) = 0.0
ELSE ! over land
tempbc = (1.0 + deltem(i,j)/(rrtem(i,j) + 1.0)) / &
(1.0 + rsttemp(i,j)*cdha(i,j)+deltem(i,j) / &
(rrtem(i,j) + 1.0) )
evaprg (i,j) = temple(i,j)
evaprtr(i,j) = veg(i,j) * evappot(i,j) * tempbc * &
(1.0 - (2.0*wetcanp(i,j))**0.5)
evaprr (i,j) = rhoa(i,j)*cdqa(i,j)*windsp(i,j)*evaprr(i,j)
END IF
lhflx(i,j) = evaprg(i,j) + evaprtr(i,j) + evaprr(i,j)
lhflx(i,j) = lhflx(i,j) * lathv ! Latent heat flux
END DO
END DO
RETURN
END SUBROUTINE leflx