c**** SLE001 E001M12 SOMTQ SLB211M9
c**** (same as frank''s soils64+2bare_soils+old runoff)
c**** change to evap calculation to prevent negative runoff
c**** soils45 but with snowmelt subroutine snmlt changed
c**** to melt snow before 1st layer ground ice.
ccc comments from new soils
c**** 8/11/97 - modified to include snow model
c**** 10/29/96 - aroot and broot back to original values; added
c**** call to cpars to change vegetation parameters for pilps.
c**** 9/7/96 - back to heat capacities/10
c**** 5/14/96 - added soils64 surface runoff changes/corrections
c**** 11/13/95 - changed aroot and broot for rainf for 1.5m root depth
c**** 10/11/95 - back to full heat capacity, to avoid cd oscillations.
c**** changes for pilps: (reversed)
c**** use soils100.com instead of soils45.com
c**** set im=36,jm=24
c**** set sdata,vdata and fdata to real*4
c**** divide canopy heat capacities by 10.d0
c**** change aroot of grass to 1.0d0, to reduce root depth.
c**** end of changes for pilps
c**** changes for pilps: (kept)
c**** modify gdtm surface flux timestep limits
c**** define new diagnostics
c**** zero out diagnostics at start of advnc
c**** end of changes for pilps
c**** modified for 2 types of bare soils
c****
c**** soils62 soils45 soils45 cdfxa 04/27/95
c**** same as soils45 but with snowmelt subroutine snmlt changed
c**** to melt snow before 1st layer ground ice.
ccc end comments from new soils
c**** also corrects evaps calculation.
c**** also includes masking effects in radation fluxes.
c**** modifies timestep for canopy fluxes.
c**** soils45 10/4/93
c**** uses bedrock as a soil texture. soil depth of 3.5m
c**** everywhere, where layers can have bedrock.
c**** requires sm693.data instead of sm691.data.
c**** sdata needs to be changed in calling program.
c**** soils44b 8/25/93
c**** uses snow conductivity of .088 w m-1 c-1 instead of .3d0
c**** soils44 8/16/93
c**** adds bedrock for heat calculations, to fill out the
c**** number of layers to ngm.
c**** soils43 6/25/93
c**** comments out call to fhlmt heat flux limits.
c**** uses ghinij to return wfc1, eliminates rewfc.
c**** soils42 6/15/93
c**** adds snow insulation
c**** soils41 5/24/93
c**** uses snow masking depth from vegetation height to determine
c**** fraction of snow that is exposed.
c**** reth must be called prior to retp.
c**** soils40 5/10/93
c**** removes snow from canopy and places it on vegetated soil.
c**** soils39 4/19/93
c**** modifications for real*8 or real*4 runs. common block
c**** ordering changed for efficient alignment. sdata,fdata,
c**** and vdata are explicitly real*4. on ibm rs/6000, should
c**** be compiled with -qdpc=e option for real*8 operation.
c**** to run real*4, change implicit statement in include file.
c**** soils38 2/9/93
c**** adds heat flux correction to handle varying coefficients
c**** of drag.
c**** soils37 1/25/93
c**** changes soil crusting parameter ku/d from .05 per hour to .1d0,
c**** to agree with morin et al.
c**** soils36 11/12/92
c**** calculates heat conductivity of soils using devries method.
c**** changes loam material heat capacity and conductivity
c**** to mineral values.
c**** soils35 10/27/92
c**** includes effect of soil crusting for infiltration by
c**** modifying hydraulic conductivity calculation of layer
c**** 1 in hydra.
c**** soils34 8/28/92
c**** uses effective leaf area index alaie for purposes of
c**** canopy conductance calculation.
c**** soils33 8/9/92
c**** changes canopy water storage capacity to .1mm per lai from 1.d0
c**** soils32 7/15/92
c**** 1) for purposes of infiltration only, reduces soil conductivity
c**** by (1-thetr*fice) instead of (1-fice).
c**** 2) betad is reduced by fraction of ice in each layer.
c**** 3) transpired water is removed by betad fraction in each layer,
c**** instead of by fraction of roots. prevents negative runoff.
c**** 4) speeds up hydra by using do loop instead of if check,
c**** by using interpolation point from bisection instead of logs,
c**** and by avoiding unnecessary calls to hydra. also elimates call
c**** to hydra in ma89ezm9.f.
c**** soils31 7/1/92
c**** 1) fixes fraction of roots when soil depth is less than root
c**** depth, thus fixing non-conservation of water.
c**** soils30 6/4/92
c**** 1) uses actual final snow depth in flux limit calculations,
c**** instead of upper and lower limits. fixes spurious drying
c**** of first layer.
c**** Added gpp, GPP terms 4/25/03 nyk
c**** Added decks parameter cond_scheme 5/1/03 nyk
!#define EVAP_VEG_GROUND
!#define RAD_VEG_GROUND
module sle001 7,2
use constant
, only : stbo,tfrz=>tf,sha,lhe,one,zero,rhow
& ,shw_kg=>shw,shi_kg=>shi,lhm
use ghycom, only : ngm, imt, nlsn
implicit none
save
private
c**** public functions:
public hl0, set_snow, advnc, evap_limits, xklh0
ccc physical constants and global model parameters
ccc converting constants from 1/kg to 1/m^3
!@var shw heat capacity of water (J/m^3 C)
real*8, parameter, public :: shw= shw_kg * rhow
!@var shi heat capacity of pure ice (J/m^3 C)
real*8, parameter, public :: shi= shi_kg * rhow
!@var fsn latent heat of melt (J/m^3)
real*8, parameter, public :: fsn= lhm * rhow
!@var elh latent heat of evaporation (J/m^3)
real*8, parameter :: elh= lhe * rhow
!@var prfr fraction (by area) of precipitation
real*8, parameter :: prfr = 0.1d0
ccc data needed for debugging
type, public :: ghy_debug_str
real*8 water(2), energy(2)
end type ghy_debug_str
type ( ghy_debug_str ), public :: ghy_debug
integer, public :: ijdebug
ccc public data
ccc main accumulators
real*8, public :: atrg,ashg,aevap,alhg,aruns,arunu,aeruns,aerunu
!@var tbcs surface temperature (C)
real*8, public :: tbcs
ccc diagnostics accumulatars
real*8, public :: aevapw,aevapd,aevapb,aepc,aepb,aepp,af0dt,af1dt
& , agpp
ccc some accumulators that are currently not computed:
& ,acna,acnc
ccc beta''s
& ,abetad,abetav,abetat,abetap,abetab,abeta
!@var zw water table (not computed ?)
real*8, public :: zw(2)
ccc private accumulators:
real*8 aedifs
ccc input fluxes
real*8, public :: pr,htpr,prs,htprs,srht,trht
ccc input bc''s
real*8, public :: ts,qs,pres,rho,vsm,ch,qm1
!@var dt earth time step (s)
real*8, public :: dt
ccc soil prognostic variables
!!! integer, parameter, public :: ngm=6, ng=ngm+1, imt=5
integer, parameter, public :: ng=ngm+1
real*8, public :: w(0:ngm,2),ht(0:ngm,2),dz(ngm)
ccc soil properties which need to be passed in:
real*8, public :: q(imt,ngm),qk(imt,ngm),sl,fv,fb
ccc topographic info which is passed in
real*8, public :: top_index, top_stdev
ccc soil internal variables wchich need to be passed from/to ghinij
real*8, public :: zb(ng),zc(ng),fr(ng),snowm !veg alaie, rs,
& ,thets(0:ng,2),thetm(0:ng,2),ws(0:ngm,2),shc(0:ng,2)
& ,thm(0:64,imt-1)
!veg & ,nm,nf,vh ! added by adf ! ,alai
!@var n the worst choice for global name: number of soil layers
!@var nth some magic number computed in hl0, used in ghinij and hydra
integer, public :: n,nth
ccc soil internal variables wchich need to be passed to driver anyway
real*8, public :: snowd(2),tp(0:ngm,2),fice(0:ngm,2)
!tp(:,:) is the temperature (C) of soil layers numbered from top
!to bottom for both bare and vegetated fractions.
!tp(:,1) - temperature for bare soil fraction
!tp(:,2) - temperature for vegetated fraction
!tp(0,2) - canopy
!tp(0,1) - not used
!tp(1,:) - upper soil layer
!tp(2,:) - second soil layer
ccc snow prognostic variables
!!! integer, parameter, public :: nlsn=3
integer, public :: nsn(2)
real*8, public :: dzsn(nlsn+1,2),wsn(nlsn,2),hsn(nlsn,2)
& ,fr_snow(2)
ccc tsn1 is private
!@var tsn1 temperature of the upper layer of snow (C)
real*8 tsn1(2)
real*8 betat,betad
real*8 gpp,dts,trans_sw
!xxx real*8, public :: cnc
ccc fractions of dry,wet,covered by snow canopy
!@var fd effective fraction of dry canopy (=0 if dew)
!@var fv effective fraction of wet canopy (=1 if dew)
!@var fd0 actual fraction of dry canopy
!@var fv0 actual fraction of wet canopy
!@var fm snow masking fraction for canopy (=1 if all snow)
real*8 fd,fw,fd0,fw0,fm
!@var theta fraction of water in the soil ( 1 = 100% )
!@var f water flux between the layers ( > 0 is up ) (m/s)
!@var fh heat flux between the layers ( > 0 is up ) (W)
!@var rnf surface runoff (m/s)
!@var rnff underground runoff (m/s)
!@var fc water flux at canopy boundaries ( > 0 is up ) (m/s)
!@var fch heat flux at canopy boundaries ( > 0 is up ) (W)
real*8 theta(0:ng,2),f(1:ng,2),fh(1:ng,2),rnf(2),rnff(ng,2)
& ,fc(0:1),fch(0:1)
ccc the following three declarations are various conductivity
ccc coefficients (or data needed to compute them)
real*8 xkh(ng,2),xkhm(ng,2),h(0:ng,2) ,xk(0:ng,2),xinfc(2)
& ,d(0:ng,2) ,xku(0:ng,2)
real*8 hlm(0:64), xklm(0:64,imt-1),dlm(0:64,imt-1)
real*8 xkus(ng,2), xkusa(2)
ccc the following are fluxes computed inside the snow model
ccc they are unit/(m^2 of snow), not multiplied by aby fractions
!@var flmlt flux of water down from snow to the ground (m/m^2)
!@var fhsng flux of heat down from snow to the ground (J/m^2)
!@var thrmsn flux of radiative heat up from snow (J/m^2)
real*8 flmlt(2),fhsng(2),thrmsn(2)
ccc similar values defined for large scale, i.e. flux to entire cell
ccc total flux is val*fr_snow + value_scale
real*8 flmlt_scale(2),fhsng_scale(2)
ccc put here some data needed for print-out during the crash
type :: debug_data_str
sequence
real*8 qc,qb
end type debug_data_str
type (debug_data_str) debug_data
ccc evaporation limiting data which one needs to pass the pbl
real*8, public :: evap_max_sat, evap_max_nsat, fr_sat
ccc potential evaporation
real*8 epb, epv
ccc evaporation fluxes: these are pure fluxes over m^2
ccc i.e. they are not multiplied by any fr_...
!@var evapb, evapbs, evapvw, evapvd, evapvs, evapvg evaporation flux (m/s)
!@+ (b=bare, v=vege, d=dry, w=wet, s=snow, g=ground)
real*8 :: evapb, evapbs, evapvw, evapvd, evapvs, evapvg
!@var devapbs_dt, devapvs_dt d(evap)/dT for snow covered soil (bare/veg)
real*8 :: devapbs_dt, devapvs_dt
c!@var evapor mean (weighted with fr_..) evaporation from bare/vege soil
ccc real*8 :: evapor(2)
!@var evapdl thanspiration from each soil layer
real*8 :: evapdl(ngm,2)
ccc sensible heat fluxes: these are pure fluxes over m^2
ccc i.e. they are not multiplied by any fr_...
!@var snsh sensible heat from soil to air (bare/vege) no snow (J/s)
!@var snshs sensible heat from snow to air (bare/vege) (J/s)
!@var dsnsh_dt d(sensible heat)/d(soil temp) : same for bare/vege/snow
real*8 :: snsh(2), snshs(2), dsnsh_dt
!!!! vars below this line were not included into GHYTPC yet !!!
!@var dripw liquid water flux from canopy (similar for bare soil) (m/s)
!@var htdripw heat of drip water (similar for bare soil) (J/s)
!@var drips snow flux from canopy (similar for bare soil) (m/s)
!@var htdrips heat of drip snow (similar for bare soil) (J/s)
real*8 dripw(2), htdripw(2), drips(2), htdrips(2)
!@var evap_tot total evap. (weighted with fr_snow,fm)(bare/veg)(m/s)
!@var snsh_tot total sens. heat(weighted with fr_snow,fm)(bare/veg)(m/s)
!@var thrm_tot total T^4 heat (weighted with fr_snow,fm)(bare/veg)(m/s)
real*8 evap_tot(2), snsh_tot(2), thrm_tot(2)
public evap_tot ! for debug only !
ccc data for tracers
!@var flux_snow water flux between snow layers (>0 down) (m/s)
!@var wsn_for_tr snow water before call to snow_adv
real*8 flux_snow(0:nlsn,2), wsn_for_tr(nlsn,2)
ccc the data below this line is not in GHYTPC yet !
ccc the following vars control if bare/vegetated fraction has to
ccc be computed (i.e. f[bv] is not zero)
integer :: i_bare, i_vege
logical :: process_bare, process_vege
C***
C*** Thread Private Common Block GHYTPC
C***
COMMON /GHYTPC/
& abeta,abetab,abetad,abetap,abetat,abetav,acna,acnc,agpp
& ,aedifs,aepb,aepc,aepp,aeruns,aerunu,aevap,aevapb
& ,aevapd,aevapw,af0dt,af1dt,alhg,aruns,arunu !veg alaie,
& ,ashg,atrg,betad,betat,ch,gpp,d,devapbs_dt,devapvs_dt
& ,drips,dripw,dsnsh_dt,dts,dz,dzsn,epb ! dt dlm
& ,epv,evap_max_nsat,evap_max_sat,evap_tot,evapb
& ,evapbs,evapdl,evapvd,evapvs,evapvw,evapvg,f !evapor,
& ,fb,fc,fch,fd,fd0,fh,fhsng,fhsng_scale,fice,flmlt,flmlt_scale
& ,fm,fr,fr_sat,fr_snow,fv,fw,fw0,h,hsn,ht !hlm
& ,htdrips,htdripw,htpr,htprs,pr,pres,prs,q,qk,qm1,qs
& ,rho,rnf,rnff,shc,sl,snowd,snowm,snsh,snsh_tot !veg rs,
& ,snshs,srht,tbcs,theta,thetm,thets,thrm_tot,thrmsn !thm
& ,top_index,top_stdev,tp,trht,ts,tsn1,vsm,w,ws,wsn,xinfc,xk
& ,xkh,xkhm,xku,xkus,xkusa,zb,zc,zw ! xklm
& ,ijdebug,n,nsn !nth
& ,flux_snow,wsn_for_tr,trans_sw
!----------------------------------------------------------------------!
& ,i_bare,i_vege,process_bare,process_vege
c not sure if it works with derived type. if not - comment the
c next line out (debug_data used only for debug output)
c & ,debug_data ! needs to go to compile on COMPAQ
!$OMP THREADPRIVATE (/GHYTPC/)
C***
C***
ccc external functions
real*8, external :: qsat,dqsatdt
contains
subroutine reth 3
!@sum computes theta(:,:), snowd(:), fm, fw, fd
c**** revises values of theta based upon w.
c**** input:
c**** w - water depth, m
c**** ws - saturated water depth, m
c**** dz - layer thickness, m
c**** thets - saturated theta
c**** snowd - snow depth, water equivalent m
c**** snowm - snow masking depth, water equivalent m
c**** output:
c**** theta - water saturation
c**** fw - fraction of wet canopy
c**** fd - fraction of dry canopy
c**** fm - fraction of snow that is exposed, or masking.
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
integer lsn,ibv,k
do ibv=i_bare,i_vege
do k=1,n
theta(k,ibv)=w(k,ibv)/dz(k)
end do
end do
c**** do canopy layer
c**** here theta is the fraction of canopy covered by water
if( process_vege .and. ws(0,2).gt.0.d0 )then
theta(0,2)=( else
theta(0,2)=0.d0
endif
theta(0,2)=min(theta(0,2),one)
c**** set up snowd variables
do ibv=i_bare,i_vege
snowd(ibv)=0.d0
do lsn=1,nsn(ibv)
ccc we compute snowd as if all snow was distributed uniformly
ccc over the cell (i.e. snowd = wsn * sn_frac / 1.d0 )
snowd(ibv) = snowd(ibv) + wsn(lsn,ibv) * fr_snow(ibv)
enddo
enddo
c**** fraction of wet canopy fw
fw=theta(0,2)
c**** determine fm from snowd depth and masking depth
fm=1.d0-exp(-snowd(2)/(snowm+1d-12))
!!! testing if setting bounds on fm will make tracers work ...
if ( fm < 1.d-3 ) fm=0.d0
! if ( fm > .99d0 ) fm=1.d0
c**** correct fraction of wet canopy by snow fraction
!fw=fw+fm*(1.d0-fw) !!! no idea what does this formula mean (IA)
fd=1.d0-fw
fw0=fw
fd0=fd
return
end subroutine reth
subroutine hydra 4
!@sum computes matr. potential h(k,ibv) and H2O condactivity xk(k,ibv)
c routine to return the equlibrium value of h in a mixed soil
c layer. the h is such that each soil texture has the same
c value of h, but differing values of theta.
c hydra also calculates the conductivity xk and diffussivity d.
c**** input:
c**** theta(k,ibv) - volumetric water concentration
c**** thetm(k,ibv) - minimum theta
c**** thets(k,ibv) - maximum theta
c**** nth - number of h0 intervals in table, a power of two.
c**** hlm(j) - table of h values, from 0 (at j=0) to hmin (at j=nth)
c**** thm(j,i) - value of relative theta at hlm(j) in texture i,
c**** ranging between thets(k,ibv) at j=0 to thetm(k,ibv) at j=nth.
c**** output:
c**** h - potential, m, including both matric and gravitational
c**** d - diffusivity, dl.
c**** xk - conductivity m/s.
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
c solve for h using bisection
c we assume that if j1.lt.j2 then hlm(j1).gt.hlm(j2)
c and thm(j1,i).gt.thm(j2,i).
c
c algdel=log(1.d0+alph0)
c
real*8 d1,d2,dl,hl,temp,thr,thr0,thr1,thr2,xk1,xkl,xklu,xku1,xku2
real*8 xkud
integer i,j,ibv,k,ith,j1,j2,jcm,jc
real*8 dz_total
xkud=2.78d-5
jcm=nint(log(float(nth))/log(2.d0))
do ibv=i_bare,i_vege
xk(n+1,ibv)=0.0d0
xku(0,ibv)=0.d0
do k=1,n
j1=0
j2=nth
thr1=thets(k,ibv)
thr2=thetm(k,ibv)
thr0=theta(k,ibv)
thr0=min(thr1,thr0)
thr0=max(thr2,thr0)
do jc=1,jcm
j=(j1+j2)/2
thr=0.d0
do i=1,imt-1
thr=thr+thm(j,i)*q(i,k)
end do
if(thr-thr0 .lt. 0.d0) then
c here thr is too small, bisect on low j end
j2=j
thr2=thr
else if (thr-thr0 .gt. 0.d0) then
c here thr is too large, bisect on high j end
j1=j
thr1=thr
else ! i.e. .eq.
c here thr is equal to thr0
hl=hlm(j)
j1=j
thr1=thr0
c the strange value for thr2 below is only for calculating temp
thr2=-10.d0
go to 500
end if
end do ! jc
c here theta is between two adjacent thr''s. interpolate.
hl=(hlm(j1)*(thr0-thr2)+hlm(j2)*(thr1-thr0))/(thr1-thr2)
500 continue
c**** only filling hl array with matric potential (gravitational to be
c**** added later)
h(k,ibv)=hl
c**** calculate diffusivity
ith=j1
temp=(thr1-thr0)/(thr1-thr2)
d1=0.d0
d2=0.d0
xku1=0.d0
xku2=0.d0
xkus(k,ibv) = 0.d0
do i=1,imt-1
d1=d1+q(i,k)*dlm(ith,i)
d2=d2+q(i,k)*dlm(ith+1,i)
xku1=xku1+q(i,k)*xklm(ith,i)
xku2=xku2+q(i,k)*xklm(ith+1,i)
xkus(k,ibv) = xkus(k,ibv) + q(i,k)*xklm(0,i)
end do
dl=(1.d0-temp)*d1+temp*d2
dl=(1.d0-fice(k,ibv))*dl
d(k,ibv)=dl
c**** calculate conductivity
xklu=(1.d0-temp)*xku1+temp*xku2
xklu=(1.d0-fice(k,ibv))*xklu
xku(k,ibv)=xklu
if(k.eq.1) then
xk1=0.d0
do i=1,imt-1
xk1=xk1+qk(i,1)*xklm(0,i)
end do
xkl=xk1
xkl=xkl/(1.d0+xkl/(-zc(1)*xkud))
xkl=(1.d0-fice(1,ibv)*theta(1,ibv)/thets(1,ibv))*xkl
xkl=max(zero,xkl)
xk(1,ibv)=xkl
else
xk(k,ibv)=sqrt(xku(k-1,ibv)*xku(k,ibv))
end if
end do ! k
end do ! ibv
ccc compute conductivity for topmodel (i.e. mean saturated conductivity)
do ibv=i_bare,i_vege
xkusa(ibv) = 0.d0
dz_total = 0.d0
do k=1,n
xkusa(ibv) = xkusa(ibv) + xkus(k,ibv)*dz(k)
dz_total = dz_total + dz(k)
enddo
xkusa(ibv) = xkusa(ibv) / dz_total
enddo
c add gravitational potential to hl
do k=1,n
do ibv=i_bare,i_vege
h(k,ibv)=h(k,ibv)+zc(k)
end do
end do
return
end subroutine hydra
subroutine hl0 1
c**** hl0 sets up a table of theta values as a function of matric
c**** potential, h. h is tabulated in a geometric series from
c**** 0 to hmin, with a first step of delh1. the theta values
c**** depend not only on the matric potential, but also on the
c**** soil texture. we solve a cubic equation to determine
c**** theta as a function of h. hl0 also outputs the conductivity
c**** and diffusivity tables.
c**** input:
c**** a - matric potential function parameters
c**** b - hydraulic conductivity function parameters
c**** p - hydraulic diffusivity function parameters
c**** sat - saturated thetas
c**** output:
c**** thm(j,i) - theta at j''th h point for texture i
c**** xklm(j,i) - conductivity at j''th h point for texture i
c**** dlm(j,i) - diffusivity at j''th h point for texture i
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
integer, parameter :: nexp=6
real*8, parameter :: c=2.3025851d0
real*8, dimension(4,imt-1), parameter :: a=reshape(
& (/ ! matric potential coefficients for
& .2514d0, 0.0136d0, -2.8319d0, 0.5958d0, ! sand
& .1481d0, 1.8726d0, 0.1025d0, -3.6416d0, ! loam
& .2484d0, 2.4842d0, 0.4583d0, -3.9470d0, ! clay
& .8781d0, -5.1816d0, 13.2385d0,-11.9501d0/), ! peat
& (/4,imt-1/))
real*8, dimension(4,imt-1), parameter :: b=reshape(
& (/ ! conductivity coefficients for
& -0.4910d0, -9.8945d0, 9.7976d0, -3.2211d0, ! sand
& -0.3238d0,-12.9013d0, 3.4247d0, 4.4929d0, ! loam
& -0.5187d0,-13.4246d0, 2.8899d0, 5.0642d0, ! clay
& -3.0848d0, 9.5497d0,-26.2868d0, 16.6930d0/), ! peat
& (/4,imt-1/))
real*8, dimension(4,imt-1), parameter :: p=reshape(
& (/ ! diffusivity coefficients for
& -0.1800d0, -7.9999d0, 5.5685d0, -1.8868d0, ! sand
& -0.1000d0,-10.0085d0, 3.6752d0, 1.2304d0, ! loam
& -0.1951d0, -9.7055d0, 2.7418d0, 2.0054d0, ! clay
& -2.1220d0, 5.9983d0,-16.9824d0, 8.7615d0/), ! peat
& (/4,imt-1/))
real*8, parameter :: sat(imt-1) = (/.394d0,.537d0,.577d0,.885d0/)
real*8 a1,a2,a3,alph0o,alpls1,arg(imt-1),delh1,delhn,dfunc,alph0
real*8 diff,func,hmin,hs,s,sxtn,testh,xtol
integer i,j,k,m,mmax
sxtn=16.d0
nth=2**nexp
hlm(0)=0.0d0
delh1=-0.00625d0
hmin=-1000.d0
delhn=delh1
c solve for alph0 in s=((1+alph0)**n-1)/alph0
s=hmin/delh1
alph0=1.d0/8.d0
10 alph0o=alph0
alph0=(s*alph0+1.d0)**(1.d0/nth)-1.d0
if(abs(alph0o-alph0).ge.1d-8) go to 10
alpls1=1.0d0+alph0
! algdel=log(1.d0+alph0) ! not used
do 100 j=1,nth
hlm(j)=hlm(j-1)+delhn
delhn=alpls1*delhn
100 continue
mmax=100
xtol=1d-6
do 200 i=1,imt-1
thm(0,i)=1.00d0
do 150 j=1,nth
hs=-exp(c*( a1=a(3,i)/a(4,i)
a2=( a3=a(1,i)/a(4,i)
testh=thm(j-1,i)
do 130 m=1,mmax
func=(testh**3)+(a1*(testh**2))+(a2*(testh))+a3
dfunc=(3*testh**2)+(2*a1*testh)+a2
diff=func/dfunc
testh=testh-diff
if(abs(diff).lt.xtol) go to 140
130 continue
print *,'max # iterations:',mmax
140 thm(j,i)=testh
150 continue
200 continue
do 280 j=0,nth
do 245 i=1,imt-1
xklm(j,i)=0.d0
arg(i)=0.d0
do 240 k=-1,2
arg(i)=arg(i)+b(k+2,i)*thm(j,i)**k
240 continue
arg(i)=min(arg(i),sxtn)
arg(i)=max(arg(i),-sxtn)
xklm(j,i)=exp(c*arg(i))
245 continue
!print *, 'xklm ', j, (xklm(j,i), i=1,imt-1)
do 265 i=1,imt-1
dlm(j,i)=0.d0
arg(i)=0.d0
do 260 k=-1,2
arg(i)=arg(i)+p(k+2,i)*thm(j,i)**k
260 continue
arg(i)=min(arg(i),sxtn)
arg(i)=max(arg(i),-sxtn)
dlm(j,i)=exp(c*arg(i))
265 continue
280 continue
do 350 j=0,nth
do 310 i=1,imt-1
thm(j,i)=thm(j,i)*sat(i)
310 continue
350 continue
return
end subroutine hl0
subroutine evap_limits( compute_evap, evap_max_out, fr_sat_out ) 2,13
!@sum computes maximal evaporation fluxes for current soil properties
!@calls cond
use vegetation, only : veg_conductance
!@var compute_evap if .true. compute evap, else just evap_max,fr_sat
!@var evap_max_out max evaporation from unsaturated soil
!@var fr_sat_out fraction of saturated soil
logical, intent(in) :: compute_evap
real*8, intent(out) :: evap_max_out, fr_sat_out
ccc local variables
!@var evap_max evap limits due to total amount of water in the soil
!@var evap_max_snow evap limits due to amount of snow
!@var evap_max_wet evap limit for saturated soil/canopy
!@var evap_max_dry evap limit for unsaturated soil/canopy
real*8 :: evap_max(2), evap_max_snow(2)
real*8 :: evap_max_wet(2), evap_max_dry(2)
!@var qb,qbs,qv,qvs specific humidity (b-bare, v-vege, s-snow)
!----------------------------------------------------------------------!
! adf
real*8 :: qb, qbs, qv, qvs, qvg
! real*8 :: qb, qbs, qvs
!----------------------------------------------------------------------!
!@var epb,epbs,epv,epvs potential evaporation (b-bare, v-vege, s-snow)
real*8 :: epbs, epvs, epvg ! , epb, epv
!@var rho3,cna,qm1dt local variable
real*8 rho3, cna, qm1dt
integer ibv, k
!@var hw the wilting point (m)
real*8, parameter :: hw = -100.d0
!@var betadl transpiration efficiency for each soil layer
real*8 betadl(ngm) ! used in evaps_limits only
real*8 pot_evap_can
real*8 cnc ! local cnc from veg_conductance, nyk
c cna is the conductance of the atmosphere
cna=ch*vsm
rho3=rho/rhow ! i.e divide by rho_water to get flux in m/s
ccc make sure that important vars are initialized (needed for ibv hack)
evap_max(:) = 0.d0
evap_max_snow(:) = 0.d0
evap_max_wet(:) = 0.d0
evap_max_dry(:) = 0.d0
betadl(:) = 0.d0
betad = 0.d0
abetad = 0.d0
acna = 0.d0
acnc = 0.d0
ccc !!! it''s a hack should call it somewhere else !!!
!call hydra
! soil moisture
do ibv=i_bare,i_vege
evap_max(ibv) = 0.
do k=1,n
evap_max(ibv) = evap_max(ibv) +
& ( enddo
enddo
! maximal evaporation from the snow fraction
! may be too restrictive with resp. to pr, but will leave for now
do ibv=i_bare,i_vege
evap_max_snow(ibv) = pr
do k=1,nsn(ibv)
evap_max_snow(ibv) = evap_max_snow(ibv) +
& wsn(k,ibv)/dt
enddo
enddo
! evaporation from bare soil
if ( process_bare ) then
ibv = 1
!!! no support for saturated soil yet, setting just in case...
evap_max_wet(ibv) = evap_max(ibv) + pr
! evap limited by diffusion and precipitation
evap_max_dry(ibv) = min( evap_max(ibv),
& 2.467d0*d(1,1)*(theta(1,1)-thetm(1,1))/dz(1) + pr )
endif
! evaporation from the canopy
if ( process_vege ) then
ibv = 2
evap_max_wet(ibv) = w(0,2)/dt !+ pr ! pr doesn''t work for snow
! soil under dry canopy
! dry canopy
!!! this needs "qs" from the previous time step
c betad is the the root beta for transpiration.
c hw is the wilting point.
c fr(k) is the fraction of roots in layer k
betad=0.d0
do 30 k=1,n
betadl(k)=(1.d0-fice(k,2))*fr(k)*max((hw-h(k,2))/hw,zero)
betad=betad+betadl(k)
30 continue
if ( betad < 1.d-12 ) betad = 0.d0 ! to avoid 0/0 divisions
abetad=betad ! return to old diagnostics
c Get canopy conductivity cnc and gpp
qv = qsat(tp(0,2)+tfrz,lhe,pres) ! for cond_scheme==2
call veg_conductance(
& cnc
& ,gpp
& ,trans_sw !nyk
& ,betad ! evaporation efficiency
& ,tp(0,2) ! canopy temperature C
& ,qv
& ,dts
& )
!!!! test
!!! trans_sw = .1d0
!!! print *,'trans_sw = ', trans_sw
!!! trans_sw = .1d0
betat=cnc/(cnc+cna+1d-12)
abetat=betat ! return to old diagnostics
acna=cna ! return to old diagnostics
acnc=cnc ! return to old diagnostics
! agpp=gpp !nyk 4/25/03. Put in subroutine accm.
pot_evap_can = betat*rho3*cna*(qsat(tp(0,2)+tfrz,lhe,pres) - qs)
evap_max_dry(ibv) = 0.d0
if ( betad > 0.d0 ) then
do k=1,n
evap_max_dry(ibv) = evap_max_dry(ibv)
& + min( pot_evap_can*betadl(k)/betad,
& ( enddo
endif
! evap_max_dry(ibv) = min( evap_max(ibv),
! & betat*rho3*cna*( qsat(tp(0,2)+tfrz,lhe,pres) - qs ) )
! may be pr should be included somehow in e_m_d
endif !process_vege
!! now we have to add the fluxes according to fractions
! bare soil
evap_max_sat = fb*fr_snow(1)*evap_max_snow(1)
evap_max_nsat = fb*(1.-fr_snow(1))*evap_max_dry(1)
! canopy
evap_max_sat = evap_max_sat +
& fv*( fr_snow(2)*fm*evap_max_snow(2) +
& (1.-fr_snow(2)*fm)*theta(0,2)*evap_max_wet(2) )
evap_max_nsat = evap_max_nsat +
& fv*( & * evap_max_dry(2) )
fr_sat = fb * fr_snow(1) +
& fv * ( fr_snow(2)*fm + (1.-fr_snow(2)*fm)*theta(0,2) )
ccc set variables for output
evap_max_out = evap_max_nsat
fr_sat_out = fr_sat
ccc not sure if all this should be in one subroutine, but otherwise
ccc one has to pass all these flux limits
if ( .not. ( evap_max_out > 0. .or. evap_max_out <= 0. ) )then
write(99,*)evap_max_out,cnc,evap_max(2)
call stop_model("GHY::evap_limits: evap_max_out = NaN",255)
endif
if ( .not. compute_evap ) return
cccccccccccccccccccccccccccccccccccccccc
ccc the rest of evap_limits does not take into account process_bare,
ccc process_vege , but it should compute ok for dummy values.
c**** qm1 has mass of water vapor in first atmosphere layer, kg m-2
qm1dt=.001d0*qm1/dt
! need this ? if(igcm.ge.0 .and. igcm.le.3) xl=eddy/(z1-zs)
c calculate bare soil, canopy and snow mixing ratios
qb = qsat(tp(1,1)+tfrz,lhe,pres)
qbs = qsat(tsn1(1)+tfrz,lhe,pres)
qv = qsat(tp(0,2)+tfrz,lhe,pres)
qvs = qsat(tsn1(2)+tfrz,lhe,pres)
c potential evaporation for bare and vegetated soil and snow
epb = rho3*cna*(qb-qs)
epbs = rho3*cna*(qbs-qs)
epv = rho3*cna*(qv-qs)
epvs = rho3*cna*(qvs-qs)
c bare soil evaporation
if ( process_bare ) then
evapb = min( epb, evap_max_dry(1) )
evapb = max( evapb, -qm1dt )
evapbs = min( epbs, evap_max_snow(1) )
evapbs = max( evapbs, -qm1dt )
! evapor(1) = fr_snow(1)*evapbs + (1.-fr_snow(1))*evapb
else
evapb = 0.d0; evapbs = 0.d0;
endif
c vegetated soil evaporation
if ( process_vege ) then
c evapvd is dry evaporation (transpiration) from canopy
c evapvw is wet evaporation from canopy (from interception)
evapvg = 0.d0 ! in case EVAP_VEG_GROUND not defined
evapvw = min( epv, evap_max_wet(2) )
evapvw = max( evapvw,-qm1dt )
evapvd = min(epv,evap_max_dry(2)) !evap_max_dry(2) depends on qs
evapvd = max( evapvd, 0.d0 )
evapvs = min( epvs, evap_max_snow(2) )
evapvs = max( evapvs, -qm1dt )
! evapor(2) = fr_snow(2)*fm*evapvs + (1.-fr_snow(2)*fm)*
! & ( theta(0,2)*evapvw + (1.-theta(0,2))*evapvd )
if ( evapvw < 0.d0 ) then ! we have dew
fw = 1.d0 ; fd = 0.d0 ! let dew fall on entire canopy
endif
else
evapvw = 0.d0; evapvd = 0.d0; evapvs = 0.d0
endif
devapbs_dt = rho3*cna*qsat(tsn1(1)+tfrz,lhe,pres)
& *dqsatdt(tsn1(1)+tfrz,lhe)
devapvs_dt = rho3*cna*qsat(tsn1(2)+tfrz,lhe,pres)
& *dqsatdt(tsn1(2)+tfrz,lhe)
c compute transpiration for separate layers (=0 for bare soil)
evapdl(1:n,1:2) = 0.d0
if ( betad > 0.d0 ) then
do k=1,n
evapdl(k,2) = evapvd*betadl(k)/betad
end do
endif
return
cccccccccccccccccccccccccccccccccccccccc
end subroutine evap_limits
subroutine sensible_heat 1
!@sum computes sensible heat for bare/vegetated soil and snow
implicit none
real*8 cna
cna=ch*vsm
snsh(1)=sha*rho*cna*(tp(1,1)-ts+tfrz) ! bare soil
snsh(2)=sha*rho*cna*(tp(0,2)-ts+tfrz) ! canopy
snshs(1) = sha*rho*cna*(tsn1(1)-ts+tfrz) ! bare soil snow
snshs(2) = sha*rho*cna*(tsn1(2)-ts+tfrz) ! canopy snow
dsnsh_dt = sha*rho*cna ! derivative is the same for all above
end subroutine sensible_heat
c subroutine qsbal
c**** finds qs that balances fluxes.
c**** obtains qs by successive approximation.
c**** calculates evaporation.
c**** input:
c**** ch - heat conductivity coefficient from ground to surface
c**** vsm - surface layer wind speed, m s-1
c**** rho - air density, kg m-3
c**** eddy - transfer coefficient from surface to first atmosphere
c**** theta - water saturation of layers and canopy
c**** tp - temperatures of layers and canopy, c
c**** pres - atmospheric pressure at ground
c**** z1 - height of first layer, m
c**** zs - height of surface layer, m
c**** dz - layer thicknesses, m
c**** snowd - snow depths, equivalent water m
c**** pr - precipitation, m s-1
c**** q1 - mixing ratio of first layer
c**** fr - fraction of roots in layer
c**** fb - fraction of bare soil
c**** fv - fraction of vegetated soil
c**** hw - wilting point, m
c**** output:
c**** qs - mixing ratio at surface layer
c**** evap - evaporation from bare and vegetated regions, m s-1
c**** evapw - evaporation from wet canopy, m s-1, including from snow
c**** evapd - evaporation from dry canopy, m s-1
c**** evaps - evaporation from snow from canopy, m s-1
c**** betad - dry canopy beta, based on roots
c****
subroutine drip_from_canopy 1,1
!@sum computes the flux of drip water (and its heat cont.) from canopy
!@+ the drip water is split into liquid water dripw() and snow drips()
!@+ for bare soil fraction the canopy is transparent
use snow_drvm, only : snow_cover_same_as_rad
real*8 ptmp,ptmps,pfac,pm,pmax
real*8 snowf,snowfs,dr,drs
integer ibv
real*8 melt_w, melt_h
c calculate snow fall. snowf is snow fall, m s-1 of water depth.
snowf=0.d0
if(htpr.lt.0.d0)snowf=min(-htpr/fsn,pr)
c snowfs is the large scale snow fall.
snowfs=0.d0
if(htprs.lt.0.d0)snowfs=min(-htprs/fsn,prs)
if ( process_vege ) then
ptmps=prs-snowfs
ptmps=ptmps-evapvw*fw
ptmp=pr-prs-(snowf-snowfs)
c use effects of subgrid scale precipitation to calculate drip
pm=1d-6
pmax=fd0*pm
drs=max(ptmps-pmax,zero)
dr=drs
if(ptmp.gt.0.d0)then
pfac=(pmax-ptmps)*prfr/ptmp
if(pfac.ge.0.d0)then
if(pfac.lt.30.d0) dr=ptmp*exp(-pfac)
else
dr=ptmp+ptmps-pmax
endif
endif
ccc make sure "dr" makes sense
dr = min( dr, pr-snowf-evapvw*fw )
dr = max( dr, pr-snowf-evapvw*fw - (ws(0,2)-w(0,2))/dts )
dr = max( dr, 0.d0 ) ! just in case (probably don''t need it)
dripw(2) = dr
htdripw(2) = shw*dr*max(tp(0,2),0.d0) !don''t allow it to freeze
! snow falls through the canopy
drips(2) = snowf
htdrips(2) = min ( htpr, 0.d0 ) ! liquid H20 is 0 C, so no heat
else ! no vegetated fraction
dripw(2) = 0.d0
htdripw(2) = 0.d0
drips(2) = 0.d0
htdrips(2) = 0.d0
endif
ccc for bare soil drip is precipitation
drips(1) = snowf
htdrips(1) = min( htpr, 0.d0 )
dripw(1) = pr - drips(1)
htdripw(1) = htpr - htdrips(1)
if ( snow_cover_same_as_rad .ne. 0 ) then
do ibv=1,2
if ( tp(1,ibv) > 0.d0 ) then
melt_w = (1.d0-fr_snow(ibv)) * drips(ibv)
melt_h = (1.d0-fr_snow(ibv)) * htdrips(ibv)
drips(ibv) = drips(ibv) - melt_w
htdrips(ibv) = htdrips(ibv) - melt_h
dripw(ibv) = dripw(ibv) + melt_w
htdripw(ibv) = htdripw(ibv) + melt_h
endif
enddo
endif
return
end subroutine drip_from_canopy
subroutine flg 1
!@sum computes the water fluxes at the surface
c bare soil
if ( process_bare ) then
f(1,1) = -flmlt(1)*fr_snow(1) - flmlt_scale(1)
& - (dripw(1)-evapb)*(1.d0-fr_snow(1))
endif
if ( process_vege ) then
c upward flux from wet canopy
fc(0) = -pr+evapvw*fw*(1.d0-fm*fr_snow(2))
! snow masking of pr is ignored since
! it is not included into drip
fc(1) = -dripw(2) - drips(2)
c vegetated soil
!!! flux down from canopy is not -f(1,2) any more !!
!!! it is = -dripw(2) - drips(2)
!!! f(1,2) is a flux up from tyhe soil
f(1,2) = -flmlt(2)*fr_snow(2) - flmlt_scale(2)
& - (dripw(2)-evapvg)*(1.d0-fr_snow(2))
endif
c compute evap_tot for accumulators
evap_tot(1) = evapb*(1.d0-fr_snow(1)) + evapbs*fr_snow(1)
evap_tot(2) = (evapvw*fw + evapvd*fd)*(1.d0-fr_snow(2)*fm)
& + evapvs*fr_snow(2)*fm + evapvg*(1.d0-fr_snow(2))
return
end subroutine flg
subroutine flhg 1
!@sum calculates the ground heat fluxes (to the surface)
c**** input:
c**** output:
c**** fh - heat fluxes from bare soil surface, and from canopy,
c**** and between canopy and vegetated soil.
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
c**** bare soil fluxes
real*8 thrm_can, thrm_soil(2)
thrm_can = stbo * (tp(0,2)+tfrz)**4
thrm_soil(1:2) = stbo * (tp(1,1:2)+tfrz)**4
!!!! test
!!! thrm_can = 0.d0
! bare soil
if ( process_bare ) then
fh(1,1) = - fhsng(1)*fr_snow(1) - fhsng_scale(1)
& + ( -htdripw(1) +
& evapb*elh + snsh(1) + thrm_soil(1) - srht - trht)
& *(1.d0-fr_snow(1))
endif
! vegetated soil
if ( process_vege ) then
fh(1,2) = - fhsng(2)*fr_snow(2) - fhsng_scale(2)
& + ( -htdripw(2)
& + evapvg*elh + thrm_soil(2) - thrm_can
& )
& *(1.d0-fr_snow(2))
! canopy
fch(0) = -htpr +
& (evapvw*elh*fw + snsh(2) + thrm_can
& -srht-trht
& + evapvd*elh*fd)
& *(1.d0-fm*fr_snow(2))
fch(1) = -(thrm_can - thrm_soil(2))
& *(1.d0-fr_snow(2)) !rad soil
& - (thrm_can - thrmsn(2))*fr_snow(2)*(1.d0-fm) !rad snow
& - htdripw(2) - htdrips(2) !heat of precip
endif
! compute thrm_tot for accumulators
thrm_tot(1) = thrm_soil(1)*(1.d0-fr_snow(1))
& + thrmsn(1)*fr_snow(1)
thrm_tot(2) = thrm_can*(1.d0-fr_snow(2)*fm)
& + thrmsn(2)*fr_snow(2)*fm
!XXXXXXXXXXXXXXXX don't forget to add trans_sw stuff to snow !
c compute total sensible heat
snsh_tot(1) = snsh(1)*(1.d0-fr_snow(1)) + snshs(1)*fr_snow(1)
snsh_tot(2) = snsh(2)*(1.d0-fr_snow(2)*fm)
& + snshs(2)*fr_snow(2)*fm
return
end subroutine flhg
subroutine runoff 1
c**** calculates surface and underground runoffs.
c**** input:
c**** xinfc - infiltration capacity, m s-1
c**** prfr - fraction of precipitation
c**** xk - conductivity, m s-1
c**** dz - layer thicknesses, m
c**** sl - slope
c**** sdstnc - interstream distance, m
c**** output:
c**** rnf - surface runoff
c**** rnff - underground runoff, m s-1
c use effects of subgrid scale rain
c use precipitation that includes smow melt
ccc surface runoff was rewritten in a more clear way 7/30/02
real*8 f_k0_exp_k !@var f_k0_exp_k coefficient for the topmodel
!@var water_down flux of water at the soil surface
!@var satfrac fraction of saturated soil
!@var prec_fr soil fraction at which precipitation is falling
real*8 water_down, satfrac, prec_fr
integer ibv,k
!@var sdstnc interstream distance (m)
real*8, parameter :: sdstnc = 100.d0
!@var rosmp used to compute saturated fraction: (w/ws)**rosmp
real*8, parameter :: rosmp = 8.
rnff(:,:) = 0.d0
rnf(:) = pr ! hack to conserve water (for ibv != 0,1)
! - should be set to 0 after testing
c**** surface runoff
do ibv=i_bare,i_vege
water_down = -f(1,ibv)
water_down = max( water_down, zero ) ! to make sure rnf > 0
! everything that falls on saturated fraction goes to runoff
satfrac = ( rnf(ibv) = satfrac * water_down
water_down = (1.d0 - satfrac) * water_down
! if we introduce large scale precipitation it should be
! applied here
!!! the following line is a hack. in a more precise approach
!!! one should treat snow-free fraction separately
prec_fr = max( prfr, fr_snow(ibv) )
if ( water_down*30.d0 > xinfc(ibv)*prec_fr ) then
rnf(ibv) = rnf(ibv) +
& water_down * exp( -xinfc(ibv)*prec_fr/water_down )
endif
enddo
c**** underground runoff
c sl is the slope, sdstnc is the interstream distance
do ibv=i_bare,i_vege
ccc this is some rough estimate for the expression f * k0 exp(zbar/f)
ccc in the topmodel expression for the runoff
! f_k0_exp = 0.d0
! do k=1,n
! f_k0_exp = f_k0_exp + xkus(k,ibv)*w(k,ibv)/ws(k,ibv)*dz(k)
! enddo
do k=1,n
rnff(k,ibv)=xku(k,ibv)*sl*dz(k)/sdstnc
!
end do
! print *,'runoff: ', sl/sdstnc, exp( -top_index )
end do
return
end subroutine runoff
subroutine fllmt 1,2
c**** places limits on the soil water fluxes
c**** input:
c**** w - water in layers, m
c**** ws - saturated water in layers, m
c**** dts - current time step size, s
c**** f - water fluxes, m s-1
c**** snk - water sink from layers, m s-1
c**** rnf - surface runoff, m s-1
c**** snowd - snow depth, equivalent water m
c**** snowf - snow fall, equivalent water m s-1
c**** output:
c**** f - limited water fluxes, m s-1
c**** snk - limited water sinks, m s-1
c**** rnf - limited surface runoff, m s-1
c**** temp variables:
c**** snowdu - the upper bound on the snow depth at end of time step
c**** snowdl - the lower bound on the snow depth at end of time step
c**** trunc - fix for truncation on ibm mainframes
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
real*8 dflux,drnf,trunc
integer k, ibv
real*8 wn
trunc=1d-6
trunc=1d-12
trunc=0.d0 ! works better since at some places thetm=thets=0
ccc prevent over/undersaturation of layers 2-n
do ibv=i_bare,i_vege
do k=n,2,-1
wn = w(k,ibv) + ( f(k+1,ibv) - f(k,ibv)
& - rnff(k,ibv) - fd*(1.-fr_snow(2)*fm)*evapdl(k,ibv) )*dts
ccc compensate oversaturation by increasing flux up
if (wn - ws(k,ibv) > trunc)
& f(k,ibv) = f(k,ibv) + (wn - ws(k,ibv) + trunc)/dts
ccc compensate undersaturation by decreasing runoff
if (wn - dz(k)*thetm(k,ibv) < trunc) then
rnff(k,ibv) = rnff(k,ibv) +
& (wn - dz(k)*thetm(k,ibv) - trunc)/dts
if ( rnff(k,ibv) < 0.d0 ) then ! have to compensate with f
f(k,ibv) = f(k,ibv) + rnff(k,ibv)
rnff(k,ibv) = 0.d0
endif
endif
end do
end do
ccc prevent over/undersaturation of first layer
do ibv=i_bare,i_vege
wn = w(1,ibv) + ( f(2,ibv) - f(1,ibv)
& - rnf(ibv) - rnff(1,ibv)
& - fd*(1.-fr_snow(2)*fm)*evapdl(1,ibv) )*dts
if ( wn - ws(1,ibv) > trunc)
& rnf(ibv) = rnf(ibv) + (wn - ws(1,ibv) + trunc)/dts
if ( wn - dz(1)*thetm(1,ibv) < trunc )
& rnf(ibv) = rnf(ibv) +
& (wn - dz(1)*thetm(1,ibv) - trunc)/dts
enddo
ccc now trying to remove negative runoff
do ibv=i_bare,i_vege
k = 1
do while ( rnf(ibv) .lt. 0.d0 .and. k .le. n )
if ( k > 1 ) then
ccc this is how much water we can take from layer k
dflux = f(k+1,ibv) + ( & - f(k,ibv) - rnff(k,ibv)
& - fd*(1.-fr_snow(2)*fm)*evapdl(k,ibv)
f(k,ibv) = f(k,ibv) - rnf(ibv)
rnf(ibv) = rnf(ibv) + min( -rnf(ibv), dflux )
endif
ccc rnff always >= 0, use it also to compensate rnf<0
if ( rnff(k,ibv) < 0.d0 )
& call stop_model('fllmt: negative underground runoff',255)
drnf = min( -rnf(ibv), rnff(k,ibv) )
rnf(ibv) = rnf(ibv) + drnf
rnff(k,ibv) = rnff(k,ibv) - drnf
k = k + 1
enddo
ccc check if rnf==0 up to machine accuracy
if ( rnf(ibv) .lt. -1d-12 ) then
print *, 'fllmt: rnf<0, ibv=',ibv,rnf(ibv)
call stop_model('fllmt: negative runoff',255)
endif
ccc if -1d-12 < rnf < 0. put it to 0 to avoid possible problems
ccc actually for ground hydrology it is not necessary
rnf(ibv) = max ( rnf(ibv), 0.d0 )
enddo
return
end subroutine fllmt
subroutine xklh 1,1
c**** evaluates the heat conductivity between layers
c**** uses the method of DeVries.
c**** input:
c**** zb - soil layer boundaries, m
c**** zc - soil layer centers, m
c**** theta - soil water saturation
c**** fice - fraction of ice in layers
c**** alami - ice heat conductivity
c**** alamw - water heat conductivity
c**** tp - temperature of layers, c
c**** shw - specific heat of water
c**** shi - specific heat of ice
c**** shc - heat capacity of soil layers
c**** dz - layer thicknesses
c**** k,ibv - soil layer
c**** dts - the current time step
c**** output:
c**** xkh(k,ibv) - heat conductivities in each of the soil layers
c**** xkhm(k,ibv) - average heat conductivity between layer k and k-1
c****
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
c dimension xsha(ng,2),xsh(ng,2),gabc(3),hcwt(imt-1)
c
c calculate with changing ga for air. ga is the depolarization
c factor for air, calculated by linear interpolation from .333d0
c at saturation to .035 at 0 water, following devries.
real*8 gaa,gabc(3),gca,xa,xb,xden,xi,xnum,xs,xw
ccc real*8, save :: ba,hcwt(imt-1),hcwta,hcwtb,hcwti,hcwtw
ccc real*8, save :: xsha(ng,2),xsh(ng,2)
real*8 :: ba,hcwt(imt-1),hcwta,hcwtb,hcwti,hcwtw
real*8 :: xsha(ng,2),xsh(ng,2)
COMMON /XKLHSAV/ BA, HCWTW, HCWTI, HCWTB, XSHA, XSH
!$OMP THREADPRIVATE (/XKLHSAV/)
integer i, j, ibv, k
c the alam''s are the heat conductivities
real*8, parameter :: alamw = .573345d0
& ,alami = 2.1762d0
& ,alama = .025d0
& ,alambr= 2.9d0
& ,alams(imt-1) = (/ 8.8d0, 2.9d0, 2.9d0, .25d0 /)
do ibv=i_bare,i_vege
do k=1,n
gaa=.298d0*theta(k,ibv)/(thets(k,ibv)+1d-6)+.035d0
gca=1.d0-2.d0*gaa
hcwta=(2.d0/(1.d0+ba*gaa)+1.d0/(1.d0+ba*gca))/3.d0
c xw,xi,xa are the volume fractions. don''t count snow in soil lyr 1
xw=w(k,ibv)*(1.d0-fice(k,ibv))/dz(k)
xi=w(k,ibv)*fice(k,ibv)/dz(k)
xa=(thets(k,ibv)-theta(k,ibv))
xb=q(imt,k)
xnum=xw*hcwtw*alamw+xi*hcwti*alami+xa*hcwta*alama+xsha(k,ibv)
& + xb*hcwtb*alambr
xden=xw*hcwtw+xi*hcwti+xa*hcwta+xsh(k,ibv)+xb*hcwtb
xkh(k,ibv)=xnum/xden
if ( xkh(k,ibv) .lt. 0.d0 )
& call stop_model('xklh: heat conductivity<0',255)
end do
end do
c get the average conductivity between layers
do ibv=i_bare,i_vege
do k=2,n
xkhm(k,ibv)=((zb(k)-zc(k-1))*xkh(k,ibv)
& + (zc(k)-zb(k))*xkh(k-1,ibv)
& )/(zc(k)-zc(k-1))
end do
end do
c****
return
entry xklh0 1
c gabc''s are the depolarization factors, or relative spheroidal axes.
gabc(1)=.125d0
gabc(2)=gabc(1)
gabc(3)=1.d0-gabc(1)-gabc(2)
c hcwt''s are the heat conductivity weighting factors
hcwtw=1.d0
hcwti=0.d0
hcwtb=1.d0
do i=1,imt-1
hcwt(i)=0.d0
end do
do j=1,3
hcwti=hcwti+1.d0/(1.d0+(alami/alamw-1.d0)*gabc(j))
do i=1,imt-1
hcwt(i)=hcwt(i)+1.d0/(1.d0+(alams(i)/alamw-1.d0)*gabc(j))
end do
end do
hcwti=hcwti/3.d0
do i=1,imt-1
hcwt(i)=hcwt(i)/3.d0
end do
do ibv=1,2 ! i_bare,i_vege
do k=1,n
xsha(k,ibv)=0.d0
xsh(k,ibv)=0.d0
do i=1,imt-1
xs=(1.d0-thm(0,i))*q(i,k)
xsha(k,ibv)=xsha(k,ibv)+xs*hcwt(i)*alams(i)
xsh(k,ibv)=xsh(k,ibv)+xs*hcwt(i)
end do
end do
end do
ba=alama/alamw-1.d0
return
end subroutine xklh
subroutine fl 1
!@sum computes water flux between the layers
c**** input:
c**** h - soil potential of layers, m
c**** xk - conductivity of layers, m s-1
c**** zc - layer centers, m
c**** output:
c**** f - fluxes between layers, m s-1
c**** xinfc - infiltration capacity, m s-1
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
c****
integer ibv,k
do ibv=i_bare,i_vege
f(n+1,ibv)=0.d0
end do
c****
do ibv=i_bare,i_vege
do k=2,n
f(k,ibv)=-xk(k,ibv)*( end do
end do
c**** put infiltration maximum into xinfc
do ibv=i_bare,i_vege
xinfc(ibv)=xk(1,ibv)*h(1,ibv)/zc(1)
end do
return
end subroutine fl
subroutine flh 1
c**** evaluates the heat flux between layers
c**** subroutine fl must be called first
c**** input:
c**** zb - soil layer boundaries, m
c**** zc - soil layer centers, m
c**** theta - soil water saturation
c**** fice - fraction of ice in layers
c**** alami - ice heat conductivity
c**** alamw - water heat conductivity
c**** tp - temperature of layers, c
c**** shw - specific heat of water
c**** output:
c**** fh - heat flux between layers
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
c****
integer ibv, k
do ibv=i_bare,i_vege
fh(n+1,ibv)=0.d0
c total heat flux is heat carried by water flow plus heat conduction
do k=2,n
fh(k,ibv)=-xkhm(k,ibv)*(tp(k-1,ibv)-tp(k,ibv))/(zc(k-1)-zc(k))
if( fh(k,ibv)=fh(k,ibv)+f(k,ibv)*tp(k,ibv)*shw
else
fh(k,ibv)=fh(k,ibv)+f(k,ibv)*tp(k-1,ibv)*shw
endif
end do
end do
return
end subroutine flh
subroutine retp 2,4
c**** evaluates the temperatures in the soil layers based on the
c**** heat values. also executes snow melt.
c**** input:
c**** w - water in soil layers, m
c**** ht - heat in soil layers
c**** fsn - heat of fusion of water
c**** shc - specific heat capacity of soil
c**** shi - specific heat capacity of ice
c**** shw - specific heat capcity of water
c**** snowd - snow depth, equivalent water m
c**** fb - fraction of bare soil
c**** fv - fraction of vegetation
c**** fm - snow vegetation masking fraction (requires reth called first)
c**** output:
c**** tp - temperature of layers, c
c**** fice - fraction of ice of layers
ccc include 'soils45.com'http://www.giss.nasa.gov/internal/
c**** soils28 common block 9/25/90
integer ibv, k, kk
tp(:,:) = 0.d0
fice(:,:) = 0.d0
do ibv=i_bare,i_vege
kk=2-ibv
do k=kk,n
! tp(k,ibv)=0.d0
!if( ! fice(k,ibv)=-ht(k,ibv)/(fsn*w(k,ibv))
!else
! fice(k,ibv)=0.d0
!endif
if( fsn*w(k,ibv)+ht(k,ibv) .lt. 0.d0 ) then ! all frozen
tp(k,ibv)=(ht(k,ibv)+w(k,ibv)*fsn)/(shc(k,ibv)+w(k,ibv)*shi)
fice(k,ibv)=1.d0
else if( ht(k,ibv) .gt. 0.d0 ) then ! all melted
tp(k,ibv)=ht(k,ibv)/(shc(k,ibv)+w(k,ibv)*shw)
!shc -- specific heat of canopy
!shw -- specific heat of water
! fice(k,ibv)=0.d0
else if( w(k,ibv) .ge. 1d-12 )then ! part frozen
fice(k,ibv)=-ht(k,ibv)/(fsn*w(k,ibv))
endif
end do
end do
ccc this is a fix for undefined tsn1 at the beginning of soil routines
ccc probably should be moved to some other place
tsn1(:) = 0.d0
do ibv=i_bare,i_vege
! tsn1(ibv) = 0.d0
if ( wsn(1,ibv) .gt. 1.d-6 .and.
& hsn(1,ibv) + wsn(1,ibv)*fsn .lt. 0.d0 ) then
tsn1(ibv) = (hsn(1,ibv) + wsn(1,ibv)*fsn)/(wsn(1,ibv)*shi)
endif
ccc the following is a hack. it is necessary only at the beginning of th
ccc run, when some temperatures are not initialized properly.
ccc should be removed when program is rewritten in a more clean way...
!!! I think it is not needed any more ...
! if ( wsn(1,ibv) .le. 1.d-6 ) then
! tsn1(ibv) = tp(2-ibv,ibv)
! endif
enddo
if(tp(1,1).gt.100.d0.or.tp(0,2).gt.100.d0)then
write(6,*)'retp tp bounds error'
write(6,*)'ijdebug',ijdebug
call reth
call hydra
call outw(1)
call stop_model(
& 'retp: tground > 100C - see soil_outw and fort.99',255)
endif
return
end subroutine retp
subroutine apply_fluxes 1,3
!@sum apply computed fluxes to adwance w and h
!@ver 1.0
c**** shw - specific heat of water
c**** tp - temperature of layers, c
integer ibv, k
ccc the canopy
w(0,2) = w(0,2) + ( fc(1) - fc(0) )*dts
ht(0,2)=ht(0,2) + ( fch(1) - fch(0) )*dts
ccc the soil
do ibv=i_bare,i_vege
ccc surface runoff
w(1,ibv) = w(1,ibv) - rnf(ibv)*dts
ht(1,ibv) = ht(1,ibv) - shw*max(tp(1,ibv),0.d0)*rnf(ibv)*dts
ccc rest of the fluxes
do k=1,n
w(k,ibv) = w(k,ibv) +
& ( f(k+1,ibv) - f(k,ibv) - rnff(k,ibv)
& - fd*(1.-fr_snow(2)*fm)*evapdl(k,ibv)
& )*dts
ht(k,ibv) = ht(k,ibv) +
& ( fh(k+1,ibv) - fh(k,ibv) -
& shw*max(tp(k,ibv),0.d0)*( rnff(k,ibv)
! & + fd*(1.-fr_snow(2)*fm)*evapdl(k,ibv) !!! hack !!!
& ) )*dts
end do
end do
ccc do we need this check ?
if ( w(0,2) < 0.d0 ) then
if (w(0,2) < -1.d-12) call stop_model('GHY:Canopy H2O < 0',255)
w(0,2) = 0.d0
endif
ccc check for under/over-saturation
do ibv=i_bare,i_vege
do k=1,n
if ( w(k,ibv) < dz(k)*thetm(k,ibv) - 1.d-14 ) then
print*,"ghy:",k,ibv,w(k,ibv),dz(k),thetm(k,ibv)
call stop_model("ghy: w < dz*thetm",255)
end if
if ( w(k,ibv) > ws(k,ibv) + 1.d-14 )
& call stop_model("ghy: w > ws",255)
w(k,ibv) = max( w(k,ibv), dz(k)*thetm(k,ibv) )
w(k,ibv) = min( w(k,ibv), ws(k,ibv) )
enddo
enddo
return
end subroutine apply_fluxes
subroutine advnc 1,32
c**** advances quantities by one time step.
c**** input:
c**** dt - time step, s
c**** dz - layer thickness, m
c**** tp - layer temperatures, c
c**** tfrz - freezing point of water, k
c**** w - soil water in layers, m
c**** snowd - snow depth, m
c**** f - water flux, m s-1
c**** snk - water sinks, m s-1
c**** ht - heat in soil layers
c**** fh - heat flux in soil layers
c**** snkh - heat sink in layers
c**** snowf - snow fall, m s-1 of equivalent water
c**** evap - evaporation, m s-1
c**** output:
c**** w - updater water in soil layers, m s-1
c**** ht - updated heat in soil layers
c**** snowd - updated snow depth, m s-1 of equivalent water
c**** rus - overall surface runoff, m s-1 replaced by aruns
c**** aruns - overall surface runoff, kg m-2
c**** aeruns - overall surface heat runoff, j m-2
c**** aerunu - underground heat runoff, j m-2
c**** uses:
c**** retp,reth,fl,flg,runoff,sink,sinkh,fllmt,flh,flhg.
c**** also uses surf with its required variables.
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
use vegetation, only: update_veg_locals
real*8 dtm,tb0,tc0,dtr,tot_w1
integer limit,nit
real*8 dum1, dum2
limit=300 ! 200 increase to avoid a few more stops
nit=0
dtr=dt
dts=dt ! make sure that dts is always initialized
ccc trying to skip fb==0 and fv==0 fractions of the cell
ccc reset main water/heat fluxes, so they are always initialized
f(:,:) = 0.d0
fh(:,:) = 0.d0
fc(:) = 0.d0
fch(:) = 0.d0
ccc normal case (both present)
i_bare = 1; i_vege = 2
process_bare = .true.; process_vege = .true.
if ( fb == 0.d0 ) then ! bare fraction is missing
i_bare = 2
process_bare = .false.
endif
if ( fv == 0.d0 ) then ! bare fraction is missing
i_vege = 1
process_vege = .false.
endif
!debug debug!
! pr = 0.d0
! htpr = 0.d0
! trpr = 0.d0
! srht = 0.d0
! trht = 0.d0
!!!
call reth
call retp
tb0=tp(1,1)
tc0=tp(0,2)
cddd print '(a,10(e12.4))', 'ghy_temp_b ',
cddd & tp(1,1),tp(2,1),tp(0,2),tp(1,2),tp(2,2)
ccc accm0 was not called here in older version - check
call accm0
do while ( dtr > 0.d0 )
nit=nit+1
if(nit.gt.limit)go to 900
call hydra
!call qsbal
!debug
! fm = 1.d0
!!!
call evap_limits( .true., dum1, dum2 )
call sensible_heat
!debug debug!
! evapb = 0.d0
! evapbs = evapbs * .1
! evapvw = 0.d0
! evapvd = 0.d0
! evapdl = 0.d0
! evapvs = 0.d0
! devapbs_dt = 0.d0
! devapvs_dt =0.d0
! dsnsh_dt = 0.d0
! if ( evapvs < 0.d0 ) then
! evapvs = 0.d0; devapvs_dt =0.d0
! endif
!!!
call xklh
call gdtm(dtm)
!print *,'dtm ', ijdebug, dtm
if ( dtm >= dtr ) then
dts = dtr
dtr = 0.d0
else
dts = min( dtm, dtr*0.5d0 )
dtr=dtr-dts
endif
!!!
! drips = 0
! dripw = 0
call drip_from_canopy
call check_water(0)
call check_energy(0)
call snow
call fl
call flg
! call check_f11
call runoff
!debug
! rnff(:,:) = 0.d0
! rnf(:) = 0.d0
!!!
!call sink
call fllmt
!debug
! rnff(:,:) = 0.d0
! rnf(:) = 0.d0
!!!
! call check_f11
!call sinkh
call flh
call flhg
c call fhlmt
! call check_f11
!call check_water(0)
call apply_fluxes
cddd print '(a,10(e12.4))', 'tr_w_1 '
cddd & , tr_w(1,:,1) - w(:,1) * 1000.d0
cddd print '(a,10(e12.4))', 'tr_w_2 '
cddd & , tr_w(1,:,2) - w(:,2) * 1000.d0
!!!
call check_water(1)
call check_energy(1)
call accm
call reth
call retp
cddd print '(a,i6,10(e12.4))', 'ghy_temp ', ijdebug,
cddd & tp(1,1),tp(2,1),tp(0,2),tp(1,2),tp(2,2)
call update_veg_locals(evap_tot(2), rho, rhow, ch, vsm,qs)
enddo
call accmf
call hydra
call wtab ! for gcm diag. only
return
900 continue
write(99,*)'limit exceeded'
write(99,*)'dtr,dtm,dts',dtr,dtm,dts
write(99,*)'tb0,tc0',tb0,tc0
call hydra
call outw(2)
call stop_model(
& 'advnc: time step too short - see soil_outw and fort.99',255)
end subroutine advnc
subroutine accm 1,2
c**** accumulates gcm diagnostics
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
c**** the following lines were originally called before retp,
c**** reth, and hydra.
real*8 qsats
real*8 cpfac,dedifs,dqdt,el0,epen,h0
integer k
ccc main fluxes which should conserve water/energy
atrg = atrg + ( thrm_tot(1)*fb + thrm_tot(2)*fv )*dts
ashg = ashg + ( snsh_tot(1)*fb + snsh_tot(2)*fv )*dts
aevap = aevap + ( evap_tot(1)*fb + evap_tot(2)*fv )*dts
alhg = elh*aevap
aruns=aruns+(fb*rnf(1)+fv*rnf(2))*dts
aeruns=aeruns+shw*( fb*max(tp(1,1),0.d0)*rnf(1)
& + fv*max(tp(1,2),0.d0)*rnf(2) )*dts
do k=1,n
arunu=arunu+(rnff(k,1)*fb+rnff(k,2)*fv)*dts
aerunu=aerunu+ shw*( max(tp(k,1),1.d0)*rnff(k,1)*fb
* + max(tp(k,2),1.d0)*rnff(k,2)*fv )*dts
end do
ccc end of main fluxes
ccc the rest of the fluxes (mostly for diagnostics)
aevapw = aevapw + evapvw*fw*(1.d0-fr_snow(2)*fm)*fv*dts
aevapd = aevapd + evapvd*fd*(1.d0-fr_snow(2)*fm)*fv*dts
aevapb = aevapb + evapb*(1.d0-fr_snow(1))*fb*dts
!!! I guess we need to accumulate evap from snow here
aepc=aepc+(epv*fv*dts)
aepb=aepb+(epb*fb*dts)
!Accumulate GPP, nyk, like evap_tot(2)
agpp = agpp + gpp*(1.d0-fr_snow(2)*fm)*fv*dts
dedifs=f(2,1)*tp(2,1)
if( aedifs=aedifs-dts*shw*dedifs*fb
dedifs=f(2,2)*tp(2,2)
if( aedifs=aedifs-dts*shw*dedifs*fv ! not used ?
af0dt=af0dt-dts*(fb*fh(1,1)+fv*fch(0)+htpr) ! E0 excludes htpr?
af1dt=af1dt-dts*(fb*fh(2,1)+fv*fh(2,2))
return
entry accmf 1
c provides accumulation units fixups, and calculates
c penman evaporation. should be called once after
c accumulations are collected.
aruns=rhow*aruns
arunu=rhow*arunu
aevapw=rhow*aevapw
aevapd=rhow*aevapd
aevapb=rhow*aevapb
aevap =rhow*aevap
aepc=rhow*aepc
aepb=rhow*aepb
af1dt=af1dt-aedifs
c**** calculation of penman value of potential evaporation, aepp
h0=fb*(snsh_tot(1)+elh*evap_tot(1))
& +fv*(snsh_tot(2)+elh*evap_tot(2))
c h0=-atrg/dt+srht+trht
ccc h0=-thrm(2)+srht+trht
el0=elh*1d-3
cpfac=sha*rho * ch*vsm
qsats=qsat(ts,lhe,pres)
dqdt = dqsatdt(ts,lhe)*qsats
epen=(dqdt*h0+cpfac*(qsats-qs))/(el0*dqdt+sha)
aepp=epen*dt
abetap=1.d0
if (aepp.gt.0.d0) abetap=(aevapw+aevapd+aevapb)/aepp
abetap=min(abetap,one)
abetap=max(abetap,zero)
ccc computing surface temperature from thermal radiation fluxes
tbcs = sqrt(sqrt( atrg/(dt*stbo) )) - tfrz
return
entry accm0 1
c zero out accumulations
atrg=0.d0 ! thermal heat from ground
ashg=0.d0 ! sensible heat from ground
aevap=0.d0 ! evaporation from ground
!alhg=0.d0 ! not accumulated : computed from aevap
aruns=0.d0
aeruns=0.d0
arunu=0.d0
aerunu=0.d0
abetad=0.d0 ! not accumulated : do we need it? YES
abetav=0.d0 ! not accumulated : do we need it?
abetat=0.d0 ! not accumulated : do we need it?
abetap=0.d0 ! not sure how it is computed : probably wrong
abetab=0.d0 ! not accumulated : do we need it?
abeta=0.d0 ! not accumulated : do we need it?
acna=0.d0 ! not accumulated : do we need it?
acnc=0.d0 ! not accumulated : do we need it?
agpp=0.d0 ! new accumulator, nyk 4/25/03
aevapw=0.d0 ! evap from wet canopy
aevapd=0.d0 ! evap from dry canopy
aevapb=0.d0 ! evap from bare soil (no snow)
aepc=0.d0 ! potential evap from canopy
aepb=0.d0 ! potential evap from bare soil (no snow)
aedifs=0.d0 ! heat transport by water
af0dt=0.d0 ! heat from ground - htpr
af1dt=0.d0 ! heat from 2-nd soil layer
! aepp=0.d0 ! not accumulated : computed in accmf
return
end subroutine accm
subroutine gdtm(dtm) 1,5
c**** calculates the maximum time step allowed by stability
c**** considerations.
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
real*8 ak1,ak2(2),betas(2),cna,xk2(2),dldz2,dqdt,rho3,sgmm
real*8, intent(out) :: dtm
real*8 dtm1,dtm2,dtm3,dtm4,t450,xk1
integer ibv,k
t450=450.d0
dqdt=dqsatdt(ts,lhe)*qsat(ts,lhe,pres)
c****
c**** first calculate timestep for water movement in soil.
sgmm=1.0d0
dldz2=0.d0
do ibv=i_bare,i_vege
do k=1,n
dldz2=max(dldz2,d(k,ibv)/dz(k)**2)
end do
end do
dtm=sgmm/(dldz2+1d-12)
if( dtm1=dtm
if ( dtm .lt. 0.d0 ) call stop_model('gdtm: dt1_ghy<0',255)
c****
c**** next calculate timestep for heat movement in soil.
do ibv=i_bare,i_vege
do k=1,n
xk1=xkh(k,ibv)
ak1=(shc(k,ibv)+((1.d0-fice(k,ibv))*shw+fice(k,ibv)*shi)
& *w(k,ibv))/dz(k)
dtm=min(dtm,.5d0*ak1*dz(k)**2/(xk1+1d-12))
end do
end do
dtm2=dtm
if ( dtm .lt. 0.d0 ) call stop_model('gdtm: dt2_ghy<0',255)
c****
c**** finally, calculate max time step for top layer bare soil
c**** and canopy interaction with surface layer.
c**** use timestep based on coefficient of drag
cna=ch*vsm
rho3=.001d0*rho
betas(1:2) = 1.d0 ! it''s an overkill but it makes the things
! simpler.
if(epb.le.0.d0)then
betas(1)=1.0d0
else
betas(1)=evapb/epb
endif
if(epv.le.0.d0)then
betas(2)=1.0d0
else
betas(2)=(evapvw*fw+evapvd*(1.d0-fw))/epv
endif
do ibv=i_bare,i_vege
k=2-ibv
xk2(ibv)=sha*rho*cna
& + betas(ibv)*rho3*cna*elh*dqdt
& + 8.d0*stbo*(tp(k,ibv)+tfrz)**3
ak2(ibv)=shc(k,ibv)+((1.d0-fice(k,ibv))*shw+fice(k,ibv)*shi)
& *w(k,ibv)
dtm=min(dtm,0.5*ak2(ibv)/(xk2(ibv)+1d-12))
if(ibv.eq.1)dtm3=dtm
if(ibv.eq.2)dtm4=dtm
c
c prevent oscillation of top snow layer
c if(isn(ibv).ne.0.or.snowd(ibv).ne.0.d0)then
c ak3(ibv)=.05d0*shi*spgsn
c dtm=min(dtm,ak3(ibv)/(xk2(ibv)+1.d-12))
c if(ibv.eq.1)dtm5=dtm
c if(ibv.eq.2)dtm6=dtm
c endif
end do
if(dtm.lt.5.d0)then
write(99,*) '*********** gdtm: ijdebug,fb,fv',ijdebug,fb,fv
write(99,*)'dtm',dtm1,dtm2,dtm3,dtm4
write(99,*)'xk2',xk2
write(99,*)'ak2',ak2
write(99,*)'snsh',snsh
write(99,*)'xlth',elh*evap_tot(1:2)
write(99,*)'dqdt',dqdt
write(99,*)'ts,tfrz',ts,tfrz
write(99,*)'dlt',tp(1,1)-ts+tfrz,tp(0,2)-ts+tfrz
call stop_model("gdtm: time step < 5 s",255)
endif
c****
return
end subroutine gdtm
subroutine outw(i) 2,3
c**** prints current state of soil for one grid box
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
use filemanager, only: openunit
integer i,k
integer, save :: ichn = 0
call wtab
if( ichn == 0 ) call openunit("soil_outw", ichn)
write(ichn,1000)
write(ichn,*)'general quantities (bare soil or vegetation)'
write(ichn,*)'dts,tag',dts,i,' after qsbal(0),retp(1),advnc(2)'
cc write(ichn,1021)
write(ichn,1045)
write(ichn,1023)'i= ',ijdebug/1000,'pr= ',pr,'ts= ',ts-tfrz,
* 'q1= ',0.d0
write(ichn,1023)'j= ',mod(ijdebug,1000),
* 'snowf= ',drips(1),'tg= ',0.d0-tfrz,'qs= ',qs
write(ichn,1043)'t1= ',0.d0-tfrz,'vg= ',0.d0,'ch= ',ch
write(ichn,1044)'vsm= ',vsm
write(ichn,1022)
write(ichn,1021)
write(ichn,1014)'bare soil fb = ',fb
1014 format(1x,a17,f4.2)
cc write(ichn,1021)
write(ichn,1025)
write(ichn,1026)
write(ichn,1027)
& snowd(1),rnf(1),evap_tot(1),xinfc(1),zw(1),debug_data%qb
write(ichn,1021)
write(ichn,1030)
write(ichn,1031)
do 100 k=1,n
write(ichn,1040)k,theta(k,1),tp(k,1),fice(k,1),rnff(k,1),f(k,1),
& h(k,1),xk(k,1),w(k,1),ws(k,1),shc(k,1),fh(k,1),ht(k,1),
& q(1,k),q(2,k),q(3,k),q(4,k)
100 continue
cc write(ichn,1021)
write(ichn,1022)
write(ichn,1021)
write(ichn,1014)'vegetation fv = ',fv
cc write(ichn,1021)
write(ichn,1035)
write(ichn,1036)
write(ichn,1037) snowd(2),rnf(2),evap_tot(2),xinfc(2),zw(2),
& debug_data%qc,evapvw,
* evapvd,dripw(2)+drips(2),fw
write(ichn,1021)
write(ichn,1030)
write(ichn,1031)
k=0
write(ichn,1049)k,theta(k,2),tp(k,2),fice(k,2),rnff(k,2),f(k,2),
& w(k,2),ws(k,2),shc(k,2),fh(k,2),ht(k,2)
do 200 k=1,n
write(ichn,1040)k,theta(k,2),tp(k,2),fice(k,2),rnff(k,2),f(k,2),
& h(k,2),xk(k,2),w(k,2),ws(k,2),shc(k,2),fh(k,2),ht(k,2),
& q(1,k),q(2,k),q(3,k),q(4,k)
200 continue
cc write(ichn,1021)
write(ichn,1022)
write(ichn,1021)
write(ichn,1055)
1055 format(1x,3x,9x,' kgm-2',2x,9x,' kgm-2',3x,9x,' kgm-2',
* 4x,9x,'1e6jm-2',3x,9x,'1e6jm-2')
write(ichn,1060) aruns,aevapw,aepc,0.d0,af0dt
1060 format(1x,3x,'aruns = ',f9.4,2x,'aevapw = ',0pf9.4,4x,'aepc = ',
* 0pf9.4,4x,'afhg = ',-6pf9.4,2x,'af0dt = ',-6pf9.4)
write(ichn,1065) arunu,aevapd,aepb,atrg,htpr*dts
1065 format(1x,3x,'arunu = ',f9.4,2x,'aevapd = ',0pf9.4,4x,'aepb = ',
* 0pf9.4,4x,'atrg = ',-6pf9.4,2x,'aphdt = ',-6pf9.4)
write(ichn,1070) aevapb,ashg,aeruns
1070 format(1x,3x,' ',9x,2x,'aevapb = ',0pf9.4,4x,' ',
* 9x,4x,'ashg = ',-6pf9.4,2x,'aerns = ',-6pf9.4)
write(ichn,1073) alhg,af1dt,aedifs
1073 format(1x,3x,' ',9x,2x,' ',9x,4x,' ',
* 9x,4x,'alhg = ',-6pf9.4,2x,'af1dt = ',-6pf9.4/
* 1x,3x,' ',9x,2x,' ',9x,4x,' ',
* 9x,4x,' ',2x,9x,'aedfs = ',-6pf9.4)
c**** more outw outputs
write(ichn,*)'thrm ',thrm_tot
write(ichn,*)'xlth ',elh*evap_tot(1:2)
write(ichn,*)'snsh ',snsh
write(ichn,*)'htpr,srht,trht ',htpr,srht,trht
return
1000 format(' ',121('='))
1001 format('1')
1010 format(1x,a20,f10.0)
1020 format(1x,a20,1pe12.4)
1021 format('0')
1022 format(' ',60('. '),'.')
1023 format(1x,a10,i10,a10,6pf10.2,2(a10,0pf8.2),a10,0pf8.4)
1043 format(1x,10x,10x,10x,10x,2(a10,0pf8.2),a10,0pf8.4)
1044 format(1x,10x,10x,38x,1(a10,f8.2))
1045 format(1x,20x,10x,' 1e-6ms-1',10x,4x,'t(c)',10x,4x,'ms-1')
1024 format(1x,6(a10,e10.2))
1015 format(1x,4(a8,f8.2))
1019 format(1x,12x,4(a8,f8.2))
1025 format(' ',5x,'snowd',7x,'rnf',6x,'evap',6x,'xinfc',
* 8x,'zw',8x,'qb')
1026 format(' ',' mh2o',2x,'1e-6ms-1',2x,'1e-6ms-1',3x,
* '1e-6ms-1',3x,' m',5x,' ')
1027 format(' ',0pf10.4,6pf10.4,6pf10.4,1x,6pf10.1,0pf10.4,
* 0pf10.4)
1030 format(' ',5x,'theta',3x,'tp',2x,'fice',4x,'runoff'
& ,8x,'fl',9x,'h',8x,'xk',6x,'w',5x,'ws',8x,
& 'shc',8x,'fh',8x,'ht',1x,'sand',1x,'loam',1x,'clay',1x,'peat')
1031 format(' ',5x,5x,2x,'(c)',2x,6x,'1e-6ms-1',2x,'1e-6ms-1',
& 3x,' m',2x,'1e-6ms-1',1x,' m',1x,' m',
& 1x,'1e6jm-3c-1',4x,' wm-2',3x,'1e6jm-2',4x,'%',4x,'%',4x,'%'
& ,4x,'%'/1x,125('-'))
1035 format(' ',5x,'snowd',7x,'rnf',6x,'evap',6x,'xinfc',
* 8x,'zw',8x,'qc',5x,'evapw',5x,'evapd',8x,'dr',8x,'fw')
1036 format(' ',' mh2o',2x,'1e-6ms-1',2x,'1e-6ms-1',3x,
* '1e-6ms-1',2x,' m',5x,' ',2x,'1e-6ms-1',2x,
* '1e-6ms-1',2x,'1e-6ms-1',' ')
1037 format(' ',0pf10.4,6pf10.4,6pf10.4,1x,6pf10.1,0pf10.4,
* 0pf10.4,6pf10.4,6pf10.4,6pf10.4,0pf10.2)
1040 format(1x,i3,f7.3,f5.1,f6.3,1p,6pf10.4,6pf10.4,0pf10.3,6pf10.4,
* 0pf7.4,0pf7.4,1x,-6pf10.4,0pf10.4,-6pf10.4,4(2pf5.1))
1049 format(1x,i3,f7.3,f5.1,f6.3,1p,6pf10.4,6pf10.4,10x,10x,
* 0pf7.4,0pf7.4,1x,-6pf10.4,0pf10.4,-6pf10.4,3(2pf5.1))
end subroutine outw
c****
subroutine wtab 2
c**** returns water table zw for ibv=1 and 2.d0
c**** input:
c**** zb - layer boundaries, m
c**** zc - soil centers, m
c**** dz - layer thicknesses, m
c**** h - soil potential of layers, m
c**** f - fluxes between layers, m s-1
c**** xk - conductivities of layers, m s-1
c**** output:
c**** zw(2) - water table for ibv=1 and 2, m
ccc include 'soils45.com'
c**** soils28 common block 9/25/90
real*8 denom,hmat,tol
integer ibv,k
tol=1d-6
do 100 ibv=1,2
c**** find non-saturated layer
do 10 k=n,1,-1
if( 10 continue
k=1
20 continue
c**** retrieve matric potential
c write(6,*)'ij,n,k,hmat,ibv,xkl,ibv',ijdebug,n,k,hmat,ibv,xk(k,ibv)
hmat=h(k,ibv)-zc(k)
c**** calculate denominator, and keep zw above zb(k+1)
if(xk(k,ibv).le.1d-20) then
denom=-2.d0*hmat/dz(k)
go to 90
end if
denom=max( 90 continue
c**** calculate water table
c write(6,*) 'denom',denom
zw(ibv)=zb(k)-sqrt(-2.d0*hmat*dz(k)/(denom+1d-20))
100 continue
return
end subroutine wtab
subroutine snow 1,2
use snow_drvm, only: snow_drv
implicit none
real*8 fmask(2), evapsn(2), devapsn_dt(2), canht
integer ibv
fmask(1) = 1.d0
fmask(2) = fm
evapsn(1) = evapbs
evapsn(2) = evapvs
devapsn_dt(1) = devapbs_dt
devapsn_dt(2) = devapvs_dt
canht = stbo*(tp(0,2)+tfrz)**4
!!! correction for canopy transmittance
ccc reset some fluxes for ibv hack
flmlt_scale(:) = 0.d0
fhsng_scale(:) = 0.d0
flmlt(:) = 0.d0
fhsng(:) = 0.d0
thrmsn(:) = 0.d0
flux_snow(:,:) = 0.d0
ccc remember initial snow water for tracers
wsn_for_tr(:,1) = wsn(:,1)*fr_snow(1)
wsn_for_tr(:,2) = wsn(:,2)*fr_snow(2)
do ibv=i_bare,i_vege
call snow_drv(
ccc input
& fmask(ibv), evapsn(ibv), snshs(ibv), srht, trht, canht,
& drips(ibv), dripw(ibv), htdrips(ibv), htdripw(ibv),
& devapsn_dt(ibv), dsnsh_dt, dts,
& tp(1,ibv), dz(1), nlsn, top_stdev,
ccc updated
& dzsn(1,ibv), wsn(1,ibv), hsn(1,ibv), nsn(ibv),
& fr_snow(ibv),
ccc output
& flmlt(ibv), fhsng(ibv), flmlt_scale(ibv),
& fhsng_scale(ibv), thrmsn(ibv), flux_snow(0,ibv)
& )
enddo
evapbs = evapsn(1)
evapvs = evapsn(2)
end subroutine snow
subroutine set_snow 1,10
!@sum set_snow extracts snow from the first soil layer and initializes
!@+ snow model prognostic variables
!@+ should be called when model restarts from the old restart file
!@+ ( which doesn''t contain new snow model (i.e. 3 layer) data )
c
c input:
c snowd(2) - landsurface snow depth
c w(k,2) - landsurface water in soil layers
c ht(k,2) - landsurface heat in soil layers
c fsn - heat of fusion
c shi - specific heat of ice
c shc(k,2) - heat capacity of soil layers
c
c output:
c dzsn(lsn,2) - snow layer thicknesses
c wsn(lsn,2) - snow layer water equivalent depths
c hsn(lsn,2) - snow layer heat contents
c tsn1(2) - snow top temperature
c isn(2) - 0 if no snow, 1 if snow
c nsn(2) - number of snow layers
c snowd(2)
c w(k,2)
c ht(k,2)
c
c calling sequence:
c
c assignment of w,ht,snowd
c call ghinij(i,j,wfc1)
c call set_snow
c note: only to be called when initializing from landsurface
c prognostic variables without the snow model.
c
use snow_model, only: snow_redistr, snow_fraction
integer ibv
c outer loop over ibv
do ibv=1,2
c initalize all cases to nsn=1
nsn(ibv)=1
ccc since we don''t know what kind of data we are dealing with,
ccc better check it
if( snowd(ibv) .gt. w(1,ibv)-dz(1)*thetm(1,ibv) ) then
write(99,*) 'snowd corrected: old=', snowd(ibv)
snowd(ibv) = w(1,ibv)-dz(1)*thetm(1,ibv) - 1.d-10
write(99,*) ' new=', snowd(ibv)
if ( snowd(ibv) .lt. -0.001d0 )
& call stop_model('set_snow: neg. snow',255)
if ( snowd(ibv) .lt. 0.d0 ) snowd(ibv) = 0.d0 ! rounding error
endif
c if there is no snow, set isn=0. set snow variables to 0.d0
if(snowd(ibv).le.0.d0)then
dzsn(1,ibv)=0.d0
wsn(1,ibv)=0.d0
hsn(1,ibv)=0.d0
tsn1(ibv)=0.d0
fr_snow(ibv) = 0.d0
else
c given snow, set isn=1.d0
c!!! dzsn(1,ibv)=snowd(ibv)/spgsn
c!!! replacing prev line considering rho_snow = 200
dzsn(1,ibv)=snowd(ibv) * 5.d0
wsn(1,ibv)=snowd(ibv)
c!!! actually have to compute fr_snow and modify dzsn ...
fr_snow(ibv) = 1.d0
c given snow, temperature of first layer can''t be positive.
c the top snow temperature is the temperatre of the first layer.
if(fsn*w(1,ibv)+ht(1,ibv).lt.0.d0)then
tsn1(ibv)=(ht(1,ibv)+w(1,ibv)*fsn)/(shc(1,ibv)+w(1,ibv)*sh