609 lines
26 KiB
Fortran
609 lines
26 KiB
Fortran
Module DryDepVelocity
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! Purpose:
|
|
! Deposition velocity (m/s)
|
|
!
|
|
! Method:
|
|
! This code simulates dry deposition velocities using the Wesely scheme.
|
|
! Details of this method can be found in:
|
|
!
|
|
! M.L Wesely. Parameterization of surface resistances to gaseous dry deposition
|
|
! in regional-scale numericl models. 1989. Atmospheric Environment vol.23 No.6
|
|
! pp. 1293-1304.
|
|
!
|
|
! In Wesely (1998) "the magnitude of the dry deposition velocity can be found
|
|
! as:
|
|
!
|
|
! |vd|=(ra+rb+rc)^-1
|
|
!
|
|
! where ra is the aerodynamic resistance (common to all gases) between a
|
|
! specific height and the surface, rb is the quasilaminar sublayer resistance
|
|
! (whose only dependence on the porperties of the gas of interest is its
|
|
! molecular diffusivity in air), and rc is the bulk surface resistance".
|
|
!
|
|
! In this subroutine both ra and rb are calculated elsewhere in CLM. Thus ra
|
|
! and rb were "globalized" in order to gain access to them for the calculation.
|
|
! "ram1" is the CLM variable used for ra. ram1 was globalized in the following
|
|
! subroutines; BareGroundFluxes.F90, Biogeophysics_lake.F90, CanopyFluxes.F90,
|
|
! and clmtype.F90.
|
|
!
|
|
! "rb" is the CLM variable used for rb in the Wesely equation above. rb was
|
|
! globalized in the following subroutines; clmtype.F90
|
|
!
|
|
! In Wesely (1989) rc is estimated for five seasonal categories and 11 landuse
|
|
! types. For each season and landuse type, Wesely compiled data into a
|
|
! look-up-table for several parameters used to calculate rc. In this subroutine
|
|
! the same values are used as found in wesely's look-up-tables, the only
|
|
! difference is that this subroutine uses a CLM generated LAI to select values
|
|
! from the look-up-table instead of seasonality. Inaddition, Wesely(1989)
|
|
! land use types are "mapped" into CLM plant function types (PFT).
|
|
!
|
|
! Subroutine written to operate at the patch level.
|
|
!
|
|
! Output:
|
|
!
|
|
! vd(n_species) !Dry deposition velocity [m s-1] for each molecule or species
|
|
!
|
|
! Author: Beth Holland and James Sulzman
|
|
!
|
|
! Modified: Francis Vitt -- 30 Mar 2007
|
|
! Modified: Maria Val Martin -- 15 Jan 2014
|
|
! Corrected major bugs in the leaf and stomatal resitances. The code is now
|
|
! coupled to LAI and Rs uses the Ball-Berry Scheme. Also, corrected minor
|
|
! bugs in rlu and rcl calculations. Added
|
|
! no vegetation removal for CO. See README for details and
|
|
! Val Martin et al., 2014 GRL for major corrections
|
|
!
|
|
!********* !!! IMPORTANT !!! ************
|
|
! STOMATAL RESISTANCE IS OPTIMIZED TO MATCH UP OBSERVATIONS
|
|
!-----------------------------------------------------------------------
|
|
|
|
use shr_kind_mod, only : r8 => shr_kind_r8
|
|
use clmtype
|
|
use abortutils, only : endrun
|
|
use clm_time_manager, only : get_nstep, get_curr_date, get_curr_time
|
|
use clm_atmlnd, only : clm_a2l
|
|
use spmdMod, only : masterproc
|
|
use seq_drydep_mod, only : n_drydep, drydep_list
|
|
use seq_drydep_mod, only : drydep_method, DD_XLND
|
|
use seq_drydep_mod, only : index_o3=>o3_ndx, index_o3a=>o3a_ndx, index_so2=>so2_ndx, index_h2=>h2_ndx, &
|
|
index_co=>co_ndx, index_ch4=>ch4_ndx, index_pan=>pan_ndx, &
|
|
index_xpan=>xpan_ndx
|
|
implicit none
|
|
save
|
|
|
|
private
|
|
|
|
public :: depvel_compute
|
|
|
|
CONTAINS
|
|
|
|
!-----------------------------------------------------------------------
|
|
! computes the dry deposition velocity of tracers
|
|
!-----------------------------------------------------------------------
|
|
subroutine depvel_compute( lbp , ubp )
|
|
use shr_const_mod , only : tmelt => shr_const_tkfrz
|
|
use seq_drydep_mod , only : seq_drydep_setHCoeff, mapping, drat, foxd, &
|
|
rcls, h2_a, h2_b, h2_c, ri, rac, rclo, rlu, &
|
|
rgss, rgso
|
|
use clm_varcon , only : istsoil, istice, istice_mec, istslak, istdlak, istwet, isturb
|
|
use clm_varctl , only : iulog
|
|
use pftvarcon , only : noveg, ndllf_evr_tmp_tree, ndllf_evr_brl_tree, &
|
|
ndllf_dcd_brl_tree, nbrdlf_evr_trp_tree, &
|
|
nbrdlf_evr_tmp_tree, nbrdlf_dcd_trp_tree, &
|
|
nbrdlf_dcd_tmp_tree, nbrdlf_dcd_brl_tree, &
|
|
nbrdlf_evr_shrub, nbrdlf_dcd_tmp_shrub, &
|
|
nbrdlf_dcd_brl_shrub, nc3_arctic_grass, &
|
|
nc3_nonarctic_grass, nc4_grass, nc3crop, &
|
|
nirrig, npcropmin, npcropmax
|
|
|
|
implicit none
|
|
|
|
!----Arguments-----------------------------------------------------
|
|
|
|
integer, intent(in) :: lbp, ubp ! pft bounds
|
|
|
|
! ------------------------ local variables ------------------------
|
|
! local pointers to implicit in arguments
|
|
integer , pointer :: plandunit(:) !pft's landunit index
|
|
integer , pointer :: ivt(:) !landunit type
|
|
integer , pointer :: itypveg(:) !vegetation type for current pft
|
|
integer , pointer :: pgridcell(:) !pft's gridcell index
|
|
real(r8), pointer :: pwtgcell(:) !weight of pft relative to corresponding gridcell
|
|
real(r8), pointer :: elai(:) !one-sided leaf area index with burying by snow
|
|
real(r8), pointer :: forc_t(:) !atmospheric temperature (Kelvin)
|
|
real(r8), pointer :: forc_q(:) !atmospheric specific humidity (kg/kg)
|
|
real(r8), pointer :: forc_psrf(:) !surface pressure (Pa)
|
|
real(r8), pointer :: latdeg(:) !latitude (degrees)
|
|
real(r8), pointer :: londeg(:) !longitude (degrees)
|
|
real(r8), pointer :: forc_rain(:) !rain rate [mm/s]
|
|
real(r8), pointer :: forc_snow(:) !snow rate [mm/s]
|
|
real(r8), pointer :: forc_lwrad(:) !direct beam radiation (visible only)
|
|
real(r8), pointer :: forc_solad(:,:) !direct beam radiation (visible only)
|
|
real(r8), pointer :: forc_solai(:,:) !direct beam radiation (visible only)
|
|
real(r8), pointer :: ram1(:) !aerodynamical resistance
|
|
real(r8), pointer :: vds(:) !aerodynamical resistance
|
|
real(r8), pointer :: rssun(:) !stomatal resistance
|
|
real(r8), pointer :: rssha(:) !shaded stomatal resistance (s/m)
|
|
real(r8), pointer :: fsun(:) !sunlit fraction of canopy
|
|
real(r8), pointer :: rb1(:) !leaf boundary layer resistance [s/m]
|
|
real(r8), pointer :: annlai(:,:) !12 months of monthly lai from input data set
|
|
real(r8), pointer :: mlaidiff(:) !difference in lai between month one and month two
|
|
real(r8), pointer :: velocity(:,:)
|
|
real(r8), pointer :: snowdp(:) ! snow height (m)
|
|
|
|
integer, pointer :: pcolumn(:) ! column index associated with each pft
|
|
integer :: c
|
|
integer , pointer :: itypelun(:) ! landunit type
|
|
|
|
real(r8), pointer :: h2osoi_vol(:,:) ! volumetric soil water (0<=h2osoi_vol<=watsat)
|
|
real(r8) :: soilw, var_soilw, fact_h2, dv_soil_h2
|
|
|
|
! new local variables
|
|
integer :: pi,g, l
|
|
integer :: ispec
|
|
integer :: length
|
|
integer :: wesveg !wesely vegegation index
|
|
integer :: clmveg !clm veg index from ivegtype
|
|
integer :: i
|
|
integer :: index_season !seasonal index based on LAI. This indexs wesely data tables
|
|
integer :: nstep !current step
|
|
integer :: indexp
|
|
|
|
real(r8) :: pg ! surface pressure
|
|
real(r8) :: tc ! temperature in celsius
|
|
real(r8) :: rs ! constant for calculating rsmx
|
|
real(r8) :: es ! saturation vapor pressur
|
|
real(r8) :: ws ! saturation mixing ratio
|
|
real(r8) :: rmx ! resistance by vegetation
|
|
real(r8) :: qs ! saturation specific humidity
|
|
real(r8) :: dewm ! multiplier for rs when dew occurs
|
|
real(r8) :: crs ! multiplier to calculate crs
|
|
real(r8) :: rdc ! part of lower canopy resistance
|
|
real(r8) :: rain ! rain fall
|
|
real(r8) :: spec_hum ! specific humidity
|
|
real(r8) :: solar_flux ! solar radiation(direct beam) W/m2
|
|
real(r8) :: lat ! latitude in degrees
|
|
real(r8) :: lon ! longitude in degrees
|
|
real(r8) :: sfc_temp ! surface temp
|
|
real(r8) :: minlai ! minimum of monthly lai
|
|
real(r8) :: maxlai ! maximum of monthly lai
|
|
real(r8) :: rds ! resistance for aerosols
|
|
|
|
!mvm 11/30/2013
|
|
real(r8) :: rlu_lai ! constant to calculate rlu over bulk canopy
|
|
real(r8) :: rs_factor ! constant to optimize stomatal resistance
|
|
|
|
logical :: has_dew
|
|
logical :: has_rain
|
|
real(r8), parameter :: rain_threshold = 1.e-7_r8 ! of the order of 1cm/day expressed in m/s
|
|
|
|
! local arrays: dependent on species only
|
|
!
|
|
|
|
real(r8), dimension(n_drydep) :: rsmx !vegetative resistance (plant mesophyll)
|
|
real(r8), dimension(n_drydep) :: rclx !lower canopy resistance
|
|
real(r8), dimension(n_drydep) :: rlux !vegetative resistance (upper canopy)
|
|
real(r8), dimension(n_drydep) :: rgsx !gournd resistance
|
|
real(r8), dimension(n_drydep) :: heff
|
|
real(r8) :: rc !combined surface resistance
|
|
real(r8) :: cts !correction to flu rcl and rgs for frost
|
|
real(r8) :: rlux_o3 !to calculate O3 leaf resistance in dew/rain conditions
|
|
|
|
! constants
|
|
real(r8), parameter :: slope = 0._r8 ! Used to calculate rdc in (lower canopy resistance)
|
|
integer, parameter :: wveg_unset = -1 ! Unset Wesley vegetation type
|
|
|
|
character(len=32), parameter :: subname = "depvel_compute"
|
|
|
|
!-------------------------------------------------------------------------------------
|
|
! jfl : mods for PAN
|
|
!-------------------------------------------------------------------------------------
|
|
real(r8) :: dv_pan
|
|
real(r8) :: c0_pan(11) = (/ 0.000_r8, 0.006_r8, 0.002_r8, 0.009_r8, 0.015_r8, &
|
|
0.006_r8, 0.000_r8, 0.000_r8, 0.000_r8, 0.002_r8, 0.002_r8 /)
|
|
real(r8) :: k_pan (11) = (/ 0.000_r8, 0.010_r8, 0.005_r8, 0.004_r8, 0.003_r8, &
|
|
0.005_r8, 0.000_r8, 0.000_r8, 0.000_r8, 0.075_r8, 0.002_r8 /)
|
|
!-----------------------------------------------------------------------
|
|
if ( n_drydep == 0 .or. drydep_method /= DD_XLND ) return
|
|
|
|
! local pointers to original implicit out arrays
|
|
|
|
! Assign local pointers to derived subtypes components (column-level)
|
|
forc_t => clm_a2l%forc_t
|
|
forc_q => clm_a2l%forc_q
|
|
forc_psrf => clm_a2l%forc_pbot
|
|
forc_rain => clm_a2l%forc_rain
|
|
|
|
latdeg => grc%latdeg
|
|
londeg => grc%londeg
|
|
ivt => pft%itype
|
|
elai => pps%elai
|
|
ram1 => pps%ram1
|
|
vds => pps%vds
|
|
fsun => pps%fsun
|
|
rssun => pps%rssun
|
|
rssha => pps%rssha
|
|
rb1 => pps%rb1
|
|
mlaidiff => pps%mlaidiff
|
|
annlai => pps%annlai
|
|
|
|
forc_solai => clm_a2l%forc_solai
|
|
forc_solad => clm_a2l%forc_solad
|
|
|
|
pwtgcell => pft%wtgcell
|
|
pgridcell => pft%gridcell
|
|
plandunit => pft%landunit
|
|
|
|
pcolumn => pft%column
|
|
itypelun => lun%itype
|
|
|
|
h2osoi_vol => cws%h2osoi_vol
|
|
|
|
velocity => pdd%drydepvel ! cm/sec
|
|
|
|
snowdp => cps%snowdp
|
|
|
|
! Assign local pointers to original implicit out arrays
|
|
!_________________________________________________________________
|
|
!
|
|
! Begin loop through pfts
|
|
pft_loop: do pi = lbp,ubp
|
|
l = plandunit(pi)
|
|
|
|
! Note: some glacier_mec pfts may have zero weight
|
|
gcell_wght: if (pwtgcell(pi)>0._r8 .or. itypelun(l)==istice_mec) then
|
|
|
|
c = pcolumn(pi)
|
|
g = pgridcell(pi)
|
|
pg = forc_psrf(g)
|
|
spec_hum = forc_q(g)
|
|
rain = forc_rain(g)
|
|
sfc_temp = forc_t(g)
|
|
lat = latdeg(g)
|
|
lon = londeg(g)
|
|
solar_flux = forc_solad(g,1)
|
|
clmveg = ivt(pi)
|
|
soilw = h2osoi_vol(c,1)
|
|
|
|
!map CLM veg type into Wesely veg type
|
|
wesveg = wveg_unset
|
|
if (clmveg == noveg ) wesveg = 8
|
|
if (clmveg == ndllf_evr_tmp_tree ) wesveg = 5
|
|
if (clmveg == ndllf_evr_brl_tree ) wesveg = 5
|
|
if (clmveg == ndllf_dcd_brl_tree ) wesveg = 5
|
|
if (clmveg == nbrdlf_evr_trp_tree ) wesveg = 4
|
|
if (clmveg == nbrdlf_evr_tmp_tree ) wesveg = 4
|
|
if (clmveg == nbrdlf_dcd_trp_tree ) wesveg = 4
|
|
if (clmveg == nbrdlf_dcd_tmp_tree ) wesveg = 4
|
|
if (clmveg == nbrdlf_dcd_brl_tree ) wesveg = 4
|
|
if (clmveg == nbrdlf_evr_shrub ) wesveg = 11
|
|
if (clmveg == nbrdlf_dcd_tmp_shrub ) wesveg = 11
|
|
if (clmveg == nbrdlf_dcd_brl_shrub ) wesveg = 11
|
|
if (clmveg == nc3_arctic_grass ) wesveg = 3
|
|
if (clmveg == nc3_nonarctic_grass ) wesveg = 3
|
|
if (clmveg == nc4_grass ) wesveg = 3
|
|
if (clmveg == nc3crop ) wesveg = 2
|
|
if (clmveg == nirrig ) wesveg = 2
|
|
if (clmveg >= npcropmin .and. clmveg <= npcropmax ) wesveg = 2
|
|
if (wesveg == wveg_unset )then
|
|
write(iulog,*) 'clmveg = ', clmveg, 'itypelun = ', itypelun(l)
|
|
call endrun( subname//': Not able to determine Wesley vegetation type')
|
|
end if
|
|
|
|
! create seasonality index used to index wesely data tables from LAI, Bascially
|
|
!if elai is between max lai from input data and half that max the index_season=1
|
|
|
|
|
|
!mail1j and mlai2j are the two monthly lai values pulled from a CLM input data set
|
|
!/fs/cgd/csm/inputdata/lnd/clm2/rawdata/mksrf_lai.nc. lai for dates in the middle
|
|
!of the month are interpolated using using these values and stored in the variable
|
|
!elai (done elsewhere). If the difference between mlai1j and mlai2j is greater
|
|
!than zero it is assumed to be fall and less than zero it is assumed to be spring.
|
|
|
|
!wesely seasonal "index_season"
|
|
! 1 - midsummer with lush vegetation
|
|
! 2 - Autumn with unharvested cropland
|
|
! 3 - Late autumn after frost, no snow
|
|
! 4 - Winter, snow on ground and subfreezing
|
|
! 5 - Transitional spring with partially green short annuals
|
|
|
|
|
|
!mlaidiff=jan-feb
|
|
minlai=minval(annlai(:,pi))
|
|
maxlai=maxval(annlai(:,pi))
|
|
|
|
index_season = -1
|
|
|
|
if ( itypelun(l) /= istsoil )then
|
|
if ( itypelun(l) == istice .or. itypelun(l) == istice_mec ) then
|
|
wesveg = 8
|
|
index_season = 4
|
|
elseif ( itypelun(l) == istdlak .or. itypelun(l) == istslak ) then
|
|
wesveg = 7
|
|
index_season = 4
|
|
elseif ( itypelun(l) == istwet ) then
|
|
wesveg = 9
|
|
index_season = 2
|
|
elseif ( itypelun(l) == isturb ) then
|
|
wesveg = 1
|
|
index_season = 2
|
|
end if
|
|
else if ( snowdp(c) > 0 ) then
|
|
index_season = 4
|
|
else if(elai(pi).gt.0.5_r8*maxlai) then
|
|
index_season = 1
|
|
endif
|
|
|
|
if (index_season<0) then
|
|
if (elai(pi).lt.(minlai+0.05*(maxlai-minlai))) then
|
|
index_season = 3
|
|
endif
|
|
endif
|
|
|
|
if (index_season<0) then
|
|
if (mlaidiff(pi).gt.0.0_r8) then
|
|
index_season = 2
|
|
elseif (mlaidiff(pi).lt.0.0_r8) then
|
|
index_season = 5
|
|
elseif (mlaidiff(pi).eq.0.0_r8) then
|
|
index_season = 3
|
|
endif
|
|
endif
|
|
|
|
if (index_season<0) then
|
|
call endrun( subname//': not able to determine season')
|
|
endif
|
|
|
|
! saturation specific humidity
|
|
!
|
|
es = 611_r8*exp(5414.77_r8*((1._r8/tmelt)-(1._r8/sfc_temp)))
|
|
ws = .622_r8*es/(pg-es)
|
|
qs = ws/(1._r8+ws)
|
|
|
|
has_dew = .false.
|
|
if( qs <= spec_hum ) then
|
|
has_dew = .true.
|
|
end if
|
|
if( sfc_temp < tmelt ) then
|
|
has_dew = .false.
|
|
end if
|
|
|
|
has_rain = rain > rain_threshold
|
|
|
|
if ( has_dew .or. has_rain ) then
|
|
dewm = 3._r8
|
|
else
|
|
dewm = 1._r8
|
|
end if
|
|
|
|
!Define tc
|
|
tc = sfc_temp - tmelt
|
|
|
|
!
|
|
! rdc (lower canopy res)
|
|
!
|
|
rdc=100._r8*(1._r8+1000._r8/(solar_flux+10._r8))/(1._r8+1000._r8*slope)
|
|
|
|
! surface resistance : depends on both land type and species
|
|
! land types are computed seperately, then resistance is computed as average of values
|
|
! following wesely rc=(1/(rs+rm) + 1/rlu +1/(rdc+rcl) + 1/(rac+rgs))**-1
|
|
|
|
!*******************************************************
|
|
call seq_drydep_setHCoeff( sfc_temp, heff(:n_drydep) )
|
|
!*********************************************************
|
|
|
|
species_loop1: do ispec=1, n_drydep
|
|
if(mapping(ispec).le.0) cycle
|
|
|
|
if(ispec.eq.index_o3.or.ispec.eq.index_o3a.or.ispec.eq.index_so2) then
|
|
rmx=0._r8
|
|
else
|
|
rmx=1._r8/((heff(ispec)/3000._r8)+(100._r8*foxd(ispec)))
|
|
endif
|
|
|
|
! correction for frost
|
|
cts = 1000._r8*exp( -tc - 4._r8 )
|
|
|
|
!ground resistance
|
|
rgsx(ispec) = 1._r8/((heff(ispec)/(1.e5_r8*(rgss(index_season,wesveg)+cts))) + &
|
|
(foxd(ispec)/(rgso(index_season,wesveg)+cts)))
|
|
|
|
!-------------------------------------------------------------------------------------
|
|
! special case for H2 and CO;; CH4 is set ot a fraction of dv(H2)
|
|
!-------------------------------------------------------------------------------------
|
|
if( ispec == index_h2 .or. ispec == index_co .or. ispec == index_ch4 ) then
|
|
|
|
if( ispec == index_co ) then
|
|
fact_h2 = 1.0_r8
|
|
elseif ( ispec == index_h2 ) then
|
|
fact_h2 = 0.5_r8
|
|
elseif ( ispec == index_ch4 ) then
|
|
fact_h2 = 50.0_r8
|
|
end if
|
|
|
|
!-------------------------------------------------------------------------------------
|
|
! no deposition on snow, ice, desert, and water
|
|
!-------------------------------------------------------------------------------------
|
|
if( wesveg == 1 .or. wesveg == 7 .or. wesveg == 8 .or. index_season == 4 ) then
|
|
rgsx(ispec) = 1.e36_r8
|
|
else
|
|
var_soilw = max( .1_r8,min( soilw,.3_r8 ) )
|
|
if( wesveg == 3 ) then
|
|
var_soilw = log( var_soilw )
|
|
end if
|
|
dv_soil_h2 = h2_c(wesveg) + var_soilw*(h2_b(wesveg) + var_soilw*h2_a(wesveg))
|
|
if( dv_soil_h2 > 0._r8 ) then
|
|
rgsx(ispec) = fact_h2/(dv_soil_h2*1.e-4_r8)
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
!-------------------------------------------------------------------------------------
|
|
! no deposition on water or no vegetation or snow (elai<=0)
|
|
!-------------------------------------------------------------------------------------
|
|
|
|
no_dep: if( wesveg == 7 .or. elai(pi).le.0_r8 ) then !mvm 11/26/2013
|
|
rclx(ispec)=1.e36_r8
|
|
rsmx(ispec)=1.e36_r8
|
|
rlux(ispec)=1.e36_r8
|
|
else
|
|
|
|
!Stomatal resistance
|
|
!MVM: adjusted rs to calculate stomata conductance over bulk canopy (CLM report pag 161)
|
|
rs=(fsun(pi)*rssun(pi)/elai(pi))+((rssha(pi)/elai(pi))*(1.-fsun(pi)))
|
|
|
|
!MVM: rs_factor=0.2 to match up Rs observations (Padro et al, 1996)
|
|
rs_factor = 0.2_r8
|
|
rsmx(ispec) = rs_factor*rs*drat(ispec)+rmx
|
|
|
|
! Leaf resistance
|
|
!MVM: adjusted rlu by LAI to get leaf resistance over bulk canopy (gao and wesely, 1995)
|
|
rlu_lai=cts+rlu(index_season,wesveg)/elai(pi)
|
|
rlux(ispec) = rlu_lai/(1.e-5_r8*heff(ispec)+foxd(ispec))
|
|
|
|
!Lower canopy resistance
|
|
rclx(ispec) = 1._r8/((heff(ispec)/(1.e5_r8*(rcls(index_season,wesveg)+cts))) + &
|
|
(foxd(ispec)/(rclo(index_season,wesveg)+cts)))
|
|
|
|
!-----------------------------------
|
|
!mvm 11/30/2013: special case for CO
|
|
!Dry deposition of CO and hydrocarbons is negligibly
|
|
!small in vegetation [Mueller and Brasseur, 1995].
|
|
!------------------------------------
|
|
if( ispec == index_co ) then
|
|
rclx(ispec)=1.e36_r8
|
|
rsmx(ispec)=1.e36_r8
|
|
rlux(ispec)=1.e36_r8
|
|
endif
|
|
|
|
!--------------------------------------------
|
|
! jfl : special case for PAN
|
|
!--------------------------------------------
|
|
if( ispec == index_pan ) then
|
|
dv_pan = c0_pan(wesveg) * (1._r8 - exp(-k_pan(wesveg)*(rs*drat(ispec))*1.e-2_r8 ))
|
|
|
|
if( dv_pan > 0._r8 .and. index_season /= 4 ) then
|
|
rsmx(ispec) = ( 1._r8/dv_pan )
|
|
end if
|
|
end if
|
|
|
|
endif no_dep
|
|
|
|
end do species_loop1
|
|
|
|
|
|
!----------------------------------------------
|
|
!Adjustment for dew and rain in leaf resitances
|
|
!---------------------------------------------
|
|
! no effect over water
|
|
no_water: if( wesveg.ne.7 ) then
|
|
!MVM: effect only on vegetated areas (elai> 0)
|
|
with_LAI: if (elai(pi).gt.0._r8) then
|
|
|
|
!
|
|
! no effect if sfc_temp < O C
|
|
!
|
|
non_freezing: if(sfc_temp.gt.tmelt) then
|
|
if( has_dew ) then
|
|
rlu_lai=cts+rlu(index_season,wesveg)/elai(pi)
|
|
rlux_o3 = 1._r8/((1._r8/3000._r8)+(1._r8/(3._r8*rlu_lai)))
|
|
|
|
if (index_o3 > 0) then
|
|
rlux(index_o3) = rlux_o3
|
|
endif
|
|
if (index_o3a > 0) then
|
|
rlux(index_o3a) = rlux_o3
|
|
endif
|
|
endif
|
|
|
|
if(has_rain) then
|
|
rlu_lai=cts+rlu(index_season,wesveg)/elai(pi)
|
|
rlux_o3 = 1._r8/((1._r8/1000._r8)+(1._r8/(3._r8*rlu_lai)))
|
|
|
|
if (index_o3 > 0) then
|
|
rlux(index_o3) = rlux_o3
|
|
endif
|
|
if (index_o3a > 0) then
|
|
rlux(index_o3a) = rlux_o3
|
|
endif
|
|
endif
|
|
|
|
species_loop2: do ispec=1,n_drydep
|
|
if(mapping(ispec).le.0) cycle
|
|
if(ispec.ne.index_o3.and.ispec.ne.index_o3a.and.ispec.ne.index_so2) then
|
|
|
|
if( has_dew .or. has_rain) then
|
|
rlu_lai=cts+rlu(index_season,wesveg)/elai(pi)
|
|
rlux(ispec)=1._r8/((1._r8/(3._r8*rlu_lai))+ &
|
|
(1.e-7_r8*heff(ispec))+(foxd(ispec)/rlux_o3))
|
|
endif
|
|
|
|
elseif(ispec.eq.index_so2) then
|
|
|
|
if( has_dew ) then
|
|
rlux(ispec) = 100._r8
|
|
endif
|
|
|
|
if(has_rain) then
|
|
rlu_lai=cts+rlu(index_season,wesveg)/elai(pi)
|
|
rlux(ispec) = 1._r8/((1._r8/5000._r8)+(1._r8/(3._r8*rlu_lai)))
|
|
endif
|
|
|
|
if( has_dew .or. has_rain ) then
|
|
!MVM:rlux=50 for SO2 in dew or rain only for *urban land* type surfaces.
|
|
if (wesveg.eq.1) then
|
|
rlux(ispec)=50._r8
|
|
endif
|
|
endif
|
|
end if
|
|
!mvm 11/30/2013: special case for CO
|
|
if( ispec.eq.index_co ) then
|
|
rlux(ispec)=1.e36_r8
|
|
endif
|
|
end do species_loop2
|
|
endif non_freezing
|
|
endif with_LAI
|
|
endif no_water
|
|
|
|
! resistance for aerosols
|
|
rds = 1._r8/vds(pi)
|
|
|
|
species_loop3: do ispec=1,n_drydep
|
|
if(mapping(ispec).le.0) cycle
|
|
|
|
!
|
|
! compute rc
|
|
!
|
|
rc = 1._r8/((1._r8/rsmx(ispec))+(1._r8/rlux(ispec)) + &
|
|
(1._r8/(rdc+rclx(ispec)))+(1._r8/(rac(index_season,wesveg)+rgsx(ispec))))
|
|
rc = max( 10._r8, rc)
|
|
!
|
|
! assume no surface resistance for SO2 over water
|
|
!
|
|
if ( drydep_list(ispec) == 'SO2' .and. wesveg == 7 ) then
|
|
rc = 0._r8
|
|
end if
|
|
|
|
select case( drydep_list(ispec) )
|
|
case ( 'SO4' )
|
|
velocity(pi,ispec) = (1._r8/(ram1(pi)+rds))*100._r8
|
|
case ( 'NH4','NH4NO3','XNH4NO3' )
|
|
velocity(pi,ispec) = (1._r8/(ram1(pi)+0.5_r8*rds))*100._r8
|
|
case ( 'Pb' )
|
|
velocity(pi,ispec) = 0.2_r8
|
|
case ( 'CB1', 'CB2', 'OC1', 'OC2', 'SOAM', 'SOAI', 'SOAT', 'SOAB', 'SOAX' )
|
|
velocity(pi,ispec) = 0.10_r8
|
|
case default
|
|
velocity(pi,ispec) = (1._r8/(ram1(pi)+rb1(pi)+rc))*100._r8
|
|
end select
|
|
end do species_loop3
|
|
endif gcell_wght
|
|
end do pft_loop
|
|
|
|
end subroutine depvel_compute
|
|
|
|
end module DryDepVelocity
|