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
c**** Added decks parameter vegCO2X_off 3/2/04 nyk
#include "rundeck_opts.h"
!#define EVAP_VEG_GROUND
!#define RAD_VEG_GROUND
module sle001 33,5
use constant
, only : stbo,tfrz=>tf,sha,lhe,one,zero,rhow
& ,shw_kg=>shw,shi_kg=>shi,lhm
use ghy_com, only : ngm, imt, nlsn, LS_NFRAC !, wsn_max
#ifdef TRACERS_WATER
* ,ntixw
use tracer_com, only : ntm,tr_wd_TYPE,nWater
#endif
implicit none
save
private
c**** public functions:
public hl0, set_snow, advnc, evap_limits, xklh
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,aflmlt,aintercep
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,ch,qm1,vs,vs0,tprime,qprime
!@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,LS_NFRAC),ht(0:ngm,LS_NFRAC),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
!@var epb,epbs,epv,epvs potential evaporation (b-bare, v-vege, s-snow)
real*8 epb, epv, epbs, epvs, epvg
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)
#ifdef TRACERS_WATER
ccc should be passed from elsewhere
!@var ntgm maximal number of tracers that can by passesd to HGY
integer, parameter, public :: ntgm = ntm
!@var ntg actual number of tracers passed to ground hydrology
integer, public :: ntg
!@var trpr flux of tracers in precipitation (?/m^2 s)
!@var trdd flux of tracers as dry deposit (?/m^2 s)
!@var tr_w amount of tracers in the soil (?)
!@var tr_wsn amount of tracers in the snow (?)
real*8, public :: trpr(ntgm), trdd(ntgm), tr_surf(ntgm),
& tr_w(ntgm,0:ngm,2), tr_wsn(ntgm,nlsn,2)
ccc tracers output:
!@var tr_evap flux of tracers to atm due to evaporation (?/m^2 s)
!@var tr_evap flux of tracers due to runoff (?/m^2 s)
real*8 tr_evap(ntgm,2),tr_rnff(ntgm,2)
!@var tr_evap flux of tracers to atm due to evaporation
!@var tr_evap flux of tracers due to runoff
real*8, public :: atr_evap(ntgm),atr_rnff(ntgm),atr_g(ntgm)
#endif
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,aflmlt,aintercep
& ,ashg,atrg,betad,betat,ch,gpp,d,devapbs_dt,devapvs_dt
& ,drips,dripw,dsnsh_dt,dts,dz,dzsn,epb,epbs,epvs,epvg ! 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,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
& ,vs,vs0,tprime,qprime
!----------------------------------------------------------------------!
& ,i_bare,i_vege,process_bare,process_vege
#ifdef TRACERS_WATER
& ,trpr,trdd, tr_surf, tr_w, tr_wsn,tr_evap,tr_rnff ! ntg
& ,atr_evap,atr_rnff,atr_g
#endif
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 6
!@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)=(w(0,2)/ws(0,2))**(2.d0/3.d0)
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 8
!@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 2
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*(a(1,i)+a(2,i)+a(3,i)+a(4,i)))
a1=a(3,i)/a(4,i)
a2=(a(2,i)-(log(-hlm(j)-hs))/c)/a(4,i)
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 ) 4,25
!@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
real*8 v_qprime ! local variable
#ifdef EVAP_VEG_GROUND
real*8 evap_max_vegsoil
#endif
c cna is the conductance of the atmosphere
cna=ch*vs
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) +
& (w(k,ibv)-dz(k)*thetm(k,ibv))/dt
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
#ifdef EVAP_VEG_GROUND
evap_max_vegsoil = min( evap_max(ibv),
& 2.467d0*d(1,2)*(theta(1,2)-thetm(1,2))/dz(1) + pr )
#endif
! 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
& )
!print *,"HGY_COND: ",ijdebug, cnc, betadl
!print *,"GHY_FORCINGS: ", ijdebug, tp(0,2)
!!!! test
!!! trans_sw = .1d0
!!! print *,'trans_sw = ', trans_sw
!!! trans_sw = .1d0
! trans_sw = min( trans_sw, .9d0 )
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.
c pot_evap_can = betat*rho3*cna*(qsat(tp(0,2)+tfrz,lhe,pres) - qs)
pot_evap_can = betat*rho3*ch*(
& vs*(qsat(tp(0,2)+tfrz,lhe,pres) - qs)
& -(vs-vs0)*qprime)
evap_max_dry(ibv) = 0.d0
if ( betad > 0.d0 .and. pot_evap_can > 0.d0 ) then
do k=1,n
evap_max_dry(ibv) = evap_max_dry(ibv)
& + min( pot_evap_can*betadl(k)/betad,
& (w(k,ibv)-dz(k)*thetm(k,ibv))/dt )
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*( (1. - fr_snow(2)*fm) * ( 1. - theta(0,2) )
& * 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)
#ifdef EVAP_VEG_GROUND
qvg = qsat(tp(1,2)+tfrz,lhe,pres)
#endif
c potential evaporation for bare and vegetated soil and snow
c epb = rho3*cna*(qb-qs)
c epbs = rho3*cna*(qbs-qs)
c epv = rho3*cna*(qv-qs)
c epvs = rho3*cna*(qvs-qs)
v_qprime=(vs-vs0)*qprime
epb = rho3*ch*( vs*(qb-qs) -v_qprime )
epbs = rho3*ch*( vs*(qbs-qs)-v_qprime )
epv = rho3*ch*( vs*(qv-qs) -v_qprime )
epvs = rho3*ch*( vs*(qvs-qs)-v_qprime )
#ifdef EVAP_VEG_GROUND
c epvg = rho3*cna*(qvg-qs) ! actually not correct !
epvg = rho3*ch*( vs*(qvg-qs)-v_qprime )
#endif
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 )
#ifdef EVAP_VEG_GROUND
evapvg = min(epvg,evap_max_vegsoil)
! keep evapv <= epv
evapvg = min( evapvg, epv - evapvd*fd - evapvw*fw )
evapvg = max( evapvg, 0.d0 )
#endif
! 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
#ifdef EVAP_VEG_GROUND
evapvg = 0.d0
#endif
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 2
!@sum computes sensible heat for bare/vegetated soil and snow
implicit none
real*8 cna,v_tprime
c cna=ch*vs
c snsh(1)=sha*rho*cna*(tp(1,1)-ts+tfrz) ! bare soil
c snsh(2)=sha*rho*cna*(tp(0,2)-ts+tfrz) ! canopy
c snshs(1) = sha*rho*cna*(tsn1(1)-ts+tfrz) ! bare soil snow
c snshs(2) = sha*rho*cna*(tsn1(2)-ts+tfrz) ! canopy snow
c dsnsh_dt = sha*rho*cna ! derivative is the same for all above
cna=ch*vs
v_tprime=(vs-vs0)*tprime
snsh(1)=sha*rho*ch*(vs*(tp(1,1)-ts+tfrz)-v_tprime)
! bare soil
snsh(2)=sha*rho*ch*(vs*(tp(0,2)-ts+tfrz)-v_tprime)
! canopy
snshs(1)=sha*rho*ch*(vs*(tsn1(1)-ts+tfrz)-v_tprime)
! bare soil snow
snshs(2)=sha*rho*ch*(vs*(tsn1(2)-ts+tfrz)-v_tprime)
! 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**** vs - 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 2,2
!@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
#define TRY_TO_MELT_FRESH_SNOW_ON_WARM_GROUND
#ifdef TRY_TO_MELT_FRESH_SNOW_ON_WARM_GROUND
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
endif
return
end subroutine drip_from_canopy
subroutine flg 2
!@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 2
!@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
#ifdef RAD_VEG_GROUND
& * (1.d0-trans_sw)
& - trans_sw*(srht + trht)
#endif
& )
& *(1.d0-fr_snow(2))
! canopy
fch(0) = -htpr +
& (evapvw*elh*fw + snsh(2) + thrm_can
#ifdef RAD_VEG_GROUND
& * (1.d0-trans_sw)
& - (1.d0 - trans_sw)*(srht + trht)
#else
& -srht-trht
#endif
& + evapvd*elh*fd)
& *(1.d0-fm*fr_snow(2))
fch(1) = -(thrm_can - thrm_soil(2))
#ifdef RAD_VEG_GROUND
& *(1.d0 - trans_sw)
#endif
& *(1.d0-fr_snow(2)) !rad soil
& - (thrm_can - thrmsn(2))*fr_snow(2)*(1.d0-fm) !rad snow
#ifdef RAD_VEG_GROUND
& *(1.d0 - trans_sw)
#endif
& - 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)
#ifdef RAD_VEG_GROUND
& *(1.d0-trans_sw)
#endif
& + thrmsn(2)*fr_snow(2)*fm
#ifdef RAD_VEG_GROUND
& + thrmsn(2)*trans_sw*fr_snow(2)*(1.d0-fm)
& + trans_sw*thrm_soil(2)*(1.d0-fr_snow(2))
#endif
!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 2
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 = (w(1,ibv)/ws(1,ibv))**rosmp
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
!/* #define do_topmodel_runoff */
#ifdef do_topmodel_runoff
if ( ws(k,ibv) > 1.d-16 ) then
f_k0_exp_k = (1.d0-fice(k,ibv))
$ * xkus(k,ibv)*w(k,ibv)/ws(k,ibv)*dz(k)
else
f_k0_exp_k = 0.d0
endif
c print *,'rnff: ', k, rnff(k,ibv), f_k0_exp_k*exp( -top_index )
c print *, xkus(k,ibv),w(k,ibv),ws(k,ibv),dz(k),top_index
rnff(k,ibv)=f_k0_exp_k*exp( -top_index )
#endif
end do
! print *,'runoff: ', sl/sdstnc, exp( -top_index )
end do
return
end subroutine runoff
subroutine fllmt 2,4
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) + (w(k,ibv)-dz(k)*thetm(k,ibv))/dts
& - 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
!!! hack to prevent taking water from empty layers
!!! this is probably not quite correct so it is disabled by default
cddd do ibv=i_bare,i_vege
cddd do k=n,1,-1
cddd if ( f(k,ibv) > 0.d0 .and. w(k,ibv) < 1.d-16 ) then
cddd print *,"GHY: corrected flux->0 ",k,ibv,f(k,ibv)
cddd f(k,ibv) = 0.d0
cddd endif
cddd enddo
cddd enddo
return
end subroutine fllmt
subroutine xklh( xklh0_flag ) 4,2
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.
integer, intent(in), optional :: xklh0_flag
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 /)
if ( present(xklh0_flag) ) goto 777
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
777 continue
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 2
!@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)*(h(k-1,ibv)-h(k,ibv))/(zc(k-1)-zc(k))
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 2
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(f(k,ibv).gt.0)then
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 4,8
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(w(k,ibv).ge.1d-12)then
! 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.120.d0.or.tp(0,2).gt.120.d0
& .or. tp(1,1)<-150.d0 .or. tp(0,2)<-150.d0 )then
write(6,*)'retp tp bounds error'
write(6,*)'ijdebug',ijdebug
call reth
call hydra
call outw(1)
call stop_model(
& 'retp: tground > 120C - see soil_outw and fort.99',255)
endif
return
end subroutine retp
subroutine apply_fluxes 2,4
!@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)write(0,*)'GHY:CanopyH2O<0 at',ijdebug,w(0,2)
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 2,66
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 accm(0)
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)
#ifdef TRACERS_WATER
call ghy_tracers
#endif
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
!!!
! if(wsn_max>0) call restrict_snow (wsn_max)
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, vs,qs)
#ifdef TRACERS_WATER
C**** finalise surface tracer concentration here
tot_w1 = fb*( w(1,1)*(1.d0-fr_snow(1))
& + wsn(1,1)*fr_snow(1) )
& + fv*( w(0,2)*(1.d0-fm*fr_snow(2))
& + wsn(1,2)*fm*fr_snow(2) )
& + fv*w(1,2)*0.01 ! hack !!! (for 0 H2O in canopy)
if ( tot_w1 > 1.d-30 ) then
atr_g(:ntg) = ! instantaneous
& ( fb*( tr_w(:ntg,1,1)*(1.d0-fr_snow(1))
& + tr_wsn(:ntg,1,1) ) !*fr_snow(1)
& + fv*( tr_w(:ntg,0,2)*(1.d0-fm*fr_snow(2))
& + tr_wsn(:ntg,1,2)*fm )
& + fv*tr_w(:ntg,1,2)*0.01 ! hack !!! (for 0 H2O in canopy)
& ) / !*fr_snow(2)
& (rhow * tot_w1) ! * dts
endif
#endif
enddo
call accm(1)
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( flag ) 6,4
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.
integer, intent(in), optional :: flag
real*8 qsats
real*8 cpfac,dedifs,dqdt,el0,epen,h0
integer k
#ifdef TRACERS_WATER
real*8 tot_w1
#endif
if ( present(flag) ) then
if ( flag == 0 ) goto 778
goto 777
endif
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),0.d0)*rnff(k,1)*fb
* + max(tp(k,2),0.d0)*rnff(k,2)*fv )*dts
end do
ccc end of main fluxes
ccc the rest of the fluxes (mostly for diagnostics)
aflmlt = aflmlt + ( fb*(flmlt(1)*fr_snow(1)+flmlt_scale(1))
$ + fv*(flmlt(2)*fr_snow(2)+flmlt_scale(2)) )*dts
! if(process_vege) aintercep = aintercep + pr - dripw(2) - drips(2)
ccc max in the following expression removes extra drip because of dew
aintercep = aintercep + fv*max(pr - dripw(2) - drips(2),0.d0)*dts
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)
aepc=aepc+( epv*(1.d0-fr_snow(2)*fm) + epvs*fr_snow(2)*fm )*fv*dts
#ifdef EVAP_VEG_GROUND
! next line is probably not correct (need to use something more
! sophisticated for epvg)
aepc=aepc+( epvg*(1.d0-fr_snow(2)) )*fv*dts
#endif
aepb=aepb+( epb*(1.d0-fr_snow(1)) + epbs*fr_snow(1) )*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(f(2,1).lt.0.d0) dedifs=f(2,1)*tp(1,1)
aedifs=aedifs-dts*shw*dedifs*fb
dedifs=f(2,2)*tp(2,2)
if(f(2,2).lt.0.d0) dedifs=f(2,2)*tp(1,2)
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))
#ifdef TRACERS_WATER
ccc accumulate tracer fluxes
atr_evap(:ntg) = atr_evap(:ntg)
& + ( tr_evap(:ntg,1)*fb + tr_evap(:ntg,2)*fv )*dts
atr_rnff(:ntg) = atr_rnff(:ntg)
& + ( tr_rnff(:ntg,1)*fb + tr_rnff(:ntg,2)*fv )*dts
C**** no point in setting this here since fm will change before end
c tot_w1 = fb*( w(1,1)*(1.d0-fr_snow(1))
c & + wsn(1,1)*fr_snow(1) )
c & + fv*( w(0,2)*(1.d0-fm*fr_snow(2)) ! *fw0 !! remove fw ?
c & + wsn(1,2)*fm*fr_snow(2) )
c if ( tot_w1 > 1.d-30 ) then
c atr_g(:ntg) = ! instantaneous ! atr_g(:ntg) +
c & ( fb*( tr_w(:ntg,1,1)*(1.d0-fr_snow(1))
c & + tr_wsn(:ntg,1,1) ) !*fr_snow(1)
c & + fv*( tr_w(:ntg,0,2)*(1.d0-fm*fr_snow(2)) ! *fw0 !! remove fw?
c & + tr_wsn(:ntg,1,2)*fm ) ) / !*fr_snow(2)
c & (rhow * tot_w1) ! * dts
c endif
cddd print '(a,100(e12.4))','tr_evap_err',
cddd & ( tr_evap(1,1)*fb + tr_evap(1,2)*fv )/1000.d0 -
cddd & ( evap_tot(1)*fb + evap_tot(2)*fv ),
cddd & ( evap_tot(1)*fb + evap_tot(2)*fv ),
cddd & evapvs, fr_snow, fv, fm
! & atr_evap(1) - aevap*1000.d0, aevap*1000.d0, evapvs, fr_snow
#endif
return
! entry accmf
777 continue
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
aflmlt=rhow*aflmlt
aintercep=rhow*aintercep
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*vs
qsats=qsat(ts,lhe,pres)
dqdt = dqsatdt(ts,lhe)*qsats
! epen=(dqdt*h0+cpfac*(qsats-qs))/(el0*dqdt+sha)
epen=(dqdt*h0+sha*rho*ch*
& ( vs*(qsats-qs)-(vs-vs0)*qprime ))/(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
778 continue
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
aflmlt=0.d0
aintercep=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
#ifdef TRACERS_WATER
ccc tracers
atr_evap(:ntg) = 0.d0
atr_rnff(:ntg) = 0.d0
atr_g(:ntg) = 0.d0
#endif
return
end subroutine accm
subroutine gdtm(dtm) 2,10
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(q(4,1).gt.0.d0)dtm=min(dtm,t450)
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*vs
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) 4,6
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)'vs= ',vs
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), 0.d0,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 4
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(w(k,ibv).lt.ws(k,ibv)*(1.d0-tol))go to 20
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(f(k,ibv)/xk(k,ibv)+1.d0,-2.d0*hmat/dz(k))
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 2,4
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
#ifdef RAD_VEG_GROUND
canht = canht*(1.d0-trans_sw) + (srht+trht)*trans_sw
#endif
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 2,20
!@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)*shi)
else
tsn1(ibv)=0.d0
endif
c use snow temperature to get the heat of the snow
hsn(1,ibv)=tsn1(ibv)*wsn(1,ibv)*shi-wsn(1,ibv)*fsn
c subtract the snow from the landsurface prognositic variables
w(1,ibv)=w(1,ibv)-wsn(1,ibv)
ht(1,ibv)=ht(1,ibv)-hsn(1,ibv)
ccc and now limit all the snow to 5cm water equivalent
if ( snowd(ibv) .gt. 0.05d0 ) then
snowd(ibv) = 0.05d0
dzsn(1,ibv)= snowd(ibv) * 5.d0
wsn(1,ibv)= snowd(ibv)
hsn(1,ibv)= tsn1(ibv)*wsn(1,ibv)*shi-wsn(1,ibv)*fsn
endif
ccc and now limit all the snow to 50cm water equivalent (for debug)
cddd if ( snowd(ibv) .gt. 0.0005d0 ) then
cddd snowd(ibv) = 0.5d0
cddd dzsn(1,ibv)= snowd(ibv) * 5.d0
cddd wsn(1,ibv)= snowd(ibv)
cddd hsn(1,ibv)= tsn1(ibv)*wsn(1,ibv)*shi-wsn(1,ibv)*fsn
cddd endif
ccc redistribute snow over the layers and recompute fr_snow
ccc (to make the data compatible with snow model)
if ( .not. ( hsn(1,ibv) > 0. .or. hsn(1,ibv) <= 0. ) )
& call stop_model("ERR in init_snow: NaN", 255)
call snow_fraction(dzsn(:,ibv), nsn(ibv), 0.d0, 0.d0,
& 1.d0, fr_snow(ibv) )
if ( fr_snow(ibv) > 0.d0 ) then
call snow_redistr(dzsn(:,ibv), wsn(:,ibv), hsn(:,ibv),
& nsn(ibv), 1.d0/fr_snow(ibv) )
if ( .not. ( hsn(1,ibv) > 0. .or. hsn(1,ibv) <= 0. ) )
& call stop_model("ERR in init_snow 2: NaN", 255)
else
snowd(ibv) = 0.d0
dzsn(1,ibv)=0.d0
wsn(1,ibv)=0.d0
hsn(1,ibv)=0.d0
tsn1(ibv)=0.d0
fr_snow(ibv) = 0.d0
endif
endif
!!! debug : check if snow is redistributed correctly
if ( dzsn(1,ibv) > 0.d0 .and. dzsn(1,ibv) < .099d0) then
call stop_model("set_snow: error in dz",255)
endif
if ( dzsn(1,ibv) > 0.d0 ) then
if ( wsn(1,ibv)/dzsn(1,ibv) < .1d0)
& call stop_model("set_snow: error in dz",255)
endif
enddo
if ( .not. ( hsn(1,1) > 0. .or. hsn(1,1) <= 0. ) )
& call stop_model("ERR in init_snow 2: NaN", 255)
if ( .not. ( hsn(1,2) > 0. .or. hsn(1,2) <= 0. ) )
& call stop_model("ERR in init_snow 2: NaN", 255)
return
end subroutine set_snow
#ifdef TRACERS_WATER
subroutine check_wc(wc) 18
real*8 wc
!!! removing test...
return
!!!! ;-(
!wc = 1000.d0
!return
if ( abs(wc-1000.d0) > 1.d-3 ) then
print *,"GHYTR: wc= ", wc, (wc-1000.d0)/1000.d0
!call stop_model("GHY: wc != 1000",255)
endif
end subroutine check_wc
subroutine ghy_tracers_drydep 2
!@sum applies dry deposit flux to tracers in upper layers of
!@+ canopy, soil and snow
ccc input from outside:
ccc trdd(ntgm) - flux of tracers as dry deposit (kg /m^2 /s)
!@var m number of tracers (=ntg)
integer m
m = ntg
if ( process_vege ) then
! +drydep
tr_w(:m,0,2) = tr_w(:m,0,2) + trdd(:m)*(1-fm*fr_snow(2))*dts
tr_wsn(:m,1,2) = tr_wsn(:m,1,2) + trdd(:m)*fm*fr_snow(2)*dts
endif
if ( process_bare ) then
! +drydep
tr_w(:m,1,1) = tr_w(:m,1,1) + trdd(:m)*(1-fr_snow(1))*dts
tr_wsn(:m,1,1) = tr_wsn(:m,1,1) + trdd(:m)*fr_snow(1)*dts
endif
end subroutine ghy_tracers_drydep
subroutine ghy_tracers 4,48
ccc will need the following data from GHY:
ccc rnff(k,ibv) - underground runoff fro layers 1:n
ccc rnf(ibv) - surface runoff
ccc evapd*betadl(k)/(betad+1d-12) - transpiration from 1:n (ibv=2 only)
ccc f(k,ibv) - flux between layers (+up) (upper edge)
ccc will need from SNOW:
ccc tr_flux(0:nsl) - flux down
ccc will do it as follows:
ccc 1. propagate tracers to canopy layer (no flux up on lower boundary)
ccc 2. propagate tracers through the snow (no flux up from soil)
ccc 3. propagate them through the soil
ccc input from outside:
ccc pr - precipitation (m/s) ~ (10^3 kg /m^2 /s)
ccc trpr(ntgm) - flux of tracers (kg /m^2 /s)
ccc trdd(ntgm) - flux of tracers as dry deposit (kg /m^2 /s)
ccc tr_surf(ntgm) - surface concentration of tracers (kg/m^3 H2O)
ccc **************************************************
ccc !!! test
ccc internal vars:
!@var wi instantaneous amount of water in soil layers (m)
!@var tr_wc ratio tracers/water in soil layers (?/m)
!@var tr_wcc ratio tracers/water in canopy (?/m)
real*8 wi(0:ngm,1:2),tr_wc(ntgm,0:ngm,1:2),tr_wcc(ntgm,1:2)
!@var wsni instantaneous amount of water in snow layers (m)
!@var tr_wsnc ratio tracers/water in snow layers (?/m)
real*8 wsni(nlsn,2),tr_wsnc(ntgm,0:nlsn,2)
!@var flux_snow_tot water flux between snow layers * fr_snow (m/s)
real*8 flux_snow_tot(0:nlsn,2)
!@var flux_dt amount of water moved in current context (m)
real*8 flux_dt, err, evap_tmp(ntgm), flux_dtm(ntgm)
integer ibv,i
!@var m = ntg number of ground tracers (to make notations shorter)
integer m,mc
real*8, parameter :: EPSF = 1.d-20 ! 1.d-11
#ifdef TRACERS_SPECIAL_O18
real*8, external :: fracvl
#endif
!@var tol error tolerance for each type of tracer
real*8 tol(ntg)
ccc for debug
real*8 tr_w_o(ntgm,0:ngm,2), tr_wsn_o(ntgm,nlsn,2)
if ( ntg < 1 ) return ! no water tracers
!m = 2 !!! testing
m = ntg
ccc set error tolerance for each time of tracer (can be increased
ccc if model stops due to round-off error)
do mc=1,m
tol(mc) = ( sum(tr_w(mc,:,:))/(size(tr_w,2)*size(tr_w,3))
& + sum(tr_wsn(mc,:,:))/(size(tr_wsn,2)*size(tr_wsn,3))
& + trpr(mc)*dts + trdd(mc)*dts + tr_surf(mc)*1.d-8*dts )
& *1.d-10 + 1.d-32
enddo
!!! for test only
! trpr(:) = .5d0 * pr
! tr_surf(:) = 1000.d0 !!!???
!tr_w = 0.d0
!tr_wsn = 0.d0
#ifdef DEBUG_TR_WITH_WATER
if ( abs(tr_surf(1)-1000.) > 1.d-5 )
& call stop_model("GHY: tr_surf != 1000",255)
#endif
cddd print *, 'starting ghy_tracers'
cddd print *, 'trpr, trpr_dt ', trpr(1), trpr(1)*dts
cddd print *, 'tr_w 1 ', tr_w(1,:,1)
cddd print *, 'tr_w 2 ', tr_w(1,:,2)
cddd print *, 'tr_wsn 1 ', tr_wsn(1,:,1)
cddd print *, 'tr_wsn 2 ', tr_wsn(1,:,2)
!!! test
! tr_wsn(1,:,1) = wsn_for_tr(:,1) * fr_snow(1) * 1000.d0
! tr_wsn(1,:,2) = wsn_for_tr(:,2) * fr_snow(2) * 1000.d0
tr_w_o(:m,:,i_bare:i_vege) = tr_w(:m,:,i_bare:i_vege)
tr_wsn_o(:m,:,i_bare:i_vege) = tr_wsn(:m,:,i_bare:i_vege)
ccc set internal vars
do ibv=i_bare,i_vege
wi(0:n,ibv) = w(0:n,ibv)
wsni(1:nlsn,ibv) = wsn_for_tr(1:nlsn,ibv) ! *fr_snow(ibv)
flux_snow_tot(0:nlsn,ibv) = flux_snow(0:nlsn,ibv)*fr_snow(ibv)
!!! hack
!!! snow doesnt suck water from ground
if ( flux_snow_tot(nlsn,ibv) < 0.d0 ) then
flux_snow_tot(0:nlsn-1,ibv) = flux_snow_tot(0:nlsn-1,ibv)
& - flux_snow_tot(nlsn,ibv)
flux_snow_tot(nlsn,ibv) = 0.d0
endif
enddo
ccc reset accumulators to 0
tr_evap(:m,1:2) = 0.d0
tr_rnff(:m,1:2) = 0.d0
ccc apply dry deposit tracers here (may need to be removed outside
ccc but will keep it here for now for compatibility with conservation tests
call ghy_tracers_drydep
ccc canopy
tr_wcc(:m,:) = 0.d0
!!!test
!!! tr_wcc(:m,:) = 1000.d0
if ( process_vege ) then
! +precip
tr_w(:m,0,2) = tr_w(:m,0,2) + trpr(:m)*dts
wi(0,2) = wi(0,2) + pr*dts
if ( wi(0,2) > 0.d0 ) tr_wcc(:m,2) = tr_w(:m,0,2)/wi(0,2)
call check_wc(tr_wcc(1,2))
! +- evap
!gg if (tr_wd_TYPE(ntixw(mc)).eq.nWATER) THEN
evap_tmp(:m) = 0.d0
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& evap_tmp(mc) = fc(0)+pr
if ( evapvw >= 0.d0 ) then ! no dew
!!!evap_tmp(:m) = fc(0)+pr
#ifdef TRACERS_SPECIAL_O18
if ( evap_tmp(1)*dts < wi(0,2) .and. tp(0,2) > 0.d0 ) then
do mc=1,m
ccc loops...
evap_tmp(mc) = evap_tmp(mc) * fracvl(tp(0,2),ntixw(mc))
enddo
endif
#endif
tr_evap(:m,2) = tr_evap(:m,2) + evap_tmp(:m)*tr_wcc(:m,2)
tr_w(:m,0,2) = tr_w(:m,0,2) - evap_tmp(:m)*tr_wcc(:m,2)*dts
else ! dew adds tr_surf to canopy
tr_evap(:m,2) = tr_evap(:m,2) + evap_tmp(:m)*tr_surf(:m)
tr_w(:m,0,2) = tr_w(:m,0,2) - evap_tmp(:m)*tr_surf(:m)*dts
wi(0,2) = wi(0,2) - (fc(0)+pr)*dts
tr_wcc(:m,2) = tr_w(:m,0,2)/wi(0,2)
endif
!gg end if
call check_wc(tr_wcc(1,2))
! -drip
tr_w(:m,0,2) = tr_w(:m,0,2) + fc(1)*tr_wcc(:m,2)*dts
endif
! trivial value for bare soil
if ( pr > 0.d0 ) tr_wcc(:m,1) = trpr(:m)/pr
call check_wc(tr_wcc(1,1))
ccc end canopy
ccc snow
!!!>> tr_wsnc(:m,:,:) = 0.d0 !!! was 0
tr_wsnc(:m,:,:) = 0.d0 !!! was 1000
! dew
!gg if (tr_wd_TYPE(ntixw(mc)).eq.nWATER) THEN
if ( process_bare .and. evapbs < 0.d0 ) then
flux_dt = - evapbs*fr_snow(1)*dts
flux_dtm(:m) = 0.d0
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& flux_dtm(mc) = flux_dt
tr_evap(:m,1) = tr_evap(:m,1) - flux_dtm(:m)/dts*tr_surf(:m)
tr_wsn(:m,1,1) = tr_wsn(:m,1,1) + flux_dtm(:m)*tr_surf(:m)
wsni(1,1) = wsni(1,1) + flux_dt
endif
if ( process_vege .and. evapvs < 0.d0 ) then
flux_dt = - evapvs*fm*fr_snow(2)*dts
flux_dtm(:m) = 0.d0
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& flux_dtm(mc) = flux_dt
tr_evap(:m,2) = tr_evap(:m,2) - flux_dtm(:m)/dts*tr_surf(:m)
tr_wsn(:m,1,2) = tr_wsn(:m,1,2) + flux_dtm(:m)*tr_surf(:m)
wsni(1,2) = wsni(1,2) + flux_dt
endif
!gg end if
do ibv=i_bare,i_vege
! init tr_wsnc
do i=1,nlsn
if ( wsni(i,ibv) > 0.d0 ) then
tr_wsnc(:m,i,ibv) = tr_wsn(:m,i,ibv)/wsni(i,ibv)
else
tr_wsnc(:m,i,ibv) = 0.d0
endif
enddo
if ( fr_snow(ibv) > 0.d0 ) then ! process snow
flux_dt = (drips(ibv) + dripw(ibv)*fr_snow(ibv))*dts
tr_wsn(:m,1,ibv) = tr_wsn(:m,1,ibv) + tr_wcc(:m,ibv)*flux_dt
wsni(1,ibv) = wsni(1,ibv) + flux_dt
if ( wsni(1,ibv) > 0.d0 )
& tr_wsnc(:m,1,ibv) = tr_wsn(:m,1,ibv)/wsni(1,ibv)
!!!
! tr_wsnc(:m,1,ibv) = 1000.d0
call check_wc(tr_wsnc(1,1,ibv))
! sweep down
do i=2,nlsn
if ( flux_snow_tot(i-1,ibv) > EPSF ) then
flux_dt = flux_snow_tot(i-1,ibv)*dts
tr_wsn(:m,i,ibv) = tr_wsn(:m,i,ibv)
& + tr_wsnc(:m,i-1,ibv)*flux_dt
wsni(i,ibv) = wsni(i,ibv) + flux_dt
if ( wsni(i,ibv) > 0.d0 )
& tr_wsnc(:m,i,ibv) = tr_wsn(:m,i,ibv)/wsni(i,ibv)
!!!
! tr_wsnc(:m,i,ibv) = 1000.d0
call check_wc(tr_wsnc(1,i,ibv))
tr_wsn(:m,i-1,ibv) = tr_wsn(:m,i-1,ibv)
& - tr_wsnc(:m,i-1,ibv)*flux_dt
wsni(i-1,ibv) = wsni(i-1,ibv) + flux_dt ! extra?
endif
enddo
! sweep up
do i=nlsn-1,1,-1
if ( flux_snow_tot(i,ibv) < -EPSF ) then
flux_dt = - flux_snow_tot(i,ibv)*dts
tr_wsn(:m,i,ibv) = tr_wsn(:m,i,ibv)
& + tr_wsnc(:m,i+1,ibv)*flux_dt
wsni(i,ibv) = wsni(i,ibv) + flux_dt
if ( wsni(i,ibv) > 0.d0 )
& tr_wsnc(:m,i,ibv) = tr_wsn(:m,i,ibv)/wsni(i,ibv)
!!!
! tr_wsnc(:m,i,ibv) = 1000.d0
call check_wc(tr_wsnc(1,i,ibv))
tr_wsn(:m,i+1,ibv) = tr_wsn(:m,i+1,ibv)
& - tr_wsnc(:m,i+1,ibv)*flux_dt
wsni(i+1,ibv) = wsni(i+1,ibv) - flux_dt ! extra?
endif
enddo
!!!flux_dt = flmlt(ibv)*fr_snow(ibv)*dts
flux_dt = flux_snow_tot(nlsn,ibv)*dts
if ( abs(flux_dt-flmlt(ibv)*fr_snow(ibv)*dts) > 1.d-12 )
& print *,"GHYTR: tracers flux_snow_tot problem",
& flux_dt-flmlt(ibv)*fr_snow(ibv)*dts
tr_w(:m,1,ibv) = tr_w(:m,1,ibv)
& + tr_wsnc(:m,nlsn,ibv)*flux_dt
wi(1,ibv) = wi(1,ibv) + flux_dt
tr_wsn(:m,nlsn,ibv) = tr_wsn(:m,nlsn,ibv)
& - tr_wsnc(:m,nlsn,ibv)*flux_dt
else ! no snow
tr_w(:m,1,ibv) = tr_w(:m,1,ibv)
& + sum( tr_wsn(:m,1:nlsn,ibv), 2 )
! & + tr_wcc(:m,ibv)*flmlt_scale(ibv)*dts
& + tr_wcc(:m,ibv)*drips(ibv)*dts
tr_wsn(:m,1:nlsn,ibv) = 0.d0
wi(1,ibv) = wi(1,ibv) + flmlt_scale(ibv)*dts
endif
! evap
!gg if (tr_wd_TYPE(ntixw(mc)).eq.nWATER) THEN
flux_dtm(:m) = 0.d0
if ( ibv == 1 ) then
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& flux_dtm(mc) = max( evapbs*fr_snow(1)*dts, 0.d0 )
else
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& flux_dtm(mc) = max( evapvs*fm*fr_snow(2)*dts, 0.d0 )
endif
!!! test
! tr_wsnc(:m,1,ibv) = 1000.d0
tr_wsn(:m,1,ibv) = tr_wsn(:m,1,ibv)
& - tr_wsnc(:m,1,ibv)*flux_dtm(:m)
tr_evap(:m,ibv)= tr_evap(:m,ibv)
& + tr_wsnc(:m,1,ibv)*flux_dtm(:m)/dts
!gg end if
enddo
!!!! for test
! tr_wsn(1,:,1) = wsn(:,1) * fr_snow(1) * 1000.d0
! tr_wsn(1,:,2) = wsn(:,2) * fr_snow(2) * 1000.d0
cddd print '(a,10(e12.4))', 'tr_wsn_1 '
cddd & , tr_wsn(1,:,1) - wsn(:,1) * fr_snow(1) * 1000.d0
cddd & , fr_snow(1), nsn(1)+0.d0 ,flux_snow_tot(0:nlsn,1)*dts
cddd print '(a,10(e12.4))', 'tr_wsn_2 '
cddd & , tr_wsn(1,:,2) - wsn(:,2) * fr_snow(2) * 1000.d0
cddd & , fr_snow(2), nsn(2)+0.d0, flux_snow_tot(0:nlsn,2)*dts
ccc soil layers
!>> tr_wc(:m,1:n,1:2) = 0.d0
!> tr_wc(:m,1:n,1:2) = 1000.d0 !!! was 0
! add dew to bare soil
!gg if (tr_wd_TYPE(ntixw(mc)).eq.nWATER) THEN
if ( process_bare .and. evapb < 0.d0 ) then
flux_dt = - evapb*(1.d0-fr_snow(1))*dts
flux_dtm(:m) = 0.d0
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& flux_dtm(mc) = flux_dt
tr_evap(:m,1) = tr_evap(:m,1) - flux_dtm(:m)/dts*tr_surf(:m)
tr_w(:m,1,1) = tr_w(:m,1,1) + flux_dtm(:m)*tr_surf(:m)
wi(1,1) = wi(1,1) + flux_dt
endif
!gg end if
do ibv=i_bare,i_vege
! initial tr_wc
do i=1,n
if ( wi(i,ibv) > 0.d0 ) then
tr_wc(:m,i,ibv) = tr_w(:m,i,ibv)/wi(i,ibv)
else ! actually 'else' should never happen
tr_wc(:m,i,ibv) = 0.d0
endif
enddo
! precip
flux_dt = dripw(ibv)*(1.d0-fr_snow(ibv))*dts
tr_w(:m,1,ibv) = tr_w(:m,1,ibv) + tr_wcc(:m,ibv)*flux_dt
wi(1,ibv) = wi(1,ibv) + flux_dt
if ( wi(1,ibv) > 0.d0 )
& tr_wc(:m,1,ibv) = tr_w(:m,1,ibv)/wi(1,ibv)
call check_wc(tr_wc(1,1,ibv))
! sweep down
do i=2,n
if ( f(i,ibv) < 0.d0 ) then
flux_dt = - f(i,ibv)*dts
tr_w(:m,i,ibv) = tr_w(:m,i,ibv) + tr_wc(:m,i-1,ibv)*flux_dt
wi(i,ibv) = wi(i,ibv) + flux_dt
if ( wi(i,ibv) > 0.d0 )
& tr_wc(:m,i,ibv) = tr_w(:m,i,ibv)/wi(i,ibv)
call check_wc(tr_wc(1,i,ibv))
tr_w(:m,i-1,ibv) = tr_w(:m,i-1,ibv)
& - tr_wc(:m,i-1,ibv)*flux_dt
endif
enddo
! sweep up
do i=n-1,1,-1
if ( f(i+1,ibv) > 0.d0 ) then
flux_dt = f(i+1,ibv)*dts
tr_w(:m,i,ibv) = tr_w(:m,i,ibv) + tr_wc(:m,i+1,ibv)*flux_dt
wi(i,ibv) = wi(i,ibv) + flux_dt
if ( wi(i,ibv) > 0.d0 )
& tr_wc(:m,i,ibv) = tr_w(:m,i,ibv)/wi(i,ibv)
call check_wc(tr_wc(1,i,ibv))
tr_w(:m,i+1,ibv) = tr_w(:m,i+1,ibv)
& - tr_wc(:m,i+1,ibv)*flux_dt
endif
enddo
enddo
ccc evap from bare soil
if ( process_bare .and. evapb >= 0.d0 ) then
evap_tmp(:m) = 0.d0
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& evap_tmp(mc) = evapb*(1.d0-fr_snow(1))
!gg if (tr_wd_TYPE(ntixw(mc)).eq.nWATER) THEN
#ifdef TRACERS_SPECIAL_O18
if ( evap_tmp(1)*dts < wi(1,1) .and. tp(1,1) > 0.d0 ) then
do mc=1,m
ccc loops...
evap_tmp(mc) = evap_tmp(mc) * fracvl(tp(1,1),ntixw(mc))
enddo
endif
#endif
tr_w(:m,1,1) = tr_w(:m,1,1) - evap_tmp(:m)*tr_wc(:m,1,1)*dts
tr_evap(:m,1) = tr_evap(:m,1) + evap_tmp(:m)*tr_wc(:m,1,1)
!gg endif
endif
ccc runoff
do ibv=i_bare,i_vege
tr_rnff(:m,ibv) = tr_rnff(:m,ibv) + tr_wc(:m,1,ibv)*rnf(ibv)
tr_w(:m,1,ibv) = tr_w(:m,1,ibv) - tr_wc(:m,1,ibv)*rnf(ibv)*dts
do i=1,n
tr_rnff(:m,ibv) = tr_rnff(:m,ibv)
& + tr_wc(:m,i,ibv)*rnff(i,ibv)
tr_w(:m,i,ibv) = tr_w(:m,i,ibv)
& - tr_wc(:m,i,ibv)*rnff(i,ibv)*dts
enddo
enddo
ccc !!! include surface runoff !!!
ccc transpiration
if ( process_vege ) then
do i=1,n
!gg if (tr_wd_TYPE(ntixw(n)).eq.nWATER) THEN
flux_dtm(:m) = 0.d0
forall(mc=1:m, tr_wd_TYPE(ntixw(mc)) == nWATER)
& flux_dtm(mc) = fd*(1.-fr_snow(2)*fm)*evapdl(i,2)*dts
tr_evap(:m,2) = tr_evap(:m,2) + tr_wc(:m,i,2)*flux_dtm(:m)/dts
tr_w(:m,i,2) = tr_w(:m,i,2) - tr_wc(:m,i,2)*flux_dtm(:m)
!gg end if
enddo
endif
cddd print *, 'end ghy_tracers'
cddd print *, 'trpr, trpr_dt ', trpr(1), trpr(1)*dts
cddd print *, 'tr_w 1 ', tr_w(1,:,1)
cddd print *, 'tr_w 2 ', tr_w(1,:,2)
cddd print *, 'tr_wsn 1 ', tr_wsn(1,:,1)
cddd print *, 'tr_wsn 2 ', tr_wsn(1,:,2)
cddd print *, 'evap ', tr_evap(1,:)*dts
cddd print *, 'runoff ', tr_rnff(1,:)*dts
do mc=1,m
do ibv=i_bare,i_vege
do i=1,n
if ( tr_w(mc,i,ibv) < 0.d0 ) then
if ( tr_w(mc,i,ibv) < -tol(mc) ) then
print *,'GHY:tr_w<0 ', tr_w(mc,i,ibv), ibv, i
call stop_model("GHY:tr_w<0",255)
endif
tr_w(mc,i,ibv) = 0.d0
endif
enddo
do i=1,nlsn
if ( tr_wsn(mc,i,ibv) < 0.d0 ) then
if ( tr_wsn(mc,i,ibv) < -tol(mc) ) then
print *,'GHY:tr_wsn<0 ', tr_wsn(mc,i,ibv),
& ibv, i, flmlt(ibv), fr_snow(ibv)
endif
tr_wsn(mc,i,ibv) = 0.d0
endif
enddo
err = trpr(mc)*dts + trdd(mc)*dts
& + sum( tr_w_o(mc,:,ibv) ) + sum( tr_wsn_o(mc,:,ibv) )
& - sum( tr_w(mc,:,ibv) ) - sum( tr_wsn(mc,:,ibv) )
& - tr_evap(mc,ibv)*dts - tr_rnff(mc,ibv)*dts
!print *,'err ', err
if ( abs( err ) > tol(mc) ) then
write(0,*)
$ 'ghy tracers not conserved at',ijdebug,' err=',err
write(99,*)
$ 'ghy tracers not conserved at',ijdebug,' err=',err
$ ,"mc= ",mc, "ibv= ",ibv
$ ,"value =",trpr(mc)*dts
& + sum(tr_w_o(mc,:,ibv) ) + sum( tr_wsn_o(mc,:,ibv))
write(99,*) "tr_w_o",tr_w_o(mc,:,ibv)
$ ,"tr_wsn_o",tr_wsn_o(mc,:,ibv)
$ ,"tr_w",tr_w(mc,:,ibv)
$ ,"tr_wsn",tr_wsn(mc,:,ibv)
$ ,"trpr(mc)*dts",trpr(mc)*dts
$ ,"trdd(mc)*dts",trdd(mc)*dts
$ ,"tr_evap*tds",tr_evap(mc,ibv)*dts
$ ,"tr_rnff",tr_rnff(mc,ibv)*dts
!!! call stop_model("ghy tracers not conserved",255)
endif
enddo
enddo
!!! hack
!!! get rid of possible garbage
do ibv=i_bare,i_vege
do i=1,nlsn
if ( wsni(i,ibv) < 1.d-14 ) then
wsni(i,ibv) = 0.d0
tr_wsn(:m,i,ibv) = 0.d0
endif
enddo
enddo
end subroutine ghy_tracers
#endif
subroutine check_water( flag ) 4
integer flag
real*8 total_water(2), error_water
real*8 old_total_water(2), old_fr_snow(2) ! save
COMMON /check_water_tp/ old_total_water, old_fr_snow
!$OMP THREADPRIVATE (/check_water_tp/)
integer k, ibv
total_water(1) = 0.
total_water(2) = w(0,2)
do ibv=1,2
do k=1,nsn(ibv)
total_water(ibv) = total_water(ibv) + wsn(k,ibv)*fr_snow(ibv)
enddo
do k=1,n
total_water(ibv) = total_water(ibv) + w(k,ibv)
enddo
enddo ! ibv
if ( flag == 0 ) then
old_total_water(:) = total_water(:)
old_fr_snow(:) = fr_snow(:)
return
endif
! bare soil
if ( process_bare ) then
error_water = (total_water(1) - old_total_water(1)) / dts
$ - pr + evap_tot(1) + sum(rnff(1:n,1)) + rnf(1)
if (abs(error_water)>1.d-15)
$ write(99,*)'bare',ijdebug,error_water
if ( abs( error_water ) > 1.d-13 ) then
c & call stop_model('GHY: water conservation problem',255)
write(0,*)'evap_tot(1)',ijdebug,evap_tot(1)
endif
endif
! vegetated soil
if ( process_vege ) then
error_water = (total_water(2) - old_total_water(2)) / dts
& - pr + evap_tot(2) + sum(rnff(1:n,2)) + rnf(2)
if (abs(error_water)>1.d-15)
$ write(99,*)'vege',ijdebug,error_water
if ( abs( error_water ) > 1.d-10 ) then
write(0,*)'evap_tot(2)',ijdebug,evap_tot(2)
c call stop_model(
c & 'GHY: water conservation problem in veg. soil',255)
endif
endif
ghy_debug%water(:) = total_water(:)
end subroutine check_water
subroutine check_energy( flag ) 4
integer flag
real*8 total_energy(2), error_energy
real*8 old_total_energy(2), old_fr_snow(2) !save
COMMON /check_energy_tp/ old_total_energy, old_fr_snow
!$OMP THREADPRIVATE (/check_energy_tp/)
integer k, ibv
total_energy(1) = 0.
total_energy(2) = ht(0,2)
do ibv=1,2
do k=1,nsn(ibv)
total_energy(ibv) = total_energy(ibv)+hsn(k,ibv)*fr_snow(ibv)
enddo
do k=1,n
total_energy(ibv) = total_energy(ibv) + ht(k,ibv)
enddo
enddo ! ibv
if ( flag == 0 ) then
old_total_energy(:) = total_energy(:)
old_fr_snow(:) = fr_snow(:)
return
endif
! bare soil
if ( process_bare ) then
error_energy = (total_energy(1) - old_total_energy(1)) / dts
$ - htpr + elh*evap_tot(1)
& + shw*( sum( rnff(1:n,1)*max(tp(1:n,1),0.d0) )
& + rnf(1)*max(tp(1,1),0.d0) )
& - srht - trht + thrm_tot(1) + snsh_tot(1)
if ( abs( error_energy ) > 1.d-5 ) then
write(0,*)'GHY:bare soil error_energy',ijdebug,error_energy
!call stop_model('GHY: energy conservation problem',255)
endif
endif
! vegetated soil
!!! if ( fr_snow(2) > 0.d0 .or. old_fr_snow(2) > 0.d0 ) return
if ( process_vege ) then
error_energy = (total_energy(2) - old_total_energy(2)) / dts
$ - htpr + elh*evap_tot(2)
& + shw*( sum( rnff(1:n,2)*max(tp(1:n,2),0.d0) )
& + rnf(2)*max(tp(1,2),0.d0) )
& - srht - trht + thrm_tot(2) + snsh_tot(2)
if ( abs( error_energy ) > 1.d-5) then
write(0,*)'GHY:veg soil error_energy',ijdebug,error_energy
!call stop_model('GHY: energy cons problem in veg. soil',255)
endif
endif
ghy_debug%energy(:) = total_energy(:)
end subroutine check_energy
subroutine restrict_snow (wsn_max) 1,2
!@sum remove the snow in excess of WSN_MAX and dump its water into
!@+ runoff, keeping associated heat in the soil
real*8 :: wsn_max,wsn_tot,d_wsn,eta,dw(2:3),dh(2:3)
integer :: ibv
do ibv=i_bare,i_vege
if ( fr_snow(ibv) <= 0.d0 ) cycle
wsn_tot = sum( wsn(1:nsn(ibv),ibv) )
if ( wsn_tot <= WSN_MAX ) cycle
! check if snow structure ok for thick snow
if ( nsn(ibv) < 3)
& call stop_model("remove_extra_snow: nsn<3",255)
d_wsn = wsn_tot - WSN_MAX
! do not remove too much at a time (for now 24mm / day)
d_wsn = min( d_wsn, 1.d-3*dts/3600.d0 )
!print *,"restricting snow: wsn_tot = ", wsn_tot,ijdebug
!print *,"before:",hsn(1:3, ibv),fr_snow(ibv) ,ht(1,ibv)
! fraction of snow to be removed:
eta = d_wsn/sum(wsn(2:3,ibv))
dw(2:3) = eta*wsn(2:3, ibv)
dh(2:3) = eta*hsn(2:3, ibv)
!print *,"dw,dh:",dw,dh
wsn(2:3, ibv) = wsn(2:3, ibv) - dw(2:3)
hsn(2:3, ibv) = hsn(2:3, ibv) - dh(2:3)
dzsn(2:3, ibv) = (1.d0-eta)*dzsn(2:3, ibv)
rnf(ibv) = rnf(ibv) + sum(dw(2:3))*fr_snow(ibv)/dts
ht(1,ibv) = ht(1,ibv) + sum(dh(2:3))*fr_snow(ibv)
& - sum(dw(2:3))*fr_snow(ibv)*max(tp(1,ibv),0.d0)*shw
!print *,"after:",hsn(1:3, ibv),fr_snow(ibv) ,ht(1,ibv)
#ifdef TRACERS_WATER
tr_rnff(1:ntg,ibv) = tr_rnff(1:ntg,ibv) +
& eta*(tr_wsn(1:ntg,2,ibv)+tr_wsn(1:ntg,3,ibv))
tr_wsn(1:ntg,2:3,ibv) = (1.d0-eta)*tr_wsn(1:ntg,2:3,ibv)
#endif
end do
end subroutine restrict_snow
end module sle001
ccccccccccc notes ccccccccccccccccc
ccc just in case, the thickness of layers is:
c n, dz = 1, 9.99999642372131348E-2
c n, dz = 2, 0.17254400253295898
c n, dz = 3, 0.2977144718170166
c n, dz = 4, 0.51368874311447144
c n, dz = 5, 0.88633960485458374
c n, dz = 6, 1.5293264389038086