MODULE module_sf_lake 3
! The lake scheme was retrieved from the Community Land Model version 4.5
! (Oleson et al. 2013) with some modifications by Gu et al. (2013). It is a
! one-dimensional mass and energy balance scheme with 20-25 model layers,
! including up to 5 snow layers on the lake ice, 10 water layers, and 10 soil
! layers on the lake bottom. The lake scheme is used with actual lake points and
! lake depth derived from the WPS, and it also can be used with user defined
! lake points and lake depth in WRF (lake_min_elev and lakedepth_default).
! The lake scheme is independent of a land surface scheme and therefore
! can be used with any land surface scheme embedded in WRF. The lake scheme
! developments and evaluations were included in Subin et al. (2012) and Gu et al. (2013)
!
! Subin et al. 2012: Improved lake model for climate simulations, J. Adv. Model.
! Earth Syst., 4, M02001. DOI:10.1029/2011MS000072;
! Gu et al. 2013: Calibration and validation of lake surface temperature simulations
! with the coupled WRF-Lake model. Climatic Change, 1-13, 10.1007/s10584-013-0978-y.
USE module_wrf_error
USE module_model_constants, ONLY : rcp
implicit none
integer, parameter :: r8 = selected_real_kind(12)
integer, parameter :: nlevsoil = 10 ! number of soil layers
integer, parameter :: nlevlake = 10 ! number of lake layers
integer, parameter :: nlevsnow = 5 ! maximum number of snow layers
integer,parameter :: lbp = 1 ! pft-index bounds
integer,parameter :: ubp = 1
integer,parameter :: lbc = 1 ! column-index bounds
integer,parameter :: ubc = 1
integer,parameter :: num_shlakec = 1 ! number of columns in lake filter
integer,parameter :: filter_shlakec(1) = 1 ! lake filter (columns)
integer,parameter :: num_shlakep = 1 ! number of pfts in lake filter
integer,parameter :: filter_shlakep(1) = 1 ! lake filter (pfts)
integer,parameter :: pcolumn(1) = 1
integer,parameter :: pgridcell(1) = 1
integer,parameter :: cgridcell(1) = 1 ! gridcell index of column
integer,parameter :: clandunit(1) = 1 ! landunit index of column
integer,parameter :: begg = 1
integer,parameter :: endg = 1
integer,parameter :: begl = 1
integer,parameter :: endl = 1
integer,parameter :: begc = 1
integer,parameter :: endc = 1
integer,parameter :: begp = 1
integer,parameter :: endp = 1
integer,parameter :: column =1
logical,parameter :: lakpoi(1) = .true.
!Initialize physical constants:
real(r8), parameter :: vkc = 0.4_r8 !von Karman constant [-]
real(r8), parameter :: pie = 3.141592653589793_r8 ! pi
real(r8), parameter :: grav = 9.80616_r8 !gravity constant [m/s2]
real(r8), parameter :: sb = 5.67e-8_r8 !stefan-boltzmann constant [W/m2/K4]
real(r8), parameter :: tfrz = 273.16_r8 !freezing temperature [K]
real(r8), parameter :: denh2o = 1.000e3_r8 !density of liquid water [kg/m3]
real(r8), parameter :: denice = 0.917e3_r8 !density of ice [kg/m3]
real(r8), parameter :: cpice = 2.11727e3_r8 !Specific heat of ice [J/kg-K]
real(r8), parameter :: cpliq = 4.188e3_r8 !Specific heat of water [J/kg-K]
real(r8), parameter :: hfus = 3.337e5_r8 !Latent heat of fusion for ice [J/kg]
real(r8), parameter :: hvap = 2.501e6_r8 !Latent heat of evap for water [J/kg]
real(r8), parameter :: hsub = 2.501e6_r8+3.337e5_r8 !Latent heat of sublimation [J/kg]
real(r8), parameter :: rair = 287.0423_r8 !gas constant for dry air [J/kg/K]
real(r8), parameter :: cpair = 1.00464e3_r8 !specific heat of dry air [J/kg/K]
real(r8), parameter :: tcrit = 2.5 !critical temperature to determine rain or snow
real(r8), parameter :: tkwat = 0.6 !thermal conductivity of water [W/m/k]
real(r8), parameter :: tkice = 2.290 !thermal conductivity of ice [W/m/k]
real(r8), parameter :: tkairc = 0.023 !thermal conductivity of air [W/m/k]
real(r8), parameter :: bdsno = 250. !bulk density snow (kg/m**3)
real(r8), public, parameter :: spval = 1.e36 !special value for missing data (ocean)
real, parameter :: depth_c = 50. ! below the level t_lake3d will be 277.0 !mchen
! These are tunable constants
real(r8), parameter :: wimp = 0.05 !Water impremeable if porosity less than wimp
real(r8), parameter :: ssi = 0.033 !Irreducible water saturation of snow
real(r8), parameter :: cnfac = 0.5 !Crank Nicholson factor between 0 and 1
! Initialize water type constants
integer,parameter :: istsoil = 1 !soil "water" type
integer, private :: i ! loop index
real(r8) :: dtime ! land model time step (sec)
real(r8) :: zlak(1:nlevlake) !lake z (layers)
real(r8) :: dzlak(1:nlevlake) !lake dz (thickness)
real(r8) :: zsoi(1:nlevsoil) !soil z (layers)
real(r8) :: dzsoi(1:nlevsoil) !soil dz (thickness)
real(r8) :: zisoi(0:nlevsoil) !soil zi (interfaces)
real(r8) :: sand(19) ! percent sand
real(r8) :: clay(19) ! percent clay
data(sand(i), i=1,19)/92.,80.,66.,20.,5.,43.,60.,&
10.,32.,51., 6.,22.,39.7,0.,100.,54.,17.,100.,92./
data(clay(i), i=1,19)/ 3., 5.,10.,15.,5.,18.,27.,&
33.,33.,41.,47.,58.,14.7,0., 0., 8.5,54., 0., 3./
! real(r8) :: dtime ! land model time step (sec)
real(r8) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
real(r8) :: tksatu(1,nlevsoil) ! thermal conductivity, saturated soil [W/m-K]
real(r8) :: tkmg(1,nlevsoil) ! thermal conductivity, soil minerals [W/m-K]
real(r8) :: tkdry(1,nlevsoil) ! thermal conductivity, dry soil (W/m/Kelvin)
real(r8) :: csol(1,nlevsoil) ! heat capacity, soil solids (J/m**3/Kelvin)
CONTAINS
SUBROUTINE Lake( t_phy ,p8w ,dz8w ,qvcurr ,& !i 1,1
u_phy ,v_phy , glw ,emiss ,&
rainbl ,dtbl ,swdown ,albedo ,&
xlat_urb2d ,z_lake3d ,dz_lake3d ,lakedepth2d ,&
watsat3d ,csol3d ,tkmg3d ,tkdry3d ,&
tksatu3d ,ivgtyp ,ht ,xland ,&
iswater, xice, xice_threshold, lake_min_elev ,&
ids ,ide ,jds ,jde ,&
kds ,kde ,ims ,ime ,&
jms ,jme ,kms ,kme ,&
its ,ite ,jts ,jte ,&
kts ,kte ,&
h2osno2d ,snowdp2d ,snl2d ,z3d ,& !h
dz3d ,zi3d ,h2osoi_vol3d ,h2osoi_liq3d ,&
h2osoi_ice3d ,t_grnd2d ,t_soisno3d ,t_lake3d ,&
savedtke12d ,lake_icefrac3d ,&
#if (EM_CORE==1)
! lakemask ,lakeflag ,&
lakemask ,&
#endif
hfx ,lh ,grdflx ,tsk ,& !o
qfx ,t2 ,th2 ,q2 )
!==============================================================================
! This subroutine was first edited by Hongping Gu and Jiming Jin for coupling
! 07/20/2010
!==============================================================================
IMPLICIT NONE
!in:
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
INTEGER , INTENT (IN) :: iswater
REAL, INTENT(IN) :: xice_threshold
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: XICE
#if (EM_CORE==1)
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: LAKEMASK
! INTEGER, INTENT(IN):: LAKEFLAG
#endif
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: t_phy
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: p8w
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: dz8w
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: qvcurr
REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: U_PHY
REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: V_PHY
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: glw
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: emiss
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: rainbl
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: swdown
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(INOUT) :: albedo
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: XLAND
REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: XLAT_URB2D
INTEGER, DIMENSION( ims:ime, jms:jme ) ,INTENT(INOUT) :: IVGTYP
REAL, INTENT(IN) :: dtbl
REAL, DIMENSION( ims:ime,1:nlevlake,jms:jme ),INTENT(IN) :: z_lake3d
REAL, DIMENSION( ims:ime,1:nlevlake,jms:jme ),INTENT(IN) :: dz_lake3d
REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: watsat3d
REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: csol3d
REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: tkmg3d
REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: tkdry3d
REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: tksatu3d
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: lakedepth2d
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: ht
REAL ,INTENT(IN) :: lake_min_elev
!out:
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: HFX
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: LH
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: GRDFLX
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: TSK
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: QFX
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: T2
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: TH2
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: Q2
!in&out:
real, dimension(ims:ime,jms:jme ) ,intent(inout) :: savedtke12d
real, dimension(ims:ime,jms:jme ) ,intent(inout) :: snowdp2d, &
h2osno2d, &
snl2d, &
t_grnd2d
real, dimension( ims:ime,1:nlevlake, jms:jme ) ,INTENT(inout) :: t_lake3d, &
lake_icefrac3d
real, dimension( ims:ime,-nlevsnow+1:nlevsoil, jms:jme ) ,INTENT(inout) :: t_soisno3d, &
h2osoi_ice3d, &
h2osoi_liq3d, &
h2osoi_vol3d, &
z3d, &
dz3d
real, dimension( ims:ime,-nlevsnow+0:nlevsoil, jms:jme ) ,INTENT(inout) :: zi3d
!local variable:
REAL :: SFCTMP,PBOT,PSFC,ZLVL,Q2K,EMISSI,LWDN,PRCP,SOLDN,SOLNET
INTEGER :: C,i,j,k
!tempory varibles in:
real(r8) :: forc_t(1) ! atmospheric temperature (Kelvin)
real(r8) :: forc_pbot(1) ! atm bottom level pressure (Pa)
real(r8) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
real(r8) :: forc_hgt(1) ! atmospheric reference height (m)
real(r8) :: forc_hgt_q(1) ! observational height of humidity [m]
real(r8) :: forc_hgt_t(1) ! observational height of temperature [m]
real(r8) :: forc_hgt_u(1) ! observational height of wind [m]
real(r8) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
real(r8) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
real(r8) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
! real(r8) :: forc_rho(1) ! density (kg/m**3)
real(r8) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
real(r8) :: prec(1) ! snow or rain rate [mm/s]
real(r8) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
real(r8) :: lat(1) ! latitude (radians)
real(r8) :: z_lake(1,nlevlake) ! layer depth for lake (m)
real(r8) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
real(r8) :: lakedepth(1) ! column lake depth (m)
logical :: do_capsnow(1) ! true => do snow capping
!in&out
real(r8) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
real(r8) :: t_grnd(1) ! ground temperature (Kelvin)
real(r8) :: h2osno(1) ! snow water (mm H2O)
real(r8) :: snowdp(1) ! snow height (m)
real(r8) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
real(r8) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
real(r8) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
real(r8) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
integer :: snl(1) ! number of snow layers
real(r8) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
real(r8) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
real(r8) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
real(r8) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
real(r8) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
!out:
real(r8) :: eflx_gnet(1) !net heat flux into ground (W/m**2)
real(r8) :: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
real(r8) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
real(r8) :: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
real(r8) :: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
real(r8) :: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
real(r8) :: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
real(r8) :: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
real(r8) :: ram1(1) ! aerodynamical resistance (s/m)
! for calculation of decay of eddy diffusivity with depth
! Change the type variable to pass back to WRF.
real(r8) :: z0mg(1) ! roughness length over ground, momentum (m(
dtime = dtbl
DO J = jts,jte
DO I = its,ite
SFCTMP = t_phy(i,1,j)
PBOT = p8w(i,2,j)
PSFC = P8w(i,1,j)
ZLVL = 0.5 * dz8w(i,1,j)
Q2K = qvcurr(i,1,j)/(1.0 + qvcurr(i,1,j))
EMISSI = EMISS(I,J)
LWDN = GLW(I,J)*EMISSI
PRCP = RAINBL(i,j)/dtbl
SOLDN = SWDOWN(I,J) ! SOLDN is total incoming solar
SOLNET = SOLDN*(1.-ALBEDO(I,J)) ! use mid-day albedo to determine net downward solar
! (no solar zenith angle correction)
! IF (XLAND(I,J).GT.1.5) THEN
! if ( xice(i,j).gt.xice_threshold) then
! ivgtyp(i,j) = iswater
! xland(i,j) = 2.
! lake_icefrac3d(i,1,j) = xice(i,j)
! endif
#if (EM_CORE==1)
if (lakemask(i,j).eq.1) THEN
#else
if (ivgtyp(i,j)==iswater.and.ht(i,j)>= lake_min_elev ) THEN
#endif
do c = 1,column
forc_t(c) = SFCTMP ! [K]
forc_pbot(c) = PBOT
forc_psrf(c) = PSFC
forc_hgt(c) = ZLVL ! [m]
forc_hgt_q(c) = ZLVL ! [m]
forc_hgt_t(c) = ZLVL ! [m]
forc_hgt_u(c) = ZLVL ! [m]
forc_q(c) = Q2K ! [kg/kg]
forc_u(c) = U_PHY(I,1,J)
forc_v(c) = V_PHY(I,1,J)
! forc_rho(c) = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
forc_lwrad(c) = LWDN ! [W/m/m]
prec(c) = PRCP ! [mm/s]
sabg(c) = SOLNET
lat(c) = XLAT_URB2D(I,J)*pie/180 ! [radian]
do_capsnow(c) = .false.
lakedepth(c) = lakedepth2d(i,j)
savedtke1(c) = savedtke12d(i,j)
snowdp(c) = snowdp2d(i,j)
h2osno(c) = h2osno2d(i,j)
snl(c) = snl2d(i,j)
t_grnd(c) = t_grnd2d(i,j)
do k = 1,nlevlake
t_lake(c,k) = t_lake3d(i,k,j)
lake_icefrac(c,k) = lake_icefrac3d(i,k,j)
z_lake(c,k) = z_lake3d(i,k,j)
dz_lake(c,k) = dz_lake3d(i,k,j)
enddo
do k = -nlevsnow+1,nlevsoil
t_soisno(c,k) = t_soisno3d(i,k,j)
h2osoi_ice(c,k) = h2osoi_ice3d(i,k,j)
h2osoi_liq(c,k) = h2osoi_liq3d(i,k,j)
h2osoi_vol(c,k) = h2osoi_vol3d(i,k,j)
z(c,k) = z3d(i,k,j)
dz(c,k) = dz3d(i,k,j)
enddo
do k = -nlevsnow+0,nlevsoil
zi(c,k) = zi3d(i,k,j)
enddo
do k = 1,nlevsoil
watsat(c,k) = watsat3d(i,k,j)
csol(c,k) = csol3d(i,k,j)
tkmg(c,k) = tkmg3d(i,k,j)
tkdry(c,k) = tkdry3d(i,k,j)
tksatu(c,k) = tksatu3d(i,k,j)
enddo
enddo
CALL LakeMain(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !I
forc_hgt_t,forc_hgt_u,forc_q, forc_u, &
forc_v,forc_lwrad,prec, sabg,lat, &
z_lake,dz_lake,lakedepth,do_capsnow, &
h2osno,snowdp,snl,z,dz,zi, & !H
h2osoi_vol,h2osoi_liq,h2osoi_ice, &
t_grnd,t_soisno,t_lake, &
savedtke1,lake_icefrac, &
eflx_lwrad_net,eflx_gnet, & !O
eflx_sh_tot,eflx_lh_tot, &
t_ref2m,q_ref2m, &
taux,tauy,ram1,z0mg)
do c = 1,column
HFX(I,J) = eflx_sh_tot(c) ![W/m/m]
LH(I,J) = eflx_lh_tot(c) !W/m/m]
GRDFLX(I,J) = eflx_gnet(c) !W/m/m]
TSK(I,J) = t_grnd(c) ![K]
T2(I,J) = t_ref2m(c)
TH2(I,J) = T2(I,J)*(1.E5/PSFC)**RCP
Q2(I,J) = q_ref2m(c)
albedo(i,j) = ( 0.6 * lake_icefrac(c,1) ) + ( (1.0-lake_icefrac(c,1)) * 0.08)
if( tsk(i,j) >= tfrz ) then
qfx(i,j) = eflx_lh_tot(c)/hvap
else
qfx(i,j) = eflx_lh_tot(c)/hsub ! heat flux (W/m^2)=>mass flux(kg/(sm^2))
endif
enddo
! Renew Lake State Varialbes:(14)
do c = 1,column
savedtke12d(i,j) = savedtke1(c)
snowdp2d(i,j) = snowdp(c)
h2osno2d(i,j) = h2osno(c)
snl2d(i,j) = snl(c)
t_grnd2d(i,j) = t_grnd(c)
do k = 1,nlevlake
t_lake3d(i,k,j) = t_lake(c,k)
lake_icefrac3d(i,k,j) = lake_icefrac(c,k)
enddo
do k = -nlevsnow+1,nlevsoil
z3d(i,k,j) = z(c,k)
dz3d(i,k,j) = dz(c,k)
t_soisno3d(i,k,j) = t_soisno(c,k)
h2osoi_liq3d(i,k,j) = h2osoi_liq(c,k)
h2osoi_ice3d(i,k,j) = h2osoi_ice(c,k)
h2osoi_vol3d(i,k,j) = h2osoi_vol(c,k)
enddo
do k = -nlevsnow+0,nlevsoil
zi3d(i,k,j) = zi(c,k)
enddo
enddo
endif
! ENDIF ! if xland = 2
ENDDO
ENDDO
END SUBROUTINE Lake
SUBROUTINE LakeMain(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !I 1,3
forc_hgt_t,forc_hgt_u,forc_q, forc_u, &
forc_v,forc_lwrad,prec, sabg,lat, &
z_lake,dz_lake,lakedepth,do_capsnow, &
h2osno,snowdp,snl,z,dz,zi, & !H
h2osoi_vol,h2osoi_liq,h2osoi_ice, &
t_grnd,t_soisno,t_lake, &
savedtke1,lake_icefrac, &
eflx_lwrad_net,eflx_gnet, & !O
eflx_sh_tot,eflx_lh_tot, &
t_ref2m,q_ref2m, &
taux,tauy,ram1,z0mg)
implicit none
!in:
real(r8),intent(in) :: forc_t(1) ! atmospheric temperature (Kelvin)
real(r8),intent(in) :: forc_pbot(1) ! atm bottom level pressure (Pa)
real(r8),intent(in) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
real(r8),intent(in) :: forc_hgt(1) ! atmospheric reference height (m)
real(r8),intent(in) :: forc_hgt_q(1) ! observational height of humidity [m]
real(r8),intent(in) :: forc_hgt_t(1) ! observational height of temperature [m]
real(r8),intent(in) :: forc_hgt_u(1) ! observational height of wind [m]
real(r8),intent(in) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
real(r8),intent(in) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
real(r8),intent(in) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
! real(r8),intent(in) :: forc_rho(1) ! density (kg/m**3)
real(r8),intent(in) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
real(r8),intent(in) :: prec(1) ! snow or rain rate [mm/s]
real(r8),intent(in) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
real(r8),intent(in) :: lat(1) ! latitude (radians)
real(r8),intent(in) :: z_lake(1,nlevlake) ! layer depth for lake (m)
real(r8),intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
real(r8), intent(in) :: lakedepth(1) ! column lake depth (m)
!!!!!!!!!!!!!!!!tep(in),hydro(in)
! real(r8), intent(in) :: watsat(1,1:nlevsoil) ! volumetric soil water at saturation (porosity)
!!!!!!!!!!!!!!!!hydro
logical , intent(in) :: do_capsnow(1) ! true => do snow capping
!in&out
real(r8),intent(inout) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
real(r8),intent(inout) :: t_grnd(1) ! ground temperature (Kelvin)
real(r8),intent(inout) :: h2osno(1) ! snow water (mm H2O)
real(r8),intent(inout) :: snowdp(1) ! snow height (m)
real(r8),intent(inout) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
real(r8),intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
real(r8),intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
real(r8),intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
integer ,intent(inout) :: snl(1) ! number of snow layers
real(r8),intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
real(r8),intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
real(r8),intent(inout) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
real(r8),intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
real(r8),intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
!out:
real(r8),intent(out) :: eflx_gnet(1) !net heat flux into ground (W/m**2)
real(r8),intent(out) :: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
real(r8),intent(out) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
real(r8),intent(out) :: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
real(r8),intent(out) :: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
real(r8),intent(out) :: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
real(r8),intent(out) :: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
real(r8),intent(out) :: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
real(r8),intent(out) :: ram1(1) ! aerodynamical resistance (s/m)
! for calculation of decay of eddy diffusivity with depth
! Change the type variable to pass back to WRF.
real(r8),intent(out) :: z0mg(1) ! roughness length over ground, momentum (m(
!local output
real(r8) :: begwb(1) ! water mass begining of the time step
real(r8) :: t_veg(1) ! vegetation temperature (Kelvin)
real(r8) :: eflx_soil_grnd(1) ! soil heat flux (W/m**2) [+ = into soil]
real(r8) :: eflx_lh_grnd(1) ! ground evaporation heat flux (W/m**2) [+ to atm]
real(r8) :: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
real(r8) :: eflx_lwrad_out(1) ! emitted infrared (longwave) radiation (W/m**2)
real(r8) :: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
real(r8) :: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
real(r8) :: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
real(r8) :: forc_snow(1) ! snow rate [mm/s]
real(r8) :: forc_rain(1) ! rain rate [mm/s]
real(r8) :: ws(1) ! surface friction velocity (m/s)
real(r8) :: ks(1) ! coefficient passed to ShalLakeTemperature
real(r8) :: qflx_snomelt(1) !snow melt (mm H2O /s) tem(out),snowwater(in)
integer :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
real(r8) :: endwb(1) ! water mass end of the time step
real(r8) :: snowage(1) ! non dimensional snow age [-]
real(r8) :: snowice(1) ! average snow ice lens
real(r8) :: snowliq(1) ! average snow liquid water
real(r8) :: t_snow(1) ! vertically averaged snow temperature
real(r8) :: qflx_drain(1) ! sub-surface runoff (mm H2O /s)
real(r8) :: qflx_surf(1) ! surface runoff (mm H2O /s)
real(r8) :: qflx_infl(1) ! infiltration (mm H2O /s)
real(r8) :: qflx_qrgwl(1) ! qflx_surf at glaciers, wetlands, lakes
real(r8) :: qcharge(1) ! aquifer recharge rate (mm/s)
real(r8) :: qflx_snowcap(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
real(r8) :: qflx_snowcap_col(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
real(r8) :: qflx_snow_grnd_pft(1) ! snow on ground after interception (mm H2O/s) [+]
real(r8) :: qflx_snow_grnd_col(1) ! snow on ground after interception (mm H2O/s) [+]
real(r8) :: qflx_rain_grnd(1) ! rain on ground after interception (mm H2O/s) [+]
real(r8) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
real(r8) :: qflx_evap_tot_col(1) !pft quantity averaged to the column (assuming one pft)
real(r8) :: soilalpha(1) !factor that reduces ground saturated specific humidity (-)
real(r8) :: zwt(1) !water table depth
real(r8) :: fcov(1) !fractional area with water table at surface
real(r8) :: rootr_column(1,1:nlevsoil) !effective fraction of roots in each soil layer
real(r8) :: qflx_evap_grnd(1) ! ground surface evaporation rate (mm H2O/s) [+]
real(r8) :: qflx_sub_snow(1) ! sublimation rate from snow pack (mm H2O /s) [+]
real(r8) :: qflx_dew_snow(1) ! surface dew added to snow pack (mm H2O /s) [+]
real(r8) :: qflx_dew_grnd(1) ! ground surface dew formation (mm H2O /s) [+]
real(r8) :: qflx_rain_grnd_col(1) !rain on ground after interception (mm H2O/s) [+]
! lat = lat*pie/180 ! [radian]
if (prec(1)> 0.) then
if ( forc_t(1) > (tfrz + tcrit)) then
forc_rain(1) = prec(1)
forc_snow(1) = 0.
! flfall(1) = 1.
else
forc_rain(1) = 0.
forc_snow(1) = prec(1)
! if ( forc_t(1) <= tfrz) then
! flfall(1) = 0.
! else if ( forc_t(1) <= tfrz+2.) then
! flfall(1) = -54.632 + 0.2 * forc_t(1)
! else
! flfall(1) = 0.4
endif
else
forc_rain(1) = 0.
forc_snow(1) = 0.
! flfall(1) = 1.
endif
CALL ShalLakeFluxes(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !i
forc_hgt_t,forc_hgt_u,forc_q, &
forc_u,forc_v,forc_lwrad,forc_snow, &
forc_rain,t_grnd,h2osno,snowdp,sabg,lat, &
dz,dz_lake,t_soisno,t_lake,snl,h2osoi_liq, &
h2osoi_ice,savedtke1, &
qflx_prec_grnd,qflx_evap_soi,qflx_evap_tot, & !o
eflx_sh_grnd,eflx_lwrad_out,eflx_lwrad_net, &
eflx_soil_grnd,eflx_sh_tot,eflx_lh_tot, &
eflx_lh_grnd,t_veg,t_ref2m,q_ref2m,taux,tauy, &
ram1,ws,ks,eflx_gnet,z0mg)
CALL ShalLakeTemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & !i
z_lake,ws,ks,snl,eflx_gnet,lakedepth, &
lake_icefrac,snowdp, & !i&o
eflx_sh_grnd,eflx_sh_tot,eflx_soil_grnd, & !o
t_lake,t_soisno,h2osoi_liq, &
h2osoi_ice,savedtke1, &
frac_iceold,qflx_snomelt,imelt)
CALL ShalLakeHydrology(dz_lake,forc_rain,forc_snow, & !i
begwb,qflx_evap_tot,forc_t,do_capsnow, &
t_grnd,qflx_evap_soi, &
qflx_snomelt,imelt,frac_iceold, & !i add by guhp
z,dz,zi,snl,h2osno,snowdp,lake_icefrac,t_lake, & !i&o
endwb,snowage,snowice,snowliq,t_snow, & !o
t_soisno,h2osoi_ice,h2osoi_liq,h2osoi_vol, &
qflx_drain,qflx_surf,qflx_infl,qflx_qrgwl, &
qcharge,qflx_prec_grnd,qflx_snowcap, &
qflx_snowcap_col,qflx_snow_grnd_pft, &
qflx_snow_grnd_col,qflx_rain_grnd, &
qflx_evap_tot_col,soilalpha,zwt,fcov, &
rootr_column,qflx_evap_grnd,qflx_sub_snow, &
qflx_dew_snow,qflx_dew_grnd,qflx_rain_grnd_col)
!==================================================================================
! !DESCRIPTION:
! Calculation of Shallow Lake Hydrology. Full hydrology of snow layers is
! done. However, there is no infiltration, and the water budget is balanced with
END SUBROUTINE LakeMain
SUBROUTINE ShalLakeFluxes(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !i 1,10
forc_hgt_t,forc_hgt_u,forc_q, &
forc_u,forc_v,forc_lwrad,forc_snow, &
forc_rain,t_grnd,h2osno,snowdp,sabg,lat, &
dz,dz_lake,t_soisno,t_lake,snl,h2osoi_liq, &
h2osoi_ice,savedtke1, &
qflx_prec_grnd,qflx_evap_soi,qflx_evap_tot, & !o
eflx_sh_grnd,eflx_lwrad_out,eflx_lwrad_net, &
eflx_soil_grnd,eflx_sh_tot,eflx_lh_tot, &
eflx_lh_grnd,t_veg,t_ref2m,q_ref2m,taux,tauy, &
ram1,ws,ks,eflx_gnet,z0mg)
!==============================================================================
! DESCRIPTION:
! Calculates lake temperatures and surface fluxes for shallow lakes.
!
! Shallow lakes have variable depth, possible snow layers above, freezing & thawing of lake water,
! and soil layers with active temperature and gas diffusion below.
!
! WARNING: This subroutine assumes lake columns have one and only one pft.
!
! REVISION HISTORY:
! Created by Zack Subin, 2009
! Reedited by Hongping Gu, 2010
!==============================================================================
! implicit none
implicit none
!in:
real(r8),intent(in) :: forc_t(1) ! atmospheric temperature (Kelvin)
real(r8),intent(in) :: forc_pbot(1) ! atmospheric pressure (Pa)
real(r8),intent(in) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
real(r8),intent(in) :: forc_hgt(1) ! atmospheric reference height (m)
real(r8),intent(in) :: forc_hgt_q(1) ! observational height of humidity [m]
real(r8),intent(in) :: forc_hgt_t(1) ! observational height of temperature [m]
real(r8),intent(in) :: forc_hgt_u(1) ! observational height of wind [m]
real(r8),intent(in) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
real(r8),intent(in) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
real(r8),intent(in) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
real(r8),intent(in) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
! real(r8),intent(in) :: forc_rho(1) ! density (kg/m**3)
real(r8),intent(in) :: forc_snow(1) ! snow rate [mm/s]
real(r8),intent(in) :: forc_rain(1) ! rain rate [mm/s]
real(r8),intent(in) :: h2osno(1) ! snow water (mm H2O)
real(r8),intent(in) :: snowdp(1) ! snow height (m)
real(r8),intent(in) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
real(r8),intent(in) :: lat(1) ! latitude (radians)
real(r8),intent(in) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
real(r8),intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
real(r8),intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
real(r8),intent(in) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
integer ,intent(in) :: snl(1) ! number of snow layers
real(r8),intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
real(r8),intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
real(r8),intent(in) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
!inout:
real(r8),intent(inout) :: t_grnd(1) ! ground temperature (Kelvin)
!out:
real(r8),intent(out):: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
real(r8),intent(out):: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
real(r8),intent(out):: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
real(r8),intent(out):: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
real(r8),intent(out):: eflx_lwrad_out(1) ! emitted infrared (longwave) radiation (W/m**2)
real(r8),intent(out):: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
real(r8),intent(out):: eflx_soil_grnd(1) ! soil heat flux (W/m**2) [+ = into soil]
real(r8),intent(out):: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
real(r8),intent(out):: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
real(r8),intent(out):: eflx_lh_grnd(1) ! ground evaporation heat flux (W/m**2) [+ to atm]
real(r8),intent(out):: t_veg(1) ! vegetation temperature (Kelvin)
real(r8),intent(out):: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
real(r8),intent(out):: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
real(r8),intent(out):: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
real(r8),intent(out):: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
real(r8),intent(out):: ram1(1) ! aerodynamical resistance (s/m)
real(r8),intent(out):: ws(1) ! surface friction velocity (m/s)
real(r8),intent(out):: ks(1) ! coefficient passed to ShalLakeTemperature
! for calculation of decay of eddy diffusivity with depth
real(r8),intent(out):: eflx_gnet(1) !net heat flux into ground (W/m**2)
! Change the type variable to pass back to WRF.
real(r8),intent(out):: z0mg(1) ! roughness length over ground, momentum (m(
!OTHER LOCAL VARIABLES:
integer , parameter :: islak = 2 ! index of lake, 1 = deep lake, 2 = shallow lake
integer , parameter :: niters = 3 ! maximum number of iterations for surface temperature
real(r8), parameter :: beta1 = 1._r8 ! coefficient of convective velocity (in computing W_*) [-]
real(r8), parameter :: emg = 0.97_r8 ! ground emissivity (0.97 for snow)
real(r8), parameter :: zii = 1000._r8! convective boundary height [m]
real(r8), parameter :: tdmax = 277._r8 ! temperature of maximum water density
real(r8) :: forc_th(1) ! atmospheric potential temperature (Kelvin)
real(r8) :: forc_vp(1) !atmospheric vapor pressure (Pa)
real(r8) :: forc_rho(1) ! density (kg/m**3)
integer :: i,fc,fp,g,c,p ! do loop or array index
integer :: fncopy ! number of values in pft filter copy
integer :: fnold ! previous number of pft filter values
integer :: fpcopy(num_shlakep) ! pft filter copy for iteration loop
integer :: iter ! iteration index
integer :: nmozsgn(lbp:ubp) ! number of times moz changes sign
integer :: jtop(lbc:ubc) ! top level for each column (no longer all 1)
! real(r8) :: dtime ! land model time step (sec)
real(r8) :: ax ! used in iteration loop for calculating t_grnd (numerator of NR solution)
real(r8) :: bx ! used in iteration loop for calculating t_grnd (denomin. of NR solution)
real(r8) :: degdT ! d(eg)/dT
real(r8) :: dqh(lbp:ubp) ! diff of humidity between ref. height and surface
real(r8) :: dth(lbp:ubp) ! diff of virtual temp. between ref. height and surface
real(r8) :: dthv ! diff of vir. poten. temp. between ref. height and surface
real(r8) :: dzsur(lbc:ubc) ! 1/2 the top layer thickness (m)
real(r8) :: eg ! water vapor pressure at temperature T [pa]
real(r8) :: htvp(lbc:ubc) ! latent heat of vapor of water (or sublimation) [j/kg]
real(r8) :: obu(lbp:ubp) ! monin-obukhov length (m)
real(r8) :: obuold(lbp:ubp) ! monin-obukhov length of previous iteration
real(r8) :: qsatg(lbc:ubc) ! saturated humidity [kg/kg]
real(r8) :: qsatgdT(lbc:ubc) ! d(qsatg)/dT
real(r8) :: qstar ! moisture scaling parameter
real(r8) :: ram(lbp:ubp) ! aerodynamical resistance [s/m]
real(r8) :: rah(lbp:ubp) ! thermal resistance [s/m]
real(r8) :: raw(lbp:ubp) ! moisture resistance [s/m]
real(r8) :: stftg3(lbp:ubp) ! derivative of fluxes w.r.t ground temperature
real(r8) :: temp1(lbp:ubp) ! relation for potential temperature profile
real(r8) :: temp12m(lbp:ubp) ! relation for potential temperature profile applied at 2-m
real(r8) :: temp2(lbp:ubp) ! relation for specific humidity profile
real(r8) :: temp22m(lbp:ubp) ! relation for specific humidity profile applied at 2-m
real(r8) :: tgbef(lbc:ubc) ! initial ground temperature
real(r8) :: thm(lbc:ubc) ! intermediate variable (forc_t+0.0098*forc_hgt_t)
real(r8) :: thv(lbc:ubc) ! virtual potential temperature (kelvin)
real(r8) :: thvstar ! virtual potential temperature scaling parameter
real(r8) :: tksur ! thermal conductivity of snow/soil (w/m/kelvin)
real(r8) :: tsur ! top layer temperature
real(r8) :: tstar ! temperature scaling parameter
real(r8) :: um(lbp:ubp) ! wind speed including the stablity effect [m/s]
real(r8) :: ur(lbp:ubp) ! wind speed at reference height [m/s]
real(r8) :: ustar(lbp:ubp) ! friction velocity [m/s]
real(r8) :: wc ! convective velocity [m/s]
real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory
real(r8) :: zldis(lbp:ubp) ! reference height "minus" zero displacement height [m]
real(r8) :: displa(lbp:ubp) ! displacement (always zero) [m]
! real(r8) :: z0mg(lbp:ubp) ! roughness length over ground, momentum [m]
real(r8) :: z0hg(lbp:ubp) ! roughness length over ground, sensible heat [m]
real(r8) :: z0qg(lbp:ubp) ! roughness length over ground, latent heat [m]
real(r8) :: beta(2) ! fraction solar rad absorbed at surface: depends on lake type
real(r8) :: u2m ! 2 m wind speed (m/s)
real(r8) :: u10(1) ! 10-m wind (m/s) (for dust model)
real(r8) :: fv(1) ! friction velocity (m/s) (for dust model)
real(r8) :: fm(lbp:ubp) ! needed for BGC only to diagnose 10m wind speed
real(r8) :: bw ! partial density of water (ice + liquid)
real(r8) :: t_grnd_temp ! Used in surface flux correction over frozen ground
real(r8) :: betaprime(lbc:ubc) ! Effective beta: 1 for snow layers, beta(islak) otherwise
character*256 :: message
! This assumes all radiation is absorbed in the top snow layer and will need
! to be changed for CLM 4.
!
! Constants for lake temperature model
!
data beta/0.4_r8, 0.4_r8/ ! (deep lake, shallow lake)
! This is the energy absorbed at the lake surface if no snow.
! data za /0.6_r8, 0.5_r8/
! data eta /0.1_r8, 0.5_r8/
!-----------------------------------------------------------------------
! dtime = get_step_size()
! Begin calculations
!dir$ concurrent
!cdir nodep
forc_th(1) = forc_t(1) * (forc_psrf(1)/ forc_pbot(1))**(rair/cpair)
forc_vp(1) = forc_q(1) * forc_pbot(1)/ (0.622 + 0.378 * forc_q(1))
forc_rho(1) = (forc_pbot(1) - 0.378 * forc_vp(1)) / (rair * forc_t(1))
do fc = 1, num_shlakec
c = filter_shlakec(fc)
g = cgridcell(c)
! Surface temperature and fluxes
! Find top layer
if (snl(c) > 0 .or. snl(c) < -5) then
WRITE(message,*) 'snl is not defined in ShalLakeFluxesMod'
CALL wrf_message(message)
CALL wrf_error_fatal("snl: out of range value")
end if
! if (snl(c) /= 0) then
! write(6,*)'snl is not equal to zero in ShalLakeFluxesMod'
! call endrun()
! end if
jtop(c) = snl(c) + 1
if (snl(c) < 0) then
betaprime(c) = 1._r8 !Assume all solar rad. absorbed at the surface of the top snow layer.
dzsur(c) = dz(c,jtop(c))/2._r8
else
betaprime(c) = beta(islak)
dzsur(c) = dz_lake(c,1)/2._r8
end if
! Originally this was 1*dz, but shouldn't it be 1/2?
! Saturated vapor pressure, specific humidity and their derivatives
! at lake surface
call QSat(t_grnd(c), forc_pbot(g), eg, degdT, qsatg(c), qsatgdT(c))
! Potential, virtual potential temperature, and wind speed at the
! reference height
thm(c) = forc_t(g) + 0.0098_r8*forc_hgt_t(g) ! intermediate variable
thv(c) = forc_th(g)*(1._r8+0.61_r8*forc_q(g)) ! virtual potential T
end do
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
g = pgridcell(p)
nmozsgn(p) = 0
obuold(p) = 0._r8
displa(p) = 0._r8
! Roughness lengths
! changed by Hongping Gu
! if (t_grnd(c) >= tfrz) then ! for unfrozen lake
! z0mg(p) = 0.01_r8
! else ! for frozen lake
! ! Is this okay even if it is snow covered? What is the roughness over
! non-veg. snow?
! z0mg(p) = 0.04_r8
! end if
if (t_grnd(c) >= tfrz) then ! for unfrozen lake
z0mg(p) = 0.001_r8 !original 0.01
else if(snl(c) == 0 ) then ! for frozen lake
! Is this okay even if it is snow covered? What is the roughness over
! non-veg. snow?
z0mg(p) = 0.005_r8 !original 0.04, now for frozen lake without snow
else ! for frozen lake with snow
z0mg(p) = 0.0024_r8
end if
z0hg(p) = z0mg(p)
z0qg(p) = z0mg(p)
! Latent heat
#if (defined PERGRO)
htvp(c) = hvap
#else
if (t_grnd(c) > tfrz) then
htvp(c) = hvap
else
htvp(c) = hsub
end if
#endif
! Zack Subin, 3/26/09: Shouldn't this be the ground temperature rather than the air temperature above?
! I'll change it for now.
! Initialize stability variables
ur(p) = max(1.0_r8,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g)))
dth(p) = thm(c)-t_grnd(c)
dqh(p) = forc_q(g)-qsatg(c)
dthv = dth(p)*(1._r8+0.61_r8*forc_q(g))+0.61_r8*forc_th(g)*dqh(p)
zldis(p) = forc_hgt_u(g) - 0._r8
! Initialize Monin-Obukhov length and wind speed
call MoninObukIni(ur(p), thv(c), dthv, zldis(p), z0mg(p), um(p), obu(p))
end do
iter = 1
fncopy = num_shlakep
fpcopy(1:num_shlakep) = filter_shlakep(1:num_shlakep)
! Begin stability iteration
ITERATION : do while (iter <= niters .and. fncopy > 0)
! Determine friction velocity, and potential temperature and humidity
! profiles of the surface boundary layer
call FrictionVelocity(pgridcell,forc_hgt,forc_hgt_u, & !i
forc_hgt_t,forc_hgt_q, & !i
lbp, ubp, fncopy, fpcopy, & !i
displa, z0mg, z0hg, z0qg, & !i
obu, iter, ur, um, & !i
ustar,temp1, temp2, temp12m, temp22m, & !o
u10,fv, & !o
fm) !i&o
!dir$ concurrent
!cdir nodep
do fp = 1, fncopy
p = fpcopy(fp)
c = pcolumn(p)
g = pgridcell(p)
tgbef(c) = t_grnd(c)
if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
tksur = savedtke1(c)
! Set this to the eddy conductivity from the last
! timestep, as the molecular conductivity will be orders of magnitude too small.
! Will have to deal with first timestep.
tsur = t_lake(c,1)
else if (snl(c) == 0) then !frozen but no snow layers
tksur = tkice
tsur = t_lake(c,1)
else
!Need to calculate thermal conductivity of the top snow layer
bw = (h2osoi_ice(c,jtop(c))+h2osoi_liq(c,jtop(c)))/dz(c,jtop(c))
tksur = tkairc + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkairc)
tsur = t_soisno(c,jtop(c))
end if
! Determine aerodynamic resistances
ram(p) = 1._r8/(ustar(p)*ustar(p)/um(p))
rah(p) = 1._r8/(temp1(p)*ustar(p))
raw(p) = 1._r8/(temp2(p)*ustar(p))
ram1(p) = ram(p) !pass value to global variable
! Get derivative of fluxes with respect to ground temperature
stftg3(p) = emg*sb*tgbef(c)*tgbef(c)*tgbef(c)
! Changed surface temperature from t_lake(c,1) to tsur.
! Also adjusted so that if there are snow layers present, all radiation is absorbed in the top layer.
ax = betaprime(c)*sabg(p) + emg*forc_lwrad(g) + 3._r8*stftg3(p)*tgbef(c) &
+ forc_rho(g)*cpair/rah(p)*thm(c) &
- htvp(c)*forc_rho(g)/raw(p)*(qsatg(c)-qsatgdT(c)*tgbef(c) - forc_q(g)) &
+ tksur*tsur/dzsur(c)
!Changed sabg(p) and to betaprime(c)*sabg(p).
bx = 4._r8*stftg3(p) + forc_rho(g)*cpair/rah(p) &
+ htvp(c)*forc_rho(g)/raw(p)*qsatgdT(c) + tksur/dzsur(c)
t_grnd(c) = ax/bx
! Update htvp
#ifndef PERGRO
if (t_grnd(c) > tfrz) then
htvp(c) = hvap
else
htvp(c) = hsub
end if
#endif
! Surface fluxes of momentum, sensible and latent heat
! using ground temperatures from previous time step
eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdT(c)*(t_grnd(c)-tgbef(c))-forc_q(g))/raw(p)
! Re-calculate saturated vapor pressure, specific humidity and their
! derivatives at lake surface
call QSat(t_grnd(c), forc_pbot(g), eg, degdT, qsatg(c), qsatgdT(c))
dth(p)=thm(c)-t_grnd(c)
dqh(p)=forc_q(g)-qsatg(c)
tstar = temp1(p)*dth(p)
qstar = temp2(p)*dqh(p)
thvstar=tstar*(1._r8+0.61_r8*forc_q(g)) + 0.61_r8*forc_th(g)*qstar
zeta=zldis(p)*vkc * grav*thvstar/(ustar(p)**2*thv(c))
if (zeta >= 0._r8) then !stable
zeta = min(2._r8,max(zeta,0.01_r8))
um(p) = max(ur(p),0.1_r8)
else !unstable
zeta = max(-100._r8,min(zeta,-0.01_r8))
wc = beta1*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
um(p) = sqrt(ur(p)*ur(p)+wc*wc)
end if
obu(p) = zldis(p)/zeta
if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1
obuold(p) = obu(p)
end do ! end of filtered pft loop
iter = iter + 1
if (iter <= niters ) then
! Rebuild copy of pft filter for next pass through the ITERATION loop
fnold = fncopy
fncopy = 0
do fp = 1, fnold
p = fpcopy(fp)
if (nmozsgn(p) < 3) then
fncopy = fncopy + 1
fpcopy(fncopy) = p
end if
end do ! end of filtered pft loop
end if
end do ITERATION ! end of stability iteration
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
g = pgridcell(p)
! If there is snow on the ground and t_grnd > tfrz: reset t_grnd = tfrz.
! Re-evaluate ground fluxes.
! h2osno > 0.5 prevents spurious fluxes.
! note that qsatg and qsatgdT should be f(tgbef) (PET: not sure what this
! comment means)
! Zack Subin, 3/27/09: Since they are now a function of whatever t_grnd was before cooling
! to freezing temperature, then this value should be used in the derivative correction term.
! Should this happen if the lake temperature is below freezing, too? I'll assume that for now.
! Also, allow convection if ground temp is colder than lake but warmer than 4C, or warmer than
! lake which is warmer than freezing but less than 4C.
!#ifndef SHLAKETEST
if ( (h2osno(c) > 0.5_r8 .or. t_lake(c,1) <= tfrz) .and. t_grnd(c) > tfrz) then
!#else
! if ( t_lake(c,1) <= tfrz .and. t_grnd(c) > tfrz) then
!#endif
t_grnd_temp = t_grnd(c)
t_grnd(c) = tfrz
eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdT(c)*(t_grnd(c)-t_grnd_temp) - forc_q(g))/raw(p)
else if ( (t_lake(c,1) > t_grnd(c) .and. t_grnd(c) > tdmax) .or. &
(t_lake(c,1) < t_grnd(c) .and. t_lake(c,1) > tfrz .and. t_grnd(c) < tdmax) ) then
! Convective mixing will occur at surface
t_grnd_temp = t_grnd(c)
t_grnd(c) = t_lake(c,1)
eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdT(c)*(t_grnd(c)-t_grnd_temp) - forc_q(g))/raw(p)
end if
! Update htvp
#ifndef PERGRO
if (t_grnd(c) > tfrz) then
htvp(c) = hvap
else
htvp(c) = hsub
end if
#endif
! Net longwave from ground to atmosphere
! eflx_lwrad_out(p) = (1._r8-emg)*forc_lwrad(g) + stftg3(p)*(-3._r8*tgbef(c)+4._r8*t_grnd(c))
! What is tgbef doing in this equation? Can't it be exact now? --Zack Subin, 4/14/09
eflx_lwrad_out(p) = (1._r8-emg)*forc_lwrad(g) + emg*sb*t_grnd(c)**4
! Ground heat flux
eflx_soil_grnd(p) = sabg(p) + forc_lwrad(g) - eflx_lwrad_out(p) - &
eflx_sh_grnd(p) - htvp(c)*qflx_evap_soi(p)
!Why is this sabg(p) and not beta*sabg(p)??
!I've kept this as the incorrect sabg so that the energy balance check will be correct.
!This is the effective energy flux into the ground including the lake solar absorption
!below the surface. The variable eflx_gnet will be used to pass the actual heat flux
!from the ground interface into the lake.
taux(p) = -forc_rho(g)*forc_u(g)/ram(p)
tauy(p) = -forc_rho(g)*forc_v(g)/ram(p)
eflx_sh_tot(p) = eflx_sh_grnd(p)
qflx_evap_tot(p) = qflx_evap_soi(p)
eflx_lh_tot(p) = htvp(c)*qflx_evap_soi(p)
eflx_lh_grnd(p) = htvp(c)*qflx_evap_soi(p)
#if (defined LAKEDEBUG)
write(message,*) 'c, sensible heat = ', c, eflx_sh_tot(p), 'latent heat = ', eflx_lh_tot(p) &
, 'ground temp = ', t_grnd(c), 'h2osno = ', h2osno(c)
CALL wrf_message(message)
if (abs(eflx_sh_tot(p)) > 1500 .or. abs(eflx_lh_tot(p)) > 1500) then
write(message,*)'WARNING: SH, LH = ', eflx_sh_tot(p), eflx_lh_tot(p)
CALL wrf_message(message)
end if
if (abs(eflx_sh_tot(p)) > 10000 .or. abs(eflx_lh_tot(p)) > 10000 &
.or. abs(t_grnd(c)-288)>200 ) CALL wrf_error_fatal ( 't_grnd is out of range' )
#endif
! 2 m height air temperature
t_ref2m(p) = thm(c) + temp1(p)*dth(p)*(1._r8/temp12m(p) - 1._r8/temp1(p))
! 2 m height specific humidity
q_ref2m(p) = forc_q(g) + temp2(p)*dqh(p)*(1._r8/temp22m(p) - 1._r8/temp2(p))
! Energy residual used for melting snow
! Effectively moved to ShalLakeTemp
! Prepare for lake layer temperature calculations below
! fin(c) = betaprime * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
! eflx_sh_tot(p) + eflx_lh_tot(p))
! NOW this is just the net ground heat flux calculated below.
eflx_gnet(p) = betaprime(c) * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
eflx_sh_tot(p) + eflx_lh_tot(p))
! This is the actual heat flux from the ground interface into the lake, not including
! the light that penetrates the surface.
! u2m = max(1.0_r8,ustar(p)/vkc*log(2._r8/z0mg(p)))
! u2 often goes below 1 m/s; it seems like the only reason for this minimum is to
! keep it from being zero in the ks equation below; 0.1 m/s is a better limit for
! stable conditions --ZS
u2m = max(0.1_r8,ustar(p)/vkc*log(2._r8/z0mg(p)))
ws(c) = 1.2e-03_r8 * u2m
ks(c) = 6.6_r8*sqrt(abs(sin(lat(g))))*(u2m**(-1.84_r8))
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! End of surface flux relevant code in original BiogeophysicsLakeMod until history loop.
! The following are needed for global average on history tape.
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
g = pgridcell(p)
! t_veg(p) = forc_t(g)
!This is an odd choice, since elsewhere t_veg = t_grnd for bare ground.
!Zack Subin, 4/09
t_veg(p) = t_grnd(c)
eflx_lwrad_net(p) = eflx_lwrad_out(p) - forc_lwrad(g)
qflx_prec_grnd(p) = forc_rain(g) + forc_snow(g)
end do
END SUBROUTINE ShalLakeFluxes
SUBROUTINE ShalLakeTemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & !i 1,13
z_lake,ws,ks,snl,eflx_gnet,lakedepth, &
lake_icefrac,snowdp, & !i&o
eflx_sh_grnd,eflx_sh_tot,eflx_soil_grnd, & !o
t_lake,t_soisno,h2osoi_liq, &
h2osoi_ice,savedtke1, &
frac_iceold,qflx_snomelt,imelt)
!=======================================================================================================
! !DESCRIPTION:
! Calculates temperatures in the 20-25 layer column of (possible) snow,
! lake water, and soil beneath lake.
! Snow and soil temperatures are determined as in SoilTemperature, except
! for appropriate boundary conditions at the top of the snow (the flux is fixed
! to be the ground heat flux calculated in ShalLakeFluxes), the bottom of the snow
! (adjacent to top lake layer), and the top of the soil (adjacent to the bottom
! lake layer). Also, the soil is assumed to be always fully saturated (ShalLakeHydrology
! will have to insure this). The whole column is solved simultaneously as one tridiagonal matrix.
! Lake temperatures are determined from the Hostetler model as before, except now:
! i) Lake water layers can freeze by any fraction and release latent heat; thermal
! and mechanical properties are adjusted for ice fraction.
! ii) Convective mixing (though not eddy diffusion) still occurs for frozen lakes.
! iii) No sunlight is absorbed in the lake if there are snow layers.
! iv) Light is allowed to reach the top soil layer (where it is assumed to be completely absorbed).
! v) Lakes have variable depth, set ultimately in surface data set but now in initShalLakeMod.
!
! Eddy + molecular diffusion:
! d ts d d ts 1 ds
! ---- = -- [(km + ke) ----] + -- --
! dt dz dz cw dz
!
! where: ts = temperature (kelvin)
! t = time (s)
! z = depth (m)
! km = molecular diffusion coefficient (m**2/s)
! ke = eddy diffusion coefficient (m**2/s)
! cw = heat capacity (j/m**3/kelvin)
! s = heat source term (w/m**2)
!
! Shallow lakes are allowed to have variable depth, set in _____.
!
! For shallow lakes: ke > 0 if unfrozen,
! and convective mixing occurs WHETHER OR NOT frozen. (See e.g. Martynov...)
!
! Use the Crank-Nicholson method to set up tridiagonal system of equations to
! solve for ts at time n+1, where the temperature equation for layer i is
! r_i = a_i [ts_i-1] n+1 + b_i [ts_i] n+1 + c_i [ts_i+1] n+1
!
! The solution conserves energy as:
!
! [For lake layers]
! cw*([ts( 1)] n+1 - [ts( 1)] n)*dz( 1)/dt + ... +
! cw*([ts(nlevlake)] n+1 - [ts(nlevlake)] n)*dz(nlevlake)/dt = fin
! But now there is phase change, so cv is not constant and there is
! latent heat.
!
! where:
! [ts] n = old temperature (kelvin)
! [ts] n+1 = new temperature (kelvin)
! fin = heat flux into lake (w/m**2)
! = betaprime*sabg + forc_lwrad - eflx_lwrad_out - eflx_sh_tot - eflx_lh_tot
! (This is now the same as the ground heat flux.)
! + phi(1) + ... + phi(nlevlake) + phi(top soil level)
! betaprime = beta(islak) for no snow layers, and 1 for snow layers.
! This assumes all radiation is absorbed in the top snow layer and will need
! to be changed for CLM 4.
!
! WARNING: This subroutine assumes lake columns have one and only one pft.
!
! Outline:
! 1!) Initialization
! 2!) Lake density
! 3!) Diffusivity
! 4!) Heat source term from solar radiation penetrating lake
! 5!) Set thermal props and find initial energy content
! 6!) Set up vectors for tridiagonal matrix solution
! 7!) Solve tridiagonal and back-substitute
! 8!) (Optional) Do first energy check using temperature change at constant heat capacity.
! 9!) Phase change
! 9.5!) (Optional) Do second energy check using temperature change and latent heat, considering changed heat capacity.
! Also do soil water balance check.
!10!) Convective mixing
!11!) Do final energy check to detect small numerical errors (especially from convection)
! and dump small imbalance into sensible heat, or pass large errors to BalanceCheckMod for abort.
!
! REVISION HISTORY:
! Created by Zack Subin, 2009.
! Reedited by Hongping Gu, 2010.
!=========================================================================================================
! use TridiagonalMod , only : Tridiagonal
implicit none
!in:
real(r8), intent(in) :: t_grnd(1) ! ground temperature (Kelvin)
real(r8), intent(inout) :: h2osno(1) ! snow water (mm H2O)
real(r8), intent(in) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
real(r8), intent(in) :: dz(1,-nlevsnow + 1:nlevsoil) ! layer thickness for snow & soil (m)
real(r8), intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
real(r8), intent(in) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
real(r8), intent(in) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
! the other z and dz variables
real(r8), intent(in) :: z_lake(1,nlevlake) ! layer depth for lake (m)
real(r8), intent(in) :: ws(1) ! surface friction velocity (m/s)
real(r8), intent(in) :: ks(1) ! coefficient passed to ShalLakeTemperature
! for calculation of decay of eddy diffusivity with depth
integer , intent(in) :: snl(1) ! negative of number of snow layers
real(r8), intent(inout) :: eflx_gnet(1) ! net heat flux into ground (W/m**2) at the surface interface
real(r8), intent(in) :: lakedepth(1) ! column lake depth (m)
! real(r8), intent(in) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
real(r8), intent(inout) :: snowdp(1) !snow height (m)
!out:
real(r8), intent(out) :: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
real(r8), intent(out) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
real(r8), intent(out) :: eflx_soil_grnd(1) ! heat flux into snow / lake (W/m**2) [+ = into soil]
! Here this includes the whole lake radiation absorbed.
#if (defined SHLAKETEST)
real(r8), intent(out) :: qmelt(1) ! snow melt [mm/s] [temporary]
#endif
real(r8), intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2) [for snow & soil layers]
real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2) [for snow & soil layers]
real(r8), intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
real(r8), intent(out) :: savedtke1(1) ! top level thermal conductivity (W/mK)
real(r8), intent(out) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
real(r8), intent(out) :: qflx_snomelt(1) !snow melt (mm H2O /s)
integer, intent(out) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
! OTHER LOCAL VARIABLES:
integer , parameter :: islak = 2 ! index of lake, 1 = deep lake, 2 = shallow lake
real(r8), parameter :: p0 = 1._r8 ! neutral value of turbulent prandtl number
integer :: i,j,fc,fp,g,c,p ! do loop or array index
! real(r8) :: dtime ! land model time step (sec)
real(r8) :: beta(2) ! fraction solar rad absorbed at surface: depends on lake type
real(r8) :: za(2) ! base of surface absorption layer (m): depends on lake type
real(r8) :: eta(2) ! light extinction coefficient (/m): depends on lake type
real(r8) :: cwat ! specific heat capacity of water (j/m**3/kelvin)
real(r8) :: cice_eff ! effective heat capacity of ice (using density of
! water because layer depth is not adjusted when freezing
real(r8) :: cfus ! effective heat of fusion per unit volume
! using water density as above
real(r8) :: km ! molecular diffusion coefficient (m**2/s)
real(r8) :: tkice_eff ! effective conductivity since layer depth is constant
real(r8) :: a(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "a" vector for tridiagonal matrix
real(r8) :: b(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "b" vector for tridiagonal matrix
real(r8) :: c1(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "c" vector for tridiagonal matrix
real(r8) :: r(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "r" vector for tridiagonal solution
real(r8) :: rhow(lbc:ubc,nlevlake) ! density of water (kg/m**3)
real(r8) :: phi(lbc:ubc,nlevlake) ! solar radiation absorbed by layer (w/m**2)
real(r8) :: kme(lbc:ubc,nlevlake) ! molecular + eddy diffusion coefficient (m**2/s)
real(r8) :: rsfin ! relative flux of solar radiation into layer
real(r8) :: rsfout ! relative flux of solar radiation out of layer
real(r8) :: phi_soil(lbc:ubc) ! solar radiation into top soil layer (W/m**2)
real(r8) :: ri ! richardson number
real(r8) :: fin(lbc:ubc) ! net heat flux into lake at ground interface (w/m**2)
real(r8) :: ocvts(lbc:ubc) ! (cwat*(t_lake[n ])*dz
real(r8) :: ncvts(lbc:ubc) ! (cwat*(t_lake[n+1])*dz
real(r8) :: ke ! eddy diffusion coefficient (m**2/s)
real(r8) :: zin ! depth at top of layer (m)
real(r8) :: zout ! depth at bottom of layer (m)
real(r8) :: drhodz ! d [rhow] /dz (kg/m**4)
real(r8) :: n2 ! brunt-vaisala frequency (/s**2)
real(r8) :: num ! used in calculating ri
real(r8) :: den ! used in calculating ri
real(r8) :: tav_froz(lbc:ubc) ! used in aver temp for convectively mixed layers (C)
real(r8) :: tav_unfr(lbc:ubc) ! "
real(r8) :: nav(lbc:ubc) ! used in aver temp for convectively mixed layers
real(r8) :: phidum ! temporary value of phi
real(r8) :: iceav(lbc:ubc) ! used in calc aver ice for convectively mixed layers
real(r8) :: qav(lbc:ubc) ! used in calc aver heat content for conv. mixed layers
integer :: jtop(lbc:ubc) ! top level for each column (no longer all 1)
real(r8) :: cv (lbc:ubc,-nlevsnow+1:nlevsoil) !heat capacity of soil/snow [J/(m2 K)]
real(r8) :: tk (lbc:ubc,-nlevsnow+1:nlevsoil) !thermal conductivity of soil/snow [W/(m K)]
!(at interface below, except for j=0)
real(r8) :: cv_lake (lbc:ubc,1:nlevlake) !heat capacity [J/(m2 K)]
real(r8) :: tk_lake (lbc:ubc,1:nlevlake) !thermal conductivity at layer node [W/(m K)]
real(r8) :: cvx (lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !heat capacity for whole column [J/(m2 K)]
real(r8) :: tkix(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !thermal conductivity at layer interfaces
!for whole column [W/(m K)]
real(r8) :: tx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! temperature of whole column [K]
real(r8) :: tktopsoillay(lbc:ubc) ! thermal conductivity [W/(m K)]
real(r8) :: fnx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !heat diffusion through the layer interface below [W/m2]
real(r8) :: phix(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !solar source term for whole column [W/m**2]
real(r8) :: zx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !interface depth (+ below surface) for whole column [m]
real(r8) :: dzm !used in computing tridiagonal matrix [m]
real(r8) :: dzp !used in computing tridiagonal matrix [m]
integer :: jprime ! j - nlevlake
real(r8) :: factx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !coefficient used in computing tridiagonal matrix
real(r8) :: t_lake_bef(lbc:ubc,1:nlevlake) !beginning lake temp for energy conservation check [K]
real(r8) :: t_soisno_bef(lbc:ubc,-nlevsnow+1:nlevsoil) !beginning soil temp for E cons. check [K]
real(r8) :: lhabs(lbc:ubc) ! total per-column latent heat abs. from phase change (J/m^2)
real(r8) :: esum1(lbc:ubc) ! temp for checking energy (J/m^2)
real(r8) :: esum2(lbc:ubc) ! ""
real(r8) :: zsum(lbc:ubc) ! temp for putting ice at the top during convection (m)
real(r8) :: wsum(lbc:ubc) ! temp for checking water (kg/m^2)
real(r8) :: wsum_end(lbc:ubc) ! temp for checking water (kg/m^2)
real(r8) :: errsoi(1) ! soil/lake energy conservation error (W/m**2)
real(r8) :: eflx_snomelt(1) !snow melt heat flux (W/m**2)
CHARACTER*256 :: message
!
! Constants for lake temperature model
!
data beta/0.4_r8, 0.4_r8/ ! (deep lake, shallow lake)
data za /0.6_r8, 0.6_r8/
! For now, keep beta and za for shallow lake the same as deep lake, until better data is found.
! It looks like eta is key and that larger values give better results for shallow lakes. Use
! empirical expression from Hakanson (below). This is still a very unconstrained parameter
! that deserves more attention.
! Some radiation will be allowed to reach the soil.
!-----------------------------------------------------------------------
! 1!) Initialization
! Determine step size
! dtime = get_step_size()
! Initialize constants
cwat = cpliq*denh2o ! water heat capacity per unit volume
cice_eff = cpice*denh2o !use water density because layer depth is not adjusted
!for freezing
cfus = hfus*denh2o ! latent heat per unit volume
tkice_eff = tkice * denice/denh2o !effective conductivity since layer depth is constant
km = tkwat/cwat ! a constant (molecular diffusivity)
! Begin calculations
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
! Initialize Ebal quantities computed below
ocvts(c) = 0._r8
ncvts(c) = 0._r8
esum1(c) = 0._r8
esum2(c) = 0._r8
end do
! Initialize set of previous time-step variables as in DriverInit,
! which is currently not called over lakes. This has to be done
! here because phase change will occur in this routine.
! Ice fraction of snow at previous time step
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j >= snl(c) + 1) then
frac_iceold(c,j) = h2osoi_ice(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))
end if
end do
end do
! Sum soil water.
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j == 1) wsum(c) = 0._r8
wsum(c) = wsum(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
end do
end do
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
! Prepare for lake layer temperature calculations below
! fin(c) = betaprime * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
! eflx_sh_tot(p) + eflx_lh_tot(p))
! fin(c) now passed from ShalLakeFluxes as eflx_gnet
fin(c) = eflx_gnet(p)
end do
! 2!) Lake density
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
rhow(c,j) = (1._r8 - lake_icefrac(c,j)) * &
1000._r8*( 1.0_r8 - 1.9549e-05_r8*(abs(t_lake(c,j)-277._r8))**1.68_r8 ) &
+ lake_icefrac(c,j)*denice
! Allow for ice fraction; assume constant ice density.
! Is this the right weighted average?
! Using this average will make sure that surface ice is treated properly during
! convective mixing.
end do
end do
! 3!) Diffusivity and implied thermal "conductivity" = diffusivity * cwat
do j = 1, nlevlake-1
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
drhodz = (rhow(c,j+1)-rhow(c,j)) / (z_lake(c,j+1)-z_lake(c,j))
n2 = grav / rhow(c,j) * drhodz
! Fixed sign error here: our z goes up going down into the lake, so no negative
! sign is needed to make this positive unlike in Hostetler. --ZS
num = 40._r8 * n2 * (vkc*z_lake(c,j))**2
den = max( (ws(c)**2) * exp(-2._r8*ks(c)*z_lake(c,j)), 1.e-10_r8 )
ri = ( -1._r8 + sqrt( max(1._r8+num/den, 0._r8) ) ) / 20._r8
if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
! ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
if( t_lake(c,1) > 277.15_r8 ) then
if (lakedepth(c) > 15.0 ) then
ke = 1.e+2_r8*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
else
ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
endif
else
if (lakedepth(c) > 15.0 ) then
if (lakedepth(c) > 150.0 ) then
ke = 1.e+5_r8*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
else
ke =1.e+4_r8*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
end if
else
ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
endif
end if
kme(c,j) = km + ke
tk_lake(c,j) = kme(c,j)*cwat
! If there is some ice in this layer (this should rarely happen because the surface
! is unfrozen and it will be unstable), still use the cwat to get out the tk b/c the eddy
! diffusivity equation assumes water.
else
kme(c,j) = km
tk_lake(c,j) = tkwat*tkice_eff / ( (1._r8-lake_icefrac(c,j))*tkice_eff &
+ tkwat*lake_icefrac(c,j) )
! Assume the resistances add as for the calculation of conductivities at layer interfaces.
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
j = nlevlake
kme(c,nlevlake) = kme(c,nlevlake-1)
if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
tk_lake(c,j) = tk_lake(c,j-1)
else
tk_lake(c,j) = tkwat*tkice_eff / ( (1._r8-lake_icefrac(c,j))*tkice_eff &
+ tkwat*lake_icefrac(c,j) )
end if
! Use in surface flux calculation for next timestep.
savedtke1(c) = kme(c,1)*cwat ! Will only be used if unfrozen
! set number of column levels for use by Tridiagonal below
jtop(c) = snl(c) + 1
end do
! 4!) Heat source term: unfrozen lakes only
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
! Set eta(:), the extinction coefficient, according to L Hakanson, Aquatic Sciences, 1995
! (regression of Secchi Depth with lake depth for small glacial basin lakes), and the
! Poole & Atkins expression for extinction coeffient of 1.7 / Secchi Depth (m).
#ifndef ETALAKE
eta(:) = 1.1925_r8*lakedepth(c)**(-0.424)
#else
eta(:) = ETALAKE
#endif
zin = z_lake(c,j) - 0.5_r8*dz_lake(c,j)
zout = z_lake(c,j) + 0.5_r8*dz_lake(c,j)
rsfin = exp( -eta(islak)*max( zin-za(islak),0._r8 ) )
rsfout = exp( -eta(islak)*max( zout-za(islak),0._r8 ) )
! Let rsfout for bottom layer go into soil.
! This looks like it should be robust even for pathological cases,
! like lakes thinner than za.
if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
phidum = (rsfin-rsfout) * sabg(p) * (1._r8-beta(islak))
if (j == nlevlake) then
phi_soil(c) = rsfout * sabg(p) * (1._r8-beta(islak))
end if
else if (j == 1 .and. snl(c) == 0) then !if frozen but no snow layers
phidum = sabg(p) * (1._r8-beta(islak))
else !radiation absorbed at surface
phidum = 0._r8
if (j == nlevlake) phi_soil(c) = 0._r8
end if
phi(c,j) = phidum
end do
end do
! 5!) Set thermal properties and check initial energy content.
! For lake
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
cv_lake(c,j) = dz_lake(c,j) * (cwat*(1._r8-lake_icefrac(c,j)) + cice_eff*lake_icefrac(c,j))
end do
end do
! For snow / soil
call SoilThermProp_Lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
tk, cv, tktopsoillay)
! Sum cv*t_lake for energy check
! Include latent heat term, and correction for changing heat capacity with phase change.
! This will need to be over all soil / lake / snow layers. Lake is below.
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
! ocvts(c) = ocvts(c) + cv_lake(c,j)*t_lake(c,j) &
ocvts(c) = ocvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
+ cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) !&
! + (cwat-cice_eff)*lake_icefrac(c)*tfrz*dz_lake(c,j) !enthalpy reconciliation term
t_lake_bef(c,j) = t_lake(c,j)
end do
end do
! Now do for soil / snow layers
do j = -nlevsnow + 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j >= jtop(c)) then
! ocvts(c) = ocvts(c) + cv(c,j)*t_soisno(c,j) &
ocvts(c) = ocvts(c) + cv(c,j)*(t_soisno(c,j)-tfrz) &
+ hfus*h2osoi_liq(c,j) !&
! + (cpliq-cpice)*h2osoi_ice(c,j)*tfrz !enthalpy reconciliation term
if (j == 1 .and. h2osno(c) > 0._r8 .and. j == jtop(c)) then
ocvts(c) = ocvts(c) - h2osno(c)*hfus
end if
t_soisno_bef(c,j) = t_soisno(c,j)
if(abs(t_soisno(c,j)-288) > 150) then
WRITE( message,* ) 'WARNING: Extreme t_soisno at c, level',c, j
CALL wrf_error_fatal ( message )
endif
end if
end do
end do
!!!!!!!!!!!!!!!!!!!
! 6!) Set up vector r and vectors a, b, c1 that define tridiagonal matrix
! Heat capacity and resistance of snow without snow layers (<1cm) is ignored during diffusion,
! but its capacity to absorb latent heat may be used during phase change.
! Set up interface depths, zx, heat capacities, cvx, solar source terms, phix, and temperatures, tx.
do j = -nlevsnow+1, nlevlake+nlevsoil
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
jprime = j - nlevlake
if (j >= jtop(c)) then
if (j < 1) then !snow layer
zx(c,j) = z(c,j)
cvx(c,j) = cv(c,j)
phix(c,j) = 0._r8
tx(c,j) = t_soisno(c,j)
else if (j <= nlevlake) then !lake layer
zx(c,j) = z_lake(c,j)
cvx(c,j) = cv_lake(c,j)
phix(c,j) = phi(c,j)
tx(c,j) = t_lake(c,j)
else !soil layer
zx(c,j) = zx(c,nlevlake) + dz_lake(c,nlevlake)/2._r8 + z(c,jprime)
cvx(c,j) = cv(c,jprime)
if (j == nlevlake + 1) then !top soil layer
phix(c,j) = phi_soil(c)
else !middle or bottom soil layer
phix(c,j) = 0._r8
end if
tx(c,j) = t_soisno(c,jprime)
end if
end if
end do
end do
! Determine interface thermal conductivities, tkix
do j = -nlevsnow+1, nlevlake+nlevsoil
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
jprime = j - nlevlake
if (j >= jtop(c)) then
if (j < 0) then !non-bottom snow layer
tkix(c,j) = tk(c,j)
else if (j == 0) then !bottom snow layer
dzp = zx(c,j+1) - zx(c,j)
tkix(c,j) = tk_lake(c,1)*tk(c,j)*dzp / &
(tk(c,j)*z_lake(c,1) + tk_lake(c,1)*(-z(c,j)) )
! tk(c,0) is the conductivity at the middle of that layer, as defined in SoilThermProp_Lake
else if (j < nlevlake) then !non-bottom lake layer
tkix(c,j) = ( tk_lake(c,j)*tk_lake(c,j+1) * (dz_lake(c,j+1)+dz_lake(c,j)) ) &
/ ( tk_lake(c,j)*dz_lake(c,j+1) + tk_lake(c,j+1)*dz_lake(c,j) )
else if (j == nlevlake) then !bottom lake layer
dzp = zx(c,j+1) - zx(c,j)
tkix(c,j) = (tktopsoillay(c)*tk_lake(c,j)*dzp / &
(tktopsoillay(c)*dz_lake(c,j)/2._r8 + tk_lake(c,j)*z(c,1) ) )
! tktopsoillay is the conductivity at the middle of that layer, as defined in SoilThermProp_Lake
else !soil layer
tkix(c,j) = tk(c,jprime)
end if
end if
end do
end do
! Determine heat diffusion through the layer interface and factor used in computing
! tridiagonal matrix and set up vector r and vectors a, b, c1 that define tridiagonal
! matrix and solve system
do j = -nlevsnow+1, nlevlake+nlevsoil
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
if (j >= jtop(c)) then
if (j < nlevlake+nlevsoil) then !top or interior layer
factx(c,j) = dtime/cvx(c,j)
fnx(c,j) = tkix(c,j)*(tx(c,j+1)-tx(c,j))/(zx(c,j+1)-zx(c,j))
else !bottom soil layer
factx(c,j) = dtime/cvx(c,j)
fnx(c,j) = 0._r8 !not used
end if
end if
enddo
end do
do j = -nlevsnow+1,nlevlake+nlevsoil
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
if (j >= jtop(c)) then
if (j == jtop(c)) then !top layer
dzp = zx(c,j+1)-zx(c,j)
a(c,j) = 0._r8
b(c,j) = 1+(1._r8-cnfac)*factx(c,j)*tkix(c,j)/dzp
c1(c,j) = -(1._r8-cnfac)*factx(c,j)*tkix(c,j)/dzp
r(c,j) = tx(c,j) + factx(c,j)*( fin(c) + phix(c,j) + cnfac*fnx(c,j) )
else if (j < nlevlake+nlevsoil) then !middle layer
dzm = (zx(c,j)-zx(c,j-1))
dzp = (zx(c,j+1)-zx(c,j))
a(c,j) = - (1._r8-cnfac)*factx(c,j)* tkix(c,j-1)/dzm
b(c,j) = 1._r8+ (1._r8-cnfac)*factx(c,j)*(tkix(c,j)/dzp + tkix(c,j-1)/dzm)
c1(c,j) = - (1._r8-cnfac)*factx(c,j)* tkix(c,j)/dzp
r(c,j) = tx(c,j) + cnfac*factx(c,j)*( fnx(c,j) - fnx(c,j-1) ) + factx(c,j)*phix(c,j)
else !bottom soil layer
dzm = (zx(c,j)-zx(c,j-1))
a(c,j) = - (1._r8-cnfac)*factx(c,j)*tkix(c,j-1)/dzm
b(c,j) = 1._r8+ (1._r8-cnfac)*factx(c,j)*tkix(c,j-1)/dzm
c1(c,j) = 0._r8
r(c,j) = tx(c,j) - cnfac*factx(c,j)*fnx(c,j-1)
end if
end if
enddo
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 7!) Solve for tdsolution
call Tridiagonal(lbc, ubc, -nlevsnow + 1, nlevlake + nlevsoil, jtop, num_shlakec, filter_shlakec, &
a, b, c1, r, tx)
! Set t_soisno and t_lake
do j = -nlevsnow+1, nlevlake + nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
jprime = j - nlevlake
! Don't do anything with invalid snow layers.
if (j >= jtop(c)) then
if (j < 1) then !snow layer
t_soisno(c,j) = tx(c,j)
else if (j <= nlevlake) then !lake layer
t_lake(c,j) = tx(c,j)
else !soil layer
t_soisno(c,jprime) = tx(c,j)
end if
end if
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!
! 8!) Sum energy content and total energy into lake for energy check. Any errors will be from the
! Tridiagonal solution.
#if (defined LAKEDEBUG)
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
esum1(c) = esum1(c) + (t_lake(c,j)-t_lake_bef(c,j))*cv_lake(c,j)
esum2(c) = esum2(c) + (t_lake(c,j)-tfrz)*cv_lake(c,j)
end do
end do
do j = -nlevsnow+1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j >= jtop(c)) then
esum1(c) = esum1(c) + (t_soisno(c,j)-t_soisno_bef(c,j))*cv(c,j)
esum2(c) = esum2(c) + (t_soisno(c,j)-tfrz)*cv(c,j)
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
! Again assuming only one pft per column
! esum1(c) = esum1(c) + lhabs(c)
errsoi(c) = esum1(c)/dtime - eflx_soil_grnd(p)
! eflx_soil_grnd includes all the solar radiation absorbed in the lake,
! unlike eflx_gnet
if(abs(errsoi(c)) > 1.e-5_r8) then
WRITE( message,* )'Primary soil energy conservation error in shlake &
column during Tridiagonal Solution,', 'error (W/m^2):', c, errsoi(c)
CALL wrf_error_fatal ( message )
end if
end do
! This has to be done before convective mixing because the heat capacities for each layer
! will get scrambled.
#endif
!!!!!!!!!!!!!!!!!!!!!!!
! 9!) Phase change
call PhaseChange_Lake (snl,h2osno,dz,dz_lake, & !i
t_soisno,h2osoi_liq,h2osoi_ice, & !i&o
lake_icefrac,t_lake, snowdp, & !i&o
qflx_snomelt,eflx_snomelt,imelt, & !o
cv, cv_lake, & !i&o
lhabs) !o
!!!!!!!!!!!!!!!!!!!!!!!
! 9.5!) Second energy check and water check. Now check energy balance before and after phase
! change, considering the possibility of changed heat capacity during phase change, by
! using initial heat capacity in the first step, final heat capacity in the second step,
! and differences from tfrz only to avoid enthalpy correction for (cpliq-cpice)*melt*tfrz.
! Also check soil water sum.
#if (defined LAKEDEBUG)
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
esum2(c) = esum2(c) - (t_lake(c,j)-tfrz)*cv_lake(c,j)
end do
end do
do j = -nlevsnow+1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j >= jtop(c)) then
esum2(c) = esum2(c) - (t_soisno(c,j)-tfrz)*cv(c,j)
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
! Again assuming only one pft per column
esum2(c) = esum2(c) - lhabs(c)
errsoi(c) = esum2(c)/dtime
if(abs(errsoi(c)) > 1.e-5_r8) then
write(message,*)'Primary soil energy conservation error in shlake column during Phase Change, error (W/m^2):', &
c, errsoi(c)
CALL wrf_error_fatal ( message )
end if
end do
! Check soil water
! Sum soil water.
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j == 1) wsum_end(c) = 0._r8
wsum_end(c) = wsum_end(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
if (j == nlevsoil) then
if (abs(wsum(c)-wsum_end(c))>1.e-7_r8) then
write(message,*)'Soil water balance error during phase change in ShalLakeTemperature.', &
'column, error (kg/m^2):', c, wsum_end(c)-wsum(c)
CALL wrf_error_fatal ( message )
end if
end if
end do
end do
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!
! 10!) Convective mixing: make sure fracice*dz is conserved, heat content c*dz*T is conserved, and
! all ice ends up at the top. Done over all lakes even if frozen.
! Either an unstable density profile or ice in a layer below an incompletely frozen layer will trigger.
!Recalculate density
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
rhow(c,j) = (1._r8 - lake_icefrac(c,j)) * &
1000._r8*( 1.0_r8 - 1.9549e-05_r8*(abs(t_lake(c,j)-277._r8))**1.68_r8 ) &
+ lake_icefrac(c,j)*denice
end do
end do
do j = 1, nlevlake-1
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
qav(c) = 0._r8
nav(c) = 0._r8
iceav(c) = 0._r8
end do
do i = 1, j+1
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (rhow(c,j) > rhow(c,j+1) .or. &
(lake_icefrac(c,j) < 1._r8 .and. lake_icefrac(c,j+1) > 0._r8) ) then
#if (defined LAKEDEBUG)
if (i==1) then
write(message,*), 'Convective Mixing in column ', c, '.'
CALL wrf_message(message)
endif
#endif
qav(c) = qav(c) + dz_lake(c,i)*(t_lake(c,i)-tfrz) * &
((1._r8 - lake_icefrac(c,i))*cwat + lake_icefrac(c,i)*cice_eff)
! tav(c) = tav(c) + t_lake(c,i)*dz_lake(c,i)
iceav(c) = iceav(c) + lake_icefrac(c,i)*dz_lake(c,i)
nav(c) = nav(c) + dz_lake(c,i)
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (rhow(c,j) > rhow(c,j+1) .or. &
(lake_icefrac(c,j) < 1._r8 .and. lake_icefrac(c,j+1) > 0._r8) ) then
qav(c) = qav(c)/nav(c)
iceav(c) = iceav(c)/nav(c)
!If the average temperature is above freezing, put the extra energy into the water.
!If it is below freezing, take it away from the ice.
if (qav(c) > 0._r8) then
tav_froz(c) = 0._r8 !Celsius
tav_unfr(c) = qav(c) / ((1._r8 - iceav(c))*cwat)
else if (qav(c) < 0._r8) then
tav_froz(c) = qav(c) / (iceav(c)*cice_eff)
tav_unfr(c) = 0._r8 !Celsius
else
tav_froz(c) = 0._r8
tav_unfr(c) = 0._r8
end if
end if
end do
do i = 1, j+1
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (nav(c) > 0._r8) then
! if(0==1) then
!Put all the ice at the top.!
!If the average temperature is above freezing, put the extra energy into the water.
!If it is below freezing, take it away from the ice.
!For the layer with both ice & water, be careful to use the average temperature
!that preserves the correct total heat content given what the heat capacity of that
!layer will actually be.
if (i == 1) zsum(c) = 0._r8
if ((zsum(c)+dz_lake(c,i))/nav(c) <= iceav(c)) then
lake_icefrac(c,i) = 1._r8
t_lake(c,i) = tav_froz(c) + tfrz
else if (zsum(c)/nav(c) < iceav(c)) then
lake_icefrac(c,i) = (iceav(c)*nav(c) - zsum(c)) / dz_lake(c,i)
! Find average value that preserves correct heat content.
t_lake(c,i) = ( lake_icefrac(c,i)*tav_froz(c)*cice_eff &
+ (1._r8 - lake_icefrac(c,i))*tav_unfr(c)*cwat ) &
/ ( lake_icefrac(c,i)*cice_eff + (1-lake_icefrac(c,i))*cwat ) + tfrz
else
lake_icefrac(c,i) = 0._r8
t_lake(c,i) = tav_unfr(c) + tfrz
end if
zsum(c) = zsum(c) + dz_lake(c,i)
rhow(c,i) = (1._r8 - lake_icefrac(c,i)) * &
1000._r8*( 1.0_r8 - 1.9549e-05_r8*(abs(t_lake(c,i)-277._r8))**1.68_r8 ) &
+ lake_icefrac(c,i)*denice
end if
end do
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!
! 11!) Re-evaluate thermal properties and sum energy content.
! For lake
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
cv_lake(c,j) = dz_lake(c,j) * (cwat*(1._r8-lake_icefrac(c,j)) + cice_eff*lake_icefrac(c,j))
#if (defined LAKEDEBUG)
write(message,*)'Lake Ice Fraction, c, level:', c, j, lake_icefrac(c,j)
CALL wrf_message(message)
#endif
end do
end do
! For snow / soil
! call SoilThermProp_Lake(lbc, ubc, num_shlakec, filter_shlakec, tk, cv, tktopsoillay)
call SoilThermProp_Lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
tk, cv, tktopsoillay)
! Do as above to sum energy content
do j = 1, nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
! ncvts(c) = ncvts(c) + cv_lake(c,j)*t_lake(c,j) &
ncvts(c) = ncvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
+ cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) !&
! + (cwat-cice_eff)*lake_icefrac(c)*tfrz*dz_lake(c,j) !enthalpy reconciliation term
fin(c) = fin(c) + phi(c,j)
end do
end do
do j = -nlevsnow + 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (j >= jtop(c)) then
! ncvts(c) = ncvts(c) + cv(c,j)*t_soisno(c,j) &
ncvts(c) = ncvts(c) + cv(c,j)*(t_soisno(c,j)-tfrz) &
+ hfus*h2osoi_liq(c,j) !&
! + (cpliq-cpice)*h2osoi_ice(c,j)*tfrz !enthalpy reconciliation term
if (j == 1 .and. h2osno(c) > 0._r8 .and. j == jtop(c)) then
ncvts(c) = ncvts(c) - h2osno(c)*hfus
end if
end if
if (j == 1) fin(c) = fin(c) + phi_soil(c)
end do
end do
! Check energy conservation.
do fp = 1, num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
errsoi(c) = (ncvts(c)-ocvts(c)) / dtime - fin(c)
#ifndef LAKEDEBUG
! if (abs(errsoi(c)) < 0.10_r8) then ! else send to Balance Check and abort
if (abs(errsoi(c)) < 10._r8) then ! else send to Balance Check and abort
#else
if (abs(errsoi(c)) < 1._r8) then
#endif
eflx_sh_tot(p) = eflx_sh_tot(p) - errsoi(c)
eflx_sh_grnd(p) = eflx_sh_grnd(p) - errsoi(c)
eflx_soil_grnd(p) = eflx_soil_grnd(p) + errsoi(c)
eflx_gnet(p) = eflx_gnet(p) + errsoi(c)
! if (abs(errsoi(c)) > 1.e-3_r8) then
if (abs(errsoi(c)) > 1.e-1_r8) then
write(message,*)'errsoi incorporated into sensible heat in ShalLakeTemperature: c, (W/m^2):', c, errsoi(c)
CALL wrf_message(message)
end if
errsoi(c) = 0._r8
#if (defined LAKEDEBUG)
else
write(message,*)'Soil Energy Balance Error at column, ', c, 'G, fintotal, column E tendency = ', &
eflx_gnet(p), fin(c), (ncvts(c)-ocvts(c)) / dtime
CALL wrf_message(message)
#endif
end if
end do
! This loop assumes only one point per column.
end subroutine ShalLakeTemperature
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!BOP
!
! ROUTINE: SoilThermProp_Lake
!
! !INTERFACE:
subroutine SoilThermProp_Lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, & 2,1
tk, cv, tktopsoillay)
!
! !DESCRIPTION:
! Calculation of thermal conductivities and heat capacities of
! snow/soil layers
! (1) The volumetric heat capacity is calculated as a linear combination
! in terms of the volumetric fraction of the constituent phases.
!
! (2) The thermal conductivity of soil is computed from the algorithm of
! Johansen (as reported by Farouki 1981), and of snow is from the
! formulation used in SNTHERM (Jordan 1991).
! The thermal conductivities at the interfaces between two neighboring
! layers (j, j+1) are derived from an assumption that the flux across
! the interface is equal to that from the node j to the interface and the
! flux from the interface to the node j+1.
!
! For lakes, the proper soil layers (not snow) should always be saturated.
!
! !USES:
implicit none
!in
integer , intent(in) :: snl(1) ! number of snow layers
! real(r8), intent(in) :: h2osno(1) ! snow water (mm H2O)
! real(r8), intent(in) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
! real(r8), intent(in) :: tksatu(1,nlevsoil) ! thermal conductivity, saturated soil [W/m-K]
! real(r8), intent(in) :: tkmg(1,nlevsoil) ! thermal conductivity, soil minerals [W/m-K]
! real(r8), intent(in) :: tkdry(1,nlevsoil) ! thermal conductivity, dry soil (W/m/Kelvin)
! real(r8), intent(in) :: csol(1,nlevsoil) ! heat capacity, soil solids (J/m**3/Kelvin)
real(r8), intent(in) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness (m)
real(r8), intent(in) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
real(r8), intent(in) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth (m)
real(r8), intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil temperature (Kelvin)
real(r8), intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
real(r8), intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
!out
real(r8), intent(out) :: cv(lbc:ubc,-nlevsnow+1:nlevsoil) ! heat capacity [J/(m2 K)]
real(r8), intent(out) :: tk(lbc:ubc,-nlevsnow+1:nlevsoil) ! thermal conductivity [W/(m K)]
real(r8), intent(out) :: tktopsoillay(lbc:ubc) ! thermal conductivity [W/(m K)]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !CALLED FROM:
! subroutine ShalLakeTemperature in this module.
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/13/02, Peter Thornton: migrated to new data structures
! 7/01/03, Mariana Vertenstein: migrated to vector code
! 4/09, Zack Subin, adjustment for ShalLake code.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !LOCAL VARIABLES:
!
! local pointers to original implicit in scalars
!
! integer , pointer :: clandunit(:) ! column's landunit
! integer , pointer :: ityplun(:) ! landunit type
!
!EOP
! OTHER LOCAL VARIABLES:
integer :: l,c,j ! indices
integer :: fc ! lake filtered column indices
real(r8) :: bw ! partial density of water (ice + liquid)
real(r8) :: dksat ! thermal conductivity for saturated soil (j/(k s m))
real(r8) :: dke ! kersten number
real(r8) :: fl ! fraction of liquid or unfrozen water to total water
real(r8) :: satw ! relative total water content of soil.
real(r8) :: thk(lbc:ubc,-nlevsnow+1:nlevsoil) ! thermal conductivity of layer
character*256 :: message
! Thermal conductivity of soil from Farouki (1981)
do j = -nlevsnow+1,nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
! Only examine levels from 1->nlevsoil
if (j >= 1) then
! l = clandunit(c)
! if (ityplun(l) /= istwet .AND. ityplun(l) /= istice) then
! This could be altered later for allowing this to be over glaciers.
! Soil should be saturated.
#if (defined LAKEDEBUG)
satw = (h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice)/(dz(c,j)*watsat(c,j))
! satw = min(1._r8, satw)
if (satw < 0.999_r8) then
write(message,*)'WARNING: soil layer unsaturated in SoilThermProp_Lake, satw, j = ', satw, j
CALL wrf_error_fatal ( message )
end if
! Could use denice because if it starts out frozen, the volume of water will go below sat.,
! since we're not yet doing excess ice.
! But take care of this in HydrologyLake.
#endif
satw = 1._r8
fl = h2osoi_liq(c,j)/(h2osoi_ice(c,j)+h2osoi_liq(c,j))
if (t_soisno(c,j) >= tfrz) then ! Unfrozen soil
dke = max(0._r8, log10(satw) + 1.0_r8)
dksat = tksatu(c,j)
else ! Frozen soil
dke = satw
dksat = tkmg(c,j)*0.249_r8**(fl*watsat(c,j))*2.29_r8**watsat(c,j)
endif
thk(c,j) = dke*dksat + (1._r8-dke)*tkdry(c,j)
! else
! thk(c,j) = tkwat
! if (t_soisno(c,j) < tfrz) thk(c,j) = tkice
! endif
endif
! Thermal conductivity of snow, which from Jordan (1991) pp. 18
! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1
if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then
bw = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/dz(c,j)
thk(c,j) = tkairc + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkairc)
end if
end do
end do
! Thermal conductivity at the layer interface
! Have to correct for the fact that bottom snow layer and top soil layer border lake.
! For the first case, the snow layer conductivity for the middle of the layer will be returned.
! Because the interfaces are below the soil layers, the conductivity for the top soil layer
! will have to be returned separately.
do j = -nlevsnow+1,nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
if (j >= snl(c)+1 .AND. j <= nlevsoil-1 .AND. j /= 0) then
tk(c,j) = thk(c,j)*thk(c,j+1)*(z(c,j+1)-z(c,j)) &
/(thk(c,j)*(z(c,j+1)-zi(c,j))+thk(c,j+1)*(zi(c,j)-z(c,j)))
else if (j == 0) then
tk(c,j) = thk(c,j)
else if (j == nlevsoil) then
tk(c,j) = 0._r8
end if
! For top soil layer.
if (j == 1) tktopsoillay(c) = thk(c,j)
end do
end do
! Soil heat capacity, from de Vires (1963)
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
! l = clandunit(c)
! if (ityplun(l) /= istwet .AND. ityplun(l) /= istice) then
cv(c,j) = csol(c,j)*(1-watsat(c,j))*dz(c,j) + &
(h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq)
! else
! cv(c,j) = (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq)
! endif
! if (j == 1) then
! if (snl(c)+1 == 1 .AND. h2osno(c) > 0._r8) then
! cv(c,j) = cv(c,j) + cpice*h2osno(c)
! end if
! end if
! Won't worry about heat capacity for thin snow on lake with no snow layers.
enddo
end do
! Snow heat capacity
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
if (snl(c)+1 < 1 .and. j >= snl(c)+1) then
cv(c,j) = cpliq*h2osoi_liq(c,j) + cpice*h2osoi_ice(c,j)
end if
end do
end do
end subroutine SoilThermProp_Lake
!-----------------------------------------------------------------------
!BOP
!
! ROUTINE: PhaseChange_Lake
!
! !INTERFACE:
subroutine PhaseChange_Lake (snl,h2osno,dz,dz_lake, & !i 1
t_soisno,h2osoi_liq,h2osoi_ice, & !i&o
lake_icefrac,t_lake, snowdp, & !i&o
qflx_snomelt,eflx_snomelt,imelt, & !o
cv, cv_lake, & !i&o
lhabs) !o
!=============================================================================================
! !DESCRIPTION:
! Calculation of the phase change within snow, soil, & lake layers:
! (1) Check the conditions for which the phase change may take place,
! i.e., the layer temperature is great than the freezing point
! and the ice mass is not equal to zero (i.e. melting),
! or the layer temperature is less than the freezing point
! and the liquid water mass is greater than the allowable supercooled
! (i.e. freezing).
! (2) Assess the amount of phase change from the energy excess (or deficit)
! after setting the layer temperature to freezing point, depending on
! how much water or ice is available.
! (3) Re-adjust the ice and liquid mass, and the layer temperature: either to
! the freezing point if enough water or ice is available to fully compensate,
! or to a remaining temperature.
! The specific heats are assumed constant. Potential cycling errors resulting from
! this assumption will be trapped at the end of ShalLakeTemperature.
! !CALLED FROM:
! subroutine ShalLakeTemperature in this module
!
! !REVISION HISTORY:
! 04/2009 Zack Subin: Initial code
!==============================================================================================
! !USES:
!
! !ARGUMENTS:
implicit none
!in:
integer , intent(in) :: snl(1) !number of snow layers
real(r8), intent(inout) :: h2osno(1) !snow water (mm H2O)
real(r8), intent(in) :: dz(1,-nlevsnow+1:nlevsoil) !layer thickness (m)
real(r8), intent(in) :: dz_lake(1,nlevlake) !lake layer thickness (m)
! Needed in case snow height is less than critical value.
!inout:
real(r8), intent(inout) :: snowdp(1) !snow height (m)
real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
real(r8), intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
real(r8), intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
!out:
real(r8), intent(out) :: qflx_snomelt(1) !snow melt (mm H2O /s)
real(r8), intent(out) :: eflx_snomelt(1) !snow melt heat flux (W/m**2)
integer, intent(out) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
!What's the sign of this? Is it just output?
real(r8), intent(inout) :: cv(lbc:ubc,-nlevsnow+1:nlevsoil) ! heat capacity [J/(m2 K)]
real(r8), intent(inout) :: cv_lake (lbc:ubc,1:nlevlake) ! heat capacity [J/(m2 K)]
real(r8), intent(out):: lhabs(lbc:ubc) ! total per-column latent heat abs. (J/m^2)
! OTHER LOCAL VARIABLES:
integer :: j,c,g !do loop index
integer :: fc !lake filtered column indices
! real(r8) :: dtime !land model time step (sec)
real(r8) :: heatavail !available energy for melting or freezing (J/m^2)
real(r8) :: heatrem !energy residual or loss after melting or freezing
real(r8) :: melt !actual melting (+) or freezing (-) [kg/m2]
real(r8), parameter :: smallnumber = 1.e-7_r8 !to prevent tiny residuals from rounding error
logical :: dophasechangeflag
!-----------------------------------------------------------------------
! dtime = get_step_size()
! Initialization
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
qflx_snomelt(c) = 0._r8
eflx_snomelt(c) = 0._r8
lhabs(c) = 0._r8
end do
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
if (j >= snl(c) + 1) imelt(c,j) = 0
end do
end do
! Check for case of snow without snow layers and top lake layer temp above freezing.
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
if (snl(c) == 0 .and. h2osno(c) > 0._r8 .and. t_lake(c,1) > tfrz) then
heatavail = (t_lake(c,1) - tfrz) * cv_lake(c,1)
melt = min(h2osno(c), heatavail/hfus)
heatrem = max(heatavail - melt*hfus, 0._r8)
!catch small negative value to keep t at tfrz
t_lake(c,1) = tfrz + heatrem/(cv_lake(c,1))
snowdp(c) = snowdp(c)*(1._r8 - melt/h2osno(c))
h2osno(c) = h2osno(c) - melt
lhabs(c) = lhabs(c) + melt*hfus
qflx_snomelt(c) = qflx_snomelt(c) + melt
! Prevent tiny residuals
if (h2osno(c) < smallnumber) h2osno(c) = 0._r8
if (snowdp(c) < smallnumber) snowdp(c) = 0._r8
end if
end do
! Lake phase change
do j = 1,nlevlake
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
dophasechangeflag = .false.
if (t_lake(c,j) > tfrz .and. lake_icefrac(c,j) > 0._r8) then ! melting
dophasechangeflag = .true.
heatavail = (t_lake(c,j) - tfrz) * cv_lake(c,j)
melt = min(lake_icefrac(c,j)*denh2o*dz_lake(c,j), heatavail/hfus)
!denh2o is used because layer thickness is not adjusted for freezing
heatrem = max(heatavail - melt*hfus, 0._r8)
!catch small negative value to keep t at tfrz
else if (t_lake(c,j) < tfrz .and. lake_icefrac(c,j) < 1._r8) then !freezing
dophasechangeflag = .true.
heatavail = (t_lake(c,j) - tfrz) * cv_lake(c,j)
melt = max(-(1._r8-lake_icefrac(c,j))*denh2o*dz_lake(c,j), heatavail/hfus)
!denh2o is used because layer thickness is not adjusted for freezing
heatrem = min(heatavail - melt*hfus, 0._r8)
!catch small positive value to keep t at tfrz
end if
! Update temperature and ice fraction.
if (dophasechangeflag) then
lake_icefrac(c,j) = lake_icefrac(c,j) - melt/(denh2o*dz_lake(c,j))
lhabs(c) = lhabs(c) + melt*hfus
! Update heat capacity
cv_lake(c,j) = cv_lake(c,j) + melt*(cpliq-cpice)
t_lake(c,j) = tfrz + heatrem/cv_lake(c,j)
! Prevent tiny residuals
if (lake_icefrac(c,j) > 1._r8 - smallnumber) lake_icefrac(c,j) = 1._r8
if (lake_icefrac(c,j) < smallnumber) lake_icefrac(c,j) = 0._r8
end if
end do
end do
! Snow & soil phase change
do j = -nlevsnow+1,nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
dophasechangeflag = .false.
if (j >= snl(c) + 1) then
if (t_soisno(c,j) > tfrz .and. h2osoi_ice(c,j) > 0._r8) then ! melting
dophasechangeflag = .true.
heatavail = (t_soisno(c,j) - tfrz) * cv(c,j)
melt = min(h2osoi_ice(c,j), heatavail/hfus)
heatrem = max(heatavail - melt*hfus, 0._r8)
!catch small negative value to keep t at tfrz
if (j <= 0) then !snow
imelt(c,j) = 1
qflx_snomelt(c) = qflx_snomelt(c) + melt
end if
else if (t_soisno(c,j) < tfrz .and. h2osoi_liq(c,j) > 0._r8) then !freezing
dophasechangeflag = .true.
heatavail = (t_soisno(c,j) - tfrz) * cv(c,j)
melt = max(-h2osoi_liq(c,j), heatavail/hfus)
heatrem = min(heatavail - melt*hfus, 0._r8)
!catch small positive value to keep t at tfrz
if (j <= 0) then !snow
imelt(c,j) = 2
qflx_snomelt(c) = qflx_snomelt(c) + melt
! Does this works for both signs of melt in SnowHydrology? I think
! qflx_snomelt(c) is just output.
end if
end if
! Update temperature and soil components.
if (dophasechangeflag) then
h2osoi_ice(c,j) = h2osoi_ice(c,j) - melt
h2osoi_liq(c,j) = h2osoi_liq(c,j) + melt
lhabs(c) = lhabs(c) + melt*hfus
! Update heat capacity
cv(c,j) = cv(c,j) + melt*(cpliq-cpice)
t_soisno(c,j) = tfrz + heatrem/cv(c,j)
! Prevent tiny residuals
if (h2osoi_ice(c,j) < smallnumber) h2osoi_ice(c,j) = 0._r8
if (h2osoi_liq(c,j) < smallnumber) h2osoi_liq(c,j) = 0._r8
end if
end if
end do
end do
! Update eflx_snomelt(c)
!dir$ concurrent
!cdir nodep
do fc = 1,num_shlakec
c = filter_shlakec(fc)
eflx_snomelt(c) = qflx_snomelt(c)*hfus
end do
!!!
end subroutine PhaseChange_Lake
subroutine ShalLakeHydrology(dz_lake,forc_rain,forc_snow, & !i 1,9
begwb,qflx_evap_tot,forc_t,do_capsnow, &
t_grnd,qflx_evap_soi, &
qflx_snomelt,imelt,frac_iceold, & !i add by guhp
z,dz,zi,snl,h2osno,snowdp,lake_icefrac,t_lake, & !i&o
endwb,snowage,snowice,snowliq,t_snow, & !o
t_soisno,h2osoi_ice,h2osoi_liq,h2osoi_vol, &
qflx_drain,qflx_surf,qflx_infl,qflx_qrgwl, &
qcharge,qflx_prec_grnd,qflx_snowcap, &
qflx_snowcap_col,qflx_snow_grnd_pft, &
qflx_snow_grnd_col,qflx_rain_grnd, &
qflx_evap_tot_col,soilalpha,zwt,fcov, &
rootr_column,qflx_evap_grnd,qflx_sub_snow, &
qflx_dew_snow,qflx_dew_grnd,qflx_rain_grnd_col)
!==================================================================================
! !DESCRIPTION:
! Calculation of Shallow Lake Hydrology. Full hydrology of snow layers is
! done. However, there is no infiltration, and the water budget is balanced with
! qflx_qrgwl. Lake water mass is kept constant. The soil is simply maintained at
! volumetric saturation if ice melting frees up pore space. Likewise, if the water
! portion alone at some point exceeds pore capacity, it is reduced. This is consistent
! with the possibility of initializing the soil layer with excess ice. The only
! real error with that is that the thermal conductivity will ignore the excess ice
! (and accompanying thickness change).
!
! If snow layers are present over an unfrozen lake, and the top layer of the lake
! is capable of absorbing the latent heat without going below freezing,
! the snow-water is runoff and the latent heat is subtracted from the lake.
!
! WARNING: This subroutine assumes lake columns have one and only one pft.
!
! Sequence is:
! ShalLakeHydrology:
! Do needed tasks from Hydrology1, Biogeophysics2, & top of Hydrology2.
! -> SnowWater: change of snow mass and snow water onto soil
! -> SnowCompaction: compaction of snow layers
! -> CombineSnowLayers: combine snow layers that are thinner than minimum
! -> DivideSnowLayers: subdivide snow layers that are thicker than maximum
! Add water to soil if melting has left it with open pore space.
! Cleanup and do water balance.
! If snow layers are found above a lake with unfrozen top layer, whose top
! layer has enough heat to melt all the snow ice without freezing, do so
! and eliminate the snow layers.
!
! !REVISION HISTORY:
! Created by Zack Subin, 2009
!
!============================================================================================
! USES:
!
implicit none
! in:
! integer , intent(in) :: clandunit(1) ! column's landunit
! integer , intent(in) :: ityplun(1) ! landunit type
! real(r8), intent(in) :: watsat(1,1:nlevsoil) ! volumetric soil water at saturation (porosity)
real(r8), intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
real(r8), intent(in) :: forc_rain(1) ! rain rate [mm/s]
real(r8), intent(in) :: forc_snow(1) ! snow rate [mm/s]
real(r8), intent(in) :: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
real(r8), intent(in) :: forc_t(1) ! atmospheric temperature (Kelvin)
#if (defined OFFLINE)
real(r8), intent(in) :: flfall(1) ! fraction of liquid water within falling precipitation
#endif
logical , intent(in) :: do_capsnow(1) ! true => do snow capping
real(r8), intent(in) :: t_grnd(1) ! ground temperature (Kelvin)
real(r8), intent(in) :: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
real(r8), intent(in) :: qflx_snomelt(1) !snow melt (mm H2O /s)
integer, intent(in) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0
!inout:
real(r8), intent(inout) :: begwb(1) ! water mass begining of the time step
! inout:
real(r8), intent(inout) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth (m)
real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness depth (m)
real(r8), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) ! interface depth (m)
integer , intent(inout) :: snl(1) ! number of snow layers
real(r8), intent(inout) :: h2osno(1) ! snow water (mm H2O)
real(r8), intent(inout) :: snowdp(1) ! snow height (m)
real(r8), intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
real(r8), intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
real(r8), intent(inout) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
! out:
real(r8), intent(out) :: endwb(1) ! water mass end of the time step
real(r8), intent(out) :: snowage(1) ! non dimensional snow age [-]
real(r8), intent(out) :: snowice(1) ! average snow ice lens
real(r8), intent(out) :: snowliq(1) ! average snow liquid water
real(r8), intent(out) :: t_snow(1) ! vertically averaged snow temperature
real(r8), intent(out) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! snow temperature (Kelvin)
real(r8), intent(out) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
real(r8), intent(out) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
real(r8), intent(out) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
real(r8), intent(out) :: qflx_drain(1) ! sub-surface runoff (mm H2O /s)
real(r8), intent(out) :: qflx_surf(1) ! surface runoff (mm H2O /s)
real(r8), intent(out) :: qflx_infl(1) ! infiltration (mm H2O /s)
real(r8), intent(out) :: qflx_qrgwl(1) ! qflx_surf at glaciers, wetlands, lakes
real(r8), intent(out) :: qcharge(1) ! aquifer recharge rate (mm/s)
real(r8), intent(out) :: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
real(r8), intent(out) :: qflx_snowcap(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
real(r8), intent(out) :: qflx_snowcap_col(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
real(r8), intent(out) :: qflx_snow_grnd_pft(1) ! snow on ground after interception (mm H2O/s) [+]
real(r8), intent(out) :: qflx_snow_grnd_col(1) ! snow on ground after interception (mm H2O/s) [+]
real(r8), intent(out) :: qflx_rain_grnd(1) ! rain on ground after interception (mm H2O/s) [+]
real(r8), intent(out) :: qflx_evap_tot_col(1) !pft quantity averaged to the column (assuming one pft)
real(r8) ,intent(out) :: soilalpha(1) !factor that reduces ground saturated specific humidity (-)
real(r8), intent(out) :: zwt(1) !water table depth
real(r8), intent(out) :: fcov(1) !fractional area with water table at surface
real(r8), intent(out) :: rootr_column(1,1:nlevsoil) !effective fraction of roots in each soil layer
real(r8), intent(out) :: qflx_evap_grnd(1) ! ground surface evaporation rate (mm H2O/s) [+]
real(r8), intent(out) :: qflx_sub_snow(1) ! sublimation rate from snow pack (mm H2O /s) [+]
real(r8), intent(out) :: qflx_dew_snow(1) ! surface dew added to snow pack (mm H2O /s) [+]
real(r8), intent(out) :: qflx_dew_grnd(1) ! ground surface dew formation (mm H2O /s) [+]
real(r8), intent(out) :: qflx_rain_grnd_col(1) !rain on ground after interception (mm H2O/s) [+]
! Block of biogeochem currently not used.
#ifndef SHLAKE
real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm)
real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b"
real(r8), pointer :: bsw2(:,:) ! Clapp and Hornberger "b" for CN code
real(r8), pointer :: psisat(:,:) ! soil water potential at saturation for CN code (MPa)
real(r8), pointer :: vwcsat(:,:) ! volumetric water content at saturation for CN code (m3/m3)
real(r8), pointer :: wf(:) ! soil water as frac. of whc for top 0.5 m
real(r8), pointer :: soilpsi(:,:) ! soil water potential in each soil layer (MPa)
real(r8) :: psi,vwc,fsat ! temporary variables for soilpsi calculation
#if (defined DGVM) || (defined CN)
real(r8) :: watdry ! temporary
real(r8) :: rwat(lbc:ubc) ! soil water wgted by depth to maximum depth of 0.5 m
real(r8) :: swat(lbc:ubc) ! same as rwat but at saturation
real(r8) :: rz(lbc:ubc) ! thickness of soil layers contributing to rwat (m)
real(r8) :: tsw ! volumetric soil water to 0.5 m
real(r8) :: stsw ! volumetric soil water to 0.5 m at saturation
#endif
#endif
! OTHER LOCAL VARIABLES:
integer :: p,fp,g,l,c,j,fc,jtop ! indices
integer :: num_shlakesnowc ! number of column snow points
integer :: filter_shlakesnowc(ubc-lbc+1) ! column filter for snow points
integer :: num_shlakenosnowc ! number of column non-snow points
integer :: filter_shlakenosnowc(ubc-lbc+1) ! column filter for non-snow points
! real(r8) :: dtime ! land model time step (sec)
integer :: newnode ! flag when new snow node is set, (1=yes, 0=no)
real(r8) :: dz_snowf ! layer thickness rate change due to precipitation [mm/s]
real(r8) :: bifall ! bulk density of newly fallen dry snow [kg/m3]
real(r8) :: fracsnow(lbp:ubp) ! frac of precipitation that is snow
real(r8) :: fracrain(lbp:ubp) ! frac of precipitation that is rain
real(r8) :: qflx_prec_grnd_snow(lbp:ubp) ! snow precipitation incident on ground [mm/s]
real(r8) :: qflx_prec_grnd_rain(lbp:ubp) ! rain precipitation incident on ground [mm/s]
real(r8) :: qflx_evap_soi_lim ! temporary evap_soi limited by top snow layer content [mm/s]
real(r8) :: h2osno_temp ! temporary h2osno [kg/m^2]
real(r8), parameter :: snow_bd = 250._r8 !constant snow bulk density (only used in special case here) [kg/m^3]
real(r8) :: sumsnowice(lbc:ubc) ! sum of snow ice if snow layers found above unfrozen lake [kg/m&2]
logical :: unfrozen(lbc:ubc) ! true if top lake layer is unfrozen with snow layers above
real(r8) :: heatrem ! used in case above [J/m^2]
real(r8) :: heatsum(lbc:ubc) ! used in case above [J/m^2]
real(r8) :: qflx_top_soil(1) !net water input into soil from top (mm/s)
character*256 :: message
#if (defined LAKEDEBUG)
real(r8) :: snow_water(lbc:ubc) ! temporary sum of snow water for Bal Check [kg/m^2]
#endif
!-----------------------------------------------------------------------
! Determine step size
! dtime = get_step_size()
! Add soil water to water balance.
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
begwb(c) = begwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Do precipitation onto ground, etc., from Hydrology1.
!dir$ concurrent
!cdir nodep
do fp = 1, num_shlakep
p = filter_shlakep(fp)
g = pgridcell(p)
! l = plandunit(p)
c = pcolumn(p)
! Precipitation onto ground (kg/(m2 s))
! ! PET, 1/18/2005: Added new terms for mass balance correction
! ! due to dynamic pft weight shifting (column-level h2ocan_loss)
! ! Because the fractionation between rain and snow is indeterminate if
! ! rain + snow = 0, I am adding this very small flux only to the rain
! ! components.
! Not relevant unless PFTs are added to lake later.
! if (frac_veg_nosno(p) == 0) then
qflx_prec_grnd_snow(p) = forc_snow(g)
qflx_prec_grnd_rain(p) = forc_rain(g) !+ h2ocan_loss(c)
! else
! qflx_prec_grnd_snow(p) = qflx_through_snow(p) + (qflx_candrip(p) * fracsnow(p))
! qflx_prec_grnd_rain(p) = qflx_through_rain(p) + (qflx_candrip(p) * fracrain(p)) + h2ocan_loss(c)
! end if
qflx_prec_grnd(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
if (do_capsnow(c)) then
qflx_snowcap(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
qflx_snow_grnd_pft(p) = 0._r8
qflx_rain_grnd(p) = 0._r8
else
qflx_snowcap(p) = 0._r8
#if (defined OFFLINE)
qflx_snow_grnd_pft(p) = qflx_prec_grnd(p)*(1._r8-flfall(g)) ! ice onto ground (mm/s)
qflx_rain_grnd(p) = qflx_prec_grnd(p)*flfall(g) ! liquid water onto ground (mm/s)
#else
qflx_snow_grnd_pft(p) = qflx_prec_grnd_snow(p) ! ice onto ground (mm/s)
qflx_rain_grnd(p) = qflx_prec_grnd_rain(p) ! liquid water onto ground (mm/s)
#endif
end if
! Assuming one PFT; needed for below
qflx_snow_grnd_col(c) = qflx_snow_grnd_pft(p)
qflx_rain_grnd_col(c) = qflx_rain_grnd(p)
end do ! (end pft loop)
! Determine snow height and snow water
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
! l = clandunit(c)
g = cgridcell(c)
! Use Alta relationship, Anderson(1976); LaChapelle(1961),
! U.S.Department of Agriculture Forest Service, Project F,
! Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
if (do_capsnow(c)) then
dz_snowf = 0._r8
else
if (forc_t(g) > tfrz + 2._r8) then
bifall=50._r8 + 1.7_r8*(17.0_r8)**1.5_r8
else if (forc_t(g) > tfrz - 15._r8) then
bifall=50._r8 + 1.7_r8*(forc_t(g) - tfrz + 15._r8)**1.5_r8
else
bifall=50._r8
end if
dz_snowf = qflx_snow_grnd_col(c)/bifall
snowdp(c) = snowdp(c) + dz_snowf*dtime
h2osno(c) = h2osno(c) + qflx_snow_grnd_col(c)*dtime ! snow water equivalent (mm)
end if
! if (itype(l)==istwet .and. t_grnd(c)>tfrz) then
! h2osno(c)=0._r8
! snowdp(c)=0._r8
! snowage(c)=0._r8
! end if
! Take care of this later in function.
! When the snow accumulation exceeds 10 mm, initialize snow layer
! Currently, the water temperature for the precipitation is simply set
! as the surface air temperature
newnode = 0 ! flag for when snow node will be initialized
if (snl(c) == 0 .and. qflx_snow_grnd_col(c) > 0.0_r8 .and. snowdp(c) >= 0.01_r8) then
newnode = 1
snl(c) = -1
dz(c,0) = snowdp(c) ! meter
z(c,0) = -0.5_r8*dz(c,0)
zi(c,-1) = -dz(c,0)
snowage(c) = 0._r8 ! snow age
t_soisno(c,0) = min(tfrz, forc_t(g)) ! K
h2osoi_ice(c,0) = h2osno(c) ! kg/m2
h2osoi_liq(c,0) = 0._r8 ! kg/m2
frac_iceold(c,0) = 1._r8
end if
! The change of ice partial density of surface node due to precipitation.
! Only ice part of snowfall is added here, the liquid part will be added
! later.
if (snl(c) < 0 .and. newnode == 0) then
h2osoi_ice(c,snl(c)+1) = h2osoi_ice(c,snl(c)+1)+dtime*qflx_snow_grnd_col(c)
dz(c,snl(c)+1) = dz(c,snl(c)+1)+dz_snowf*dtime
end if
end do
! Calculate sublimation and dew, adapted from HydrologyLake and Biogeophysics2.
!dir$ concurrent
!cdir nodep
do fp = 1,num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
jtop = snl(c)+1
! Use column variables here
qflx_evap_grnd(c) = 0._r8
qflx_sub_snow(c) = 0._r8
qflx_dew_snow(c) = 0._r8
qflx_dew_grnd(c) = 0._r8
if (jtop <= 0) then ! snow layers
j = jtop
! Assign ground evaporation to sublimation from soil ice or to dew
! on snow or ground
if (qflx_evap_soi(p) >= 0._r8) then
! for evaporation partitioning between liquid evap and ice sublimation,
! use the ratio of liquid to (liquid+ice) in the top layer to determine split
! Since we're not limiting evap over lakes, but still can't remove more from top
! snow layer than there is there, create temp. limited evap_soi.
qflx_evap_soi_lim = min(qflx_evap_soi(p), (h2osoi_liq(c,j)+h2osoi_ice(c,j))/dtime)
if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
qflx_evap_grnd(c) = max(qflx_evap_soi_lim*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
else
qflx_evap_grnd(c) = 0._r8
end if
qflx_sub_snow(c) = qflx_evap_soi_lim - qflx_evap_grnd(c)
else
if (t_grnd(c) < tfrz) then
qflx_dew_snow(c) = abs(qflx_evap_soi(p))
else
qflx_dew_grnd(c) = abs(qflx_evap_soi(p))
end if
end if
! Update the pft-level qflx_snowcap
! This was moved in from Hydrology2 to keep all pft-level
! calculations out of Hydrology2
if (do_capsnow(c)) qflx_snowcap(p) = qflx_snowcap(p) + qflx_dew_snow(c) + qflx_dew_grnd(c)
else ! No snow layers: do as in HydrologyLake but with actual clmtype variables
if (qflx_evap_soi(p) >= 0._r8) then
! Sublimation: do not allow for more sublimation than there is snow
! after melt. Remaining surface evaporation used for infiltration.
qflx_sub_snow(c) = min(qflx_evap_soi(p), h2osno(c)/dtime)
qflx_evap_grnd(c) = qflx_evap_soi(p) - qflx_sub_snow(c)
else
if (t_grnd(c) < tfrz-0.1_r8) then
qflx_dew_snow(c) = abs(qflx_evap_soi(p))
else
qflx_dew_grnd(c) = abs(qflx_evap_soi(p))
end if
end if
! Update snow pack for dew & sub.
h2osno_temp = h2osno(c)
if (do_capsnow(c)) then
h2osno(c) = h2osno(c) - qflx_sub_snow(c)*dtime
qflx_snowcap(p) = qflx_snowcap(p) + qflx_dew_snow(c) + qflx_dew_grnd(c)
else
h2osno(c) = h2osno(c) + (-qflx_sub_snow(c)+qflx_dew_snow(c))*dtime
end if
if (h2osno_temp > 0._r8) then
snowdp(c) = snowdp(c) * h2osno(c) / h2osno_temp
else
snowdp(c) = h2osno(c)/snow_bd !Assume a constant snow bulk density = 250.
end if
#if (defined PERGRO)
if (abs(h2osno(c)) < 1.e-10_r8) h2osno(c) = 0._r8
#else
h2osno(c) = max(h2osno(c), 0._r8)
#endif
end if
qflx_snowcap_col(c) = qflx_snowcap(p)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Determine initial snow/no-snow filters (will be modified possibly by
! routines CombineSnowLayers and DivideSnowLayers below
call BuildSnowFilter(lbc, ubc, num_shlakec, filter_shlakec,snl, & !i
num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc) !o
! Determine the change of snow mass and the snow water onto soil
call SnowWater(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, & !i
num_shlakenosnowc, filter_shlakenosnowc, & !i
snl,do_capsnow,qflx_snomelt,qflx_rain_grnd, & !i
qflx_sub_snow,qflx_evap_grnd, & !i
qflx_dew_snow,qflx_dew_grnd,dz, & !i
h2osoi_ice,h2osoi_liq, & !i&o
qflx_top_soil) !o
! Determine soil hydrology
! Here this consists only of making sure that soil is saturated even as it melts and 10%
! of pore space opens up. Conversely, if excess ice is melting and the liquid water exceeds the
! saturation value, then remove water.
do j = 1,nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (h2osoi_vol(c,j) < watsat(c,j)) then
h2osoi_liq(c,j) = (watsat(c,j)*dz(c,j) - h2osoi_ice(c,j)/denice)*denh2o
! h2osoi_vol will be updated below, and this water addition will come from qflx_qrgwl
else if (h2osoi_liq(c,j) > watsat(c,j)*denh2o*dz(c,j)) then
h2osoi_liq(c,j) = watsat(c,j)*denh2o*dz(c,j)
end if
end do
end do
!!!!!!!!!!
! if (.not. is_perpetual()) then
if (1==1) then
! Natural compaction and metamorphosis.
call SnowCompaction(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, &!i
snl,imelt,frac_iceold,t_soisno, &!i
h2osoi_ice,h2osoi_liq, &!i
dz) !&o
! Combine thin snow elements
call CombineSnowLayers(lbc, ubc, & !i
num_shlakesnowc, filter_shlakesnowc, & !i&o
snl,h2osno,snowdp,dz,zi, & !i&o
t_soisno,h2osoi_ice,h2osoi_liq, & !i&o
z) !o
! Divide thick snow elements
call DivideSnowLayers(lbc, ubc, & !i
num_shlakesnowc, filter_shlakesnowc, & !i&o
snl,dz,zi,t_soisno, & !i&o
h2osoi_ice,h2osoi_liq, & !i&o
z) !o
else
do fc = 1, num_shlakesnowc
c = filter_shlakesnowc(fc)
h2osno(c) = 0._r8
end do
do j = -nlevsnow+1,0
do fc = 1, num_shlakesnowc
c = filter_shlakesnowc(fc)
if (j >= snl(c)+1) then
h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
end if
end do
end do
end if
! Check for snow layers above lake with unfrozen top layer. Mechanically,
! the snow will fall into the lake and melt or turn to ice. If the top layer has
! sufficient heat to melt the snow without freezing, then that will be done.
! Otherwise, the top layer will undergo freezing, but only if the top layer will
! not freeze completely. Otherwise, let the snow layers persist and melt by diffusion.
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (t_lake(c,1) > tfrz .and. lake_icefrac(c,1) == 0._r8 .and. snl(c) < 0) then
unfrozen(c) = .true.
else
unfrozen(c) = .false.
end if
end do
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (unfrozen(c)) then
if (j == -nlevsnow+1) then
sumsnowice(c) = 0._r8
heatsum(c) = 0._r8
end if
if (j >= snl(c)+1) then
sumsnowice(c) = sumsnowice(c) + h2osoi_ice(c,j)
heatsum(c) = heatsum(c) + h2osoi_ice(c,j)*cpice*(tfrz - t_soisno(c,j)) &
+ h2osoi_liq(c,j)*cpliq*(tfrz - t_soisno(c,j))
end if
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (unfrozen(c)) then
heatsum(c) = heatsum(c) + sumsnowice(c)*hfus
heatrem = (t_lake(c,1) - tfrz)*cpliq*denh2o*dz_lake(c,1) - heatsum(c)
if (heatrem + denh2o*dz_lake(c,1)*hfus > 0._r8) then
! Remove snow and subtract the latent heat from the top layer.
h2osno(c) = 0._r8
snl(c) = 0
! The rest of the bookkeeping for the removed snow will be done below.
#if (defined LAKEDEBUG)
write(message,*)'Snow layers removed above unfrozen lake for column, snowice:', &
c, sumsnowice(c)
CALL wrf_message(message)
#endif
if (heatrem > 0._r8) then ! simply subtract the heat from the layer
t_lake(c,1) = t_lake(c,1) - heatrem/(cpliq*denh2o*dz_lake(c,1))
else !freeze part of the layer
t_lake(c,1) = tfrz
lake_icefrac(c,1) = -heatrem/(denh2o*dz_lake(c,1)*hfus)
end if
end if
end if
end do
!!!!!!!!!!!!
! Set snow age to zero if no snow
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakesnowc
c = filter_shlakesnowc(fc)
if (snl(c) == 0) then
snowage(c) = 0._r8
end if
end do
! Set empty snow layers to zero
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakesnowc
c = filter_shlakesnowc(fc)
if (j <= snl(c) .and. snl(c) > -nlevsnow) then
h2osoi_ice(c,j) = 0._r8
h2osoi_liq(c,j) = 0._r8
t_soisno(c,j) = 0._r8
dz(c,j) = 0._r8
z(c,j) = 0._r8
zi(c,j-1) = 0._r8
end if
end do
end do
! Build new snow filter
call BuildSnowFilter(lbc, ubc, num_shlakec, filter_shlakec, snl,& !i
num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc) !o
! Vertically average t_soisno and sum of h2osoi_liq and h2osoi_ice
! over all snow layers for history output
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakesnowc
c = filter_shlakesnowc(fc)
t_snow(c) = 0._r8
snowice(c) = 0._r8
snowliq(c) = 0._r8
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakenosnowc
c = filter_shlakenosnowc(fc)
t_snow(c) = spval
snowice(c) = spval
snowliq(c) = spval
end do
do j = -nlevsnow+1, 0
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakesnowc
c = filter_shlakesnowc(fc)
if (j >= snl(c)+1) then
t_snow(c) = t_snow(c) + t_soisno(c,j)
snowice(c) = snowice(c) + h2osoi_ice(c,j)
snowliq(c) = snowliq(c) + h2osoi_liq(c,j)
end if
end do
end do
! Determine ending water balance and volumetric soil water
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
if (snl(c) < 0) t_snow(c) = t_snow(c)/abs(snl(c))
endwb(c) = h2osno(c)
end do
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
endwb(c) = endwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice)
end do
end do
#if (defined LAKEDEBUG)
! Check to make sure snow water adds up correctly.
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1, num_shlakec
c = filter_shlakec(fc)
jtop = snl(c)+1
if(j == jtop) snow_water(c) = 0._r8
if(j >= jtop) then
snow_water(c) = snow_water(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
if(j == 0 .and. abs(snow_water(c)-h2osno(c))>1.e-7_r8) then
write(message,*)'h2osno does not equal sum of snow layers in ShalLakeHydrology:', &
'column, h2osno, sum of snow layers =', c, h2osno(c), snow_water(c)
CALL wrf_error_fatal ( message )
end if
end if
end do
end do
#endif
!!!!!!!!!!!!!
! Do history variables and set special landunit runoff (adapted from end of HydrologyLake)
!dir$ concurrent
!cdir nodep
do fp = 1,num_shlakep
p = filter_shlakep(fp)
c = pcolumn(p)
g = pgridcell(p)
qflx_infl(c) = 0._r8
qflx_surf(c) = 0._r8
qflx_drain(c) = 0._r8
rootr_column(c,:) = spval
soilalpha(c) = spval
zwt(c) = spval
fcov(c) = spval
qcharge(c) = spval
! h2osoi_vol(c,:) = spval
! Insure water balance using qflx_qrgwl
qflx_qrgwl(c) = forc_rain(g) + forc_snow(g) - qflx_evap_tot(p) - (endwb(c)-begwb(c))/dtime
#if (defined LAKEDEBUG)
write(message,*)'c, rain, snow, evap, endwb, begwb, qflx_qrgwl:', &
c, forc_rain(g), forc_snow(g), qflx_evap_tot(p), endwb(c), begwb(c), qflx_qrgwl(c)
CALL wrf_message(message)
#endif
! The pft average must be done here for output to history tape
qflx_evap_tot_col(c) = qflx_evap_tot(p)
end do
!!!!!!!!!!!!!
!For now, bracket off the remaining biogeochem code. May need to bring it back
!to do soil carbon and methane beneath lakes.
#if (defined CN)
#ifndef SHLAKE
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do fc = 1, num_soilc
c = filter_soilc(fc)
if (h2osoi_liq(c,j) > 0._r8) then
vwc = h2osoi_liq(c,j)/(dz(c,j)*denh2o)
! the following limit set to catch very small values of
! fractional saturation that can crash the calculation of psi
fsat = max(vwc/vwcsat(c,j), 0.001_r8)
psi = psisat(c,j) * (fsat)**bsw2(c,j)
soilpsi(c,j) = min(max(psi,-15.0_r8),0._r8)
else
soilpsi(c,j) = -15.0_r8
end if
end do
end do
#endif
#endif
#if (defined DGVM) || (defined CN)
#ifndef SHLAKE
! Available soil water up to a depth of 0.5 m.
! Potentially available soil water (=whc) up to a depth of 0.5 m.
! Water content as fraction of whc up to a depth of 0.5 m.
!dir$ concurrent
!cdir nodep
do c = lbc,ubc
l = clandunit(c)
if (ityplun(l) == istsoil) then
rwat(c) = 0._r8
swat(c) = 0._r8
rz(c) = 0._r8
end if
end do
do j = 1, nlevsoil
!dir$ concurrent
!cdir nodep
do c = lbc,ubc
l = clandunit(c)
if (ityplun(l) == istsoil) then
if (z(c,j)+0.5_r8*dz(c,j) <= 0.5_r8) then
watdry = watsat(c,j) * (316230._r8/sucsat(c,j)) ** (-1._r8/bsw(c,j))
rwat(c) = rwat(c) + (h2osoi_vol(c,j)-watdry) * dz(c,j)
swat(c) = swat(c) + (watsat(c,j) -watdry) * dz(c,j)
rz(c) = rz(c) + dz(c,j)
end if
end if
end do
end do
!dir$ concurrent
!cdir nodep
do c = lbc,ubc
l = clandunit(c)
if (ityplun(l) == istsoil) then
if (rz(c) /= 0._r8) then
tsw = rwat(c)/rz(c)
stsw = swat(c)/rz(c)
else
watdry = watsat(c,1) * (316230._r8/sucsat(c,1)) ** (-1._r8/bsw(c,1))
tsw = h2osoi_vol(c,1) - watdry
stsw = watsat(c,1) - watdry
end if
wf(c) = tsw/stsw
else
wf(c) = 1.0_r8
end if
end do
#endif
#endif
end subroutine ShalLakeHydrology
subroutine QSat (T, p, es, esdT, qs, qsdT) 29,3
!
! !DESCRIPTION:
! Computes saturation mixing ratio and the change in saturation
! mixing ratio with respect to temperature.
! Reference: Polynomial approximations from:
! Piotr J. Flatau, et al.,1992: Polynomial fits to saturation
! vapor pressure. Journal of Applied Meteorology, 31, 1507-1513.
!
! !USES:
!
! !ARGUMENTS:
implicit none
real(r8), intent(in) :: T ! temperature (K)
real(r8), intent(in) :: p ! surface atmospheric pressure (pa)
real(r8), intent(out) :: es ! vapor pressure (pa)
real(r8), intent(out) :: esdT ! d(es)/d(T)
real(r8), intent(out) :: qs ! humidity (kg/kg)
real(r8), intent(out) :: qsdT ! d(qs)/d(T)
!
! !CALLED FROM:
! subroutine Biogeophysics1 in module Biogeophysics1Mod
! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod
! subroutine CanopyFluxesMod CanopyFluxesMod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
!
!EOP
!
! !LOCAL VARIABLES:
!
real(r8) :: T_limit
real(r8) :: td,vp,vp1,vp2
!
! For water vapor (temperature range 0C-100C)
!
real(r8), parameter :: a0 = 6.11213476
real(r8), parameter :: a1 = 0.444007856
real(r8), parameter :: a2 = 0.143064234e-01
real(r8), parameter :: a3 = 0.264461437e-03
real(r8), parameter :: a4 = 0.305903558e-05
real(r8), parameter :: a5 = 0.196237241e-07
real(r8), parameter :: a6 = 0.892344772e-10
real(r8), parameter :: a7 = -0.373208410e-12
real(r8), parameter :: a8 = 0.209339997e-15
!
! For derivative:water vapor
!
real(r8), parameter :: b0 = 0.444017302
real(r8), parameter :: b1 = 0.286064092e-01
real(r8), parameter :: b2 = 0.794683137e-03
real(r8), parameter :: b3 = 0.121211669e-04
real(r8), parameter :: b4 = 0.103354611e-06
real(r8), parameter :: b5 = 0.404125005e-09
real(r8), parameter :: b6 = -0.788037859e-12
real(r8), parameter :: b7 = -0.114596802e-13
real(r8), parameter :: b8 = 0.381294516e-16
!
! For ice (temperature range -75C-0C)
!
real(r8), parameter :: c0 = 6.11123516
real(r8), parameter :: c1 = 0.503109514
real(r8), parameter :: c2 = 0.188369801e-01
real(r8), parameter :: c3 = 0.420547422e-03
real(r8), parameter :: c4 = 0.614396778e-05
real(r8), parameter :: c5 = 0.602780717e-07
real(r8), parameter :: c6 = 0.387940929e-09
real(r8), parameter :: c7 = 0.149436277e-11
real(r8), parameter :: c8 = 0.262655803e-14
!
! For derivative:ice
!
real(r8), parameter :: d0 = 0.503277922
real(r8), parameter :: d1 = 0.377289173e-01
real(r8), parameter :: d2 = 0.126801703e-02
real(r8), parameter :: d3 = 0.249468427e-04
real(r8), parameter :: d4 = 0.313703411e-06
real(r8), parameter :: d5 = 0.257180651e-08
real(r8), parameter :: d6 = 0.133268878e-10
real(r8), parameter :: d7 = 0.394116744e-13
real(r8), parameter :: d8 = 0.498070196e-16
!-----------------------------------------------------------------------
T_limit = T - tfrz
if (T_limit > 100.0) T_limit=100.0
if (T_limit < -75.0) T_limit=-75.0
td = T_limit
if (td >= 0.0) then
es = a0 + td*(a1 + td*(a2 + td*(a3 + td*(a4 &
+ td*(a5 + td*(a6 + td*(a7 + td*a8)))))))
esdT = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4 &
+ td*(b5 + td*(b6 + td*(b7 + td*b8)))))))
else
es = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 &
+ td*(c5 + td*(c6 + td*(c7 + td*c8)))))))
esdT = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4 &
+ td*(d5 + td*(d6 + td*(d7 + td*d8)))))))
endif
es = es * 100. ! pa
esdT = esdT * 100. ! pa/K
vp = 1.0 / (p - 0.378*es)
vp1 = 0.622 * vp
vp2 = vp1 * vp
qs = es * vp1 ! kg/kg
qsdT = esdT * vp2 * p ! 1 / K
end subroutine QSat
subroutine Tridiagonal (lbc, ubc, lbj, ubj, jtop, numf, filter, & 4,1
a, b, c, r, u)
!
! !DESCRIPTION:
! Tridiagonal matrix solution
!
! !USES:
! use shr_kind_mod, only: r8 => shr_kind_r8
!
! !ARGUMENTS:
implicit none
integer , intent(in) :: lbc, ubc ! lbinning and ubing column indices
integer , intent(in) :: lbj, ubj ! lbinning and ubing level indices
integer , intent(in) :: jtop(lbc:ubc) ! top level for each column
integer , intent(in) :: numf ! filter dimension
integer , intent(in) :: filter(1:numf) ! filter
real(r8), intent(in) :: a(lbc:ubc, lbj:ubj) ! "a" left off diagonal of tridiagonal matrix
real(r8), intent(in) :: b(lbc:ubc, lbj:ubj) ! "b" diagonal column for tridiagonal matrix
real(r8), intent(in) :: c(lbc:ubc, lbj:ubj) ! "c" right off diagonal tridiagonal matrix
real(r8), intent(in) :: r(lbc:ubc, lbj:ubj) ! "r" forcing term of tridiagonal matrix
real(r8), intent(inout) :: u(lbc:ubc, lbj:ubj) ! solution
!
! !CALLED FROM:
! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod
! subroutine SoilTemperature in module SoilTemperatureMod
! subroutine SoilWater in module HydrologyMod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 1 July 2003: Mariana Vertenstein; modified for vectorization
!
!EOP
!
! !OTHER LOCAL VARIABLES:
!
integer :: j,ci,fc !indices
real(r8) :: gam(lbc:ubc,lbj:ubj) !temporary
real(r8) :: bet(lbc:ubc) !temporary
!-----------------------------------------------------------------------
! Solve the matrix
!dir$ concurrent
!cdir nodep
do fc = 1,numf
ci = filter(fc)
bet(ci) = b(ci,jtop(ci))
end do
do j = lbj, ubj
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1,numf
ci = filter(fc)
if (j >= jtop(ci)) then
if (j == jtop(ci)) then
u(ci,j) = r(ci,j) / bet(ci)
else
gam(ci,j) = c(ci,j-1) / bet(ci)
bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j)
u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci)
end if
end if
end do
end do
!Cray X1 unroll directive used here as work-around for compiler issue 2003/10/20
!dir$ unroll 0
do j = ubj-1,lbj,-1
!dir$ prefervector
!dir$ concurrent
!cdir nodep
do fc = 1,numf
ci = filter(fc)
if (j >= jtop(ci)) then
u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1)
end if
end do
end do
end subroutine Tridiagonal
subroutine SnowWater(lbc, ubc, num_snowc, filter_snowc, & !i 3,9
num_nosnowc, filter_nosnowc, & !i
snl,do_capsnow,qflx_snomelt,qflx_rain_grnd, & !i
qflx_sub_snow,qflx_evap_grnd, & !i
qflx_dew_snow,qflx_dew_grnd,dz, & !i
h2osoi_ice,h2osoi_liq, & !i&o
qflx_top_soil) !o
!===============================================================================
! !DESCRIPTION:
! Evaluate the change of snow mass and the snow water onto soil.
! Water flow within snow is computed by an explicit and non-physical
! based scheme, which permits a part of liquid water over the holding
! capacity (a tentative value is used, i.e. equal to 0.033*porosity) to
! percolate into the underlying layer. Except for cases where the
! porosity of one of the two neighboring layers is less than 0.05, zero
! flow is assumed. The water flow out of the bottom of the snow pack will
! participate as the input of the soil water and runoff. This subroutine
! uses a filter for columns containing snow which must be constructed prior
! to being called.
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 15 November 2000: Mariana Vertenstein
! 2/26/02, Peter Thornton: Migrated to new data structures.
!=============================================================================
! !USES:
! use clmtype
implicit none
!in:
integer, intent(in) :: lbc, ubc ! column bounds
integer, intent(in) :: num_snowc ! number of snow points in column filter
integer, intent(in) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
integer, intent(in) :: num_nosnowc ! number of non-snow points in column filter
integer, intent(in) :: filter_nosnowc(ubc-lbc+1) ! column filter for non-snow points
integer , intent(in) :: snl(1) !number of snow layers
logical , intent(in) :: do_capsnow(1) !true => do snow capping
real(r8), intent(in) :: qflx_snomelt(1) !snow melt (mm H2O /s)
real(r8), intent(in) :: qflx_rain_grnd(1) !rain on ground after interception (mm H2O/s) [+]
real(r8), intent(in) :: qflx_sub_snow(1) !sublimation rate from snow pack (mm H2O /s) [+]
real(r8), intent(in) :: qflx_evap_grnd(1) !ground surface evaporation rate (mm H2O/s) [+]
real(r8), intent(in) :: qflx_dew_snow(1) !surface dew added to snow pack (mm H2O /s) [+]
real(r8), intent(in) :: qflx_dew_grnd(1) !ground surface dew formation (mm H2O /s) [+]
real(r8), intent(in) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
!inout:
real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
!out:
real(r8), intent(out) :: qflx_top_soil(1) !net water input into soil from top (mm/s)
! OTHER LOCAL VARIABLES:
integer :: c, j, fc !do loop/array indices
real(r8) :: qin(lbc:ubc) !water flow into the elmement (mm/s)
real(r8) :: qout(lbc:ubc) !water flow out of the elmement (mm/s)
real(r8) :: wgdif !ice mass after minus sublimation
real(r8) :: vol_liq(lbc:ubc,-nlevsnow+1:0) !partial volume of liquid water in layer
real(r8) :: vol_ice(lbc:ubc,-nlevsnow+1:0) !partial volume of ice lens in layer
real(r8) :: eff_porosity(lbc:ubc,-nlevsnow+1:0) !effective porosity = porosity - vol_ice
!-----------------------------------------------------------------------
! Renew the mass of ice lens (h2osoi_ice) and liquid (h2osoi_liq) in the
! surface snow layer resulting from sublimation (frost) / evaporation (condense)
!dir$ concurrent
!cdir nodep
do fc = 1,num_snowc
c = filter_snowc(fc)
if (do_capsnow(c)) then
wgdif = h2osoi_ice(c,snl(c)+1) - qflx_sub_snow(c)*dtime
h2osoi_ice(c,snl(c)+1) = wgdif
if (wgdif < 0.) then
h2osoi_ice(c,snl(c)+1) = 0.
h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
end if
h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) - qflx_evap_grnd(c) * dtime
else
wgdif = h2osoi_ice(c,snl(c)+1) + (qflx_dew_snow(c) - qflx_sub_snow(c)) * dtime
h2osoi_ice(c,snl(c)+1) = wgdif
if (wgdif < 0.) then
h2osoi_ice(c,snl(c)+1) = 0.
h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
end if
h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + &
(qflx_rain_grnd(c) + qflx_dew_grnd(c) - qflx_evap_grnd(c)) * dtime
end if
h2osoi_liq(c,snl(c)+1) = max(0._r8, h2osoi_liq(c,snl(c)+1))
end do
! Porosity and partial volume
do j = -nlevsnow+1, 0
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
vol_ice(c,j) = min(1._r8, h2osoi_ice(c,j)/(dz(c,j)*denice))
eff_porosity(c,j) = 1. - vol_ice(c,j)
vol_liq(c,j) = min(eff_porosity(c,j),h2osoi_liq(c,j)/(dz(c,j)*denh2o))
end if
end do
end do
! Capillary forces within snow are usually two or more orders of magnitude
! less than those of gravity. Only gravity terms are considered.
! the genernal expression for water flow is "K * ss**3", however,
! no effective parameterization for "K". Thus, a very simple consideration
! (not physically based) is introduced:
! when the liquid water of layer exceeds the layer's holding
! capacity, the excess meltwater adds to the underlying neighbor layer.
qin(:) = 0._r8
do j = -nlevsnow+1, 0
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
h2osoi_liq(c,j) = h2osoi_liq(c,j) + qin(c)
if (j <= -1) then
! No runoff over snow surface, just ponding on surface
if (eff_porosity(c,j) < wimp .OR. eff_porosity(c,j+1) < wimp) then
qout(c) = 0._r8
else
qout(c) = max(0._r8,(vol_liq(c,j)-ssi*eff_porosity(c,j))*dz(c,j))
qout(c) = min(qout(c),(1.-vol_ice(c,j+1)-vol_liq(c,j+1))*dz(c,j+1))
end if
else
qout(c) = max(0._r8,(vol_liq(c,j) - ssi*eff_porosity(c,j))*dz(c,j))
end if
qout(c) = qout(c)*1000.
h2osoi_liq(c,j) = h2osoi_liq(c,j) - qout(c)
qin(c) = qout(c)
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
! Qout from snow bottom
qflx_top_soil(c) = qout(c) / dtime
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_nosnowc
c = filter_nosnowc(fc)
qflx_top_soil(c) = qflx_rain_grnd(c) + qflx_snomelt(c)
end do
end subroutine SnowWater
subroutine SnowCompaction(lbc, ubc, num_snowc, filter_snowc, &!i 2,3
snl,imelt,frac_iceold,t_soisno, &!i
h2osoi_ice,h2osoi_liq, &!i
dz) !i&o
!================================================================================
! !DESCRIPTION:
! Determine the change in snow layer thickness due to compaction and
! settling.
! Three metamorphisms of changing snow characteristics are implemented,
! i.e., destructive, overburden, and melt. The treatments of the former
! two are from SNTHERM.89 and SNTHERM.99 (1991, 1999). The contribution
! due to melt metamorphism is simply taken as a ratio of snow ice
! fraction after the melting versus before the melting.
!
! CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/28/02, Peter Thornton: Migrated to new data structures
!==============================================================================
! USES:
! use clmtype
!
! !ARGUMENTS:
implicit none
!in:
integer, intent(in) :: lbc, ubc ! column bounds
integer, intent(in) :: num_snowc ! number of column snow points in column filter
integer, intent(in) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
integer, intent(in) :: snl(1) !number of snow layers
integer, intent(in) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0
real(r8), intent(in) :: frac_iceold(1,-nlevsnow+1:nlevsoil) !fraction of ice relative to the tot water
real(r8), intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
real(r8), intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
real(r8), intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
!inout:
real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
! OTHER LOCAL VARIABLES:
integer :: j, c, fc ! indices
real(r8), parameter :: c2 = 23.e-3 ! [m3/kg]
real(r8), parameter :: c3 = 2.777e-6 ! [1/s]
real(r8), parameter :: c4 = 0.04 ! [1/K]
real(r8), parameter :: c5 = 2.0 !
real(r8), parameter :: dm = 100.0 ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
real(r8), parameter :: eta0 = 9.e+5 ! The Viscosity Coefficient Eta0 [kg-s/m2]
real(r8) :: burden(lbc:ubc) ! pressure of overlying snow [kg/m2]
real(r8) :: ddz1 ! Rate of settling of snowpack due to destructive metamorphism.
real(r8) :: ddz2 ! Rate of compaction of snowpack due to overburden.
real(r8) :: ddz3 ! Rate of compaction of snowpack due to melt [1/s]
real(r8) :: dexpf ! expf=exp(-c4*(273.15-t_soisno)).
real(r8) :: fi ! Fraction of ice relative to the total water content at current time step
real(r8) :: td ! t_soisno - tfrz [K]
real(r8) :: pdzdtc ! Nodal rate of change in fractional-thickness due to compaction [fraction/s]
real(r8) :: void ! void (1 - vol_ice - vol_liq)
real(r8) :: wx ! water mass (ice+liquid) [kg/m2]
real(r8) :: bi ! partial density of ice [kg/m3]
!-----------------------------------------------------------------------
! Begin calculation - note that the following column loops are only invoked if snl(c) < 0
burden(:) = 0._r8
do j = -nlevsnow+1, 0
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
wx = h2osoi_ice(c,j) + h2osoi_liq(c,j)
void = 1. - (h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o) / dz(c,j)
! Allow compaction only for non-saturated node and higher ice lens node.
if (void > 0.001 .and. h2osoi_ice(c,j) > .1) then
bi = h2osoi_ice(c,j) / dz(c,j)
fi = h2osoi_ice(c,j) / wx
td = tfrz-t_soisno(c,j)
dexpf = exp(-c4*td)
! Settling as a result of destructive metamorphism
ddz1 = -c3*dexpf
if (bi > dm) ddz1 = ddz1*exp(-46.0e-3*(bi-dm))
! Liquid water term
if (h2osoi_liq(c,j) > 0.01*dz(c,j)) ddz1=ddz1*c5
! Compaction due to overburden
ddz2 = -burden(c)*exp(-0.08*td - c2*bi)/eta0
! Compaction occurring during melt
if (imelt(c,j) == 1) then
ddz3 = - 1./dtime * max(0._r8,(frac_iceold(c,j) - fi)/frac_iceold(c,j))
else
ddz3 = 0._r8
end if
! Time rate of fractional change in dz (units of s-1)
pdzdtc = ddz1 + ddz2 + ddz3
! The change in dz due to compaction
dz(c,j) = dz(c,j) * (1.+pdzdtc*dtime)
end if
! Pressure of overlying snow
burden(c) = burden(c) + wx
end if
end do
end do
end subroutine SnowCompaction
subroutine CombineSnowLayers(lbc, ubc, & !i 2,5
num_snowc, filter_snowc, & !i&o
snl,h2osno,snowdp,dz,zi, & !i&o
t_soisno,h2osoi_ice,h2osoi_liq, & !i&o
z) !o
!==========================================================================
! !DESCRIPTION:
! Combine snow layers that are less than a minimum thickness or mass
! If the snow element thickness or mass is less than a prescribed minimum,
! then it is combined with a neighboring element. The subroutine
! clm\_combo.f90 then executes the combination of mass and energy.
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/28/02, Peter Thornton: Migrated to new data structures.
!=========================================================================
! !USES:
! use clmtype
!
! !ARGUMENTS:
implicit none
!in:
integer, intent(in) :: lbc, ubc ! column bounds
! integer, intent(in) :: clandunit(1) !landunit index for each column
! integer, intent(in) :: ityplun(1) !landunit type
!inout:
integer, intent(inout) :: num_snowc ! number of column snow points in column filter
integer, intent(inout) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
integer , intent(inout) :: snl(1) !number of snow layers
real(r8), intent(inout) :: h2osno(1) !snow water (mm H2O)
real(r8), intent(inout) :: snowdp(1) !snow height (m)
real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
real(r8), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) !interface level below a "z" level (m)
real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
!out:
real(r8), intent(out) :: z(1,-nlevsnow+1:nlevsoil) !layer thickness (m)
!
!EOP
!
! !OTHER LOCAL VARIABLES:
!
integer :: c, fc ! column indices
integer :: i,k ! loop indices
integer :: j,l ! node indices
integer :: msn_old(lbc:ubc) ! number of top snow layer
integer :: mssi(lbc:ubc) ! node index
integer :: neibor ! adjacent node selected for combination
real(r8):: zwice(lbc:ubc) ! total ice mass in snow
real(r8):: zwliq (lbc:ubc) ! total liquid water in snow
real(r8):: dzmin(5) ! minimum of top snow layer
data dzmin /0.010, 0.015, 0.025, 0.055, 0.115/
!-----------------------------------------------------------------------
! Check the mass of ice lens of snow, when the total is less than a small value,
! combine it with the underlying neighbor.
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
msn_old(c) = snl(c)
end do
! The following loop is NOT VECTORIZED
do fc = 1, num_snowc
c = filter_snowc(fc)
! l = clandunit(c)
do j = msn_old(c)+1,0
if (h2osoi_ice(c,j) <= .1) then
! if (ityplun(l) == istsoil) then
! h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
! h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
! else if (ityplun(l) /= istsoil .and. j /= 0) then
h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
! end if
! shift all elements above this down one.
if (j > snl(c)+1 .and. snl(c) < -1) then
do i = j, snl(c)+2, -1
t_soisno(c,i) = t_soisno(c,i-1)
h2osoi_liq(c,i) = h2osoi_liq(c,i-1)
h2osoi_ice(c,i) = h2osoi_ice(c,i-1)
dz(c,i) = dz(c,i-1)
end do
end if
snl(c) = snl(c) + 1
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
h2osno(c) = 0._r8
snowdp(c) = 0._r8
zwice(c) = 0._r8
zwliq(c) = 0._r8
end do
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
snowdp(c) = snowdp(c) + dz(c,j)
zwice(c) = zwice(c) + h2osoi_ice(c,j)
zwliq(c) = zwliq(c) + h2osoi_liq(c,j)
end if
end do
end do
! Check the snow depth - all snow gone
! The liquid water assumes ponding on soil surface.
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
! l = clandunit(c)
if (snowdp(c) < 0.01 .and. snowdp(c) > 0.) then
snl(c) = 0
h2osno(c) = zwice(c)
if (h2osno(c) <= 0.) snowdp(c) = 0._r8
! if (ityplun(l) == istsoil) h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) !change by guhp
end if
end do
! Check the snow depth - snow layers combined
! The following loop IS NOT VECTORIZED
do fc = 1, num_snowc
c = filter_snowc(fc)
! Two or more layers
if (snl(c) < -1) then
msn_old(c) = snl(c)
mssi(c) = 1
do i = msn_old(c)+1,0
if (dz(c,i) < dzmin(mssi(c))) then
if (i == snl(c)+1) then
! If top node is removed, combine with bottom neighbor.
neibor = i + 1
else if (i == 0) then
! If the bottom neighbor is not snow, combine with the top neighbor.
neibor = i - 1
else
! If none of the above special cases apply, combine with the thinnest neighbor
neibor = i + 1
if ((dz(c,i-1)+dz(c,i)) < (dz(c,i+1)+dz(c,i))) neibor = i-1
end if
! Node l and j are combined and stored as node j.
if (neibor > i) then
j = neibor
l = i
else
j = i
l = neibor
end if
call Combo (dz(c,j), h2osoi_liq(c,j), h2osoi_ice(c,j), &
t_soisno(c,j), dz(c,l), h2osoi_liq(c,l), h2osoi_ice(c,l), t_soisno(c,l) )
! Now shift all elements above this down one.
if (j-1 > snl(c)+1) then
do k = j-1, snl(c)+2, -1
t_soisno(c,k) = t_soisno(c,k-1)
h2osoi_ice(c,k) = h2osoi_ice(c,k-1)
h2osoi_liq(c,k) = h2osoi_liq(c,k-1)
dz(c,k) = dz(c,k-1)
end do
end if
! Decrease the number of snow layers
snl(c) = snl(c) + 1
if (snl(c) >= -1) EXIT
else
! The layer thickness is greater than the prescribed minimum value
mssi(c) = mssi(c) + 1
end if
end do
end if
end do
! Reset the node depth and the depth of layer interface
do j = 0, -nlevsnow+1, -1
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c) + 1) then
z(c,j) = zi(c,j) - 0.5*dz(c,j)
zi(c,j-1) = zi(c,j) - dz(c,j)
end if
end do
end do
end subroutine CombineSnowLayers
subroutine DivideSnowLayers(lbc, ubc, & !i 2,10
num_snowc, filter_snowc, & !i&o
snl,dz,zi,t_soisno, & !i&o
h2osoi_ice,h2osoi_liq, & !i&o
z) !o
!============================================================================
! !DESCRIPTION:
! Subdivides snow layers if they exceed their prescribed maximum thickness.
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/28/02, Peter Thornton: Migrated to new data structures.
!============================================================================
! !USES:
! use clmtype
!
! !ARGUMENTS:
implicit none
!in:
integer, intent(in) :: lbc, ubc ! column bounds
!inout:
integer, intent(inout) :: num_snowc ! number of column snow points in column filter
integer, intent(inout) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
integer , intent(inout) :: snl(1) !number of snow layers
real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
real(r8), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) !interface level below a "z" level (m)
real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
!out:
real(r8), intent(out) :: z(1,-nlevsnow+1:nlevsoil) !layer thickness (m)
! OTHER LOCAL VARIABLES:
integer :: j, c, fc ! indices
real(r8) :: drr ! thickness of the combined [m]
integer :: msno ! number of snow layer 1 (top) to msno (bottom)
real(r8) :: dzsno(lbc:ubc,nlevsnow) ! Snow layer thickness [m]
real(r8) :: swice(lbc:ubc,nlevsnow) ! Partial volume of ice [m3/m3]
real(r8) :: swliq(lbc:ubc,nlevsnow) ! Partial volume of liquid water [m3/m3]
real(r8) :: tsno(lbc:ubc ,nlevsnow) ! Nodel temperature [K]
real(r8) :: zwice ! temporary
real(r8) :: zwliq ! temporary
real(r8) :: propor ! temporary
!-----------------------------------------------------------------------
! Begin calculation - note that the following column loops are only invoked
! for snow-covered columns
do j = 1,nlevsnow
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j <= abs(snl(c))) then
dzsno(c,j) = dz(c,j+snl(c))
swice(c,j) = h2osoi_ice(c,j+snl(c))
swliq(c,j) = h2osoi_liq(c,j+snl(c))
tsno(c,j) = t_soisno(c,j+snl(c))
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
msno = abs(snl(c))
if (msno == 1) then
! Specify a new snow layer
if (dzsno(c,1) > 0.03) then
msno = 2
dzsno(c,1) = dzsno(c,1)/2.
swice(c,1) = swice(c,1)/2.
swliq(c,1) = swliq(c,1)/2.
dzsno(c,2) = dzsno(c,1)
swice(c,2) = swice(c,1)
swliq(c,2) = swliq(c,1)
tsno(c,2) = tsno(c,1)
end if
end if
if (msno > 1) then
if (dzsno(c,1) > 0.02) then
drr = dzsno(c,1) - 0.02
propor = drr/dzsno(c,1)
zwice = propor*swice(c,1)
zwliq = propor*swliq(c,1)
propor = 0.02/dzsno(c,1)
swice(c,1) = propor*swice(c,1)
swliq(c,1) = propor*swliq(c,1)
dzsno(c,1) = 0.02
call Combo (dzsno(c,2), swliq(c,2), swice(c,2), tsno(c,2), drr, &
zwliq, zwice, tsno(c,1))
! Subdivide a new layer
if (msno <= 2 .and. dzsno(c,2) > 0.07) then
msno = 3
dzsno(c,2) = dzsno(c,2)/2.
swice(c,2) = swice(c,2)/2.
swliq(c,2) = swliq(c,2)/2.
dzsno(c,3) = dzsno(c,2)
swice(c,3) = swice(c,2)
swliq(c,3) = swliq(c,2)
tsno(c,3) = tsno(c,2)
end if
end if
end if
if (msno > 2) then
if (dzsno(c,2) > 0.05) then
drr = dzsno(c,2) - 0.05
propor = drr/dzsno(c,2)
zwice = propor*swice(c,2)
zwliq = propor*swliq(c,2)
propor = 0.05/dzsno(c,2)
swice(c,2) = propor*swice(c,2)
swliq(c,2) = propor*swliq(c,2)
dzsno(c,2) = 0.05
call Combo (dzsno(c,3), swliq(c,3), swice(c,3), tsno(c,3), drr, &
zwliq, zwice, tsno(c,2))
! Subdivided a new layer
if (msno <= 3 .and. dzsno(c,3) > 0.18) then
msno = 4
dzsno(c,3) = dzsno(c,3)/2.
swice(c,3) = swice(c,3)/2.
swliq(c,3) = swliq(c,3)/2.
dzsno(c,4) = dzsno(c,3)
swice(c,4) = swice(c,3)
swliq(c,4) = swliq(c,3)
tsno(c,4) = tsno(c,3)
end if
end if
end if
if (msno > 3) then
if (dzsno(c,3) > 0.11) then
drr = dzsno(c,3) - 0.11
propor = drr/dzsno(c,3)
zwice = propor*swice(c,3)
zwliq = propor*swliq(c,3)
propor = 0.11/dzsno(c,3)
swice(c,3) = propor*swice(c,3)
swliq(c,3) = propor*swliq(c,3)
dzsno(c,3) = 0.11
call Combo (dzsno(c,4), swliq(c,4), swice(c,4), tsno(c,4), drr, &
zwliq, zwice, tsno(c,3))
! Subdivided a new layer
if (msno <= 4 .and. dzsno(c,4) > 0.41) then
msno = 5
dzsno(c,4) = dzsno(c,4)/2.
swice(c,4) = swice(c,4)/2.
swliq(c,4) = swliq(c,4)/2.
dzsno(c,5) = dzsno(c,4)
swice(c,5) = swice(c,4)
swliq(c,5) = swliq(c,4)
tsno(c,5) = tsno(c,4)
end if
end if
end if
if (msno > 4) then
if (dzsno(c,4) > 0.23) then
drr = dzsno(c,4) - 0.23
propor = drr/dzsno(c,4)
zwice = propor*swice(c,4)
zwliq = propor*swliq(c,4)
propor = 0.23/dzsno(c,4)
swice(c,4) = propor*swice(c,4)
swliq(c,4) = propor*swliq(c,4)
dzsno(c,4) = 0.23
call Combo (dzsno(c,5), swliq(c,5), swice(c,5), tsno(c,5), drr, &
zwliq, zwice, tsno(c,4))
end if
end if
snl(c) = -msno
end do
do j = -nlevsnow+1,0
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
dz(c,j) = dzsno(c,j-snl(c))
h2osoi_ice(c,j) = swice(c,j-snl(c))
h2osoi_liq(c,j) = swliq(c,j-snl(c))
t_soisno(c,j) = tsno(c,j-snl(c))
end if
end do
end do
do j = 0, -nlevsnow+1, -1
!dir$ concurrent
!cdir nodep
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
z(c,j) = zi(c,j) - 0.5*dz(c,j)
zi(c,j-1) = zi(c,j) - dz(c,j)
end if
end do
end do
end subroutine DivideSnowLayers
subroutine Combo(dz, wliq, wice, t, dz2, wliq2, wice2, t2) 17,1
!
! !DESCRIPTION:
! Combines two elements and returns the following combined
! variables: dz, t, wliq, wice.
! The combined temperature is based on the equation:
! the sum of the enthalpies of the two elements =
! that of the combined element.
!
! !USES:
!
! !ARGUMENTS:
implicit none
real(r8), intent(in) :: dz2 ! nodal thickness of 2 elements being combined [m]
real(r8), intent(in) :: wliq2 ! liquid water of element 2 [kg/m2]
real(r8), intent(in) :: wice2 ! ice of element 2 [kg/m2]
real(r8), intent(in) :: t2 ! nodal temperature of element 2 [K]
real(r8), intent(inout) :: dz ! nodal thickness of 1 elements being combined [m]
real(r8), intent(inout) :: wliq ! liquid water of element 1
real(r8), intent(inout) :: wice ! ice of element 1 [kg/m2]
real(r8), intent(inout) :: t ! nodel temperature of elment 1 [K]
!
! !CALLED FROM:
! subroutine CombineSnowLayers in this module
! subroutine DivideSnowLayers in this module
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
!
!EOP
!
! !LOCAL VARIABLES:
!
real(r8) :: dzc ! Total thickness of nodes 1 and 2 (dzc=dz+dz2).
real(r8) :: wliqc ! Combined liquid water [kg/m2]
real(r8) :: wicec ! Combined ice [kg/m2]
real(r8) :: tc ! Combined node temperature [K]
real(r8) :: h ! enthalpy of element 1 [J/m2]
real(r8) :: h2 ! enthalpy of element 2 [J/m2]
real(r8) :: hc ! temporary
!-----------------------------------------------------------------------
dzc = dz+dz2
wicec = (wice+wice2)
wliqc = (wliq+wliq2)
h = (cpice*wice+cpliq*wliq) * (t-tfrz)+hfus*wliq
h2= (cpice*wice2+cpliq*wliq2) * (t2-tfrz)+hfus*wliq2
hc = h + h2
if(hc < 0.)then
tc = tfrz + hc/(cpice*wicec + cpliq*wliqc)
else if (hc.le.hfus*wliqc) then
tc = tfrz
else
tc = tfrz + (hc - hfus*wliqc) / (cpice*wicec + cpliq*wliqc)
end if
dz = dzc
wice = wicec
wliq = wliqc
t = tc
end subroutine Combo
subroutine BuildSnowFilter(lbc, ubc, num_nolakec, filter_nolakec,snl, & !i 4,1
num_snowc, filter_snowc, & !o
num_nosnowc, filter_nosnowc) !o
!
! !DESCRIPTION:
! Constructs snow filter for use in vectorized loops for snow hydrology.
!
! !USES:
! use clmtype
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbc, ubc ! column bounds
integer, intent(in) :: num_nolakec ! number of column non-lake points in column filter
integer, intent(in) :: filter_nolakec(ubc-lbc+1) ! column filter for non-lake points
integer, intent(in) :: snl(1) ! number of snow layers
integer, intent(out) :: num_snowc ! number of column snow points in column filter
integer, intent(out) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
integer, intent(out) :: num_nosnowc ! number of column non-snow points in column filter
integer, intent(out) :: filter_nosnowc(ubc-lbc+1) ! column filter for non-snow points
!
! !CALLED FROM:
! subroutine Hydrology2 in Hydrology2Mod
! subroutine CombineSnowLayers in this module
!
! !REVISION HISTORY:
! 2003 July 31: Forrest Hoffman
!
! !LOCAL VARIABLES:
! local pointers to implicit in arguments
!
!EOP
!
! !OTHER LOCAL VARIABLES:
integer :: fc, c
!-----------------------------------------------------------------------
! Build snow/no-snow filters for other subroutines
num_snowc = 0
num_nosnowc = 0
do fc = 1, num_nolakec
c = filter_nolakec(fc)
if (snl(c) < 0) then
num_snowc = num_snowc + 1
filter_snowc(num_snowc) = c
else
num_nosnowc = num_nosnowc + 1
filter_nosnowc(num_nosnowc) = c
end if
end do
end subroutine BuildSnowFilter
subroutine FrictionVelocity(pgridcell,forc_hgt,forc_hgt_u, & !i 4,2
forc_hgt_t,forc_hgt_q, & !i
lbp, ubp, fn, filterp, & !i
displa, z0m, z0h, z0q, & !i
obu, iter, ur, um, & !i
ustar,temp1, temp2, temp12m, temp22m, & !o
u10,fv, & !o
fm) !i&o
!=============================================================================
! !DESCRIPTION:
! Calculation of the friction velocity, relation for potential
! temperature and humidity profiles of surface boundary layer.
! The scheme is based on the work of Zeng et al. (1998):
! Intercomparison of bulk aerodynamic algorithms for the computation
! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
! Vol. 11, 2628-2644.
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 12/19/01, Peter Thornton
! Added arguments to eliminate passing clm derived type into this function.
! Created by Mariana Vertenstein
!============================================================================
! !USES:
! use clmtype
!!use clm_atmlnd, only : clm_a2l
!
! !ARGUMENTS:
implicit none
!in:
integer , intent(in) :: pgridcell(1) ! pft's gridcell index
real(r8), intent(in) :: forc_hgt(1) ! atmospheric reference height (m)
real(r8), intent(in) :: forc_hgt_u(1) ! observational height of wind [m]
real(r8), intent(in) :: forc_hgt_t(1) ! observational height of temperature [m]
real(r8), intent(in) :: forc_hgt_q(1) ! observational height of humidity [m]
integer , intent(in) :: lbp, ubp ! pft array bounds
integer , intent(in) :: fn ! number of filtered pft elements
integer , intent(in) :: filterp(fn) ! pft filter
real(r8), intent(in) :: displa(lbp:ubp) ! displacement height (m)
real(r8), intent(in) :: z0m(lbp:ubp) ! roughness length over vegetation, momentum [m]
real(r8), intent(in) :: z0h(lbp:ubp) ! roughness length over vegetation, sensible heat [m]
real(r8), intent(in) :: z0q(lbp:ubp) ! roughness length over vegetation, latent heat [m]
real(r8), intent(in) :: obu(lbp:ubp) ! monin-obukhov length (m)
integer, intent(in) :: iter ! iteration number
real(r8), intent(in) :: ur(lbp:ubp) ! wind speed at reference height [m/s]
real(r8), intent(in) :: um(lbp:ubp) ! wind speed including the stablity effect [m/s]
!out:
real(r8), intent(out) :: ustar(lbp:ubp) ! friction velocity [m/s]
real(r8), intent(out) :: temp1(lbp:ubp) ! relation for potential temperature profile
real(r8), intent(out) :: temp12m(lbp:ubp) ! relation for potential temperature profile applied at 2-m
real(r8), intent(out) :: temp2(lbp:ubp) ! relation for specific humidity profile
real(r8), intent(out) :: temp22m(lbp:ubp) ! relation for specific humidity profile applied at 2-m
real(r8), intent(out) :: u10(1) ! 10-m wind (m/s) (for dust model)
real(r8), intent(out) :: fv(1) ! friction velocity (m/s) (for dust model)
!inout:
real(r8), intent(inout) :: fm(lbp:ubp) ! needed for DGVM only to diagnose 10m wind
! OTHER LOCAL VARIABLES:
real(r8), parameter :: zetam = 1.574_r8 ! transition point of flux-gradient relation (wind profile)
real(r8), parameter :: zetat = 0.465_r8 ! transition point of flux-gradient relation (temp. profile)
integer :: f ! pft-filter index
integer :: p ! pft index
integer :: g ! gridcell index
real(r8):: zldis(lbp:ubp) ! reference height "minus" zero displacement heght [m]
real(r8):: zeta(lbp:ubp) ! dimensionless height used in Monin-Obukhov theory
#if (defined DGVM) || (defined DUST)
real(r8) :: tmp1,tmp2,tmp3,tmp4 ! Used to diagnose the 10 meter wind
real(r8) :: fmnew ! Used to diagnose the 10 meter wind
real(r8) :: fm10 ! Used to diagnose the 10 meter wind
real(r8) :: zeta10 ! Used to diagnose the 10 meter wind
#endif
!------------------------------------------------------------------------------
! Adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions.
#if (!defined PERGRO)
!dir$ concurrent
!cdir nodep
do f = 1, fn
p = filterp(f)
g = pgridcell(p)
! Wind profile
zldis(p) = forc_hgt_u(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetam) then
ustar(p) = vkc*um(p)/(log(-zetam*obu(p)/z0m(p))&
- StabilityFunc1(-zetam) &
+ StabilityFunc1(z0m(p)/obu(p)) &
+ 1.14_r8*((-zeta(p))**0.333_r8-(zetam)**0.333_r8))
else if (zeta(p) < 0._r8) then
ustar(p) = vkc*um(p)/(log(zldis(p)/z0m(p))&
- StabilityFunc1(zeta(p))&
+ StabilityFunc1(z0m(p)/obu(p)))
else if (zeta(p) <= 1._r8) then
ustar(p) = vkc*um(p)/(log(zldis(p)/z0m(p)) + 5._r8*zeta(p) -5._r8*z0m(p)/obu(p))
else
ustar(p) = vkc*um(p)/(log(obu(p)/z0m(p))+5._r8-5._r8*z0m(p)/obu(p) &
+(5._r8*log(zeta(p))+zeta(p)-1._r8))
end if
! Temperature profile
zldis(p) = forc_hgt_t(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp1(p) = vkc/(log(-zetat*obu(p)/z0h(p))&
- StabilityFunc2(-zetat) &
+ StabilityFunc2(z0h(p)/obu(p)) &
+ 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
else if (zeta(p) < 0._r8) then
temp1(p) = vkc/(log(zldis(p)/z0h(p)) &
- StabilityFunc2(zeta(p)) &
+ StabilityFunc2(z0h(p)/obu(p)))
else if (zeta(p) <= 1._r8) then
temp1(p) = vkc/(log(zldis(p)/z0h(p)) + 5._r8*zeta(p) - 5._r8*z0h(p)/obu(p))
else
temp1(p) = vkc/(log(obu(p)/z0h(p)) + 5._r8 - 5._r8*z0h(p)/obu(p) &
+ (5._r8*log(zeta(p))+zeta(p)-1._r8))
end if
! Humidity profile
if (forc_hgt_q(g) == forc_hgt_t(g) .and. z0q(p) == z0h(p)) then
temp2(p) = temp1(p)
else
zldis(p) = forc_hgt_q(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp2(p) = vkc/(log(-zetat*obu(p)/z0q(p)) &
- StabilityFunc2(-zetat) &
+ StabilityFunc2(z0q(p)/obu(p)) &
+ 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
else if (zeta(p) < 0._r8) then
temp2(p) = vkc/(log(zldis(p)/z0q(p)) &
- StabilityFunc2(zeta(p)) &
+ StabilityFunc2(z0q(p)/obu(p)))
else if (zeta(p) <= 1._r8) then
temp2(p) = vkc/(log(zldis(p)/z0q(p)) + 5._r8*zeta(p)-5._r8*z0q(p)/obu(p))
else
temp2(p) = vkc/(log(obu(p)/z0q(p)) + 5._r8 - 5._r8*z0q(p)/obu(p) &
+ (5._r8*log(zeta(p))+zeta(p)-1._r8))
end if
endif
! Temperature profile applied at 2-m
zldis(p) = 2.0_r8 + z0h(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp12m(p) = vkc/(log(-zetat*obu(p)/z0h(p))&
- StabilityFunc2(-zetat) &
+ StabilityFunc2(z0h(p)/obu(p)) &
+ 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
else if (zeta(p) < 0._r8) then
temp12m(p) = vkc/(log(zldis(p)/z0h(p)) &
- StabilityFunc2(zeta(p)) &
+ StabilityFunc2(z0h(p)/obu(p)))
else if (zeta(p) <= 1._r8) then
temp12m(p) = vkc/(log(zldis(p)/z0h(p)) + 5._r8*zeta(p) - 5._r8*z0h(p)/obu(p))
else
temp12m(p) = vkc/(log(obu(p)/z0h(p)) + 5._r8 - 5._r8*z0h(p)/obu(p) &
+ (5._r8*log(zeta(p))+zeta(p)-1._r8))
end if
! Humidity profile applied at 2-m
if (z0q(p) == z0h(p)) then
temp22m(p) = temp12m(p)
else
zldis(p) = 2.0_r8 + z0q(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp22m(p) = vkc/(log(-zetat*obu(p)/z0q(p)) - &
StabilityFunc2(-zetat) + StabilityFunc2(z0q(p)/obu(p)) &
+ 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
else if (zeta(p) < 0._r8) then
temp22m(p) = vkc/(log(zldis(p)/z0q(p)) - &
StabilityFunc2(zeta(p))+StabilityFunc2(z0q(p)/obu(p)))
else if (zeta(p) <= 1._r8) then
temp22m(p) = vkc/(log(zldis(p)/z0q(p)) + 5._r8*zeta(p)-5._r8*z0q(p)/obu(p))
else
temp22m(p) = vkc/(log(obu(p)/z0q(p)) + 5._r8 - 5._r8*z0q(p)/obu(p) &
+ (5._r8*log(zeta(p))+zeta(p)-1._r8))
end if
end if
#if (defined DGVM) || (defined DUST)
! diagnose 10-m wind for dust model (dstmbl.F)
! Notes from C. Zender's dst.F:
! According to Bon96 p. 62, the displacement height d (here displa) is
! 0.0 <= d <= 0.34 m in dust source regions (i.e., regions w/o trees).
! Therefore d <= 0.034*z1 and may safely be neglected.
! Code from LSM routine SurfaceTemperature was used to obtain u10
zldis(p) = forc_hgt_u(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (min(zeta(p), 1._r8) < 0._r8) then
tmp1 = (1._r8 - 16._r8*min(zeta(p),1._r8))**0.25_r8
tmp2 = log((1._r8+tmp1*tmp1)/2._r8)
tmp3 = log((1._r8+tmp1)/2._r8)
fmnew = 2._r8*tmp3 + tmp2 - 2._r8*atan(tmp1) + 1.5707963_r8
else
fmnew = -5._r8*min(zeta(p),1._r8)
endif
if (iter == 1) then
fm(p) = fmnew
else
fm(p) = 0.5_r8 * (fm(p)+fmnew)
end if
zeta10 = min(10._r8/obu(p), 1._r8)
if (zeta(p) == 0._r8) zeta10 = 0._r8
if (zeta10 < 0._r8) then
tmp1 = (1.0_r8 - 16.0_r8 * zeta10)**0.25_r8
tmp2 = log((1.0_r8 + tmp1*tmp1)/2.0_r8)
tmp3 = log((1.0_r8 + tmp1)/2.0_r8)
fm10 = 2.0_r8*tmp3 + tmp2 - 2.0_r8*atan(tmp1) + 1.5707963_r8
else ! not stable
fm10 = -5.0_r8 * zeta10
end if
tmp4 = log(forc_hgt(g) / 10._r8)
u10(p) = ur(p) - ustar(p)/vkc * (tmp4 - fm(p) + fm10)
fv(p) = ustar(p)
#endif
end do
#endif
#if (defined PERGRO)
!===============================================================================
! The following only applies when PERGRO is defined
!===============================================================================
!dir$ concurrent
!cdir nodep
do f = 1, fn
p = filterp(f)
g = pgridcell(p)
zldis(p) = forc_hgt_u(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetam) then ! zeta < -1
ustar(p) = vkc * um(p) / log(-zetam*obu(p)/z0m(p))
else if (zeta(p) < 0._r8) then ! -1 <= zeta < 0
ustar(p) = vkc * um(p) / log(zldis(p)/z0m(p))
else if (zeta(p) <= 1._r8) then ! 0 <= ztea <= 1
ustar(p)=vkc * um(p)/log(zldis(p)/z0m(p))
else ! 1 < zeta, phi=5+zeta
ustar(p)=vkc * um(p)/log(obu(p)/z0m(p))
endif
zldis(p) = forc_hgt_t(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp1(p)=vkc/log(-zetat*obu(p)/z0h(p))
else if (zeta(p) < 0._r8) then
temp1(p)=vkc/log(zldis(p)/z0h(p))
else if (zeta(p) <= 1._r8) then
temp1(p)=vkc/log(zldis(p)/z0h(p))
else
temp1(p)=vkc/log(obu(p)/z0h(p))
end if
zldis(p) = forc_hgt_q(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp2(p)=vkc/log(-zetat*obu(p)/z0q(p))
else if (zeta(p) < 0._r8) then
temp2(p)=vkc/log(zldis(p)/z0q(p))
else if (zeta(p) <= 1._r8) then
temp2(p)=vkc/log(zldis(p)/z0q(p))
else
temp2(p)=vkc/log(obu(p)/z0q(p))
end if
zldis(p) = 2.0_r8 + z0h(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp12m(p)=vkc/log(-zetat*obu(p)/z0h(p))
else if (zeta(p) < 0._r8) then
temp12m(p)=vkc/log(zldis(p)/z0h(p))
else if (zeta(p) <= 1._r8) then
temp12m(p)=vkc/log(zldis(p)/z0h(p))
else
temp12m(p)=vkc/log(obu(p)/z0h(p))
end if
zldis(p) = 2.0_r8 + z0q(p)
zeta(p) = zldis(p)/obu(p)
if (zeta(p) < -zetat) then
temp22m(p)=vkc/log(-zetat*obu(p)/z0q(p))
else if (zeta(p) < 0._r8) then
temp22m(p)=vkc/log(zldis(p)/z0q(p))
else if (zeta(p) <= 1._r8) then
temp22m(p)=vkc/log(zldis(p)/z0q(p))
else
temp22m(p)=vkc/log(obu(p)/z0q(p))
end if
#if (defined DGVM) || (defined DUST)
! diagnose 10-m wind for dust model (dstmbl.F)
! Notes from C. Zender's dst.F:
! According to Bon96 p. 62, the displacement height d (here displa) is
! 0.0 <= d <= 0.34 m in dust source regions (i.e., regions w/o trees).
! Therefore d <= 0.034*z1 and may safely be neglected.
! Code from LSM routine SurfaceTemperature was used to obtain u10
zldis(p) = forc_hgt_u(g)-displa(p)
zeta(p) = zldis(p)/obu(p)
if (min(zeta(p), 1._r8) < 0._r8) then
tmp1 = (1._r8 - 16._r8*min(zeta(p),1._r8))**0.25_r8
tmp2 = log((1._r8+tmp1*tmp1)/2._r8)
tmp3 = log((1._r8+tmp1)/2._r8)
fmnew = 2._r8*tmp3 + tmp2 - 2._r8*atan(tmp1) + 1.5707963_r8
else
fmnew = -5._r8*min(zeta(p),1._r8)
endif
if (iter == 1) then
fm(p) = fmnew
else
fm(p) = 0.5_r8 * (fm(p)+fmnew)
end if
zeta10 = min(10._r8/obu(p), 1._r8)
if (zeta(p) == 0._r8) zeta10 = 0._r8
if (zeta10 < 0._r8) then
tmp1 = (1.0_r8 - 16.0_r8 * zeta10)**0.25_r8
tmp2 = log((1.0_r8 + tmp1*tmp1)/2.0_r8)
tmp3 = log((1.0_r8 + tmp1)/2.0_r8)
fm10 = 2.0_r8*tmp3 + tmp2 - 2.0_r8*atan(tmp1) + 1.5707963_r8
else ! not stable
fm10 = -5.0_r8 * zeta10
end if
tmp4 = log(forc_hgt(g) / 10._r8)
u10(p) = ur(p) - ustar(p)/vkc * (tmp4 - fm(p) + fm10)
fv(p) = ustar(p)
#endif
end do
#endif
end subroutine FrictionVelocity
! !IROUTINE: StabilityFunc
!
! !INTERFACE:
real(r8) function StabilityFunc1(zeta),1
!
! !DESCRIPTION:
! Stability function for rib < 0.
!
! !USES:
! use shr_const_mod, only: SHR_CONST_PI
!Zack Subin, 7/8/08
!
! !ARGUMENTS:
implicit none
real(r8), intent(in) :: zeta ! dimensionless height used in Monin-Obukhov theory
!
! !CALLED FROM:
! subroutine FrictionVelocity in this module
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
!
!EOP
!
! !LOCAL VARIABLES:
real(r8) :: chik, chik2
!------------------------------------------------------------------------------
chik2 = sqrt(1._r8-16._r8*zeta)
chik = sqrt(chik2)
StabilityFunc1 = 2._r8*log((1._r8+chik)*0.5_r8) &
!Changed to pie, Zack Subin, 7/9/08
+ log((1._r8+chik2)*0.5_r8)-2._r8*atan(chik)+pie*0.5_r8
end function StabilityFunc1
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: StabilityFunc2
!
! !INTERFACE:
real(r8) function StabilityFunc2(zeta),1
!
! !DESCRIPTION:
! Stability function for rib < 0.
!
! !USES:
!Removed by Zack Subin, 7/9/08
! use shr_const_mod, only: SHR_CONST_PI
!
! !ARGUMENTS:
implicit none
real(r8), intent(in) :: zeta ! dimensionless height used in Monin-Obukhov theory
!
! !CALLED FROM:
! subroutine FrictionVelocity in this module
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
!
!EOP
!
! !LOCAL VARIABLES:
real(r8) :: chik2
!------------------------------------------------------------------------------
chik2 = sqrt(1._r8-16._r8*zeta)
StabilityFunc2 = 2._r8*log((1._r8+chik2)*0.5_r8)
end function StabilityFunc2
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: MoninObukIni
!
! !INTERFACE:
subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu) 4,1
!
! !DESCRIPTION:
! Initialization of the Monin-Obukhov length.
! The scheme is based on the work of Zeng et al. (1998):
! Intercomparison of bulk aerodynamic algorithms for the computation
! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
! Vol. 11, 2628-2644.
!
! !USES:
!
! !ARGUMENTS:
implicit none
real(r8), intent(in) :: ur ! wind speed at reference height [m/s]
real(r8), intent(in) :: thv ! virtual potential temperature (kelvin)
real(r8), intent(in) :: dthv ! diff of vir. poten. temp. between ref. height and surface
real(r8), intent(in) :: zldis ! reference height "minus" zero displacement heght [m]
real(r8), intent(in) :: z0m ! roughness length, momentum [m]
real(r8), intent(out) :: um ! wind speed including the stability effect [m/s]
real(r8), intent(out) :: obu ! monin-obukhov length (m)
!
! !CALLED FROM:
! subroutine BareGroundFluxes in module BareGroundFluxesMod.F90
! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod.F90
! subroutine CanopyFluxes in module CanopyFluxesMod.F90
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
!
!EOP
!
! !LOCAL VARIABLES:
!
real(r8) :: wc ! convective velocity [m/s]
real(r8) :: rib ! bulk Richardson number
real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory
real(r8) :: ustar ! friction velocity [m/s]
!-----------------------------------------------------------------------
! Initial values of u* and convective velocity
ustar=0.06_r8
wc=0.5_r8
if (dthv >= 0._r8) then
um=max(ur,0.1_r8)
else
um=sqrt(ur*ur+wc*wc)
endif
rib=grav*zldis*dthv/(thv*um*um)
#if (defined PERGRO)
rib = 0._r8
#endif
if (rib >= 0._r8) then ! neutral or stable
zeta = rib*log(zldis/z0m)/(1._r8-5._r8*min(rib,0.19_r8))
zeta = min(2._r8,max(zeta,0.01_r8 ))
else ! unstable
zeta=rib*log(zldis/z0m)
zeta = max(-100._r8,min(zeta,-0.01_r8 ))
endif
obu=zldis/zeta
end subroutine MoninObukIni
subroutine LakeDebug( str ) ,1
IMPLICIT NONE
CHARACTER*(*), str
CALL wrf_debug( 0 , TRIM(str) )
end subroutine LakeDebug
SUBROUTINE lakeini(IVGTYP, ISLTYP, HT, SNOW, & !i 1,3
lake_min_elev, restart, lakedepth_default, lake_depth, &
lakedepth2d, savedtke12d, snowdp2d, h2osno2d, & !o
snl2d, t_grnd2d, t_lake3d, lake_icefrac3d, &
z_lake3d, dz_lake3d, t_soisno3d, h2osoi_ice3d, &
h2osoi_liq3d, h2osoi_vol3d, z3d, dz3d, &
zi3d, watsat3d, csol3d, tkmg3d, &
iswater, xice, xice_threshold, xland, tsk, &
#if (EM_CORE == 1)
lakemask, lakeflag, &
#endif
lake_depth_flag, use_lakedepth, &
tkdry3d, tksatu3d, lake, its, ite, jts, jte, &
ims,ime, jms,jme)
!==============================================================================
! This subroutine was first edited by Hongping Gu for coupling
! 07/20/2010
!==============================================================================
USE module_wrf_error
implicit none
INTEGER , INTENT (IN) :: iswater
REAL, INTENT(IN) :: xice_threshold
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: XICE
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN):: TSK
REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(INOUT) :: XLAND
#if (EM_CORE == 1)
REAL, DIMENSION( ims:ime , jms:jme ) :: LAKEMASK
INTEGER , INTENT (IN) :: lakeflag
#endif
INTEGER , INTENT (INOUT) :: lake_depth_flag
INTEGER , INTENT (IN) :: use_lakedepth
LOGICAL , INTENT(IN) :: restart
INTEGER, INTENT(IN ) :: ims,ime, jms,jme
INTEGER, INTENT(IN ) :: its,ite, jts,jte
INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: IVGTYP, &
ISLTYP
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: HT
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SNOW
real, intent(in) :: lakedepth_default,lake_min_elev
real, dimension(ims:ime,jms:jme ),intent(out) :: lakedepth2d, &
savedtke12d
real, dimension(ims:ime,jms:jme ),intent(out) :: snowdp2d, &
h2osno2d, &
snl2d, &
t_grnd2d
real, dimension( ims:ime,1:nlevlake, jms:jme ),INTENT(out) :: t_lake3d, &
lake_icefrac3d, &
z_lake3d, &
dz_lake3d
real, dimension( ims:ime,-nlevsnow+1:nlevsoil, jms:jme ),INTENT(out) :: t_soisno3d, &
h2osoi_ice3d, &
h2osoi_liq3d, &
h2osoi_vol3d, &
z3d, &
dz3d
real, dimension( ims:ime,1:nlevsoil, jms:jme ),INTENT(out) :: watsat3d, &
csol3d, &
tkmg3d, &
tkdry3d, &
tksatu3d
real, dimension( ims:ime,-nlevsnow+0:nlevsoil, jms:jme ),INTENT(out) :: zi3d
LOGICAL, DIMENSION( ims:ime, jms:jme ),intent(out) :: lake
REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lake_depth
real, dimension( ims:ime,1:nlevsoil, jms:jme ) :: bsw3d, &
bsw23d, &
psisat3d, &
vwcsat3d, &
watdry3d, &
watopt3d, &
hksat3d, &
sucsat3d, &
clay3d, &
sand3d
integer :: n,i,j,k,ib,lev,bottom ! indices
real(r8),dimension(ims:ime,jms:jme ) :: bd2d ! bulk density of dry soil material [kg/m^3]
real(r8),dimension(ims:ime,jms:jme ) :: tkm2d ! mineral conductivity
real(r8),dimension(ims:ime,jms:jme ) :: xksat2d ! maximum hydraulic conductivity of soil [mm/s]
real(r8),dimension(ims:ime,jms:jme ) :: depthratio2d ! ratio of lake depth to standard deep lake depth
real(r8),dimension(ims:ime,jms:jme ) :: clay2d ! temporary
real(r8),dimension(ims:ime,jms:jme ) :: sand2d ! temporary
real(r8) :: scalez = 0.025_r8 ! Soil layer thickness discretization (m)
logical,parameter :: arbinit = .true.
real,parameter :: defval = -999.0
integer :: isl
integer :: numb_lak ! for debug
character*256 :: message
IF ( RESTART ) RETURN
DO j = jts,jte
DO i = its,ite
snowdp2d(i,j) = snow(i,j)*0.005 ! SNOW in kg/m^2 and snowdp in m
h2osno2d(i,j) = snow(i,j) ! mm
ENDDO
ENDDO
! initialize all the grid with default value
DO j = jts,jte
DO i = its,ite
lakedepth2d(i,j) = defval
snl2d(i,j) = defval
do k = -nlevsnow+1,nlevsoil
h2osoi_liq3d(i,k,j) = defval
h2osoi_ice3d(i,k,j) = defval
t_soisno3d(i,k,j) = defval
z3d(i,k,j) = defval
dz3d(i,k,j) = defval
enddo
do k = 1,nlevlake
t_lake3d(i,k,j) = defval
lake_icefrac3d(i,k,j) = defval
z_lake3d(i,k,j) = defval
dz_lake3d(i,k,j) = defval
enddo
ENDDO
ENDDO
! judge whether the grid is lake grid
numb_lak = 0
do i=its,ite
do j=jts,jte
#if (EM_CORE==1)
IF (lakeflag.eq.0) THEN
if(ht(i,j)>=lake_min_elev) then
if ( xice(i,j).gt.xice_threshold) then !mchen
ivgtyp(i,j) = iswater
xland(i,j) = 2.
lake_icefrac3d(i,1,j) = xice(i,j)
xice(i,j)=0.0
endif
endif
if(ivgtyp(i,j)==iswater.and.ht(i,j)>=lake_min_elev) then
lake(i,j) = .true.
lakemask(i,j) = 1
numb_lak = numb_lak + 1
else
lake(i,j) = .false.
lakemask(i,j) = 0
end if
ELSE
if(lakemask(i,j).eq.1) then
lake(i,j) = .true.
numb_lak = numb_lak + 1
if ( xice(i,j).gt.xice_threshold) then !mchen
ivgtyp(i,j) = iswater
xland(i,j) = 2.
lake_icefrac3d(i,1,j) = xice(i,j)
xice(i,j)=0.0
endif
else
lake(i,j) = .false.
endif
ENDIF ! end if lakeflag=0
#else
if(ht(i,j)>=lake_min_elev) then
if ( xice(i,j).gt.xice_threshold) then !mchen
ivgtyp(i,j) = iswater
xland(i,j) = 2.
lake_icefrac3d(i,1,j) = xice(i,j)
xice(i,j)=0.0
endif
endif
if(ivgtyp(i,j)==iswater.and.ht(i,j)>=lake_min_elev) then
lake(i,j) = .true.
numb_lak = numb_lak + 1
else
lake(i,j) = .false.
end if
#endif
end do
end do
write(message,*) "the total number of lake grid is :", numb_lak
CALL wrf_message(message)
! CALL LakeDebug(msg)
! initialize lake grid
DO j = jts,jte
DO i = its,ite
if ( lake(i,j) ) then
! t_soisno3d(i,:,j) = tsk(i,j)
! t_lake3d(i,:,j) = tsk(i,j)
! t_grnd2d(i,j) = tsk(i,j)
z3d(i,:,j) = 0.0
dz3d(i,:,j) = 0.0
zi3d(i,:,j) = 0.0
h2osoi_liq3d(i,:,j) = 0.0
h2osoi_ice3d(i,:,j) = 0.0
lake_icefrac3d(i,:,j) = 0.0
h2osoi_vol3d(i,:,j) = 0.0
snl2d(i,j) = 0.0
if ( use_lakedepth.eq.1 .and.lake_depth_flag.eq.0 ) then !mchen
call wrf_error_fatal ( 'STOP: You need lake-depth information. Rerun WPS or set use_lakedepth = 0')
end if
if ( use_lakedepth.eq.0 .and.lake_depth_flag.eq.1 ) then !mchen
lake_depth_flag = 0
end if
if ( lake_depth_flag.eq.1 ) then
if (lake_depth(i,j) > 0.0) then
lakedepth2d(i,j) = lake_depth(i,j)
else
if ( lakedepth_default > 0.0 ) then
lakedepth2d(i,j) = lakedepth_default
else
lakedepth2d(i,j) = spval
endif
endif
else
if ( lakedepth_default > 0.0 ) then
lakedepth2d(i,j) = lakedepth_default
else
lakedepth2d(i,j) = spval
endif
endif
endif
ENDDO
ENDDO
#ifndef EXTRALAKELAYERS
! dzlak(1) = 0.1_r8
! dzlak(2) = 1._r8
! dzlak(3) = 2._r8
! dzlak(4) = 3._r8
! dzlak(5) = 4._r8
! dzlak(6) = 5._r8
! dzlak(7) = 7._r8
! dzlak(8) = 7._r8
! dzlak(9) = 10.45_r8
! dzlak(10)= 10.45_r8
!
! zlak(1) = 0.05_r8
! zlak(2) = 0.6_r8
! zlak(3) = 2.1_r8
! zlak(4) = 4.6_r8
! zlak(5) = 8.1_r8
! zlak(6) = 12.6_r8
! zlak(7) = 18.6_r8
! zlak(8) = 25.6_r8
! zlak(9) = 34.325_r8
! zlak(10)= 44.775_r8
dzlak(1) = 0.1_r8
dzlak(2) = 0.1_r8
dzlak(3) = 0.1_r8
dzlak(4) = 0.1_r8
dzlak(5) = 0.1_r8
dzlak(6) = 0.1_r8
dzlak(7) = 0.1_r8
dzlak(8) = 0.1_r8
dzlak(9) = 0.1_r8
dzlak(10)= 0.1_r8
zlak(1) = 0.05_r8
zlak(2) = 0.15_r8
zlak(3) = 0.25_r8
zlak(4) = 0.35_r8
zlak(5) = 0.45_r8
zlak(6) = 0.55_r8
zlak(7) = 0.65_r8
zlak(8) = 0.75_r8
zlak(9) = 0.85_r8
zlak(10)= 0.95_r8
#else
dzlak(1) =0.1_r8
dzlak(2) =0.25_r8
dzlak(3) =0.25_r8
dzlak(4) =0.25_r8
dzlak(5) =0.25_r8
dzlak(6) =0.5_r8
dzlak(7) =0.5_r8
dzlak(8) =0.5_r8
dzlak(9) =0.5_r8
dzlak(10) =0.75_r8
dzlak(11) =0.75_r8
dzlak(12) =0.75_r8
dzlak(13) =0.75_r8
dzlak(14) =2_r8
dzlak(15) =2_r8
dzlak(16) =2.5_r8
dzlak(17) =2.5_r8
dzlak(18) =3.5_r8
dzlak(19) =3.5_r8
dzlak(20) =3.5_r8
dzlak(21) =3.5_r8
dzlak(22) =5.225_r8
dzlak(23) =5.225_r8
dzlak(24) =5.225_r8
dzlak(25) =5.225_r8
zlak(1) = dzlak(1)/2._r8
do k = 2,nlevlake
zlak(k) = zlak(k-1) + (dzlak(k-1)+dzlak(k))/2._r8
end do
#endif
! "0" refers to soil surface and "nlevsoil" refers to the bottom of model soil
do j = 1, nlevsoil
zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths
enddo
dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2)) !thickness b/n two interfaces
do j = 2,nlevsoil-1
dzsoi(j)= 0.5_r8*(zsoi(j+1)-zsoi(j-1))
enddo
dzsoi(nlevsoil) = zsoi(nlevsoil)-zsoi(nlevsoil-1)
zisoi(0) = 0._r8
do j = 1, nlevsoil-1
zisoi(j) = 0.5_r8*(zsoi(j)+zsoi(j+1)) !interface depths
enddo
zisoi(nlevsoil) = zsoi(nlevsoil) + 0.5_r8*dzsoi(nlevsoil)
!!!!!!!!!!!!!!!!!!begin to initialize lake variables!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DO j = jts,jte
DO i = its,ite
if ( lake(i,j) ) then
! Soil hydraulic and thermal properties
isl = ISLTYP(i,j)
if (isl == 14 ) isl = isl + 1
do k = 1,nlevsoil
sand3d(i,k,j) = sand(isl)
clay3d(i,k,j) = clay(isl)
enddo
do k = 1,nlevsoil
clay2d(i,j) = clay3d(i,k,j)
sand2d(i,j) = sand3d(i,k,j)
watsat3d(i,k,j) = 0.489_r8 - 0.00126_r8*sand2d(i,j)
bd2d(i,j) = (1._r8-watsat3d(i,k,j))*2.7e3_r8
xksat2d(i,j) = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand2d(i,j)) ) ! mm/s
tkm2d(i,j) = (8.80_r8*sand2d(i,j)+2.92_r8*clay2d(i,j))/(sand2d(i,j)+clay2d(i,j)) ! W/(m K)
bsw3d(i,k,j) = 2.91_r8 + 0.159_r8*clay2d(i,j)
bsw23d(i,k,j) = -(3.10_r8 + 0.157_r8*clay2d(i,j) - 0.003_r8*sand2d(i,j))
psisat3d(i,k,j) = -(exp((1.54_r8 - 0.0095_r8*sand2d(i,j) + 0.0063_r8*(100.0_r8-sand2d(i,j) &
-clay2d(i,j)))*log(10.0_r8))*9.8e-5_r8)
vwcsat3d(i,k,j) = (50.5_r8 - 0.142_r8*sand2d(i,j) - 0.037_r8*clay2d(i,j))/100.0_r8
hksat3d(i,k,j) = xksat2d(i,j)
sucsat3d(i,k,j) = 10._r8 * ( 10._r8**(1.88_r8-0.0131_r8*sand2d(i,j)) )
tkmg3d(i,k,j) = tkm2d(i,j) ** (1._r8- watsat3d(i,k,j))
tksatu3d(i,k,j) = tkmg3d(i,k,j)*0.57_r8**watsat3d(i,k,j)
tkdry3d(i,k,j) = (0.135_r8*bd2d(i,j) + 64.7_r8) / (2.7e3_r8 - 0.947_r8*bd2d(i,j))
csol3d(i,k,j) = (2.128_r8*sand2d(i,j)+2.385_r8*clay2d(i,j)) / (sand2d(i,j)+clay2d(i,j))*1.e6_r8 ! J/(m3 K)
watdry3d(i,k,j) = watsat3d(i,k,j) * (316230._r8/sucsat3d(i,k,j)) ** (-1._r8/bsw3d(i,k,j))
watopt3d(i,k,j) = watsat3d(i,k,j) * (158490._r8/sucsat3d(i,k,j)) ** (-1._r8/bsw3d(i,k,j))
end do
if (lakedepth2d(i,j) == spval) then
lakedepth2d(i,j) = zlak(nlevlake) + 0.5_r8*dzlak(nlevlake)
z_lake3d(i,1:nlevlake,j) = zlak(1:nlevlake)
dz_lake3d(i,1:nlevlake,j) = dzlak(1:nlevlake)
else
depthratio2d(i,j) = lakedepth2d(i,j) / (zlak(nlevlake) + 0.5_r8*dzlak(nlevlake))
z_lake3d(i,1,j) = zlak(1)
dz_lake3d(i,1,j) = dzlak(1)
dz_lake3d(i,2:nlevlake,j) = dzlak(2:nlevlake)*depthratio2d(i,j)
z_lake3d(i,2:nlevlake,j) = zlak(2:nlevlake)*depthratio2d(i,j) + dz_lake3d(i,1,j)*(1._r8 - depthratio2d(i,j))
end if
! initial t_lake3d here
t_soisno3d(i,1,j) = tsk(i,j)
t_lake3d(i,1,j) = tsk(i,j)
t_grnd2d(i,j) = 277.0
do k = 2, nlevlake
if(z_lake3d(i,k,j).le.depth_c) then
t_soisno3d(i,k,j)=tsk(i,j)+(277.0-tsk(i,j))/depth_c*z_lake3d(i,k,j)
t_lake3d(i,k,j)=tsk(i,j)+(277.0-tsk(i,j))/depth_c*z_lake3d(i,k,j)
else
t_soisno3d(i,k,j) = 277.0
t_lake3d(i,k,j) = 277.0
end if
enddo
!end initial t_lake3d here
z3d(i,1:nlevsoil,j) = zsoi(1:nlevsoil)
zi3d(i,0:nlevsoil,j) = zisoi(0:nlevsoil)
dz3d(i,1:nlevsoil,j) = dzsoi(1:nlevsoil)
savedtke12d(i,j) = tkwat ! Initialize for first timestep.
if (snowdp2d(i,j) < 0.01_r8) then
snl2d(i,j) = 0
dz3d(i,-nlevsnow+1:0,j) = 0._r8
z3d (i,-nlevsnow+1:0,j) = 0._r8
zi3d(i,-nlevsnow+0:0,j) = 0._r8
else
if ((snowdp2d(i,j) >= 0.01_r8) .and. (snowdp2d(i,j) <= 0.03_r8)) then
snl2d(i,j) = -1
dz3d(i,0,j) = snowdp2d(i,j)
else if ((snowdp2d(i,j) > 0.03_r8) .and. (snowdp2d(i,j) <= 0.04_r8)) then
snl2d(i,j) = -2
dz3d(i,-1,j) = snowdp2d(i,j)/2._r8
dz3d(i, 0,j) = dz3d(i,-1,j)
else if ((snowdp2d(i,j) > 0.04_r8) .and. (snowdp2d(i,j) <= 0.07_r8)) then
snl2d(i,j) = -2
dz3d(i,-1,j) = 0.02_r8
dz3d(i, 0,j) = snowdp2d(i,j) - dz3d(i,-1,j)
else if ((snowdp2d(i,j) > 0.07_r8) .and. (snowdp2d(i,j) <= 0.12_r8)) then
snl2d(i,j) = -3
dz3d(i,-2,j) = 0.02_r8
dz3d(i,-1,j) = (snowdp2d(i,j) - 0.02_r8)/2._r8
dz3d(i, 0,j) = dz3d(i,-1,j)
else if ((snowdp2d(i,j) > 0.12_r8) .and. (snowdp2d(i,j) <= 0.18_r8)) then
snl2d(i,j) = -3
dz3d(i,-2,j) = 0.02_r8
dz3d(i,-1,j) = 0.05_r8
dz3d(i, 0,j) = snowdp2d(i,j) - dz3d(i,-2,j) - dz3d(i,-1,j)
else if ((snowdp2d(i,j) > 0.18_r8) .and. (snowdp2d(i,j) <= 0.29_r8)) then
snl2d(i,j) = -4
dz3d(i,-3,j) = 0.02_r8
dz3d(i,-2,j) = 0.05_r8
dz3d(i,-1,j) = (snowdp2d(i,j) - dz3d(i,-3,j) - dz3d(i,-2,j))/2._r8
dz3d(i, 0,j) = dz3d(i,-1,j)
else if ((snowdp2d(i,j) > 0.29_r8) .and. (snowdp2d(i,j) <= 0.41_r8)) then
snl2d(i,j) = -4
dz3d(i,-3,j) = 0.02_r8
dz3d(i,-2,j) = 0.05_r8
dz3d(i,-1,j) = 0.11_r8
dz3d(i, 0,j) = snowdp2d(i,j) - dz3d(i,-3,j) - dz3d(i,-2,j) - dz3d(i,-1,j)
else if ((snowdp2d(i,j) > 0.41_r8) .and. (snowdp2d(i,j) <= 0.64_r8)) then
snl2d(i,j) = -5
dz3d(i,-4,j) = 0.02_r8
dz3d(i,-3,j) = 0.05_r8
dz3d(i,-2,j) = 0.11_r8
dz3d(i,-1,j) = (snowdp2d(i,j) - dz3d(i,-4,j) - dz3d(i,-3,j) - dz3d(i,-2,j))/2._r8
dz3d(i, 0,j) = dz3d(i,-1,j)
else if (snowdp2d(i,j) > 0.64_r8) then
snl2d(i,j) = -5
dz3d(i,-4,j) = 0.02_r8
dz3d(i,-3,j) = 0.05_r8
dz3d(i,-2,j) = 0.11_r8
dz3d(i,-1,j) = 0.23_r8
dz3d(i, 0,j)=snowdp2d(i,j)-dz3d(i,-4,j)-dz3d(i,-3,j)-dz3d(i,-2,j)-dz3d(i,-1,j)
endif
end if
do k = 0, snl2d(i,j)+1, -1
z3d(i,k,j) = zi3d(i,k,j) - 0.5_r8*dz3d(i,k,j)
zi3d(i,k-1,j) = zi3d(i,k,j) - dz3d(i,k,j)
end do
! 3:subroutine makearbinit
if (snl2d(i,j) < 0) then
do k = snl2d(i,j)+1, 0
! Be careful because there may be new snow layers with bad temperatures like 0 even if
! coming from init. con. file.
if(arbinit .or. t_soisno3d(i,k,j) > 300 .or. t_soisno3d(i,k,j) < 200) t_soisno3d(i,k,j) = 250._r8
enddo
end if
do k = 1, nlevsoil
if(arbinit .or. t_soisno3d(i,k,j) > 1000 .or. t_soisno3d(i,k,j) < 0) t_soisno3d(i,k,j) = t_lake3d(i,nlevlake,j)
end do
do k = 1, nlevlake
if(arbinit .or. lake_icefrac3d(i,k,j) > 1._r8 .or. lake_icefrac3d(i,k,j) < 0._r8) then
if(t_lake3d(i,k,j) >= tfrz) then
lake_icefrac3d(i,k,j) = 0._r8
else
lake_icefrac3d(i,k,j) = 1._r8
end if
end if
end do
do k = 1,nlevsoil
if (arbinit .or. h2osoi_vol3d(i,k,j) > 10._r8 .or. h2osoi_vol3d(i,k,j) < 0._r8) h2osoi_vol3d(i,k,j) = 1.0_r8
h2osoi_vol3d(i,k,j) = min(h2osoi_vol3d(i,k,j),watsat3d(i,k,j))
! soil layers
if (t_soisno3d(i,k,j) <= tfrz) then
h2osoi_ice3d(i,k,j) = dz3d(i,k,j)*denice*h2osoi_vol3d(i,k,j)
h2osoi_liq3d(i,k,j) = 0._r8
else
h2osoi_ice3d(i,k,j) = 0._r8
h2osoi_liq3d(i,k,j) = dz3d(i,k,j)*denh2o*h2osoi_vol3d(i,k,j)
endif
enddo
do k = -nlevsnow+1, 0
if (k > snl2d(i,j)) then
h2osoi_ice3d(i,k,j) = dz3d(i,k,j)*bdsno
h2osoi_liq3d(i,k,j) = 0._r8
end if
end do
end if !lake(i,j)
ENDDO
ENDDO
END SUBROUTINE lakeini
END MODULE module_sf_lake