module SurfaceResistanceMod #include "shr_assert.h" !----------------------------------------------------------------------- ! !DESCRIPTION: ! Module holding routines for calculation of surface resistances of the different tracers ! transported with BeTR. The surface here refers to water and soil, not including canopy ! ! !USES: use shr_kind_mod , only: r8 => shr_kind_r8 use shr_const_mod , only: SHR_CONST_TKFRZ use clm_varctl , only: iulog use SoilStateType , only: soilstate_type use WaterStateType, only: waterstate_type use TemperatureType , only : temperature_type implicit none save private integer :: soil_resis_method !choose the method for soil resistance calculation integer, parameter :: leepielke_1992 = 0 ! integer, parameter :: sl_14 = 1 ! ! !PUBLIC MEMBER FUNCTIONS: public :: calc_soilevap_resis public :: do_soilevap_beta, do_soil_resistance_sl14 ! public :: init_soil_resistance public :: soil_resistance_readNL character(len=*), parameter, private :: sourcefile = & __FILE__ ! ! !REVISION HISTORY: ! 6/25/2013 Created by Jinyun Tang !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !!$ subroutine init_soil_resistance() !!$ ! !!$ !DESCRIPTIONS !!$ ! initialize method for soil resis calculation !!$ ! !!$ ! !USES: !!$ use abortutils , only : endrun !!$ use fileutils , only : getavu, relavu !!$ use spmdMod , only : mpicom, masterproc !!$ use shr_mpi_mod , only : shr_mpi_bcast !!$ use clm_varctl , only : iulog, use_bedrock !!$ use controlMod , only : NLFilename !!$ use clm_nlUtilsMod , only : find_nlgroup_name !!$ !!$ ! !ARGUMENTS: !!$ !------------------------------------------------------------------------------ !!$ implicit none !!$ integer :: nu_nml ! unit for namelist file !!$ integer :: nml_error ! namelist i/o error flag !!$ character(*), parameter :: subName = "('init_soil_resistance')" !!$ !!$ !----------------------------------------------------------------------- !!$ !!$! MUST agree with name in namelist and read statement !!$ namelist /soil_resis_inparm/ soil_resis_method !!$ !!$ ! Default values for namelist !!$ !!$ soil_resis_method = sl_14 !!$ !!$ ! Read soil_resis namelist !!$ if (masterproc) then !!$ nu_nml = getavu() !!$ open( nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) !!$ call find_nlgroup_name(nu_nml, 'soil_resis_inparm', status=nml_error) !!$ if (nml_error == 0) then !!$ read(nu_nml, nml=soil_resis_inparm,iostat=nml_error) !!$ if (nml_error /= 0) then !!$ call endrun(subname // ':: ERROR reading soil_resis namelist') !!$ end if !!$ end if !!$ close(nu_nml) !!$ call relavu( nu_nml ) !!$ !!$ endif !!$ !!$ call shr_mpi_bcast(soil_resis_method, mpicom) !!$ !!$ if (masterproc) then !!$ write(iulog,*) ' ' !!$ write(iulog,*) 'soil_resis settings:' !!$ write(iulog,*) ' soil_resis_method = ',soil_resis_method !!$ endif !!$!scs !!$! soil_resis_method = leepielke_1992 !!$! soil_resis_method = sl_14 !!$!scs !!$ !!$ end subroutine init_soil_resistance !----------------------------------------------------------------------- subroutine soil_resistance_readNL(NLFilename) ! !DESCRIPTIONS ! Read the namelist for soil resistance method ! ! !USES: use abortutils , only : endrun use fileutils , only : getavu, relavu use spmdMod , only : mpicom, masterproc use shr_mpi_mod , only : shr_mpi_bcast use clm_varctl , only : iulog use clm_nlUtilsMod , only : find_nlgroup_name ! !ARGUMENTS: !------------------------------------------------------------------------------ implicit none character(len=*), intent(IN) :: NLFilename ! Namelist filename integer :: nu_nml ! unit for namelist file integer :: nml_error ! namelist i/o error flag character(*), parameter :: subName = "('init_soil_resistance')" !----------------------------------------------------------------------- ! MUST agree with name in namelist and read statement namelist /soil_resis_inparm/ soil_resis_method ! Default values for namelist soil_resis_method = sl_14 ! Read soil_resis namelist if (masterproc) then nu_nml = getavu() open( nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) call find_nlgroup_name(nu_nml, 'soil_resis_inparm', status=nml_error) if (nml_error == 0) then read(nu_nml, nml=soil_resis_inparm,iostat=nml_error) if (nml_error /= 0) then call endrun(subname // ':: ERROR reading soil_resis namelist') end if else call endrun(subname // ':: ERROR reading soil_resis namelist') end if close(nu_nml) call relavu( nu_nml ) endif call shr_mpi_bcast(soil_resis_method, mpicom) if (masterproc) then write(iulog,*) ' ' write(iulog,*) 'soil_resis settings:' write(iulog,*) ' soil_resis_method = ',soil_resis_method endif end subroutine soil_resistance_readNL !------------------------------------------------------------------------------ subroutine calc_soilevap_resis(bounds, num_nolakec, filter_nolakec, & soilstate_inst, waterstate_inst, temperature_inst) ! ! DESCRIPTIONS ! compute the resis factor for soil evaporation calculation ! use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_PI use decompMod , only : bounds_type use ColumnType , only : col use LandunitType , only : lun use abortutils , only : endrun ! ! !ARGUMENTS: implicit none type(bounds_type) , intent(in) :: bounds ! bounds integer , intent(in) :: num_nolakec integer , intent(in) :: filter_nolakec(:) type(soilstate_type) , intent(inout) :: soilstate_inst type(waterstate_type) , intent(in) :: waterstate_inst type(temperature_type), intent(in) :: temperature_inst character(len=32) :: subname = 'calc_soilevap_resis' ! subroutine name associate( & soilbeta => soilstate_inst%soilbeta_col , & ! Output: [real(r8) (:)] factor that reduces ground evaporation dsl => soilstate_inst%dsl_col , & ! Output: [real(r8) (:)] soil dry surface layer thickness soilresis => soilstate_inst%soilresis_col & ! Output: [real(r8) (:)] soil evaporative resistance ) !select the right method and do the calculation select case (soil_resis_method) case (leepielke_1992) call calc_beta_leepielke1992(bounds, num_nolakec, filter_nolakec, & soilstate_inst, waterstate_inst, soilbeta(bounds%begc:bounds%endc)) case (sl_14) call calc_soil_resistance_sl14(bounds, num_nolakec, filter_nolakec, & soilstate_inst, waterstate_inst,temperature_inst, & dsl(bounds%begc:bounds%endc), soilresis(bounds%begc:bounds%endc)) case default call endrun(subname // ':: a soilevap resis function must be specified!') end select end associate end subroutine calc_soilevap_resis !------------------------------------------------------------------------------ subroutine calc_beta_leepielke1992(bounds, num_nolakec, filter_nolakec, & soilstate_inst, waterstate_inst, soilbeta) ! ! DESCRIPTION ! compute the lee-pielke beta factor to scal actual soil evaporation from potential evaporation ! ! USES use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_PI use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use decompMod , only : bounds_type use clm_varcon , only : denh2o, denice use landunit_varcon , only : istice_mec, istwet, istsoil, istcrop use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall use column_varcon , only : icol_road_imperv, icol_road_perv use ColumnType , only : col use LandunitType , only : lun ! implicit none type(bounds_type) , intent(in) :: bounds ! bounds integer , intent(in) :: num_nolakec integer , intent(in) :: filter_nolakec(:) type(soilstate_type) , intent(in) :: soilstate_inst type(waterstate_type) , intent(in) :: waterstate_inst real(r8) , intent(inout) :: soilbeta(bounds%begc:bounds%endc) !local variables real(r8) :: fac, fac_fc, wx !temporary variables integer :: c, l, fc !indices SHR_ASSERT_ALL((ubound(soilbeta) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) associate( & watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:)] volumetric soil water at saturation (porosity) watfc => soilstate_inst%watfc_col , & ! Input: [real(r8) (:,:)] volumetric soil water at field capacity h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:)] ice lens (kg/m2) h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:)] liquid water (kg/m2) frac_sno => waterstate_inst%frac_sno_col , & ! Input: [real(r8) (:)] fraction of ground covered by snow (0 to 1) frac_h2osfc => waterstate_inst%frac_h2osfc_col & ! Input: [real(r8) (:)] fraction of ground covered by surface water (0 to 1) ) do fc = 1,num_nolakec c = filter_nolakec(fc) l = col%landunit(c) if (lun%itype(l)/=istwet .AND. lun%itype(l)/=istice_mec) then if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then wx = (h2osoi_liq(c,1)/denh2o+h2osoi_ice(c,1)/denice)/col%dz(c,1) fac = min(1._r8, wx/watsat(c,1)) fac = max( fac, 0.01_r8 ) !! Lee and Pielke 1992 beta, added by K.Sakaguchi if (wx < watfc(c,1) ) then !when water content of ths top layer is less than that at F.C. fac_fc = min(1._r8, wx/watfc(c,1)) !eqn5.66 but divided by theta at field capacity fac_fc = max( fac_fc, 0.01_r8 ) ! modify soil beta by snow cover. soilbeta for snow surface is one soilbeta(c) = (1._r8-frac_sno(c)-frac_h2osfc(c)) & *0.25_r8*(1._r8 - cos(SHR_CONST_PI*fac_fc))**2._r8 & + frac_sno(c)+ frac_h2osfc(c) else !when water content of ths top layer is more than that at F.C. soilbeta(c) = 1._r8 end if else if (col%itype(c) == icol_road_perv) then soilbeta(c) = 0._r8 else if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall) then soilbeta(c) = 0._r8 else if (col%itype(c) == icol_roof .or. col%itype(c) == icol_road_imperv) then soilbeta(c) = 0._r8 endif else soilbeta(c) = 1._r8 endif enddo end associate end subroutine calc_beta_leepielke1992 !------------------------------------------------------------------------------ function do_soilevap_beta()result(lres) ! !DESCRIPTION ! return true if the moisture stress for soil evaporation is computed as beta factor ! otherwise false implicit none logical :: lres if(soil_resis_method==leepielke_1992)then lres=.true. else lres=.false. endif return end function do_soilevap_beta !------------------------------------------------------------------------------ subroutine calc_soil_resistance_sl14(bounds, num_nolakec, filter_nolakec, & soilstate_inst, waterstate_inst, temperature_inst, dsl, soilresis) ! ! DESCRIPTION ! compute the lee-pielke beta factor to scal actual soil evaporation from potential evaporation ! ! USES use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_PI use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use decompMod , only : bounds_type use clm_varcon , only : denh2o, denice use landunit_varcon , only : istice_mec, istwet, istsoil, istcrop use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall use column_varcon , only : icol_road_imperv, icol_road_perv use ColumnType , only : col use LandunitType , only : lun ! implicit none type(bounds_type) , intent(in) :: bounds ! bounds integer , intent(in) :: num_nolakec integer , intent(in) :: filter_nolakec(:) type(soilstate_type) , intent(in) :: soilstate_inst type(waterstate_type) , intent(in) :: waterstate_inst type(temperature_type), intent(in) :: temperature_inst real(r8) , intent(inout) :: dsl(bounds%begc:bounds%endc) real(r8) , intent(inout) :: soilresis(bounds%begc:bounds%endc) !local variables real(r8) :: aird, eps, dg, d0, vwc_liq real(r8) :: eff_por_top integer :: c, l, fc !indices SHR_ASSERT_ALL((ubound(dsl) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) SHR_ASSERT_ALL((ubound(soilresis) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:)] volumetric soil water at saturation (porosity) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) ! eff_porosity => soilstate_inst%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:)] ice lens (kg/m2) h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Input: [real(r8) (:,:)] liquid water (kg/m2) ) do fc = 1,num_nolakec c = filter_nolakec(fc) l = col%landunit(c) if (lun%itype(l)/=istwet .AND. lun%itype(l)/=istice_mec) then if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then vwc_liq = max(h2osoi_liq(c,1),1.0e-6_r8)/(dz(c,1)*denh2o) ! eff_porosity not calculated til SoilHydrology eff_por_top = max(0.01_r8,watsat(c,1)-min(watsat(c,1), h2osoi_ice(c,1)/(dz(c,1)*denice))) ! calculate diffusivity and air free pore space aird = watsat(c,1)*(sucsat(c,1)/1.e7_r8)**(1./bsw(c,1)) d0 = 2.12e-5*(t_soisno(c,1)/273.15)**1.75 ![Bitelli et al., JH, 08] eps = watsat(c,1) - aird dg = eps*d0*(eps/watsat(c,1))**(3._r8/max(3._r8,bsw(c,1))) ! dsl(c) = dzmm(c,1)*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) & ! try arbitrary scaling (not top layer thickness) ! dsl(c) = 15._r8*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) & dsl(c) = 15._r8*max(0.001_r8,(0.8*eff_por_top - vwc_liq)) & ! /max(0.001_r8,(watsat(c,1)- aird)) /max(0.001_r8,(0.8*watsat(c,1)- aird)) dsl(c)=max(dsl(c),0._r8) dsl(c)=min(dsl(c),200._r8) soilresis(c) = dsl(c)/(dg*eps*1.e3) + 20._r8 soilresis(c) = min(1.e6_r8,soilresis(c)) else if (col%itype(c) == icol_road_perv) then soilresis(c) = 1.e6_r8 else if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall) then soilresis(c) = 1.e6_r8 else if (col%itype(c) == icol_roof .or. col%itype(c) == icol_road_imperv) then soilresis(c) = 1.e6_r8 endif else soilresis(c) = 0._r8 endif enddo end associate end subroutine calc_soil_resistance_sl14 !------------------------------------------------------------------------------ function do_soil_resistance_sl14()result(lres) ! !DESCRIPTION ! return true if the soil evaporative resistance is computed using a DSL ! otherwise false implicit none logical :: lres if(soil_resis_method==sl_14)then lres=.true. else lres=.false. endif return end function do_soil_resistance_sl14 end module SurfaceResistanceMod