#include <misc.h>
#include <preproc.h>
module SnowHydrologyMod 1,2
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: SnowHydrologyMod
!
! !DESCRIPTION:
! Calculate snow hydrology.
!
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
use clm_varpar , only : nlevsno
!
! !PUBLIC TYPES:
implicit none
save
!
! !PUBLIC MEMBER FUNCTIONS:
public :: SnowWater ! Change of snow mass and the snow water onto soil
public :: SnowCompaction ! Change in snow layer thickness due to compaction
public :: CombineSnowLayers ! Combine snow layers less than a min thickness
public :: DivideSnowLayers ! Subdivide snow layers if they exceed maximum thickness
public :: BuildSnowFilter ! Construct snow/no-snow filters
!
! !PRIVATE MEMBER FUNCTIONS:
private :: Combo ! Returns the combined variables: dz, t, wliq, wice.
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SnowWater
!
! !INTERFACE:
subroutine SnowWater(lbc, ubc, num_snowc, filter_snowc, & 1,6
num_nosnowc, filter_nosnowc)
!
! !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.
!
! !USES:
use clmtype
use clm_varcon , only : denh2o, denice, wimp, ssi
use clm_time_manager, only : get_step_size
use clm_atmlnd , only : clm_a2l
use SNICARMod , only : scvng_fct_mlt_bcphi, scvng_fct_mlt_bcpho, &
scvng_fct_mlt_ocphi, scvng_fct_mlt_ocpho, &
scvng_fct_mlt_dst1, scvng_fct_mlt_dst2, &
scvng_fct_mlt_dst3, scvng_fct_mlt_dst4
!
! !ARGUMENTS:
implicit none
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
!
! !CALLED FROM:
!
! !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.
! 03/28/08, Mark Flanner: Added aerosol deposition and flushing with meltwater
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
integer , pointer :: snl(:) !number of snow layers
logical , pointer :: do_capsnow(:) !true => do snow capping
real(r8), pointer :: qflx_snomelt(:) !snow melt (mm H2O /s)
real(r8), pointer :: qflx_rain_grnd(:) !rain on ground after interception (mm H2O/s) [+]
real(r8), pointer :: qflx_sub_snow(:) !sublimation rate from snow pack (mm H2O /s) [+]
real(r8), pointer :: qflx_evap_grnd(:) !ground surface evaporation rate (mm H2O/s) [+]
real(r8), pointer :: qflx_dew_snow(:) !surface dew added to snow pack (mm H2O /s) [+]
real(r8), pointer :: qflx_dew_grnd(:) !ground surface dew formation (mm H2O /s) [+]
real(r8), pointer :: dz(:,:) !layer depth (m)
!
! local pointers to implicit out arguments
!
real(r8), pointer :: qflx_top_soil(:) !net water input into soil from top (mm/s)
!
! local pointers to implicit inout arguments
!
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
integer , pointer :: cgridcell(:) ! columns's gridcell (col)
real(r8), pointer :: mss_bcphi(:,:) ! hydrophillic BC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_bcpho(:,:) ! hydrophobic BC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_ocphi(:,:) ! hydrophillic OC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_ocpho(:,:) ! hydrophobic OC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst1(:,:) ! mass of dust species 1 in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst2(:,:) ! mass of dust species 2 in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst3(:,:) ! mass of dust species 3 in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst4(:,:) ! mass of dust species 4 in snow (col,lyr) [kg]
real(r8), pointer :: flx_bc_dep_dry(:) ! dry BC deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_bc_dep_wet(:) ! wet BC deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_bc_dep(:) ! total BC deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_bc_dep_pho(:) ! hydrophobic BC deposition (col) [kg m-1 s-1]
real(r8), pointer :: flx_bc_dep_phi(:) ! hydrophillic BC deposition (col) [kg m-1 s-1]
real(r8), pointer :: flx_oc_dep_dry(:) ! dry OC deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_oc_dep_wet(:) ! wet OC deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_oc_dep(:) ! total OC deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_oc_dep_pho(:) ! hydrophobic OC deposition (col) [kg m-1 s-1]
real(r8), pointer :: flx_oc_dep_phi(:) ! hydrophillic OC deposition (col) [kg m-1 s-1]
real(r8), pointer :: flx_dst_dep_dry1(:) ! dry dust (species 1) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_wet1(:) ! wet dust (species 1) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_dry2(:) ! dry dust (species 2) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_wet2(:) ! wet dust (species 2) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_dry3(:) ! dry dust (species 3) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_wet3(:) ! wet dust (species 3) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_dry4(:) ! dry dust (species 4) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep_wet4(:) ! wet dust (species 4) deposition (col) [kg m-2 s-1]
real(r8), pointer :: flx_dst_dep(:) ! total dust deposition (col) [kg m-2 s-1]
real(r8), pointer :: forc_aer(:,:) ! aerosol deposition from atmosphere model (grd,aer) [kg m-1 s-1]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
integer :: c, j, fc !do loop/array indices
real(r8) :: dtime !land model time step (sec)
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,-nlevsno+1:0) !partial volume of liquid water in layer
real(r8) :: vol_ice(lbc:ubc,-nlevsno+1:0) !partial volume of ice lens in layer
real(r8) :: eff_porosity(lbc:ubc,-nlevsno+1:0) !effective porosity = porosity - vol_ice
integer :: g ! gridcell loop index
real(r8) :: qin_bc_phi(lbc:ubc) ! flux of hydrophilic BC into layer [kg]
real(r8) :: qout_bc_phi(lbc:ubc) ! flux of hydrophilic BC out of layer [kg]
real(r8) :: qin_bc_pho(lbc:ubc) ! flux of hydrophobic BC into layer [kg]
real(r8) :: qout_bc_pho(lbc:ubc) ! flux of hydrophobic BC out of layer [kg]
real(r8) :: qin_oc_phi(lbc:ubc) ! flux of hydrophilic OC into layer [kg]
real(r8) :: qout_oc_phi(lbc:ubc) ! flux of hydrophilic OC out of layer [kg]
real(r8) :: qin_oc_pho(lbc:ubc) ! flux of hydrophobic OC into layer [kg]
real(r8) :: qout_oc_pho(lbc:ubc) ! flux of hydrophobic OC out of layer [kg]
real(r8) :: qin_dst1(lbc:ubc) ! flux of dust species 1 into layer [kg]
real(r8) :: qout_dst1(lbc:ubc) ! flux of dust species 1 out of layer [kg]
real(r8) :: qin_dst2(lbc:ubc) ! flux of dust species 2 into layer [kg]
real(r8) :: qout_dst2(lbc:ubc) ! flux of dust species 2 out of layer [kg]
real(r8) :: qin_dst3(lbc:ubc) ! flux of dust species 3 into layer [kg]
real(r8) :: qout_dst3(lbc:ubc) ! flux of dust species 3 out of layer [kg]
real(r8) :: qin_dst4(lbc:ubc) ! flux of dust species 4 into layer [kg]
real(r8) :: qout_dst4(lbc:ubc) ! flux of dust species 4 out of layer [kg]
real(r8) :: mss_liqice ! mass of liquid+ice in a layer
!-----------------------------------------------------------------------
! Assign local pointers to derived subtype components (column-level)
snl => clm3%g%l%c%cps%snl
do_capsnow => clm3%g%l%c%cps%do_capsnow
qflx_snomelt => clm3%g%l%c%cwf%qflx_snomelt
qflx_rain_grnd => clm3%g%l%c%cwf%pwf_a%qflx_rain_grnd
qflx_sub_snow => clm3%g%l%c%cwf%pwf_a%qflx_sub_snow
qflx_evap_grnd => clm3%g%l%c%cwf%pwf_a%qflx_evap_grnd
qflx_dew_snow => clm3%g%l%c%cwf%pwf_a%qflx_dew_snow
qflx_dew_grnd => clm3%g%l%c%cwf%pwf_a%qflx_dew_grnd
qflx_top_soil => clm3%g%l%c%cwf%qflx_top_soil
dz => clm3%g%l%c%cps%dz
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
cgridcell => clm3%g%l%c%gridcell
mss_bcphi => clm3%g%l%c%cps%mss_bcphi
mss_bcpho => clm3%g%l%c%cps%mss_bcpho
mss_ocphi => clm3%g%l%c%cps%mss_ocphi
mss_ocpho => clm3%g%l%c%cps%mss_ocpho
mss_dst1 => clm3%g%l%c%cps%mss_dst1
mss_dst2 => clm3%g%l%c%cps%mss_dst2
mss_dst3 => clm3%g%l%c%cps%mss_dst3
mss_dst4 => clm3%g%l%c%cps%mss_dst4
flx_bc_dep => clm3%g%l%c%cwf%flx_bc_dep
flx_bc_dep_wet => clm3%g%l%c%cwf%flx_bc_dep_wet
flx_bc_dep_dry => clm3%g%l%c%cwf%flx_bc_dep_dry
flx_bc_dep_phi => clm3%g%l%c%cwf%flx_bc_dep_phi
flx_bc_dep_pho => clm3%g%l%c%cwf%flx_bc_dep_pho
flx_oc_dep => clm3%g%l%c%cwf%flx_oc_dep
flx_oc_dep_wet => clm3%g%l%c%cwf%flx_oc_dep_wet
flx_oc_dep_dry => clm3%g%l%c%cwf%flx_oc_dep_dry
flx_oc_dep_phi => clm3%g%l%c%cwf%flx_oc_dep_phi
flx_oc_dep_pho => clm3%g%l%c%cwf%flx_oc_dep_pho
flx_dst_dep => clm3%g%l%c%cwf%flx_dst_dep
flx_dst_dep_wet1 => clm3%g%l%c%cwf%flx_dst_dep_wet1
flx_dst_dep_dry1 => clm3%g%l%c%cwf%flx_dst_dep_dry1
flx_dst_dep_wet2 => clm3%g%l%c%cwf%flx_dst_dep_wet2
flx_dst_dep_dry2 => clm3%g%l%c%cwf%flx_dst_dep_dry2
flx_dst_dep_wet3 => clm3%g%l%c%cwf%flx_dst_dep_wet3
flx_dst_dep_dry3 => clm3%g%l%c%cwf%flx_dst_dep_dry3
flx_dst_dep_wet4 => clm3%g%l%c%cwf%flx_dst_dep_wet4
flx_dst_dep_dry4 => clm3%g%l%c%cwf%flx_dst_dep_dry4
forc_aer => clm_a2l%forc_aer
! Determine model time step
dtime = get_step_size()
! Renew the mass of ice lens (h2osoi_ice) and liquid (h2osoi_liq) in the
! surface snow layer resulting from sublimation (frost) / evaporation (condense)
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._r8) then
h2osoi_ice(c,snl(c)+1) = 0._r8
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._r8) then
h2osoi_ice(c,snl(c)+1) = 0._r8
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 = -nlevsno+1, 0
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._r8 - 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.
! Also compute aerosol fluxes through snowpack in this loop:
! 1) compute aerosol mass in each layer
! 2) add aerosol mass flux from above layer to mass of this layer
! 3) qout_xxx is mass flux of aerosol species xxx out bottom of
! layer in water flow, proportional to (current) concentration
! of aerosol in layer multiplied by a scavenging ratio.
! 4) update mass of aerosol in top layer, accordingly
! 5) update mass concentration of aerosol accordingly
qin(:) = 0._r8
qin_bc_phi(:) = 0._r8
qin_bc_pho(:) = 0._r8
qin_oc_phi(:) = 0._r8
qin_oc_pho(:) = 0._r8
qin_dst1(:) = 0._r8
qin_dst2(:) = 0._r8
qin_dst3(:) = 0._r8
qin_dst4(:) = 0._r8
do j = -nlevsno+1, 0
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)
mss_bcphi(c,j) = mss_bcphi(c,j) + qin_bc_phi(c)
mss_bcpho(c,j) = mss_bcpho(c,j) + qin_bc_pho(c)
mss_ocphi(c,j) = mss_ocphi(c,j) + qin_oc_phi(c)
mss_ocpho(c,j) = mss_ocpho(c,j) + qin_oc_pho(c)
mss_dst1(c,j) = mss_dst1(c,j) + qin_dst1(c)
mss_dst2(c,j) = mss_dst2(c,j) + qin_dst2(c)
mss_dst3(c,j) = mss_dst3(c,j) + qin_dst3(c)
mss_dst4(c,j) = mss_dst4(c,j) + qin_dst4(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._r8-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._r8
h2osoi_liq(c,j) = h2osoi_liq(c,j) - qout(c)
qin(c) = qout(c)
! mass of ice+water: in extremely rare circumstances, this can
! be zero, even though there is a snow layer defined. In
! this case, set the mass to a very small value to
! prevent division by zero.
mss_liqice = h2osoi_liq(c,j)+h2osoi_ice(c,j)
if (mss_liqice < 1E-30_r8) then
mss_liqice = 1E-30_r8
endif
! BCPHI:
! 1. flux with meltwater:
qout_bc_phi(c) = qout(c)*scvng_fct_mlt_bcphi*(mss_bcphi(c,j)/mss_liqice)
if (qout_bc_phi(c) > mss_bcphi(c,j)) then
qout_bc_phi(c) = mss_bcphi(c,j)
endif
mss_bcphi(c,j) = mss_bcphi(c,j) - qout_bc_phi(c)
qin_bc_phi(c) = qout_bc_phi(c)
! BCPHO:
! 1. flux with meltwater:
qout_bc_pho(c) = qout(c)*scvng_fct_mlt_bcpho*(mss_bcpho(c,j)/mss_liqice)
if (qout_bc_pho(c) > mss_bcpho(c,j)) then
qout_bc_pho(c) = mss_bcpho(c,j)
endif
mss_bcpho(c,j) = mss_bcpho(c,j) - qout_bc_pho(c)
qin_bc_pho(c) = qout_bc_pho(c)
! OCPHI:
! 1. flux with meltwater:
qout_oc_phi(c) = qout(c)*scvng_fct_mlt_ocphi*(mss_ocphi(c,j)/mss_liqice)
if (qout_oc_phi(c) > mss_ocphi(c,j)) then
qout_oc_phi(c) = mss_ocphi(c,j)
endif
mss_ocphi(c,j) = mss_ocphi(c,j) - qout_oc_phi(c)
qin_oc_phi(c) = qout_oc_phi(c)
! OCPHO:
! 1. flux with meltwater:
qout_oc_pho(c) = qout(c)*scvng_fct_mlt_ocpho*(mss_ocpho(c,j)/mss_liqice)
if (qout_oc_pho(c) > mss_ocpho(c,j)) then
qout_oc_pho(c) = mss_ocpho(c,j)
endif
mss_ocpho(c,j) = mss_ocpho(c,j) - qout_oc_pho(c)
qin_oc_pho(c) = qout_oc_pho(c)
! DUST 1:
! 1. flux with meltwater:
qout_dst1(c) = qout(c)*scvng_fct_mlt_dst1*(mss_dst1(c,j)/mss_liqice)
if (qout_dst1(c) > mss_dst1(c,j)) then
qout_dst1(c) = mss_dst1(c,j)
endif
mss_dst1(c,j) = mss_dst1(c,j) - qout_dst1(c)
qin_dst1(c) = qout_dst1(c)
! DUST 2:
! 1. flux with meltwater:
qout_dst2(c) = qout(c)*scvng_fct_mlt_dst2*(mss_dst2(c,j)/mss_liqice)
if (qout_dst2(c) > mss_dst2(c,j)) then
qout_dst2(c) = mss_dst2(c,j)
endif
mss_dst2(c,j) = mss_dst2(c,j) - qout_dst2(c)
qin_dst2(c) = qout_dst2(c)
! DUST 3:
! 1. flux with meltwater:
qout_dst3(c) = qout(c)*scvng_fct_mlt_dst3*(mss_dst3(c,j)/mss_liqice)
if (qout_dst3(c) > mss_dst3(c,j)) then
qout_dst3(c) = mss_dst3(c,j)
endif
mss_dst3(c,j) = mss_dst3(c,j) - qout_dst3(c)
qin_dst3(c) = qout_dst3(c)
! DUST 4:
! 1. flux with meltwater:
qout_dst4(c) = qout(c)*scvng_fct_mlt_dst4*(mss_dst4(c,j)/mss_liqice)
if (qout_dst4(c) > mss_dst4(c,j)) then
qout_dst4(c) = mss_dst4(c,j)
endif
mss_dst4(c,j) = mss_dst4(c,j) - qout_dst4(c)
qin_dst4(c) = qout_dst4(c)
end if
end do
end do
! Adjust layer thickness for any water+ice content changes in excess of previous
! layer thickness. Strictly speaking, only necessary for top snow layer, but doing
! it for all snow layers will catch problems with older initial files.
! Layer interfaces (zi) and node depths (z) do not need adjustment here because they
! are adjusted in CombineSnowLayers and are not used up to that point.
do j = -nlevsno+1, 0
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
dz(c,j) = max(dz(c,j),h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice)
end if
end do
end do
do fc = 1, num_snowc
c = filter_snowc(fc)
! Qout from snow bottom
qflx_top_soil(c) = qout(c) / dtime
end do
do fc = 1, num_nosnowc
c = filter_nosnowc(fc)
qflx_top_soil(c) = qflx_rain_grnd(c) + qflx_snomelt(c)
end do
! set aerosol deposition fluxes from forcing array
! The forcing array is either set from an external file
! or from fluxes received from the atmosphere model
do c = lbc,ubc
g = cgridcell(c)
flx_bc_dep_dry(c) = forc_aer(g,1) + forc_aer(g,2)
flx_bc_dep_wet(c) = forc_aer(g,3)
flx_bc_dep_phi(c) = forc_aer(g,1) + forc_aer(g,3)
flx_bc_dep_pho(c) = forc_aer(g,2)
flx_bc_dep(c) = forc_aer(g,1) + forc_aer(g,2) + forc_aer(g,3)
flx_oc_dep_dry(c) = forc_aer(g,4) + forc_aer(g,5)
flx_oc_dep_wet(c) = forc_aer(g,6)
flx_oc_dep_phi(c) = forc_aer(g,4) + forc_aer(g,6)
flx_oc_dep_pho(c) = forc_aer(g,5)
flx_oc_dep(c) = forc_aer(g,4) + forc_aer(g,5) + forc_aer(g,6)
flx_dst_dep_wet1(c) = forc_aer(g,7)
flx_dst_dep_dry1(c) = forc_aer(g,8)
flx_dst_dep_wet2(c) = forc_aer(g,9)
flx_dst_dep_dry2(c) = forc_aer(g,10)
flx_dst_dep_wet3(c) = forc_aer(g,11)
flx_dst_dep_dry3(c) = forc_aer(g,12)
flx_dst_dep_wet4(c) = forc_aer(g,13)
flx_dst_dep_dry4(c) = forc_aer(g,14)
flx_dst_dep(c) = forc_aer(g,7) + forc_aer(g,8) + forc_aer(g,9) + &
forc_aer(g,10) + forc_aer(g,11) + forc_aer(g,12) + &
forc_aer(g,13) + forc_aer(g,14)
end do
! aerosol deposition fluxes into top layer
! This is done after the inter-layer fluxes so that some aerosol
! is in the top layer after deposition, and is not immediately
! washed out before radiative calculations are done
do fc = 1, num_snowc
c = filter_snowc(fc)
mss_bcphi(c,snl(c)+1) = mss_bcphi(c,snl(c)+1) + (flx_bc_dep_phi(c)*dtime)
mss_bcpho(c,snl(c)+1) = mss_bcpho(c,snl(c)+1) + (flx_bc_dep_pho(c)*dtime)
mss_ocphi(c,snl(c)+1) = mss_ocphi(c,snl(c)+1) + (flx_oc_dep_phi(c)*dtime)
mss_ocpho(c,snl(c)+1) = mss_ocpho(c,snl(c)+1) + (flx_oc_dep_pho(c)*dtime)
mss_dst1(c,snl(c)+1) = mss_dst1(c,snl(c)+1) + (flx_dst_dep_dry1(c) + flx_dst_dep_wet1(c))*dtime
mss_dst2(c,snl(c)+1) = mss_dst2(c,snl(c)+1) + (flx_dst_dep_dry2(c) + flx_dst_dep_wet2(c))*dtime
mss_dst3(c,snl(c)+1) = mss_dst3(c,snl(c)+1) + (flx_dst_dep_dry3(c) + flx_dst_dep_wet3(c))*dtime
mss_dst4(c,snl(c)+1) = mss_dst4(c,snl(c)+1) + (flx_dst_dep_dry4(c) + flx_dst_dep_wet4(c))*dtime
end do
end subroutine SnowWater
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SnowCompaction
!
! !INTERFACE:
subroutine SnowCompaction(lbc, ubc, num_snowc, filter_snowc) 1,4
!
! !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.
!
! !USES:
use clmtype
use clm_time_manager, only : get_step_size
use clm_varcon , only : denice, denh2o, tfrz, istice_mec
!
! !ARGUMENTS:
implicit none
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
!
! !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
! 2/29/08, David Lawrence: Revised snow overburden to be include 0.5 weight of current layer
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in scalars
!
integer, pointer :: snl(:) !number of snow layers
!
! local pointers to implicit in arguments
!
integer, pointer :: imelt(:,:) !flag for melting (=1), freezing (=2), Not=0
real(r8), pointer :: frac_iceold(:,:) !fraction of ice relative to the tot water
real(r8), pointer :: t_soisno(:,:) !soil temperature (Kelvin)
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
!
! local pointers to implicit inout arguments
!
real(r8), pointer :: dz(:,:) !layer depth (m)
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
integer :: j, l, c, fc ! indices
real(r8):: dtime ! land model time step (sec)
real(r8), parameter :: c2 = 23.e-3_r8 ! [m3/kg]
real(r8), parameter :: c3 = 2.777e-6_r8 ! [1/s]
real(r8), parameter :: c4 = 0.04_r8 ! [1/K]
real(r8), parameter :: c5 = 2.0_r8 !
real(r8), parameter :: dm = 100.0_r8 ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
real(r8), parameter :: eta0 = 9.e+5_r8 ! 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]
integer, pointer :: clandunit(:) !landunit index for each column
integer, pointer :: ltype(:) !landunit type
!-----------------------------------------------------------------------
! Assign local pointers to derived subtypes (column-level)
snl => clm3%g%l%c%cps%snl
dz => clm3%g%l%c%cps%dz
imelt => clm3%g%l%c%cps%imelt
frac_iceold => clm3%g%l%c%cps%frac_iceold
t_soisno => clm3%g%l%c%ces%t_soisno
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
clandunit => clm3%g%l%c%landunit
ltype => clm3%g%l%itype
! Get time step
dtime = get_step_size()
! Begin calculation - note that the following column loops are only invoked if snl(c) < 0
burden(:) = 0._r8
do j = -nlevsno+1, 0
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._r8 - (h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o) / dz(c,j)
! If void is negative, then increase dz such that void = 0.
! This should be done for any landunit, but for now is done only for glacier_mec 1andunits.
l = clandunit(c)
if (ltype(l)==istice_mec .and. void < 0._r8) then
dz(c,j) = h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o
void = 0._r8
endif
! Allow compaction only for non-saturated node and higher ice lens node.
if (void > 0.001_r8 .and. h2osoi_ice(c,j) > .1_r8) 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_r8*(bi-dm))
! Liquid water term
if (h2osoi_liq(c,j) > 0.01_r8*dz(c,j)) ddz1=ddz1*c5
! Compaction due to overburden
ddz2 = -(burden(c)+wx/2._r8)*exp(-0.08_r8*td - c2*bi)/eta0
! Compaction occurring during melt
if (imelt(c,j) == 1) then
ddz3 = - 1._r8/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
! Limit compaction to no less than fully saturated layer thickness
dz(c,j) = max(dz(c,j) * (1._r8+pdzdtc*dtime),h2osoi_ice(c,j)/denice &
+ h2osoi_liq(c,j)/denh2o)
end if
! Pressure of overlying snow
burden(c) = burden(c) + wx
end if
end do
end do
end subroutine SnowCompaction
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: CombineSnowLayers
!
! !INTERFACE:
subroutine CombineSnowLayers(lbc, ubc, num_snowc, filter_snowc) 1,3
!
! !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.
!
! !USES:
use clmtype
use clm_varcon, only : istsoil, isturb
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbc, ubc ! column bounds
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
!
! !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.
! 03/28/08, Mark Flanner: Added aerosol masses and snow grain radius
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
integer, pointer :: clandunit(:) !landunit index for each column
integer, pointer :: ltype(:) !landunit type
!
! local pointers to implicit inout arguments
!
integer , pointer :: snl(:) !number of snow layers
real(r8), pointer :: h2osno(:) !snow water (mm H2O)
real(r8), pointer :: snowdp(:) !snow height (m)
real(r8), pointer :: dz(:,:) !layer depth (m)
real(r8), pointer :: zi(:,:) !interface level below a "z" level (m)
real(r8), pointer :: t_soisno(:,:) !soil temperature (Kelvin)
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
!
! local pointers to implicit out arguments
!
real(r8), pointer :: z(:,:) ! layer thickness (m)
real(r8), pointer :: mss_bcphi(:,:) ! hydrophilic BC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_bcpho(:,:) ! hydrophobic BC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_ocphi(:,:) ! hydrophilic OC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_ocpho(:,:) ! hydrophobic OC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst1(:,:) ! dust species 1 mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst2(:,:) ! dust species 2 mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst3(:,:) ! dust species 3 mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst4(:,:) ! dust species 4 mass in snow (col,lyr) [kg]
real(r8), pointer :: snw_rds(:,:) ! effective snow grain radius (col,lyr) [microns, m^-6]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
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_r8, 0.015_r8, 0.025_r8, 0.055_r8, 0.115_r8/
!-----------------------------------------------------------------------
! Assign local pointers to derived subtypes (landunit-level)
ltype => clm3%g%l%itype
! Assign local pointers to derived subtypes (column-level)
clandunit => clm3%g%l%c%landunit
snl => clm3%g%l%c%cps%snl
snowdp => clm3%g%l%c%cps%snowdp
h2osno => clm3%g%l%c%cws%h2osno
dz => clm3%g%l%c%cps%dz
zi => clm3%g%l%c%cps%zi
z => clm3%g%l%c%cps%z
t_soisno => clm3%g%l%c%ces%t_soisno
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
mss_bcphi => clm3%g%l%c%cps%mss_bcphi
mss_bcpho => clm3%g%l%c%cps%mss_bcpho
mss_ocphi => clm3%g%l%c%cps%mss_ocphi
mss_ocpho => clm3%g%l%c%cps%mss_ocpho
mss_dst1 => clm3%g%l%c%cps%mss_dst1
mss_dst2 => clm3%g%l%c%cps%mss_dst2
mss_dst3 => clm3%g%l%c%cps%mss_dst3
mss_dst4 => clm3%g%l%c%cps%mss_dst4
snw_rds => clm3%g%l%c%cps%snw_rds
! Check the mass of ice lens of snow, when the total is less than a small value,
! combine it with the underlying neighbor.
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_r8) then
if (ltype(l) == istsoil .or. ltype(l)==isturb) 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)
if (j /= 0) dz(c,j+1) = dz(c,j+1) + dz(c,j)
! NOTE: Temperature, and similarly snw_rds, of the
! underlying snow layer are NOT adjusted in this case.
! Because the layer being eliminated has a small mass,
! this should not make a large difference, but it
! would be more thorough to do so.
if (j /= 0) then
mss_bcphi(c,j+1) = mss_bcphi(c,j+1) + mss_bcphi(c,j)
mss_bcpho(c,j+1) = mss_bcpho(c,j+1) + mss_bcpho(c,j)
mss_ocphi(c,j+1) = mss_ocphi(c,j+1) + mss_ocphi(c,j)
mss_ocpho(c,j+1) = mss_ocpho(c,j+1) + mss_ocpho(c,j)
mss_dst1(c,j+1) = mss_dst1(c,j+1) + mss_dst1(c,j)
mss_dst2(c,j+1) = mss_dst2(c,j+1) + mss_dst2(c,j)
mss_dst3(c,j+1) = mss_dst3(c,j+1) + mss_dst3(c,j)
mss_dst4(c,j+1) = mss_dst4(c,j+1) + mss_dst4(c,j)
end if
else if (ltype(l) /= istsoil .and. ltype(l) /= isturb .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)
dz(c,j+1) = dz(c,j+1) + dz(c,j)
mss_bcphi(c,j+1) = mss_bcphi(c,j+1) + mss_bcphi(c,j)
mss_bcpho(c,j+1) = mss_bcpho(c,j+1) + mss_bcpho(c,j)
mss_ocphi(c,j+1) = mss_ocphi(c,j+1) + mss_ocphi(c,j)
mss_ocpho(c,j+1) = mss_ocpho(c,j+1) + mss_ocpho(c,j)
mss_dst1(c,j+1) = mss_dst1(c,j+1) + mss_dst1(c,j)
mss_dst2(c,j+1) = mss_dst2(c,j+1) + mss_dst2(c,j)
mss_dst3(c,j+1) = mss_dst3(c,j+1) + mss_dst3(c,j)
mss_dst4(c,j+1) = mss_dst4(c,j+1) + mss_dst4(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)
mss_bcphi(c,i) = mss_bcphi(c,i-1)
mss_bcpho(c,i) = mss_bcpho(c,i-1)
mss_ocphi(c,i) = mss_ocphi(c,i-1)
mss_ocpho(c,i) = mss_ocpho(c,i-1)
mss_dst1(c,i) = mss_dst1(c,i-1)
mss_dst2(c,i) = mss_dst2(c,i-1)
mss_dst3(c,i) = mss_dst3(c,i-1)
mss_dst4(c,i) = mss_dst4(c,i-1)
snw_rds(c,i) = snw_rds(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
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 = -nlevsno+1,0
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.
do fc = 1, num_snowc
c = filter_snowc(fc)
l = clandunit(c)
if (snowdp(c) < 0.01_r8 .and. snowdp(c) > 0._r8) then
snl(c) = 0
h2osno(c) = zwice(c)
mss_bcphi(c,:) = 0._r8
mss_bcpho(c,:) = 0._r8
mss_ocphi(c,:) = 0._r8
mss_ocpho(c,:) = 0._r8
mss_dst1(c,:) = 0._r8
mss_dst2(c,:) = 0._r8
mss_dst3(c,:) = 0._r8
mss_dst4(c,:) = 0._r8
if (h2osno(c) <= 0._r8) snowdp(c) = 0._r8
if (ltype(l) == istsoil .or. ltype(l) == isturb) then
h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c)
end if
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
! this should be included in 'Combo' for consistency,
! but functionally it is the same to do it here
mss_bcphi(c,j)=mss_bcphi(c,j)+mss_bcphi(c,l)
mss_bcpho(c,j)=mss_bcpho(c,j)+mss_bcpho(c,l)
mss_ocphi(c,j)=mss_ocphi(c,j)+mss_ocphi(c,l)
mss_ocpho(c,j)=mss_ocpho(c,j)+mss_ocpho(c,l)
mss_dst1(c,j)=mss_dst1(c,j)+mss_dst1(c,l)
mss_dst2(c,j)=mss_dst2(c,j)+mss_dst2(c,l)
mss_dst3(c,j)=mss_dst3(c,j)+mss_dst3(c,l)
mss_dst4(c,j)=mss_dst4(c,j)+mss_dst4(c,l)
! mass-weighted combination of effective grain size:
snw_rds(c,j) = (snw_rds(c,j)*(h2osoi_liq(c,j)+h2osoi_ice(c,j)) + &
snw_rds(c,l)*(h2osoi_liq(c,l)+h2osoi_ice(c,l))) / &
(h2osoi_liq(c,j)+h2osoi_ice(c,j)+h2osoi_liq(c,l)+h2osoi_ice(c,l))
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)
mss_bcphi(c,k) = mss_bcphi(c,k-1)
mss_bcpho(c,k) = mss_bcpho(c,k-1)
mss_ocphi(c,k) = mss_ocphi(c,k-1)
mss_ocpho(c,k) = mss_ocpho(c,k-1)
mss_dst1(c,k) = mss_dst1(c,k-1)
mss_dst2(c,k) = mss_dst2(c,k-1)
mss_dst3(c,k) = mss_dst3(c,k-1)
mss_dst4(c,k) = mss_dst4(c,k-1)
snw_rds(c,k) = snw_rds(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, -nlevsno+1, -1
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c) + 1) then
z(c,j) = zi(c,j) - 0.5_r8*dz(c,j)
zi(c,j-1) = zi(c,j) - dz(c,j)
end if
end do
end do
end subroutine CombineSnowLayers
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: DivideSnowLayers
!
! !INTERFACE:
subroutine DivideSnowLayers(lbc, ubc, num_snowc, filter_snowc) 1,6
!
! !DESCRIPTION:
! Subdivides snow layers if they exceed their prescribed maximum thickness.
!
! !USES:
use clmtype
use clm_varcon, only : tfrz
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbc, ubc ! column bounds
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
!
! !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.
! 2/29/08, David Lawrence: Snowpack T profile maintained during layer splitting
! 03/28/08, Mark Flanner: Added aerosol masses and snow grain radius
!
! !LOCAL VARIABLES:
!
! local pointers to implicit inout arguments
!
integer , pointer :: snl(:) !number of snow layers
real(r8), pointer :: dz(:,:) !layer depth (m)
real(r8), pointer :: zi(:,:) !interface level below a "z" level (m)
real(r8), pointer :: t_soisno(:,:) !soil temperature (Kelvin)
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
!
! local pointers to implicit out arguments
!
real(r8), pointer :: z(:,:) ! layer thickness (m)
real(r8), pointer :: mss_bcphi(:,:) ! hydrophilic BC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_bcpho(:,:) ! hydrophobic BC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_ocphi(:,:) ! hydrophilic OC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_ocpho(:,:) ! hydrophobic OC mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst1(:,:) ! dust species 1 mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst2(:,:) ! dust species 2 mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst3(:,:) ! dust species 3 mass in snow (col,lyr) [kg]
real(r8), pointer :: mss_dst4(:,:) ! dust species 4 mass in snow (col,lyr) [kg]
real(r8), pointer :: snw_rds(:,:) ! effective snow grain radius (col,lyr) [microns, m^-6]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
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,nlevsno) ! Snow layer thickness [m]
real(r8) :: swice(lbc:ubc,nlevsno) ! Partial volume of ice [m3/m3]
real(r8) :: swliq(lbc:ubc,nlevsno) ! Partial volume of liquid water [m3/m3]
real(r8) :: tsno(lbc:ubc ,nlevsno) ! Nodel temperature [K]
real(r8) :: zwice ! temporary
real(r8) :: zwliq ! temporary
real(r8) :: propor ! temporary
real(r8) :: dtdz ! temporary
! temporary variables mimicking the structure of other layer division variables
real(r8) :: mbc_phi(lbc:ubc,nlevsno) ! mass of BC in each snow layer
real(r8) :: zmbc_phi ! temporary
real(r8) :: mbc_pho(lbc:ubc,nlevsno) ! mass of BC in each snow layer
real(r8) :: zmbc_pho ! temporary
real(r8) :: moc_phi(lbc:ubc,nlevsno) ! mass of OC in each snow layer
real(r8) :: zmoc_phi ! temporary
real(r8) :: moc_pho(lbc:ubc,nlevsno) ! mass of OC in each snow layer
real(r8) :: zmoc_pho ! temporary
real(r8) :: mdst1(lbc:ubc,nlevsno) ! mass of dust 1 in each snow layer
real(r8) :: zmdst1 ! temporary
real(r8) :: mdst2(lbc:ubc,nlevsno) ! mass of dust 2 in each snow layer
real(r8) :: zmdst2 ! temporary
real(r8) :: mdst3(lbc:ubc,nlevsno) ! mass of dust 3 in each snow layer
real(r8) :: zmdst3 ! temporary
real(r8) :: mdst4(lbc:ubc,nlevsno) ! mass of dust 4 in each snow layer
real(r8) :: zmdst4 ! temporary
real(r8) :: rds(lbc:ubc,nlevsno)
!-----------------------------------------------------------------------
! Assign local pointers to derived subtype components (column-level)
snl => clm3%g%l%c%cps%snl
dz => clm3%g%l%c%cps%dz
zi => clm3%g%l%c%cps%zi
z => clm3%g%l%c%cps%z
t_soisno => clm3%g%l%c%ces%t_soisno
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
mss_bcphi => clm3%g%l%c%cps%mss_bcphi
mss_bcpho => clm3%g%l%c%cps%mss_bcpho
mss_ocphi => clm3%g%l%c%cps%mss_ocphi
mss_ocpho => clm3%g%l%c%cps%mss_ocpho
mss_dst1 => clm3%g%l%c%cps%mss_dst1
mss_dst2 => clm3%g%l%c%cps%mss_dst2
mss_dst3 => clm3%g%l%c%cps%mss_dst3
mss_dst4 => clm3%g%l%c%cps%mss_dst4
snw_rds => clm3%g%l%c%cps%snw_rds
! Begin calculation - note that the following column loops are only invoked
! for snow-covered columns
do j = 1,nlevsno
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))
mbc_phi(c,j) = mss_bcphi(c,j+snl(c))
mbc_pho(c,j) = mss_bcpho(c,j+snl(c))
moc_phi(c,j) = mss_ocphi(c,j+snl(c))
moc_pho(c,j) = mss_ocpho(c,j+snl(c))
mdst1(c,j) = mss_dst1(c,j+snl(c))
mdst2(c,j) = mss_dst2(c,j+snl(c))
mdst3(c,j) = mss_dst3(c,j+snl(c))
mdst4(c,j) = mss_dst4(c,j+snl(c))
rds(c,j) = snw_rds(c,j+snl(c))
end if
end do
end do
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_r8) then
msno = 2
dzsno(c,1) = dzsno(c,1)/2._r8
swice(c,1) = swice(c,1)/2._r8
swliq(c,1) = swliq(c,1)/2._r8
dzsno(c,2) = dzsno(c,1)
swice(c,2) = swice(c,1)
swliq(c,2) = swliq(c,1)
tsno(c,2) = tsno(c,1)
mbc_phi(c,1) = mbc_phi(c,1)/2._r8
mbc_phi(c,2) = mbc_phi(c,1)
mbc_pho(c,1) = mbc_pho(c,1)/2._r8
mbc_pho(c,2) = mbc_pho(c,1)
moc_phi(c,1) = moc_phi(c,1)/2._r8
moc_phi(c,2) = moc_phi(c,1)
moc_pho(c,1) = moc_pho(c,1)/2._r8
moc_pho(c,2) = moc_pho(c,1)
mdst1(c,1) = mdst1(c,1)/2._r8
mdst1(c,2) = mdst1(c,1)
mdst2(c,1) = mdst2(c,1)/2._r8
mdst2(c,2) = mdst2(c,1)
mdst3(c,1) = mdst3(c,1)/2._r8
mdst3(c,2) = mdst3(c,1)
mdst4(c,1) = mdst4(c,1)/2._r8
mdst4(c,2) = mdst4(c,1)
rds(c,2) = rds(c,1)
end if
end if
if (msno > 1) then
if (dzsno(c,1) > 0.02_r8) then
drr = dzsno(c,1) - 0.02_r8
propor = drr/dzsno(c,1)
zwice = propor*swice(c,1)
zwliq = propor*swliq(c,1)
zmbc_phi = propor*mbc_phi(c,1)
zmbc_pho = propor*mbc_pho(c,1)
zmoc_phi = propor*moc_phi(c,1)
zmoc_pho = propor*moc_pho(c,1)
zmdst1 = propor*mdst1(c,1)
zmdst2 = propor*mdst2(c,1)
zmdst3 = propor*mdst3(c,1)
zmdst4 = propor*mdst4(c,1)
propor = 0.02_r8/dzsno(c,1)
swice(c,1) = propor*swice(c,1)
swliq(c,1) = propor*swliq(c,1)
mbc_phi(c,1) = propor*mbc_phi(c,1)
mbc_pho(c,1) = propor*mbc_pho(c,1)
moc_phi(c,1) = propor*moc_phi(c,1)
moc_pho(c,1) = propor*moc_pho(c,1)
mdst1(c,1) = propor*mdst1(c,1)
mdst2(c,1) = propor*mdst2(c,1)
mdst3(c,1) = propor*mdst3(c,1)
mdst4(c,1) = propor*mdst4(c,1)
dzsno(c,1) = 0.02_r8
mbc_phi(c,2) = mbc_phi(c,2)+zmbc_phi ! (combo)
mbc_pho(c,2) = mbc_pho(c,2)+zmbc_pho ! (combo)
moc_phi(c,2) = moc_phi(c,2)+zmoc_phi ! (combo)
moc_pho(c,2) = moc_pho(c,2)+zmoc_pho ! (combo)
mdst1(c,2) = mdst1(c,2)+zmdst1 ! (combo)
mdst2(c,2) = mdst2(c,2)+zmdst2 ! (combo)
mdst3(c,2) = mdst3(c,2)+zmdst3 ! (combo)
mdst4(c,2) = mdst4(c,2)+zmdst4 ! (combo)
rds(c,2) = rds(c,1) ! (combo)
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_r8) then
msno = 3
dtdz = (tsno(c,1) - tsno(c,2))/((dzsno(c,1)+dzsno(c,2))/2._r8)
dzsno(c,2) = dzsno(c,2)/2._r8
swice(c,2) = swice(c,2)/2._r8
swliq(c,2) = swliq(c,2)/2._r8
dzsno(c,3) = dzsno(c,2)
swice(c,3) = swice(c,2)
swliq(c,3) = swliq(c,2)
tsno(c,3) = tsno(c,2) - dtdz*dzsno(c,2)/2._r8
if (tsno(c,3) >= tfrz) then
tsno(c,3) = tsno(c,2)
else
tsno(c,2) = tsno(c,2) + dtdz*dzsno(c,2)/2._r8
endif
mbc_phi(c,2) = mbc_phi(c,2)/2._r8
mbc_phi(c,3) = mbc_phi(c,2)
mbc_pho(c,2) = mbc_pho(c,2)/2._r8
mbc_pho(c,3) = mbc_pho(c,2)
moc_phi(c,2) = moc_phi(c,2)/2._r8
moc_phi(c,3) = moc_phi(c,2)
moc_pho(c,2) = moc_pho(c,2)/2._r8
moc_pho(c,3) = moc_pho(c,2)
mdst1(c,2) = mdst1(c,2)/2._r8
mdst1(c,3) = mdst1(c,2)
mdst2(c,2) = mdst2(c,2)/2._r8
mdst2(c,3) = mdst2(c,2)
mdst3(c,2) = mdst3(c,2)/2._r8
mdst3(c,3) = mdst3(c,2)
mdst4(c,2) = mdst4(c,2)/2._r8
mdst4(c,3) = mdst4(c,2)
rds(c,3) = rds(c,2)
end if
end if
end if
if (msno > 2) then
if (dzsno(c,2) > 0.05_r8) then
drr = dzsno(c,2) - 0.05_r8
propor = drr/dzsno(c,2)
zwice = propor*swice(c,2)
zwliq = propor*swliq(c,2)
zmbc_phi = propor*mbc_phi(c,2)
zmbc_pho = propor*mbc_pho(c,2)
zmoc_phi = propor*moc_phi(c,2)
zmoc_pho = propor*moc_pho(c,2)
zmdst1 = propor*mdst1(c,2)
zmdst2 = propor*mdst2(c,2)
zmdst3 = propor*mdst3(c,2)
zmdst4 = propor*mdst4(c,2)
propor = 0.05_r8/dzsno(c,2)
swice(c,2) = propor*swice(c,2)
swliq(c,2) = propor*swliq(c,2)
mbc_phi(c,2) = propor*mbc_phi(c,2)
mbc_pho(c,2) = propor*mbc_pho(c,2)
moc_phi(c,2) = propor*moc_phi(c,2)
moc_pho(c,2) = propor*moc_pho(c,2)
mdst1(c,2) = propor*mdst1(c,2)
mdst2(c,2) = propor*mdst2(c,2)
mdst3(c,2) = propor*mdst3(c,2)
mdst4(c,2) = propor*mdst4(c,2)
dzsno(c,2) = 0.05_r8
mbc_phi(c,3) = mbc_phi(c,3)+zmbc_phi ! (combo)
mbc_pho(c,3) = mbc_pho(c,3)+zmbc_pho ! (combo)
moc_phi(c,3) = moc_phi(c,3)+zmoc_phi ! (combo)
moc_pho(c,3) = moc_pho(c,3)+zmoc_pho ! (combo)
mdst1(c,3) = mdst1(c,3)+zmdst1 ! (combo)
mdst2(c,3) = mdst2(c,3)+zmdst2 ! (combo)
mdst3(c,3) = mdst3(c,3)+zmdst3 ! (combo)
mdst4(c,3) = mdst4(c,3)+zmdst4 ! (combo)
rds(c,3) = rds(c,2) ! (combo)
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_r8) then
msno = 4
dtdz = (tsno(c,2) - tsno(c,3))/((dzsno(c,2)+dzsno(c,3))/2._r8)
dzsno(c,3) = dzsno(c,3)/2._r8
swice(c,3) = swice(c,3)/2._r8
swliq(c,3) = swliq(c,3)/2._r8
dzsno(c,4) = dzsno(c,3)
swice(c,4) = swice(c,3)
swliq(c,4) = swliq(c,3)
tsno(c,4) = tsno(c,3) - dtdz*dzsno(c,3)/2._r8
if (tsno(c,4) >= tfrz) then
tsno(c,4) = tsno(c,3)
else
tsno(c,3) = tsno(c,3) + dtdz*dzsno(c,3)/2._r8
endif
mbc_phi(c,3) = mbc_phi(c,3)/2._r8
mbc_phi(c,4) = mbc_phi(c,3)
mbc_pho(c,3) = mbc_pho(c,3)/2._r8
mbc_pho(c,4) = mbc_pho(c,3)
moc_phi(c,3) = moc_phi(c,3)/2._r8
moc_phi(c,4) = moc_phi(c,3)
moc_pho(c,3) = moc_pho(c,3)/2._r8
moc_pho(c,4) = moc_pho(c,3)
mdst1(c,3) = mdst1(c,3)/2._r8
mdst1(c,4) = mdst1(c,3)
mdst2(c,3) = mdst2(c,3)/2._r8
mdst2(c,4) = mdst2(c,3)
mdst3(c,3) = mdst3(c,3)/2._r8
mdst3(c,4) = mdst3(c,3)
mdst4(c,3) = mdst4(c,3)/2._r8
mdst4(c,4) = mdst4(c,3)
rds(c,4) = rds(c,3)
end if
end if
end if
if (msno > 3) then
if (dzsno(c,3) > 0.11_r8) then
drr = dzsno(c,3) - 0.11_r8
propor = drr/dzsno(c,3)
zwice = propor*swice(c,3)
zwliq = propor*swliq(c,3)
zmbc_phi = propor*mbc_phi(c,3)
zmbc_pho = propor*mbc_pho(c,3)
zmoc_phi = propor*moc_phi(c,3)
zmoc_pho = propor*moc_pho(c,3)
zmdst1 = propor*mdst1(c,3)
zmdst2 = propor*mdst2(c,3)
zmdst3 = propor*mdst3(c,3)
zmdst4 = propor*mdst4(c,3)
propor = 0.11_r8/dzsno(c,3)
swice(c,3) = propor*swice(c,3)
swliq(c,3) = propor*swliq(c,3)
mbc_phi(c,3) = propor*mbc_phi(c,3)
mbc_pho(c,3) = propor*mbc_pho(c,3)
moc_phi(c,3) = propor*moc_phi(c,3)
moc_pho(c,3) = propor*moc_pho(c,3)
mdst1(c,3) = propor*mdst1(c,3)
mdst2(c,3) = propor*mdst2(c,3)
mdst3(c,3) = propor*mdst3(c,3)
mdst4(c,3) = propor*mdst4(c,3)
dzsno(c,3) = 0.11_r8
mbc_phi(c,4) = mbc_phi(c,4)+zmbc_phi ! (combo)
mbc_pho(c,4) = mbc_pho(c,4)+zmbc_pho ! (combo)
moc_phi(c,4) = moc_phi(c,4)+zmoc_phi ! (combo)
moc_pho(c,4) = moc_pho(c,4)+zmoc_pho ! (combo)
mdst1(c,4) = mdst1(c,4)+zmdst1 ! (combo)
mdst2(c,4) = mdst2(c,4)+zmdst2 ! (combo)
mdst3(c,4) = mdst3(c,4)+zmdst3 ! (combo)
mdst4(c,4) = mdst4(c,4)+zmdst4 ! (combo)
rds(c,4) = rds(c,3) ! (combo)
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_r8) then
msno = 5
dtdz = (tsno(c,3) - tsno(c,4))/((dzsno(c,3)+dzsno(c,4))/2._r8)
dzsno(c,4) = dzsno(c,4)/2._r8
swice(c,4) = swice(c,4)/2._r8
swliq(c,4) = swliq(c,4)/2._r8
dzsno(c,5) = dzsno(c,4)
swice(c,5) = swice(c,4)
swliq(c,5) = swliq(c,4)
tsno(c,5) = tsno(c,4) - dtdz*dzsno(c,4)/2._r8
if (tsno(c,5) >= tfrz) then
tsno(c,5) = tsno(c,4)
else
tsno(c,4) = tsno(c,4) + dtdz*dzsno(c,4)/2._r8
endif
mbc_phi(c,4) = mbc_phi(c,4)/2._r8
mbc_phi(c,5) = mbc_phi(c,4)
mbc_pho(c,4) = mbc_pho(c,4)/2._r8
mbc_pho(c,5) = mbc_pho(c,4)
moc_phi(c,4) = moc_phi(c,4)/2._r8
moc_phi(c,5) = moc_phi(c,4)
moc_pho(c,4) = moc_pho(c,4)/2._r8
moc_pho(c,5) = moc_pho(c,4)
mdst1(c,4) = mdst1(c,4)/2._r8
mdst1(c,5) = mdst1(c,4)
mdst2(c,4) = mdst2(c,4)/2._r8
mdst2(c,5) = mdst2(c,4)
mdst3(c,4) = mdst3(c,4)/2._r8
mdst3(c,5) = mdst3(c,4)
mdst4(c,4) = mdst4(c,4)/2._r8
mdst4(c,5) = mdst4(c,4)
rds(c,5) = rds(c,4)
end if
end if
end if
if (msno > 4) then
if (dzsno(c,4) > 0.23_r8) then
drr = dzsno(c,4) - 0.23_r8
propor = drr/dzsno(c,4)
zwice = propor*swice(c,4)
zwliq = propor*swliq(c,4)
zmbc_phi = propor*mbc_phi(c,4)
zmbc_pho = propor*mbc_pho(c,4)
zmoc_phi = propor*moc_phi(c,4)
zmoc_pho = propor*moc_pho(c,4)
zmdst1 = propor*mdst1(c,4)
zmdst2 = propor*mdst2(c,4)
zmdst3 = propor*mdst3(c,4)
zmdst4 = propor*mdst4(c,4)
propor = 0.23_r8/dzsno(c,4)
swice(c,4) = propor*swice(c,4)
swliq(c,4) = propor*swliq(c,4)
mbc_phi(c,4) = propor*mbc_phi(c,4)
mbc_pho(c,4) = propor*mbc_pho(c,4)
moc_phi(c,4) = propor*moc_phi(c,4)
moc_pho(c,4) = propor*moc_pho(c,4)
mdst1(c,4) = propor*mdst1(c,4)
mdst2(c,4) = propor*mdst2(c,4)
mdst3(c,4) = propor*mdst3(c,4)
mdst4(c,4) = propor*mdst4(c,4)
dzsno(c,4) = 0.23_r8
mbc_phi(c,5) = mbc_phi(c,5)+zmbc_phi ! (combo)
mbc_pho(c,5) = mbc_pho(c,5)+zmbc_pho ! (combo)
moc_phi(c,5) = moc_phi(c,5)+zmoc_phi ! (combo)
moc_pho(c,5) = moc_pho(c,5)+zmoc_pho ! (combo)
mdst1(c,5) = mdst1(c,5)+zmdst1 ! (combo)
mdst2(c,5) = mdst2(c,5)+zmdst2 ! (combo)
mdst3(c,5) = mdst3(c,5)+zmdst3 ! (combo)
mdst4(c,5) = mdst4(c,5)+zmdst4 ! (combo)
rds(c,5) = rds(c,4) ! (combo)
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 = -nlevsno+1,0
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))
mss_bcphi(c,j) = mbc_phi(c,j-snl(c))
mss_bcpho(c,j) = mbc_pho(c,j-snl(c))
mss_ocphi(c,j) = moc_phi(c,j-snl(c))
mss_ocpho(c,j) = moc_pho(c,j-snl(c))
mss_dst1(c,j) = mdst1(c,j-snl(c))
mss_dst2(c,j) = mdst2(c,j-snl(c))
mss_dst3(c,j) = mdst3(c,j-snl(c))
mss_dst4(c,j) = mdst4(c,j-snl(c))
snw_rds(c,j) = rds(c,j-snl(c))
end if
end do
end do
do j = 0, -nlevsno+1, -1
do fc = 1, num_snowc
c = filter_snowc(fc)
if (j >= snl(c)+1) then
z(c,j) = zi(c,j) - 0.5_r8*dz(c,j)
zi(c,j-1) = zi(c,j) - dz(c,j)
end if
end do
end do
end subroutine DivideSnowLayers
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Combo
!
! !INTERFACE:
subroutine Combo(dz, wliq, wice, t, dz2, wliq2, wice2, t2) 5,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:
use clm_varcon, only : cpice, cpliq, tfrz, hfus
!
! !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
!
!
! !LOCAL VARIABLES:
!EOP
!
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
tc = tfrz + (hc - hfus*wliqc) / (cpice*wicec + cpliq*wliqc)
dz = dzc
wice = wicec
wliq = wliqc
t = tc
end subroutine Combo
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: BuildSnowFilter
!
! !INTERFACE:
subroutine BuildSnowFilter(lbc, ubc, num_nolakec, filter_nolakec, & 2,1
num_snowc, filter_snowc, &
num_nosnowc, filter_nosnowc)
!
! !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(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
integer , pointer :: snl(:) ! number of snow layers
!
!
! !OTHER LOCAL VARIABLES:
!EOP
integer :: fc, c
!-----------------------------------------------------------------------
! Assign local pointers to derived subtype components (column-level)
snl => clm3%g%l%c%cps%snl
! 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
end module SnowHydrologyMod