SUBROUTINE SFC_LAND(IM,KM,PS,U1,V1,T1,Q1,
Clu_Rev6: add TPRCP and SRFLAG
!** & SHELEG,TSKIN,QSURF,
& SHELEG,TSKIN,QSURF,TPRCP,SRFLAG,
& SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,
& DLWFLX,SLRAD,SNOWMT,DELT,Z0RL,TG3,
& GFLUX,ZSOIL,
& CM, CH, RHSCNPY,RHSMC,AIM,BIM,CIM,
& RCL,PRSL1,PRSLKI,SLIMSK,
& DRAIN,EVAP,HFLX,EP,DDVEL,
+ CMM,CHH,Z1,EVBS,EVCW,TRANS,SBSNO,
+ SNOWC,STM,SNOHF,
& tsurf,flag_iter, flag_guess)
!
USE MACHINE , ONLY : kind_phys
USE FUNCPHYS, ONLY : fpvs
USE PHYSCONS, grav => con_g, SBC => con_sbc, HVAP => con_HVAP
&, CP => con_CP, HFUS => con_HFUS, JCAL => con_JCAL
&, EPS => con_eps, EPSM1 => con_epsm1, t0c => con_t0c
&, RVRDM1 => con_FVirt, RD => con_RD
implicit none
!
! include 'constant.h'
!
integer IM, km
!
real(kind=kind_phys), parameter :: cpinv=1.0/cp, HVAPI=1.0/HVAP
real(kind=kind_phys) DELT
INTEGER SOILTYP(IM), VEGTYPE(IM)
real(kind=kind_phys) PS(IM), U1(IM), V1(IM),
& T1(IM), Q1(IM), SHELEG(IM),
& TSKIN(IM), QSURF(IM), SMC(IM,KM),
& STC(IM,KM), DM(IM), SIGMAF(IM),
& CANOPY(IM), DLWFLX(IM), SLRAD(IM),
& SNOWMT(IM), Z0RL(IM), TG3(IM),
& GFLUX(IM),
& ZSOIL(IM,KM), CM(IM), CH(IM),
& RHSCNPY(IM), RHSMC(IM,KM),
& AIM(IM,KM), BIM(IM,KM), CIM(IM,KM),
& RCL(IM), PRSL1(IM), PRSLKI(IM),
& SLIMSK(IM), DRAIN(IM), EVAP(IM),
& HFLX(IM), RNET(IM), EP(IM),
& WIND(IM), DDVEL(IM)
Clu_Rev6: add TPRCP and SRFLAG
&, TPRCP(IM), SRFLAG(IM)
Cwei added 10/24/2006
&, CHH(IM),CMM(IM),EVBS(IM),EVCW(IM)
&, TRANS(IM),SBSNO(IM),SNOWC(IM),STM(IM),SNOHF(IM)
!
! land-related prognostic fields
!
real(kind=kind_phys) SHELEG_OLD(IM),
+ TPRCP_OLD(IM), SRFLAG_OLD(IM),
+ TSKIN_OLD(IM), CANOPY_OLD(IM),
+ STC_OLD(IM,KM)
logical flag_iter(im), flag_guess(im)
!
! Locals
!
integer k,i
!
real(kind=kind_phys) CANFAC(IM),
& DDZ(IM), DDZ2(IM), DELTA(IM),
& DEW(IM), DF1(IM), DFT0(IM),
& DFT2(IM), DFT1(IM),
& DMDZ(IM), DMDZ2(IM), DTDZ1(IM),
& DTDZ2(IM), DTV(IM), EC(IM),
& EDIR(IM), ETPFAC(IM),
& FACTSNW(IM), FH2(IM),
& FX(IM), GX(IM),
& HCPCT(IM), HL1(IM), HL12(IM),
& HLINF(IM), PARTLND(IM), PH(IM),
& PH2(IM), PM(IM), PM10(IM),
& PSURF(IM), Q0(IM), QS1(IM),
& QSS(IM), RAT(IM), RCAP(IM),
& RCH(IM), RHO(IM), RS(IM),
& RSMALL(IM), SLWD(IM), SMCZ(IM),
& SNET(IM), SNOEVP(IM), SNOWD(IM),
& T1O(IM), T2MO(IM), TERM1(IM),
& TERM2(IM), THETA1(IM), THV1(IM),
& TREF(IM), TSURF(IM), TV1(IM),
& TVS(IM), TSEA(IM), TWILT(IM),
& XX(IM), XRCL(IM), YY(IM),
& Z0(IM), Z0MAX(IM), Z1(IM),
& ZTMAX(IM), ZZ(IM), PS1(IM)
!
real(kind=kind_phys) a0, a0p, a1, a1p, aa, aa0,
& aa1, adtv, alpha, arnu, b1, b1p,
& b2, b2p, bb, bb0, bb1, bb2,
& bfact, ca, cc, cc1, cc2, cfactr,
& ch2o, charnock, cice, convrad, cq, csoil,
& ctfil1,ctfil2, delt2, df2, dfsnow,
& elocp, eth, ff, FMS,
& fhs, funcdf, funckt,g, hl0, hl0inf,
& hl110, hlt, hltinf,OLINF, rcq, rcs,
& rct, restar, rhoh2o,rnu, RSI,
& rss, scanop, sig2k, sigma, smcdry,
& t12, t14, tflx, tgice, topt,
& val, vis, zbot, snomin, tem
!
cc
PARAMETER (CHARNOCK=.014,CA=.4)!C CA IS THE VON KARMAN CONSTANT
PARAMETER (G=grav,sigma=sbc)
PARAMETER (ALPHA=5.,A0=-3.975,A1=12.32,B1=-7.755,B2=6.041)
PARAMETER (A0P=-7.941,A1P=24.75,B1P=-8.705,B2P=7.899,VIS=1.4E-5)
PARAMETER (AA1=-1.076,BB1=.7045,CC1=-.05808)
PARAMETER (BB2=-.1954,CC2=.009999)
PARAMETER (ELOCP=HVAP/CP,DFSNOW=.31,CH2O=4.2E6,CSOIL=1.26E6)
PARAMETER (SCANOP=.5,CFACTR=.5,ZBOT=-3.,TGICE=271.2)
PARAMETER (CICE=1880.*917.,topt=298.)
PARAMETER (RHOH2O=1000.,CONVRAD=JCAL*1.E4/60.)
PARAMETER (CTFIL1=.5,CTFIL2=1.-CTFIL1)
PARAMETER (RNU=1.51E-5,ARNU=.135*RNU)
parameter (snomin=1.0e-9)
!
LOGICAL FLAG(IM), FLAGSNW(IM)
real(kind=kind_phys) KT1(IM), KT2(IM), KTSOIL,
& ET(IM,KM),
& STSOIL(IM,KM), AI(IM,KM), BI(IM,KM),
& CI(IM,KM), RHSTC(IM,KM)
real(kind=kind_phys) rsmax(13), rgl(13), rsmin(13), hs(13),
& smmax(9), smdry(9), smref(9), smwlt(9)
c
c the 13 vegetation types are:
c
c 1 ... broadleave-evergreen trees (tropical forest)
c 2 ... broadleave-deciduous trees
c 3 ... broadleave and needle leave trees (mixed forest)
c 4 ... needleleave-evergreen trees
c 5 ... needleleave-deciduous trees (larch)
c 6 ... broadleave trees with groundcover (savanna)
c 7 ... groundcover only (perenial)
c 8 ... broadleave shrubs with perenial groundcover
c 9 ... broadleave shrubs with bare soil
c 10 ... dwarf trees and shrubs with ground cover (trunda)
c 11 ... bare soil
c 12 ... cultivations (use parameters from type 7)
c 13 ... glacial
c
data rsmax/13*5000./
data rsmin/150.,100.,125.,150.,100.,70.,40.,
& 300.,400.,150.,999.,40.,999./
data rgl/5*30.,65.,4*100.,999.,100.,999./
data hs/41.69,54.53,51.93,47.35,47.35,54.53,36.35,
& 3*42.00,999.,36.35,999./
data smmax/.421,.464,.468,.434,.406,.465,.404,.439,.421/
data smdry/.07,.14,.22,.08,.18,.16,.12,.10,.07/
data smref/.283,.387,.412,.312,.338,.382,.315,.329,.283/
data smwlt/.029,.119,.139,.047,.010,.103,.069,.066,.029/
!
save rsmax, rsmin, rgl, hs, smmax, smdry, smref, smwlt
!
DELT2 = DELT * 2.
C
C ESTIMATE SIGMA ** K AT 2 M
C
SIG2K = 1. - 4. * G * 2. / (CP * 280.)
C
C FLAG for land
C
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ. 1.
ENDDO
C
C save land-related prognostic fields for guess run
C
do i=1, im
if(FLAG(I) .AND. flag_guess(i)) then
sheleg_old(i) = sheleg(i)
tskin_old(i) = tskin(i)
canopy_old(i) = canopy(i)
tprcp_old(i) = tprcp(i)
srflag_old(i) = srflag(i)
do k=1, km
stc_old(i,k) = stc(i,k)
enddo
endif
enddo
C
CWei_FIX_3: you need to remove snow-rain detection here
C For OSU LSM, the snow-rain detection is done inside gbphys routine !!
C
Clu_Rev6: snow-rain detection
C
C DO I=1,IM
C IF(FLAG(I).AND. flag_guess(i)) THEN
C IF(SRFLAG(I) .EQ. 1.) THEN
C SHELEG(i) = SHELEG(i) + 1.E3*TPRCP(i)
C TPRCP(i) = 0.
C ENDIF
C ENDIF
C ENDDO
C
C INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIE
C PSURF IS IN PASCALS
C WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
C RHO IS DENSITY, QS1 IS SAT. HUM. AT LEVEL1 AND QSS IS SAT. HUM. AT
C SURFACE
C CONVERT SLRAD TO THE CIVILIZED UNIT FROM LANGLEY MINUTE-1 K-4
C SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM
C
!!
! qs1 = fpvs(t1)
! qss = fpvs(tskin)
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
XRCL(I) = SQRT(RCL(I))
PSURF(I) = 1000. * PS(I)
PS1(I) = 1000. * PRSL1(I)
! SLWD(I) = SLRAD(I) * CONVRAD
SLWD(I) = SLRAD(I)
c
c DLWFLX has been given a negative sign for downward longwave
c snet is the net shortwave flux
c
SNET(I) = -SLWD(I) - DLWFLX(I)
WIND(I) = XRCL(I) * SQRT(U1(I) * U1(I) + V1(I) * V1(I))
& + MAX(0.0, MIN(DDVEL(I), 30.0))
WIND(I) = MAX(WIND(I),1.)
Q0(I) = MAX(Q1(I),1.E-8)
! TSURF(I) = TSKIN(I)
TSEA(I) = tsurf(I) !!! Cwei_q2m_iter
THETA1(I) = T1(I) * PRSLKI(I)
TV1(I) = T1(I) * (1. + RVRDM1 * Q0(I))
THV1(I) = THETA1(I) * (1. + RVRDM1 * Q0(I))
TVS(I) = TSEA(I) * (1. + RVRDM1 * Q0(I))
RHO(I) = PS1(I) / (RD * TV1(I))
cjfe QS1(I) = 1000. * FPVS(T1(I))
qs1(i) = fpvs(t1(i))
QS1(I) = EPS * QS1(I) / (PS1(I) + EPSM1 * QS1(I))
QS1(I) = MAX(QS1(I), 1.E-8)
Q0(I) = min(QS1(I),Q0(I))
cjfe QSS(I) = 1000. * FPVS(TSEA(I))
qss(i) = fpvs(tskin(i)) !!! change to tsurf?
QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
c RS = PLANTR
RS(I) = 0.
if(VEGTYPE(I).gt.0.) RS(I) = rsmin(VEGTYPE(I))
Z0(I) = .01 * Z0RL(i)
CANOPY(I)= MAX(CANOPY(I),0.)
DM(I) = 1.
FACTSNW(I) = 10.
Clu IF(SLIMSK(I).EQ.2.) FACTSNW(I) = 3.
C
C SNOW DEPTH IN WATER EQUIVALENT IS CONVERTED FROM MM TO M UNIT
C
SNOWD(I) = SHELEG(I) / 1000.
FLAGSNW(I) = .FALSE.
C
C WHEN SNOW DEPTH IS LESS THAN 1 MM, A PATCHY SNOW IS ASSUMED AND
C SOIL IS ALLOWED TO INTERACT WITH THE ATMOSPHERE.
C WE SHOULD EVENTUALLY MOVE TO A LINEAR COMBINATION OF SOIL AND
C SNOW UNDER THE CONDITION OF PATCHY SNOW.
C
IF(SNOWD(I).GT..001.OR.SLIMSK(I).EQ.2.) RS(I) = 0.
IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE.
C##DG IF(LAT.EQ.LATD) THEN
C##DG PRINT *, ' WIND,TV1,TVS,Q1,QS1,SNOW,SLIMSK=',
C##DG& WIND,TV1,TVS,Q1,QS1,SNOWD,SLIMSK
C##DG PRINT *, ' SNET, SLWD =', SNET, SLWD(I)
C##DG ENDIF
ENDIF
ENDDO
!!
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
ZSOIL(I,1) = -.10
DO K = 2, KM
ZSOIL(I,K) = ZSOIL(I,K-1)
& + (-2. - ZSOIL(I,1)) / (KM - 1)
ENDDO
CWei [+5] use the same soil layer structure as Noah if running with 4-layer
if(km.gt.2)then
ZSOIL(I,2) = -.40
ZSOIL(I,3) = -1.0
ZSOIL(I,4) = -2.0
endif
ENDIF
ENDDO
!!
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
Z1(I) = -RD * TV1(I) * LOG(PS1(I)/PSURF(I)) / G
DRAIN(I) = 0.
ENDIF
ENDDO
!!
DO K = 1, KM
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
ET(I,K) = 0.
RHSMC(I,K) = 0.
AIM(I,K) = 0.
BIM(I,K) = 1.
CIM(I,K) = 0.
STSOIL(I,K) = STC(I,K)
ENDIF
ENDDO
ENDDO
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
EDIR(I) = 0.
EC(I) = 0.
EVAP(I) = 0.
HFLX(I) = 0.
EP(I) = 0.
SNOWMT(I) = 0.
GFLUX(I) = 0.
RHSCNPY(I) = 0.
FX(I) = 0.
ETPFAC(I) = 0.
CANFAC(I) = 0.
Cwei added 10/24/2006
EVBS(I)=0
EVCW(I)=0
TRANS(I)=0
SBSNO(I)=0
SNOWC(I)=0
SNOHF(I)=0
ENDIF
ENDDO
C
C RCP = RHO CP CH V
C
DO I = 1, IM
IF(flag_iter(i).and. FLAG(I)) THEN
RCH(I) = RHO(I) * CP * CH(I) * WIND(I)
Cwei added 10/24/2006
CMM(I)=CM(I)* WIND(I)
CHH(I)=RHO(I)*CH(I)* WIND(I)
ENDIF
ENDDO
C
C COMPUTE SOIL/SNOW/ICE HEAT FLUX IN PREPARATION FOR SURFACE ENERGY
C BALANCE CALCULATION
C
DO I = 1, IM
Clu GFLUX(I) = 0.
IF(flag_iter(i).and. FLAG(I)) THEN
SMCZ(I) = .5 * (SMC(I,1) + .20)
DFT0(I) = KTSOIL(SMCZ(I),SOILTYP(I))
Clu ELSEIF(SLIMSK(I).EQ.2.) THEN
C DF FOR ICE IS TAKEN FROM MAYKUT AND UNTERSTEINER
C DF IS IN SI UNIT OF W K-1 M-1
Clu DFT0(I) = 2.2
ENDIF
ENDDO
!!
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
C IF(SNOWD(I).GT..001) THEN
IF(FLAGSNW(I)) THEN
C
C WHEN SNOW COVERED, GROUND HEAT FLUX COMES FROM SNOW
C
TFLX = MIN(T1(I), TSEA(I))
GFLUX(I) = -DFSNOW * (TFLX - STSOIL(I,1))
& / (FACTSNW(I) * MAX(SNOWD(I),.001))
ELSE
GFLUX(I) = DFT0(I) * (STSOIL(I,1) - TSEA(I))
& / (-.5 * ZSOIL(I,1))
ENDIF
GFLUX(I) = MAX(GFLUX(I),-200.)
GFLUX(I) = MIN(GFLUX(I),+200.)
ENDIF
ENDDO
DO I = 1, IM
IF(flag_iter(i).and. FLAG(I)) THEN
PARTLND(I) = 1.
IF(SNOWD(I).GT.0..AND.SNOWD(I).LE..001) THEN
PARTLND(I) = 1. - SNOWD(I) / .001
ENDIF
ENDIF
ENDDO
DO I = 1, IM
IF(flag_iter(i).and. FLAG(I)) THEN
SNOEVP(I) = 0.
if(SNOWD(I).gt..001) PARTLND(I) = 0.
ENDIF
ENDDO
C
C COMPUTE POTENTIAL EVAPORATION FOR LAND AND SEA ICE
C
DO I = 1, IM
IF(flag_iter(i).and. FLAG(I)) THEN
T12 = T1(I) * T1(I)
T14 = T12 * T12
C
C RCAP = FNET - SIGMA T**4 + GFLX - RHO CP CH V (T1-THETA1)
C
RCAP(I) = -SLWD(I) - SIGMA * T14 + GFLUX(I)
& - RCH(I) * (T1(I) - THETA1(I))
C
C RSMALL = 4 SIGMA T**3 / RCH(I) + 1
C
RSMALL(I) = 4. * SIGMA * T1(I) * T12 / RCH(I) + 1.
C
C DELTA = L / CP * DQS/DT
C
DELTA(I) = ELOCP * EPS * HVAP * QS1(I) / (RD * T12)
C
C POTENTIAL EVAPOTRANSPIRATION ( WATTS / M**2 ) AND
C POTENTIAL EVAPORATION
C
TERM1(I) = ELOCP * RSMALL(I) * RCH(I)*(QS1(I)-Q0(I))
TERM2(I) = RCAP(I) * DELTA(I)
EP(I) = (ELOCP * RSMALL(I) * RCH(I) * (QS1(I) - Q0(I))
& + RCAP(I) * DELTA(I))
EP(I) = EP(I) / (RSMALL(I) + DELTA(I))
ENDIF
ENDDO
C
C ACTUAL EVAPORATION OVER LAND IN THREE PARTS : EDIR, ET, AND EC
C
C DIRECT EVAPORATION FROM SOIL, THE UNIT GOES FROM M S-1 TO KG M-2 S-1
C
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ.1..AND.EP(I).GT.0.
ENDDO
DO I = 1, IM
IF(flag_iter(i))THEN
IF(FLAG(I)) THEN
DF1(I) = FUNCDF(SMC(I,1),SOILTYP(I))
KT1(I) = FUNCKT(SMC(I,1),SOILTYP(I))
endif
if(FLAG(I).and.STC(I,1).lt.t0c) then
DF1(I) = 0.
KT1(I) = 0.
endif
IF(FLAG(I)) THEN
c TREF = .75 * THSAT(SOILTYP(I))
TREF(I) = smref(SOILTYP(I))
c TWILT = TWLT(SOILTYP(I))
TWILT(I) = smwlt(SOILTYP(I))
smcdry = smdry(SOILTYP(I))
c FX(I) = -2. * DF1(I) * (SMC(I,1) - .23) / ZSOIL(I,1)
c & - KT1(I)
FX(I) = -2. * DF1(I) * (SMC(I,1) - smcdry) / ZSOIL(I,1)
& - KT1(I)
FX(I) = MIN(FX(I), EP(I)/HVAP)
FX(I) = MAX(FX(I),0.)
C
C SIGMAF IS THE FRACTION OF AREA COVERED BY VEGETATION
C
EDIR(I) = FX(I) * (1. - SIGMAF(I)) * PARTLND(I)
ENDIF
endif
ENDDO
c
c calculate stomatal resistance
c
DO I = 1, IM
if(flag_iter(i).and.FLAG(I)) then
c
c resistance due to PAR. We use net solar flux as proxy at the present time
c
ff = .55 * 2. * SNET(I) / rgl(VEGTYPE(I))
rcs = (ff + RS(I)/rsmax(VEGTYPE(I))) / (1. + ff)
rcs = max(rcs,.0001)
rct = 1.
rcq = 1.
c
c resistance due to thermal effect
c
c rct = 1. - .0016 * (topt - theta1) ** 2
c rct = max(rct,.0001)
c
c resistance due to humidity
c
c rcq = 1. / (1. + hs(VEGTYPE(I)) * (QS1(I) - Q0(I)))
c rcq = max(rcq,.0001)
c
c compute resistance without the effect of soil moisture
c
RS(I) = RS(I) / (rcs * rct * rcq)
endif
ENDDO
C
C TRANSPIRATION FROM ALL LEVELS OF THE SOIL
C
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
CANFAC(I) = (CANOPY(I) / SCANOP) ** CFACTR
ETPFAC(I) = SIGMAF(I)
& * (1. - CANFAC(I)) / HVAP
GX(I) = (SMC(I,1) - TWILT(I)) / (TREF(I) - TWILT(I))
GX(I) = MAX(GX(I),0.)
GX(I) = MIN(GX(I),1.)
c
c resistance due to soil moisture deficit
c
rss = GX(I) * (ZSOIL(I,1) / ZSOIL(I,km))
rss = max(rss,.0001)
RSI = RS(I) / rss
c
c transpiration a la Monteith
c
eth = (TERM1(I) + TERM2(I)) /
& (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
ET(I,1) = ETPFAC(I) * eth
& * PARTLND(I)
ENDIF
ENDDO
!!
DO K = 2, KM
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
GX(I) = (SMC(I,K) - TWILT(I)) / (TREF(I) - TWILT(I))
GX(I) = MAX(GX(I),0.)
GX(I) = MIN(GX(I),1.)
c
c resistance due to soil moisture deficit
c
rss = GX(I) * ((ZSOIL(I,k) - ZSOIL(I,k-1))/ZSOIL(I,km))
rss = max(rss,1.e-6)
RSI = RS(I) / rss
c
c transpiration a la Monteith
c
eth = (TERM1(I) + TERM2(I)) /
& (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
ET(I,K) = eth
& * ETPFAC(I) * PARTLND(I)
ENDIF
ENDDO
ENDDO
!!
400 CONTINUE
C
C CANOPY RE-EVAPORATION
C
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
EC(I) = SIGMAF(I) * CANFAC(I) * EP(I) / HVAP
EC(I) = EC(I) * PARTLND(I)
EC(I) = min(EC(I),CANOPY(I)/delt)
ENDIF
ENDDO
C
C SUM UP TOTAL EVAPORATION
C
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
EVAP(I) = EDIR(I) + EC(I)
ENDIF
ENDDO
!!
DO K = 1, KM
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
EVAP(I) = EVAP(I) + ET(I,K)
ENDIF
ENDDO
ENDDO
!!
C
C RETURN EVAP UNIT FROM KG M-2 S-1 TO WATTS M-2
C
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
EVAP(I) = MIN(EVAP(I)*HVAP,EP(I))
ENDIF
ENDDO
C##DG IF(LAT.EQ.LATD) THEN
C##DG PRINT *, 'FX(I), SIGMAF, EDIR(I), ETPFAC=', FX(I)*HVAP,SIGMAF,
C##DG& EDIR(I)*HVAP,ETPFAC*HVAP
C##DG PRINT *, ' ET =', (ET(K)*HVAP,K=1,KM)
C##DG PRINT *, ' CANFAC(I), EC(I), EVAP', CANFAC(I),EC(I)*HVAP,EVAP
C##DG ENDIF
C
C TREAT DOWNWARD MOISTURE FLUX SITUATION
C (EVAP WAS PRESET TO ZERO SO NO UPDATE NEEDED)
C DEW IS CONVERTED FROM KG M-2 TO M TO CONFORM TO PRECIP UNIT
C
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ.1..AND.EP(I).LE.0.
DEW(I) = 0.
ENDDO
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
DEW(I) = -EP(I) * DELT / (HVAP * RHOH2O)
EVAP(I) = EP(I)
DEW(I) = DEW(I) * PARTLND(I)
EVAP(I) = EVAP(I) * PARTLND(I)
DM(I) = 1.
ENDIF
ENDDO
C
C SNOW COVERED LAND
C
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ.1..AND.SNOWD(I).GT.0.
ENDDO
C
C CHANGE OF SNOW DEPTH DUE TO EVAPORATION OR SUBLIMATION
C
C CONVERT EVAP FROM KG M-2 S-1 TO M S-1 TO DETERMINE THE REDUCTION OF S
C
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
BFACT = SNOWD(I) / (DELT * EP(I) / (HVAP * RHOH2O))
BFACT = MIN(BFACT,1.)
C
C THE EVAPORATION OF SNOW
C
IF(EP(I).LE.0.) BFACT = 1.
IF(SNOWD(I).LE..001) THEN
c EVAP = (SNOWD(I)/.001)*BFACT*EP(I) + EVAP
! SNOEVP(I) = bfact * EP(I) * (1. - PARTLND(I))
! EVAP = EVAP + SNOEVP(I)
SNOEVP(I) = bfact * EP(I)
! EVAP = EVAP + SNOEVP(I) * (1. - PARTLND(I))
EVAP(I)=EVAP(I)+SNOEVP(I)*(1.-PARTLND(I))
ELSE
c EVAP(I) = BFACT * EP(I)
SNOEVP(I) = bfact * EP(I)
EVAP(I) = SNOEVP(I)
ENDIF
TSEA(I) = T1(I) +
& (RCAP(I) - GFLUX(I) - DFSNOW * (T1(I) - STSOIL(I,1))
& /(FACTSNW(I) * MAX(SNOWD(I),.001))
c & + THETA1 - T1
c & - BFACT * EP(I)) / (RSMALL(I) * RCH(I)
& - SNOEVP(I)) / (RSMALL(I) * RCH(I)
& + DFSNOW / (FACTSNW(I)* MAX(SNOWD(I),.001)))
c SNOWD(I) = SNOWD(I) - BFACT * EP(I) * DELT / (RHOH2O * HVAP)
SNOWD(I) = SNOWD(I) - SNOEVP(I) * delt / (rhoh2o * hvap)
SNOWD(I) = MAX(SNOWD(I),0.)
ENDIF
ENDDO
C
C SNOW MELT (M)
C
500 CONTINUE
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ.1.
& .AND.SNOWD(I).GT..0
ENDDO
DO I = 1, IM
IF(flag_iter(i))THEN
IF(FLAG(I).AND.TSEA(I).GT.T0C) THEN
SNOWMT(I) = RCH(I) * RSMALL(I) * DELT
& * (TSEA(I) - T0C) / (RHOH2O * HFUS)
SNOWMT(I) = min(SNOWMT(I),SNOWD(I))
SNOWD(I) = SNOWD(I) - SNOWMT(I)
SNOWD(I) = MAX(SNOWD(I),0.)
TSEA(I) = MAX(T0C,TSEA(I)
& -HFUS*SNOWMT(I)*RHOH2O/(RCH(I)*RSMALL(I)*DELT))
ENDIF
ENDIF
ENDDO
c
c We need to re-evaluate evaporation because of snow melt
c the skin temperature is now bounded to 0 deg C
c
! qss = fpvs(TSEA)
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ. 1.
IF(flag_iter(i).and.FLAG(I))THEN
! IF (SNOWD(I) .GT. 0.0) THEN
IF (SNOWD(I) .GT. snomin) THEN
cjfe QSS(I) = 1000. * FPVS(TSEA(I))
qss(i) = fpvs(TSEA(i))
QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
EVAP(I) = elocp * RCH(I) * (QSS(I) - Q0(I))
ENDIF
ENDIF
ENDDO
C
C PREPARE TENDENCY TERMS FOR THE SOIL MOISTURE FIELD WITHOUT PRECIPITAT
C THE UNIT OF MOISTURE FLUX NEEDS TO BECOME M S-1 FOR SOIL MOISTURE
C HENCE THE FACTOR OF RHOH2O
C
DO I = 1, IM
IF(flag_iter(i))THEN
if(FLAG(I)) then
DF1(I) = FUNCDF(SMCZ(I),SOILTYP(I))
KT1(I) = FUNCKT(SMCZ(I),SOILTYP(I))
endif
if(FLAG(I).and.STC(I,1).lt.t0c) then
DF1(I) = 0.
KT1(I) = 0.
endif
IF(FLAG(I)) THEN
RHSCNPY(I) = -EC(I) + SIGMAF(I) * RHOH2O * DEW(I) / DELT
SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
DMDZ(I) = (SMC(I,1) - SMC(I,2)) / (-.5 * ZSOIL(I,2))
RHSMC(I,1) = (DF1(I) * DMDZ(I) + KT1(I)
& + (EDIR(I) + ET(I,1))) / (ZSOIL(I,1) * RHOH2O)
RHSMC(I,1) = RHSMC(I,1) - (1. - SIGMAF(I)) * DEW(I) /
& ( ZSOIL(I,1) * delt)
DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
C
C AIM, BIM, AND CIM ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
C IMPLICIT UPDATE OF THE SOIL MOISTURE
C
AIM(I,1) = 0.
BIM(I,1) = DF1(I) * DDZ(I) / (-ZSOIL(I,1) * RHOH2O)
CIM(I,1) = -BIM(I,1)
ENDIF
ENDIF
ENDDO
!!
DO K = 2, KM
IF(K.LT.KM) THEN
DO I=1,IM
IF(flag_iter(i))THEN
IF(FLAG(I)) THEN
SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
DF2 = FUNCDF(SMCZ(I),SOILTYP(I))
KT2(I) = FUNCKT(SMCZ(I),SOILTYP(I))
ENDIF
IF(FLAG(I).and.STC(I,k).lt.t0c) THEN
df2 = 0.
KT2(I) = 0.
ENDIF
IF(FLAG(I)) THEN
DMDZ2(I) = (SMC(I,K) - SMC(I,K+1))
& / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
RHSMC(I,K) = (DF2 * DMDZ2(I) + KT2(I)
& - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K))
& / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
CIM(I,K) = -DF2 * DDZ2(I)
& / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
ENDIF
ENDIF
ENDDO
ELSE
DO I = 1, IM
IF(flag_iter(i))THEN
IF(FLAG(I)) THEN
KT2(I) = FUNCKT(SMC(I,K),SOILTYP(I))
ENDIF
if(FLAG(I).and.STC(I,k).lt.t0c) KT2(I) = 0.
IF(FLAG(I)) THEN
RHSMC(I,K) = (KT2(I)
& - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K))
& / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
DRAIN(I) = KT2(I)
CIM(I,K) = 0.
ENDIF
ENDIF
ENDDO
ENDIF
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
AIM(I,K) = -DF1(I) * DDZ(I)
& / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
BIM(I,K) = -(AIM(I,K) + CIM(I,K))
DF1(I) = DF2
KT1(I) = KT2(I)
DMDZ(I) = DMDZ2(I)
DDZ(I) = DDZ2(I)
ENDIF
ENDDO
ENDDO
!!
600 CONTINUE
C
C UPDATE SOIL TEMPERATURE
C
DO I=1,IM
Clu FLAG(I) = SLIMSK(I).NE.0.
FLAG(I) = SLIMSK(I).EQ.1.
ENDDO
C
C SURFACE TEMPERATURE IS PART OF THE UPDATE WHEN SNOW IS ABSENT
C
DO I=1,IM
C IF(FLAG(I).AND.SNOWD(I).LE..001) THEN
IF(flag_iter(i))THEN
IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
YY(I) = T1(I) +
c & (RCAP(I)-GFLUX(I) + THETA1 - T1(I)
& (RCAP(I)-GFLUX(I)
& - EVAP(I)) / (RSMALL(I) * RCH(I))
ZZ(I) = 1. + DFT0(I) / (-.5 * ZSOIL(I,1) * RCH(I) * RSMALL(I))
XX(I) = DFT0(I) * (STSOIL(I,1) - YY(I)) /
& (.5 * ZSOIL(I,1) * ZZ(I))
ENDIF
C IF(FLAG(I).AND.SNOWD(I).GT..001) THEN
IF(FLAG(I).AND.FLAGSNW(I)) THEN
YY(I) = STSOIL(I,1)
C
C HEAT FLUX FROM SNOW IS EXPLICIT IN TIME
C
ZZ(I) = 1.
XX(I) = DFSNOW * (STSOIL(I,1) - TSEA(I))
& / (-FACTSNW(I) * MAX(SNOWD(I),.001))
ENDIF
ENDIF
ENDDO
C
C COMPUTE THE FORCING AND THE IMPLICIT MATRIX ELEMENTS FOR UPDATE
C
C CH2O IS THE HEAT CAPACITY OF WATER AND CSOIL IS THE HEAT CAPACITY OF
C
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
DTDZ1(I) = (STSOIL(I,1) - STSOIL(I,2)) / (-.5 * ZSOIL(I,2))
DFT1(I) = KTSOIL(SMCZ(I),SOILTYP(I))
HCPCT(I) = SMC(I,1) * CH2O + (1. - SMC(I,1)) * CSOIL
DFT2(I) = DFT1(I)
DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
C
C AI, BI, AND CI ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
C IMPLICIT UPDATE OF THE SOIL TEMPERATURE
C
AI(I,1) = 0.
BI(I,1) = DFT1(I) * DDZ(I) / (-ZSOIL(I,1) * HCPCT(I))
CI(I,1) = -BI(I,1)
BI(I,1) = BI(I,1)
& + DFT0(I) / (.5 * ZSOIL(I,1) **2 * HCPCT(I) * ZZ(I))
C SS = DFT0(I) * (STSOIL(I,1) - YY(I))
C & / (.5 * ZSOIL(I,1) * ZZ(I))
C RHSTC(1) = (DFT1(I) * DTDZ1(I) - SS)
RHSTC(I,1) = (DFT1(I) * DTDZ1(I) - XX(I))
& / (ZSOIL(I,1) * HCPCT(I))
ENDIF
ENDDO
!!
DO K = 2, KM
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
HCPCT(I) = SMC(I,K) * CH2O + (1. - SMC(I,K)) * CSOIL
ENDIF
ENDDO
IF(K.LT.KM) THEN
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
DTDZ2(I) = (STSOIL(I,K) - STSOIL(I,K+1))
& / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
DFT2(I) = KTSOIL(SMCZ(I),SOILTYP(I))
DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
CI(I,K) = -DFT2(I) * DDZ2(I)
& / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
ENDIF
ENDDO
ELSE
C
C AT THE BOTTOM, CLIMATOLOGY IS ASSUMED AT 2M DEPTH FOR LAND AND
C FREEZING TEMPERATURE IS ASSUMED FOR SEA ICE AT Z(KM)
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
DTDZ2(I) = (STSOIL(I,K) - TG3(I))
& / (.5 * (ZSOIL(I,K-1) + ZSOIL(I,K)) - ZBOT)
DFT2(I) = KTSOIL(SMC(I,K),SOILTYP(I))
CI(I,K) = 0.
ENDIF
ENDDO
ENDIF
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
RHSTC(I,K) = (DFT2(I) * DTDZ2(I) - DFT1(I) * DTDZ1(I))
& / ((ZSOIL(I,K) - ZSOIL(I,K-1)) * HCPCT(I))
AI(I,K) = -DFT1(I) * DDZ(I)
& / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
BI(I,K) = -(AI(I,K) + CI(I,K))
DFT1(I) = DFT2(I)
DTDZ1(I) = DTDZ2(I)
DDZ(I) = DDZ2(I)
ENDIF
ENDDO
ENDDO
!!
700 CONTINUE
C
C SOLVE THE TRI-DIAGONAL MATRIX
C
DO K = 1, KM
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
RHSTC(I,K) = RHSTC(I,K) * DELT2
AI(I,K) = AI(I,K) * DELT2
BI(I,K) = 1. + BI(I,K) * DELT2
CI(I,K) = CI(I,K) * DELT2
ENDIF
ENDDO
ENDDO
C FORWARD ELIMINATION
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
CI(I,1) = -CI(I,1) / BI(I,1)
RHSTC(I,1) = RHSTC(I,1) / BI(I,1)
ENDIF
ENDDO
!!
DO K = 2, KM
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
CC = 1. / (BI(I,K) + AI(I,K) * CI(I,K-1))
CI(I,K) = -CI(I,K) * CC
RHSTC(I,K) = (RHSTC(I,K) - AI(I,K) * RHSTC(I,K-1)) * CC
ENDIF
ENDDO
ENDDO
!!
C BACKWARD SUBSTITUTTION
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
CI(I,KM) = RHSTC(I,KM)
ENDIF
ENDDO
!!
DO K = KM-1, 1
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
CI(I,K) = CI(I,K) * CI(I,K+1) + RHSTC(I,K)
ENDIF
ENDDO
ENDDO
C
C UPDATE SOIL AND ICE TEMPERATURE
C
DO K = 1, KM
DO I=1,IM
IF(flag_iter(i).and.FLAG(I)) THEN
STSOIL(I,K) = STSOIL(I,K) + CI(I,K)
ENDIF
ENDDO
ENDDO
C
C UPDATE SURFACE TEMPERATURE FOR SNOW FREE SURFACES
C
DO I=1,IM
IF(flag_iter(i))THEN
IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
TSEA(I) = (YY(I) + (ZZ(I) - 1.) * STSOIL(I,1)) / ZZ(I)
ENDIF
ENDIF
ENDDO
!!
C
C TIME FILTER FOR SOIL AND SKIN TEMPERATURE
C
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
TSURF(I) = CTFIL1 * TSEA(I) + CTFIL2 * TSURF(I)
ENDIF
ENDDO
DO K = 1, KM
DO I=1,IM
IF(flag_iter(i).and. FLAG(I)) THEN
STC(I,K) = CTFIL1 * STSOIL(I,K) + CTFIL2 * STC(I,K)
ENDIF
ENDDO
ENDDO
C
C GFLUX CALCULATION
C
c DO I=1,IM
c FLAG(I) = SLIMSK(I).EQ.1.
c & .AND.FLAGSNW(I)
c ENDDO
DO I = 1, IM
IF(flag_iter(i).and.FLAG(I)) THEN
if(FLAGSNW(I))then
GFLUX(I) = -DFSNOW * (TSURF(I) - STC(I,1))
& / (FACTSNW(I) * MAX(SNOWD(I),.001))
else
GFLUX(I) = DFT0(I) * (STC(I,1) - TSURF(I))
& / (-.5 * ZSOIL(I,1))
endif
ENDIF
ENDDO
c DO I = 1, IM
c FLAG(I) = SLIMSK(I).EQ.1.
c IF(flag_iter(i))THEN
c IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
c GFLUX(I) = DFT0(I) * (STC(I,1) - TSURF(I))
c & / (-.5 * ZSOIL(I,1))
c ENDIF
c ENDIF
c ENDDO
C
C CALCULATE SENSIBLE HEAT FLUX
C
DO I = 1, IM
IF(flag_iter(i).and. FLAG(I)) THEN
HFLX(I) = RCH(I) * (TSURF(I) - THETA1(I))
ENDIF
ENDDO
C
C THE REST OF THE OUTPUT
C
DO I = 1, IM
IF(flag_iter(i).and. FLAG(I)) THEN
QSURF(I) = Q1(I) + EVAP(I) / (ELOCP * RCH(I))
DM(I) = 1.
C
Cwei added 10/24/2006
EVBS(I)=EDIR(I)
EVCW(I)=EC(I)
SBSNO(I)=SNOEVP(I)
SNOWC(I)=1-PARTLND(I)
STM(I)=-1.0*SMC(I,1)*ZSOIL(I,1)
SNOHF(I)=DFSNOW * (T1(I) - STSOIL(I,1))
TRANS(I)=ET(I,1)
DO K=2,KM
STM(I)=STM(I)+SMC(I,K)*(ZSOIL(I,K-1)-ZSOIL(I,K))
TRANS(I)=TRANS(I)+ET(I,K)
ENDDO
C CONVERT SNOW DEPTH BACK TO MM OF WATER EQUIVALENT
C
SHELEG(I) = SNOWD(I) * 1000.
ENDIF
ENDDO
!
do i=1,im
IF(flag_iter(i).and. FLAG(I)) THEN
tem = 1.0 / rho(i)
hflx(i) = hflx(i) * tem * cpinv
evap(i) = evap(i) * tem * hvapi
ENDIF
enddo
Clu_q2m_iter [+17L]: restore land-related prognostic fields for guess run
do i=1, im
IF(FLAG(I)) THEN
if(flag_guess(i)) then
sheleg(i) = sheleg_old(i)
tskin(i) = tskin_old(i)
canopy(i) = canopy_old(i)
tprcp(i) = tprcp_old(i)
srflag(i) = srflag_old(i)
do k=1, km
stc(i,k) = stc_old(i,k)
enddo
else
tskin(i) = tsurf(i)
endif
ENDIF
enddo
!
C##DG IF(LAT.EQ.LATD) THEN
CC RBAL = -SLWD-SIGMA*TSKIN**4+GFLUX
CC & -EVAP - HFLX
C##DG PRINT 6000,HFLX,EVAP,GFLUX,
C##DG& STC(1), STC(2),TSKIN,RNET,SLWD
C##DG PRINT *, ' T1 =', T1
6000 FORMAT(8(F8.2,','))
CC PRINT *, ' EP, ETP,T2M(I) =', EP, ETP,T2M(I)
CC PRINT *, ' FH, FH2 =', FH, FH2
CC PRINT *, ' PH, PH2 =', PH, PH2
CC PRINT *, ' CH, RCH =', CH, RCH
CC PRINT *, ' TERM1, TERM2 =', TERM1, TERM2
CC PRINT *, ' RS(I), PLANTR =', RS(I), PLANTR
C##DG ENDIF
RETURN
END