!!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_sf_gfs 2
CONTAINS
!-------------------------------------------------------------------
SUBROUTINE SF_GFS(U3D,V3D,T3D,QV3D,P3D, & 3,5
CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
ZNT,UST,PSIM,PSIH, &
XLAND,HFX,QFX,LH,TSK,FLHC,FLQC, &
QGH,QSFC,U10,V10, &
GZ1OZ0,WSPD,BR,ISFFLX, &
EP1,EP2,KARMAN,itimestep, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!-------------------------------------------------------------------
USE MODULE_GFS_MACHINE
, ONLY : kind_phys
USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
!-- U3D 3D u-velocity interpolated to theta points (m/s)
!-- V3D 3D v-velocity interpolated to theta points (m/s)
!-- T3D temperature (K)
!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
!-- P3D 3D pressure (Pa)
!-- CP heat capacity at constant pressure for dry air (J/kg/K)
!-- ROVCP R/CP
!-- R gas constant for dry air (J/kg/K)
!-- XLV latent heat of vaporization for water (J/kg)
!-- PSFC surface pressure (Pa)
!-- ZNT roughness length (m)
!-- UST u* in similarity theory (m/s)
!-- PSIM similarity stability function for momentum
!-- PSIH similarity stability function for heat
!-- XLAND land mask (1 for land, 2 for water)
!-- HFX upward heat flux at the surface (W/m^2)
!-- QFX upward moisture flux at the surface (kg/m^2/s)
!-- LH net upward latent heat flux at surface (W/m^2)
!-- TSK surface temperature (K)
!-- FLHC exchange coefficient for heat (m/s)
!-- FLQC exchange coefficient for moisture (m/s)
!-- QGH lowest-level saturated mixing ratio
!-- U10 diagnostic 10m u wind
!-- V10 diagnostic 10m v wind
!-- GZ1OZ0 log(z/z0) where z0 is roughness length
!-- WSPD wind speed at lowest model level (m/s)
!-- BR bulk Richardson number in surface layer
!-- ISFFLX isfflx=1 for surface heat and moisture fluxes
!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
!-- KARMAN Von Karman constant
!-- ids start index for i in domain
!-- ide end index for i in domain
!-- jds start index for j in domain
!-- jde end index for j in domain
!-- kds start index for k in domain
!-- kde end index for k in domain
!-- ims start index for i in memory
!-- ime end index for i in memory
!-- jms start index for j in memory
!-- jme end index for j in memory
!-- kms start index for k in memory
!-- kme end index for k in memory
!-- its start index for i in tile
!-- ite end index for i in tile
!-- jts start index for j in tile
!-- jte end index for j in tile
!-- kts start index for k in tile
!-- kte end index for k in tile
!-------------------------------------------------------------------
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
ISFFLX,itimestep
REAL, INTENT(IN) :: &
CP, &
EP1, &
EP2, &
KARMAN, &
R, &
ROVCP, &
XLV
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
P3D, &
QV3D, &
T3D, &
U3D, &
V3D
REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
TSK, &
PSFC, &
XLAND
REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
UST, &
ZNT
REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
BR, &
CHS, &
CHS2, &
CPM, &
CQS2, &
FLHC, &
FLQC, &
GZ1OZ0, &
HFX, &
LH, &
PSIM, &
PSIH, &
QFX, &
QGH, &
QSFC, &
U10, &
V10, &
WSPD
!--------------------------- LOCAL VARS ------------------------------
REAL :: ESAT
REAL (kind=kind_phys) :: &
RHOX
REAL (kind=kind_phys), DIMENSION(its:ite) :: &
CH, &
CM, &
DDVEL, &
DRAIN, &
EP, &
EVAP, &
FH, &
FH2, &
FM, &
HFLX, &
PH, &
PM, &
PRSL1, &
PRSLKI, &
PS, &
Q1, &
Q2M, &
QSS, &
QSURF, &
RB, &
RCL, &
RHO1, &
SLIMSK, &
STRESS, &
T1, &
T2M, &
THGB, &
THX, &
TSKIN, &
SHELEG, &
U1, &
U10M, &
USTAR, &
V1, &
V10M, &
WIND, &
Z0RL, &
Z1
INTEGER :: &
I, &
IM, &
J, &
K, &
KM
if(itimestep.eq.0) then
CALL GFUNCPHYS
endif
IM=ITE-ITS+1
KM=KTE-KTS+1
DO J=jts,jte
DO i=its,ite
DDVEL(I)=0.
RCL(i)=1.
PRSL1(i)=P3D(i,kts,j)*.001
PS(i)=PSFC(i,j)*.001
Q1(I) = QV3D(i,kts,j)
! QSURF(I)=QSFC(I,J)
QSURF(I)=0.
SHELEG(I)=0.
SLIMSK(i)=ABS(XLAND(i,j)-2.)
TSKIN(i)=TSK(i,j)
T1(I) = T3D(i,kts,j)
U1(I) = U3D(i,kts,j)
USTAR(I) = UST(i,j)
V1(I) = V3D(i,kts,j)
Z0RL(I) = ZNT(i,j)*100.
ENDDO
DO i=its,ite
PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
Q1(I)=Q1(I)/(1.+Q1(I))
ENDDO
CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
SHELEG,TSKIN,QSURF, &
!WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, &
!WRF SLRAD,SNOWMT,DELT, &
Z0RL, &
!WRF TG3,GFLUX,F10M, &
U10M,V10M,T2M,Q2M, &
!WRF ZSOIL, &
CM,CH,RB, &
!WRF RHSCNPY,RHSMC,AIM,BIM,CIM, &
RCL,PRSL1,PRSLKI,SLIMSK, &
DRAIN,EVAP,HFLX,STRESS,EP, &
FM,FH,USTAR,WIND,DDVEL, &
PM,PH,FH2,QSS,Z1 )
DO i=its,ite
U10(i,j)=U10M(i)
V10(i,j)=V10M(i)
BR(i,j)=RB(i)
CHS(I,J)=CH(I)*WIND(I)
CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
esat = fpvs(t1(i))
QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
QSFC(I,J)=qss(i)
PSIH(i,j)=PH(i)
PSIM(i,j)=PM(i)
UST(i,j)=ustar(i)
WSPD(i,j)=WIND(i)
ZNT(i,j)=Z0RL(i)*.01
ENDDO
DO i=its,ite
FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
FLQC(i,j)=RHO1(I)*CHS(I,J)
GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
CQS2(i,j)=CHS2(I,J)
ENDDO
IF (ISFFLX.EQ.0) THEN
DO i=its,ite
HFX(i,j)=0.
LH(i,j)=0.
QFX(i,j)=0.
ENDDO
ELSE
DO i=its,ite
IF(XLAND(I,J)-1.5.GT.0.)THEN
HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
HFX(I,J)=AMAX1(HFX(I,J),-250.)
ENDIF
QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
QFX(I,J)=AMAX1(QFX(I,J),0.)
LH(I,J)=XLV*QFX(I,J)
ENDDO
ENDIF
ENDDO
END SUBROUTINE SF_GFS
!-------------------------------------------------------------------
SUBROUTINE PROGTM(IM,KM,PS,U1,V1,T1,Q1, & 1,17
& SHELEG,TSKIN,QSURF, &
!WRF & SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY, &
!WRF & DLWFLX,SLRAD,SNOWMT,DELT, &
& Z0RL, &
!WRF & TG3,GFLUX,F10M, &
& U10M,V10M,T2M,Q2M, &
!WRF & ZSOIL, &
& CM, CH, RB, &
!WRF & RHSCNPY,RHSMC,AIM,BIM,CIM, &
& RCL,PRSL1,PRSLKI,SLIMSK, &
& DRAIN,EVAP,HFLX,STRESS,EP, &
& FM,FH,USTAR,WIND,DDVEL, &
& PM,PH,FH2,QSS,Z1 )
!
USE MODULE_GFS_MACHINE, ONLY : kind_phys
USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
USE MODULE_GFS_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), F10M(IM), U10M(IM), &
& V10M(IM), T2M(IM), Q2M(IM), &
& ZSOIL(IM,KM), CM(IM), CH(IM), &
& RB(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), &
& FM(IM), FH(IM), USTAR(IM), &
& WIND(IM), DDVEL(IM), STRESS(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), FM10(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), TSURFO(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, &
!WRF & fhs, funcdf, funckt,g, hl0, hl0inf, &
& fhs, 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
!
!
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)
!WRF real(kind=kind_phys) KT1(IM), KT2(IM), KTSOIL, &
real(kind=kind_phys) KT1(IM), KT2(IM), &
& 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)
!
! the 13 vegetation types are:
!
! 1 ... broadleave-evergreen trees (tropical forest)
! 2 ... broadleave-deciduous trees
! 3 ... broadleave and needle leave trees (mixed forest)
! 4 ... needleleave-evergreen trees
! 5 ... needleleave-deciduous trees (larch)
! 6 ... broadleave trees with groundcover (savanna)
! 7 ... groundcover only (perenial)
! 8 ... broadleave shrubs with perenial groundcover
! 9 ... broadleave shrubs with bare soil
! 10 ... dwarf trees and shrubs with ground cover (trunda)
! 11 ... bare soil
! 12 ... cultivations (use parameters from type 7)
! 13 ... glacial
!
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
!
!WRF DELT2 = DELT * 2.
!
! ESTIMATE SIGMA ** K AT 2 M
!
SIG2K = 1. - 4. * G * 2. / (CP * 280.)
!
! INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIE
! PSURF IS IN PASCALS
! WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
! RHO IS DENSITY, QS1 IS SAT. HUM. AT LEVEL1 AND QSS IS SAT. HUM. AT
! SURFACE
! CONVERT SLRAD TO THE CIVILIZED UNIT FROM LANGLEY MINUTE-1 K-4
! SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM
!
!!
! qs1 = fpvs(t1)
! qss = fpvs(tskin)
DO I=1,IM
XRCL(I) = SQRT(RCL(I))
PSURF(I) = 1000. * PS(I)
PS1(I) = 1000. * PRSL1(I)
! SLWD(I) = SLRAD(I) * CONVRAD
!WRF SLWD(I) = SLRAD(I)
!
! DLWFLX has been given a negative sign for downward longwave
! snet is the net shortwave flux
!
!WRF SNET(I) = -SLWD(I) - DLWFLX(I)
WIND(I) = XRCL(I) * SQRT(U1(I) * U1(I) + V1(I) * V1(I)) &
& + MAX(0.0_kind_phys, MIN(DDVEL(I), 30.0_kind_phys))
WIND(I) = MAX(WIND(I),1._kind_phys)
Q0(I) = MAX(Q1(I),1.E-8_kind_phys)
TSURF(I) = TSKIN(I)
THETA1(I) = T1(I) * PRSLKI(I)
TV1(I) = T1(I) * (1. + RVRDM1 * Q0(I))
THV1(I) = THETA1(I) * (1. + RVRDM1 * Q0(I))
TVS(I) = TSURF(I) * (1. + RVRDM1 * Q0(I))
RHO(I) = PS1(I) / (RD * TV1(I))
!jfe 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_kind_phys)
Q0(I) = min(QS1(I),Q0(I))
!jfe QSS(I) = 1000. * FPVS(TSURF(I))
qss(i) = fpvs(tskin(i))
QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
! RS = PLANTR
RS(I) = 0.
!WRF if(VEGTYPE(I).gt.0.) RS(I) = rsmin(VEGTYPE(I))
Z0(I) = .01 * Z0RL(i)
!WRF CANOPY(I)= MAX(CANOPY(I),0._kind_phys)
DM(I) = 1.
!WRF
GOTO 1111
!WRF
FACTSNW(I) = 10.
IF(SLIMSK(I).EQ.2.) FACTSNW(I) = 3.
!
! SNOW DEPTH IN WATER EQUIVALENT IS CONVERTED FROM MM TO M UNIT
!
SNOWD(I) = SHELEG(I) / 1000.
FLAGSNW(I) = .FALSE.
!
! WHEN SNOW DEPTH IS LESS THAN 1 MM, A PATCHY SNOW IS ASSUMED AND
! SOIL IS ALLOWED TO INTERACT WITH THE ATMOSPHERE.
! WE SHOULD EVENTUALLY MOVE TO A LINEAR COMBINATION OF SOIL AND
! SNOW UNDER THE CONDITION OF PATCHY SNOW.
!
IF(SNOWD(I).GT..001.OR.SLIMSK(I).EQ.2.) RS(I) = 0.
IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE.
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' WIND,TV1,TVS,Q1,QS1,SNOW,SLIMSK=',
!##DG& WIND,TV1,TVS,Q1,QS1,SNOWD,SLIMSK
!##DG PRINT *, ' SNET, SLWD =', SNET, SLWD(I)
!##DG ENDIF
IF(SLIMSK(I).EQ.0.) THEN
ZSOIL(I,1) = 0.
ELSEIF(SLIMSK(I).EQ.1.) THEN
ZSOIL(I,1) = -.10
ELSE
ZSOIL(I,1) = -3. / KM
ENDIF
!WRF
1111 CONTINUE
!WRF
ENDDO
!!
!WRF
GOTO 2222
!WRF
DO K = 2, KM
DO I=1,IM
IF(SLIMSK(I).EQ.0.) THEN
ZSOIL(I,K) = 0.
ELSEIF(SLIMSK(I).EQ.1.) THEN
ZSOIL(I,K) = ZSOIL(I,K-1) &
& + (-2. - ZSOIL(I,1)) / (KM - 1)
ELSE
ZSOIL(I,K) = - 3. * FLOAT(K) / FLOAT(KM)
ENDIF
ENDDO
ENDDO
!WRF
2222 CONTINUE
!WRF
!!
DO I=1,IM
Z1(I) = -RD * TV1(I) * LOG(PS1(I)/PSURF(I)) / G
DRAIN(I) = 0.
ENDDO
!!
DO K = 1, KM
DO I=1,IM
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)
ENDDO
ENDDO
DO I=1,IM
EDIR(I) = 0.
EC(I) = 0.
EVAP(I) = 0.
EP(I) = 0.
SNOWMT(I) = 0.
GFLUX(I) = 0.
RHSCNPY(I) = 0.
FX(I) = 0.
ETPFAC(I) = 0.
CANFAC(I) = 0.
ENDDO
!
! COMPUTE STABILITY DEPENDENT EXCHANGE COEFFICIENTS
!
! THIS PORTION OF THE CODE IS PRESENTLY SUPPRESSED
!
DO I=1,IM
IF(SLIMSK(I).EQ.0.) THEN
USTAR(I) = SQRT(G * Z0(I) / CHARNOCK)
ENDIF
!
! COMPUTE STABILITY INDICES (RB AND HLINF)
!
Z0MAX(I) = MIN(Z0(I),0.1 * Z1(I))
ZTMAX(I) = Z0MAX(I)
IF(SLIMSK(I).EQ.0.) THEN
RESTAR = USTAR(I) * Z0MAX(I) / VIS
RESTAR = MAX(RESTAR,.000001_kind_phys)
! RESTAR = ALOG(RESTAR)
! RESTAR = MIN(RESTAR,5.)
! RESTAR = MAX(RESTAR,-5.)
! RAT(I) = AA1 + BB1 * RESTAR + CC1 * RESTAR ** 2
! RAT(I) = RAT(I) / (1. + BB2 * RESTAR
! & + CC2 * RESTAR ** 2)
! Rat taken from Zeng, Zhao and Dickinson 1997
RAT(I) = 2.67 * restar ** .25 - 2.57
RAT(I) = min(RAT(I),7._kind_phys)
ZTMAX(I) = Z0MAX(I) * EXP(-RAT(I))
ENDIF
ENDDO
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' z0max, ztmax, restar, RAT(I) =',
!##DG & z0max, ztmax, restar, RAT(I)
!##DG ENDIF
DO I = 1, IM
DTV(I) = THV1(I) - TVS(I)
ADTV = ABS(DTV(I))
ADTV = MAX(ADTV,.001_kind_phys)
DTV(I) = SIGN(1._kind_phys,DTV(I)) * ADTV
RB(I) = G * DTV(I) * Z1(I) / (.5 * (THV1(I) + TVS(I)) &
& * WIND(I) * WIND(I))
RB(I) = MAX(RB(I),-5000._kind_phys)
! FM(I) = LOG((Z0MAX(I)+Z1(I)) / Z0MAX(I))
! FH(I) = LOG((ZTMAX(I)+Z1(I)) / ZTMAX(I))
FM(I) = LOG((Z1(I)) / Z0MAX(I))
FH(I) = LOG((Z1(I)) / ZTMAX(I))
HLINF(I) = RB(I) * FM(I) * FM(I) / FH(I)
FM10(I) = LOG((Z0MAX(I)+10.) / Z0MAX(I))
FH2(I) = LOG((ZTMAX(I)+2.) / ZTMAX(I))
ENDDO
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' DTV, RB(I), FM(I), FH(I), HLINF =',
!##DG & dtv, rb, FM(I), FH(I), hlinf
!##DG ENDIF
!
! STABLE CASE
!
DO I = 1, IM
IF(DTV(I).GE.0.) THEN
HL1(I) = HLINF(I)
ENDIF
IF(DTV(I).GE.0..AND.HLINF(I).GT..25) THEN
HL0INF = Z0MAX(I) * HLINF(I) / Z1(I)
HLTINF = ZTMAX(I) * HLINF(I) / Z1(I)
AA = SQRT(1. + 4. * ALPHA * HLINF(I))
AA0 = SQRT(1. + 4. * ALPHA * HL0INF)
BB = AA
BB0 = SQRT(1. + 4. * ALPHA * HLTINF)
PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
FMS = FM(I) - PM(I)
FHS = FH(I) - PH(I)
HL1(I) = FMS * FMS * RB(I) / FHS
ENDIF
ENDDO
!
! SECOND ITERATION
!
DO I = 1, IM
IF(DTV(I).GE.0.) THEN
HL0 = Z0MAX(I) * HL1(I) / Z1(I)
HLT = ZTMAX(I) * HL1(I) / Z1(I)
AA = SQRT(1. + 4. * ALPHA * HL1(I))
AA0 = SQRT(1. + 4. * ALPHA * HL0)
BB = AA
BB0 = SQRT(1. + 4. * ALPHA * HLT)
PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
HL110 = HL1(I) * 10. / Z1(I)
AA = SQRT(1. + 4. * ALPHA * HL110)
PM10(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
HL12(I) = HL1(I) * 2. / Z1(I)
! AA = SQRT(1. + 4. * ALPHA * HL12(I))
BB = SQRT(1. + 4. * ALPHA * HL12(I))
PH2(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
ENDIF
ENDDO
!!
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' HL1(I), PM, PH =',
!##DG & HL1(I), pm, ph
!##DG ENDIF
!
! UNSTABLE CASE
!
!
! CHECK FOR UNPHYSICAL OBUKHOV LENGTH
!
DO I=1,IM
IF(DTV(I).LT.0.) THEN
OLINF = Z1(I) / HLINF(I)
IF(ABS(OLINF).LE.50. * Z0MAX(I)) THEN
HLINF(I) = -Z1(I) / (50. * Z0MAX(I))
ENDIF
ENDIF
ENDDO
!
! GET PM AND PH
!
DO I = 1, IM
IF(DTV(I).LT.0..AND.HLINF(I).GE.-.5) THEN
HL1(I) = HLINF(I)
PM(I) = (A0 + A1 * HL1(I)) * HL1(I) &
& / (1. + B1 * HL1(I) + B2 * HL1(I) * HL1(I))
PH(I) = (A0P + A1P * HL1(I)) * HL1(I) &
& / (1. + B1P * HL1(I) + B2P * HL1(I) * HL1(I))
HL110 = HL1(I) * 10. / Z1(I)
PM10(I) = (A0 + A1 * HL110) * HL110 &
& / (1. + B1 * HL110 + B2 * HL110 * HL110)
HL12(I) = HL1(I) * 2. / Z1(I)
PH2(I) = (A0P + A1P * HL12(I)) * HL12(I) &
& / (1. + B1P * HL12(I) + B2P * HL12(I) * HL12(I))
ENDIF
IF(DTV(I).LT.0.AND.HLINF(I).LT.-.5) THEN
HL1(I) = -HLINF(I)
PM(I) = LOG(HL1(I)) + 2. * HL1(I) ** (-.25) - .8776
PH(I) = LOG(HL1(I)) + .5 * HL1(I) ** (-.5) + 1.386
HL110 = HL1(I) * 10. / Z1(I)
PM10(I) = LOG(HL110) + 2. * HL110 ** (-.25) - .8776
HL12(I) = HL1(I) * 2. / Z1(I)
PH2(I) = LOG(HL12(I)) + .5 * HL12(I) ** (-.5) + 1.386
ENDIF
ENDDO
!
! FINISH THE EXCHANGE COEFFICIENT COMPUTATION TO PROVIDE FM AND FH
!
DO I = 1, IM
FM(I) = FM(I) - PM(I)
FH(I) = FH(I) - PH(I)
FM10(I) = FM10(I) - PM10(I)
FH2(I) = FH2(I) - PH2(I)
CM(I) = CA * CA / (FM(I) * FM(I))
CH(I) = CA * CA / (FM(I) * FH(I))
CQ = CH(I)
STRESS(I) = CM(I) * WIND(I) * WIND(I)
USTAR(I) = SQRT(STRESS(I))
! USTAR(I) = SQRT(CM(I) * WIND(I) * WIND(I))
ENDDO
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' FM, FH, CM, CH(I), USTAR =',
!##DG & FM, FH, CM, ch, USTAR
!##DG ENDIF
!
! UPDATE Z0 OVER OCEAN
!
DO I = 1, IM
IF(SLIMSK(I).EQ.0.) THEN
Z0(I) = (CHARNOCK / G) * USTAR(I) ** 2
! NEW IMPLEMENTATION OF Z0
! CC = USTAR(I) * Z0 / RNU
! PP = CC / (1. + CC)
! FF = G * ARNU / (CHARNOCK * USTAR(I) ** 3)
! Z0 = ARNU / (USTAR(I) * FF ** PP)
Z0(I) = MIN(Z0(I),.1_kind_phys)
Z0(I) = MAX(Z0(I),1.E-7_kind_phys)
Z0RL(I) = 100. * Z0(I)
ENDIF
ENDDO
GOTO 5555
!
! RCP = RHO CP CH V
!
DO I = 1, IM
RCH(I) = RHO(I) * CP * CH(I) * WIND(I)
ENDDO
!
! SENSIBLE AND LATENT HEAT FLUX OVER OPEN WATER
!
DO I = 1, IM
IF(SLIMSK(I).EQ.0.) THEN
EVAP(I) = ELOCP * RCH(I) * (QSS(I) - Q1(I))
DM(I) = 1.
QSURF(I) = QSS(I)
ENDIF
ENDDO
!
! COMPUTE SOIL/SNOW/ICE HEAT FLUX IN PREPARATION FOR SURFACE ENERGY
! BALANCE CALCULATION
!
DO I = 1, IM
GFLUX(I) = 0.
IF(SLIMSK(I).EQ.1.) THEN
SMCZ(I) = .5 * (SMC(I,1) + .20)
DFT0(I) = KTSOIL(SMCZ(I),SOILTYP(I))
ELSEIF(SLIMSK(I).EQ.2.) THEN
! DF FOR ICE IS TAKEN FROM MAYKUT AND UNTERSTEINER
! DF IS IN SI UNIT OF W K-1 M-1
DFT0(I) = 2.2
ENDIF
ENDDO
!!
DO I=1,IM
IF(SLIMSK(I).NE.0.) THEN
! IF(SNOWD(I).GT..001) THEN
IF(FLAGSNW(I)) THEN
!
! WHEN SNOW COVERED, GROUND HEAT FLUX COMES FROM SNOW
!
TFLX = MIN(T1(I), TSURF(I))
GFLUX(I) = -DFSNOW * (TFLX - STSOIL(I,1)) &
& / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
ELSE
GFLUX(I) = DFT0(I) * (STSOIL(I,1) - TSURF(I)) &
& / (-.5 * ZSOIL(I,1))
ENDIF
GFLUX(I) = MAX(GFLUX(I),-200._kind_phys)
GFLUX(I) = MIN(GFLUX(I),+200._kind_phys)
ENDIF
ENDDO
DO I = 1, IM
FLAG(I) = SLIMSK(I).NE.0.
PARTLND(I) = 1.
IF(SNOWD(I).GT.0..AND.SNOWD(I).LE..001) THEN
PARTLND(I) = 1. - SNOWD(I) / .001
ENDIF
ENDDO
DO I = 1, IM
SNOEVP(I) = 0.
if(SNOWD(I).gt..001) PARTLND(I) = 0.
ENDDO
!
! COMPUTE POTENTIAL EVAPORATION FOR LAND AND SEA ICE
!
DO I = 1, IM
IF(FLAG(I)) THEN
T12 = T1(I) * T1(I)
T14 = T12 * T12
!
! RCAP = FNET - SIGMA T**4 + GFLX - RHO CP CH V (T1-THETA1)
!
RCAP(I) = -SLWD(I) - SIGMA * T14 + GFLUX(I) &
& - RCH(I) * (T1(I) - THETA1(I))
!
! RSMALL = 4 SIGMA T**3 / RCH(I) + 1
!
RSMALL(I) = 4. * SIGMA * T1(I) * T12 / RCH(I) + 1.
!
! DELTA = L / CP * DQS/DT
!
DELTA(I) = ELOCP * EPS * HVAP * QS1(I) / (RD * T12)
!
! POTENTIAL EVAPOTRANSPIRATION ( WATTS / M**2 ) AND
! POTENTIAL EVAPORATION
!
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
!
! ACTUAL EVAPORATION OVER LAND IN THREE PARTS : EDIR, ET, AND EC
!
! DIRECT EVAPORATION FROM SOIL, THE UNIT GOES FROM M S-1 TO KG M-2 S-1
!
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ.1..AND.EP(I).GT.0.
ENDDO
DO I = 1, IM
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
! TREF = .75 * THSAT(SOILTYP(I))
TREF(I) = smref(SOILTYP(I))
! TWILT = TWLT(SOILTYP(I))
TWILT(I) = smwlt(SOILTYP(I))
smcdry = smdry(SOILTYP(I))
! FX(I) = -2. * DF1(I) * (SMC(I,1) - .23) / ZSOIL(I,1)
! & - 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._kind_phys)
!
! SIGMAF IS THE FRACTION OF AREA COVERED BY VEGETATION
!
EDIR(I) = FX(I) * (1. - SIGMAF(I)) * PARTLND(I)
ENDIF
ENDDO
!
! calculate stomatal resistance
!
DO I = 1, IM
if(FLAG(I)) then
!
! resistance due to PAR. We use net solar flux as proxy at the present time
!
ff = .55 * 2. * SNET(I) / rgl(VEGTYPE(I))
rcs = (ff + RS(I)/rsmax(VEGTYPE(I))) / (1. + ff)
rcs = max(rcs,.0001_kind_phys)
rct = 1.
rcq = 1.
!
! resistance due to thermal effect
!
! rct = 1. - .0016 * (topt - theta1) ** 2
! rct = max(rct,.0001)
!
! resistance due to humidity
!
! rcq = 1. / (1. + hs(VEGTYPE(I)) * (QS1(I) - Q0(I)))
! rcq = max(rcq,.0001)
!
! compute resistance without the effect of soil moisture
!
RS(I) = RS(I) / (rcs * rct * rcq)
endif
ENDDO
!
! TRANSPIRATION FROM ALL LEVELS OF THE SOIL
!
DO I = 1, IM
IF(FLAG(I)) THEN
CANFAC(I) = (CANOPY(I) / SCANOP) ** CFACTR
endif
IF(FLAG(I)) THEN
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._kind_phys)
GX(I) = MIN(GX(I),1._kind_phys)
!
! resistance due to soil moisture deficit
!
rss = GX(I) * (ZSOIL(I,1) / ZSOIL(I,km))
rss = max(rss,.0001_kind_phys)
RSI = RS(I) / rss
!
! transpiration a la Monteith
!
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(I)) THEN
GX(I) = (SMC(I,K) - TWILT(I)) / (TREF(I) - TWILT(I))
GX(I) = MAX(GX(I),0._kind_phys)
GX(I) = MIN(GX(I),1._kind_phys)
!
! resistance due to soil moisture deficit
!
rss = GX(I) * ((ZSOIL(I,k) - ZSOIL(I,k-1))/ZSOIL(I,km))
rss = max(rss,1.e-6_kind_phys)
RSI = RS(I) / rss
!
! transpiration a la Monteith
!
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
!
! CANOPY RE-EVAPORATION
!
DO I=1,IM
IF(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
!
! SUM UP TOTAL EVAPORATION
!
DO I = 1, IM
IF(FLAG(I)) THEN
EVAP(I) = EDIR(I) + EC(I)
ENDIF
ENDDO
!!
DO K = 1, KM
DO I=1,IM
IF(FLAG(I)) THEN
EVAP(I) = EVAP(I) + ET(I,K)
ENDIF
ENDDO
ENDDO
!!
!
! RETURN EVAP UNIT FROM KG M-2 S-1 TO WATTS M-2
!
DO I=1,IM
IF(FLAG(I)) THEN
EVAP(I) = MIN(EVAP(I)*HVAP,EP(I))
ENDIF
ENDDO
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, 'FX(I), SIGMAF, EDIR(I), ETPFAC=', FX(I)*HVAP,SIGMAF,
!##DG& EDIR(I)*HVAP,ETPFAC*HVAP
!##DG PRINT *, ' ET =', (ET(K)*HVAP,K=1,KM)
!##DG PRINT *, ' CANFAC(I), EC(I), EVAP', CANFAC(I),EC(I)*HVAP,EVAP
!##DG ENDIF
!
! EVAPORATION OVER BARE SEA ICE
!
DO I = 1, IM
! IF(SLIMSK(I).EQ.2.AND.SNOWD(I).LE..001) THEN
IF(SLIMSK(I).EQ.2.) THEN
EVAP(I) = PARTLND(I) * EP(I)
ENDIF
ENDDO
!
! TREAT DOWNWARD MOISTURE FLUX SITUATION
! (EVAP WAS PRESET TO ZERO SO NO UPDATE NEEDED)
! DEW IS CONVERTED FROM KG M-2 TO M TO CONFORM TO PRECIP UNIT
!
DO I = 1, IM
FLAG(I) = SLIMSK(I).NE.0..AND.EP(I).LE.0.
DEW(I) = 0.
ENDDO
DO I = 1, IM
IF(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
!
! SNOW COVERED LAND AND SEA ICE
!
DO I = 1, IM
FLAG(I) = SLIMSK(I).NE.0..AND.SNOWD(I).GT.0.
ENDDO
!
! CHANGE OF SNOW DEPTH DUE TO EVAPORATION OR SUBLIMATION
!
! CONVERT EVAP FROM KG M-2 S-1 TO M S-1 TO DETERMINE THE REDUCTION OF S
!
DO I = 1, IM
IF(FLAG(I)) THEN
BFACT = SNOWD(I) / (DELT * EP(I) / (HVAP * RHOH2O))
BFACT = MIN(BFACT,1._kind_phys)
!
! THE EVAPORATION OF SNOW
!
IF(EP(I).LE.0.) BFACT = 1.
IF(SNOWD(I).LE..001) THEN
! 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
! EVAP(I) = BFACT * EP(I)
SNOEVP(I) = bfact * EP(I)
EVAP(I) = SNOEVP(I)
ENDIF
TSURF(I) = T1(I) + &
& (RCAP(I) - GFLUX(I) - DFSNOW * (T1(I) - STSOIL(I,1)) &
& /(FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys)) &
! & + THETA1 - T1 &
! & - BFACT * EP(I)) / (RSMALL(I) * RCH(I) &
& - SNOEVP(I)) / (RSMALL(I) * RCH(I) &
& + DFSNOW / (FACTSNW(I)* MAX(SNOWD(I),.001_kind_phys)))
! 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._kind_phys)
ENDIF
ENDDO
!
! SNOW MELT (M)
!
500 CONTINUE
DO I = 1, IM
FLAG(I) = SLIMSK(I).NE.0. &
& .AND.SNOWD(I).GT..0
ENDDO
DO I = 1, IM
IF(FLAG(I).AND.TSURF(I).GT.T0C) THEN
SNOWMT(I) = RCH(I) * RSMALL(I) * DELT &
& * (TSURF(I) - T0C) / (RHOH2O * HFUS)
SNOWMT(I) = min(SNOWMT(I),SNOWD(I))
SNOWD(I) = SNOWD(I) - SNOWMT(I)
SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
TSURF(I) = MAX(T0C,TSURF(I) &
& -HFUS*SNOWMT(I)*RHOH2O/(RCH(I)*RSMALL(I)*DELT))
ENDIF
ENDDO
!
! We need to re-evaluate evaporation because of snow melt
! the skin temperature is now bounded to 0 deg C
!
! qss = fpvs(tsurf)
DO I = 1, IM
! IF (SNOWD(I) .GT. 0.0) THEN
IF (SNOWD(I) .GT. snomin) THEN
!jfe QSS(I) = 1000. * FPVS(TSURF(I))
qss(i) = fpvs(tsurf(i))
QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
EVAP(I) = elocp * RCH(I) * (QSS(I) - Q0(I))
ENDIF
ENDDO
!
! PREPARE TENDENCY TERMS FOR THE SOIL MOISTURE FIELD WITHOUT PRECIPITAT
! THE UNIT OF MOISTURE FLUX NEEDS TO BECOME M S-1 FOR SOIL MOISTURE
! HENCE THE FACTOR OF RHOH2O
!
DO I = 1, IM
FLAG(I) = SLIMSK(I).EQ.1.
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))
!
! AIM, BIM, AND CIM ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
! IMPLICIT UPDATE OF THE SOIL MOISTURE
!
AIM(I,1) = 0.
BIM(I,1) = DF1(I) * DDZ(I) / (-ZSOIL(I,1) * RHOH2O)
CIM(I,1) = -BIM(I,1)
ENDIF
ENDDO
!!
DO K = 2, KM
IF(K.LT.KM) THEN
DO I=1,IM
IF(FLAG(I)) THEN
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
ENDDO
ELSE
DO I = 1, IM
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
ENDDO
ENDIF
DO I = 1, IM
IF(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
!
! UPDATE SOIL TEMPERATURE AND SEA ICE TEMPERATURE
!
DO I=1,IM
FLAG(I) = SLIMSK(I).NE.0.
ENDDO
!
! SURFACE TEMPERATURE IS PART OF THE UPDATE WHEN SNOW IS ABSENT
!
DO I=1,IM
! IF(FLAG(I).AND.SNOWD(I).LE..001) THEN
IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
YY(I) = T1(I) + &
! & (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
! IF(FLAG(I).AND.SNOWD(I).GT..001) THEN
IF(FLAG(I).AND.FLAGSNW(I)) THEN
YY(I) = STSOIL(I,1)
!
! HEAT FLUX FROM SNOW IS EXPLICIT IN TIME
!
ZZ(I) = 1.
XX(I) = DFSNOW * (STSOIL(I,1) - TSURF(I)) &
& / (-FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
ENDIF
ENDDO
!
! COMPUTE THE FORCING AND THE IMPLICIT MATRIX ELEMENTS FOR UPDATE
!
! CH2O IS THE HEAT CAPACITY OF WATER AND CSOIL IS THE HEAT CAPACITY OF
!
DO I = 1, IM
IF(FLAG(I)) THEN
SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
DTDZ1(I) = (STSOIL(I,1) - STSOIL(I,2)) / (-.5 * ZSOIL(I,2))
IF(SLIMSK(I).EQ.1.) THEN
DFT1(I) = KTSOIL(SMCZ(I),SOILTYP(I))
HCPCT(I) = SMC(I,1) * CH2O + (1. - SMC(I,1)) * CSOIL
ELSE
DFT1(I) = DFT0(I)
HCPCT(I) = CICE
ENDIF
DFT2(I) = DFT1(I)
DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
!
! AI, BI, AND CI ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
! IMPLICIT UPDATE OF THE SOIL TEMPERATURE
!
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))
! SS = DFT0(I) * (STSOIL(I,1) - YY(I)) &
! & / (.5 * ZSOIL(I,1) * ZZ(I))
! 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(SLIMSK(I).EQ.1.) THEN
HCPCT(I) = SMC(I,K) * CH2O + (1. - SMC(I,K)) * CSOIL
ELSEIF(SLIMSK(I).EQ.2.) THEN
HCPCT(I) = CICE
ENDIF
ENDDO
IF(K.LT.KM) THEN
DO I = 1, IM
IF(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))
IF(SLIMSK(I).EQ.1.) THEN
DFT2(I) = KTSOIL(SMCZ(I),SOILTYP(I))
ENDIF
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
!
! AT THE BOTTOM, CLIMATOLOGY IS ASSUMED AT 2M DEPTH FOR LAND AND
! FREEZING TEMPERATURE IS ASSUMED FOR SEA ICE AT Z(KM)
DO I = 1, IM
IF(SLIMSK(I).EQ.1.) 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
IF(SLIMSK(I).EQ.2.) THEN
DTDZ2(I) = (STSOIL(I,K) - TGICE) &
& / (.5 * ZSOIL(I,K-1) - .5 * ZSOIL(I,K))
DFT2(I) = DFT1(I)
CI(I,K) = 0.
ENDIF
ENDDO
ENDIF
DO I = 1, IM
IF(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
!
! SOLVE THE TRI-DIAGONAL MATRIX
!
DO K = 1, KM
DO I=1,IM
IF(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
! FORWARD ELIMINATION
DO I=1,IM
IF(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(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
!!
! BACKWARD SUBSTITUTTION
DO I=1,IM
IF(FLAG(I)) THEN
CI(I,KM) = RHSTC(I,KM)
ENDIF
ENDDO
!!
DO K = KM-1, 1
DO I=1,IM
IF(FLAG(I)) THEN
CI(I,K) = CI(I,K) * CI(I,K+1) + RHSTC(I,K)
ENDIF
ENDDO
ENDDO
!
! UPDATE SOIL AND ICE TEMPERATURE
!
DO K = 1, KM
DO I=1,IM
IF(FLAG(I)) THEN
STSOIL(I,K) = STSOIL(I,K) + CI(I,K)
ENDIF
ENDDO
ENDDO
!
! UPDATE SURFACE TEMPERATURE FOR SNOW FREE SURFACES
!
DO I=1,IM
! IF(SLIMSK(I).NE.0..AND.SNOWD(I).LE..001) THEN
IF(SLIMSK(I).NE.0..AND..NOT.FLAGSNW(I)) THEN
TSURF(I) = (YY(I) + (ZZ(I) - 1.) * STSOIL(I,1)) / ZZ(I)
ENDIF
! IF(SLIMSK(I).EQ.2..AND.SNOWD(I).LE..001) THEN
IF(SLIMSK(I).EQ.2..AND..NOT.FLAGSNW(I)) THEN
TSURF(I) = MIN(TSURF(I),T0C)
ENDIF
ENDDO
!!
DO K = 1, KM
DO I=1,IM
IF(SLIMSK(I).EQ.2) THEN
STSOIL(I,K) = MIN(STSOIL(I,K),T0C)
ENDIF
ENDDO
ENDDO
!
! TIME FILTER FOR SOIL AND SKIN TEMPERATURE
!
DO I=1,IM
IF(SLIMSK(I).NE.0.) THEN
TSKIN(I) = CTFIL1 * TSURF(I) + CTFIL2 * TSKIN(I)
ENDIF
ENDDO
DO K = 1, KM
DO I=1,IM
IF(SLIMSK(I).NE.0.) THEN
STC(I,K) = CTFIL1 * STSOIL(I,K) + CTFIL2 * STC(I,K)
ENDIF
ENDDO
ENDDO
!
! GFLUX CALCULATION
!
DO I=1,IM
FLAG(I) = SLIMSK(I).NE.0. &
! & .AND.SNOWD(I).GT..001 &
& .AND.FLAGSNW(I)
ENDDO
DO I = 1, IM
IF(FLAG(I)) THEN
GFLUX(I) = -DFSNOW * (TSKIN(I) - STC(I,1)) &
& / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
ENDIF
ENDDO
DO I = 1, IM
! IF(SLIMSK(I).NE.0..AND.SNOWD(I).LE..001) THEN
IF( SLIMSK(I).NE.0..AND..NOT.FLAGSNW(I)) THEN
GFLUX(I) = DFT0(I) * (STC(I,1) - TSKIN(I)) &
& / (-.5 * ZSOIL(I,1))
ENDIF
ENDDO
5555 CONTINUE
!
! CALCULATE SENSIBLE HEAT FLUX
!
!WRF DO I = 1, IM
!WRF HFLX(I) = RCH(I) * (TSKIN(I) - THETA1(I))
!WRF ENDDO
!
! THE REST OF THE OUTPUT
!
!WRF DO I = 1, IM
!WRF QSURF(I) = Q1(I) + EVAP(I) / (ELOCP * RCH(I))
!WRF DM(I) = 1.
!
! CONVERT SNOW DEPTH BACK TO MM OF WATER EQUIVALENT
!
!WRF SHELEG(I) = SNOWD(I) * 1000.
!WRF ENDDO
!
DO I = 1, IM
F10M(I) = FM10(I) / FM(I)
F10M(I) = min(F10M(I),1._kind_phys)
U10M(I) = F10M(I) * XRCL(I) * U1(I)
V10M(I) = F10M(I) * XRCL(I) * V1(I)
!WRF T2M(I) = TSKIN(I) * (1. - FH2(I) / FH(I)) &
!WRF & + THETA1(I) * FH2(I) / FH(I)
!WRF T2M(I) = T2M(I) * SIG2K
! Q2M(I) = QSURF(I) * (1. - FH2(I) / FH(I)) &
! & + Q1(I) * FH2(I) / FH(I)
! T2M(I) = T1
! Q2M(I) = Q1
!WRF IF(EVAP(I).GE.0.) THEN
!
! IN CASE OF EVAPORATION, USE THE INFERRED QSURF TO DEDUCE Q2M
!
!WRF Q2M(I) = QSURF(I) * (1. - FH2(I) / FH(I)) &
!WRF & + Q1(I) * FH2(I) / FH(I)
!WRF ELSE
!
! FOR DEW FORMATION SITUATION, USE SATURATED Q AT TSKIN
!
!jfe QSS(I) = 1000. * FPVS(TSKIN(I))
!WRF qss(I) = fpvs(tskin(I))
!WRF QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
!WRF Q2M(I) = QSS(I) * (1. - FH2(I) / FH(I)) &
!WRF & + Q1(I) * FH2(I) / FH(I)
!WRF ENDIF
!jfe QSS(I) = 1000. * FPVS(T2M(I))
!WRF QSS(I) = fpvs(t2m(I))
! QSS(I) = 1000. * T2MO(I)
!WRF QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
!WRF Q2M(I) = MIN(Q2M(I),QSS(I))
ENDDO
!!
! DO I = 1, IM
! RNET(I) = -SLWD(I) - SIGMA * TSKIN(I) **4
! ENDDO
!!
!
!WRF do i=1,im
!WRF tem = 1.0 / rho(i)
!WRF hflx(i) = hflx(i) * tem * cpinv
!WRF evap(i) = evap(i) * tem * hvapi
!WRF enddo
!
!##DG IF(LAT.EQ.LATD) THEN
!C RBAL = -SLWD-SIGMA*TSKIN**4+GFLUX
!C & -EVAP - HFLX
!##DG PRINT 6000,HFLX,EVAP,GFLUX,
!##DG& STC(1), STC(2),TSKIN,RNET,SLWD
!##DG PRINT *, ' T1 =', T1
6000 FORMAT(8(F8.2,','))
!C PRINT *, ' EP, ETP,T2M(I) =', EP, ETP,T2M(I)
!C PRINT *, ' FH, FH2 =', FH, FH2
!C PRINT *, ' PH, PH2 =', PH, PH2
!C PRINT *, ' CH, RCH =', CH, RCH
!C PRINT *, ' TERM1, TERM2 =', TERM1, TERM2
!C PRINT *, ' RS(I), PLANTR =', RS(I), PLANTR
!##DG ENDIF
RETURN
END SUBROUTINE PROGTM
!
! PROGT2 IS THE SECOND PART OF THE SOIL MODEL THAT IS EXECUTED
! AFTER PRECIPITATION FOR THE TIME STEP HAS BEEN CALCULATED
!
!FPP$ NOCONCUR R
!FPP$ EXPAND(FUNCDF,FUNCKT,THSAT)
SUBROUTINE PROGT2(IM,KM,RHSCNPY, &,2
& RHSMC,AI,BI,CI,SMC,SLIMSK, &
& CANOPY,PRECIP,RUNOFF,SNOWMT, &
& ZSOIL,SOILTYP,SIGMAF,DELT,me)
!c
USE MODULE_GFS_MACHINE , ONLY : kind_phys
implicit none
integer km, IM, me
real(kind=kind_phys) delt
real(kind=kind_phys) RHSCNPY(IM), RHSMC(IM,KM), AI(IM,KM), &
& BI(IM,KM), CI(IM,KM), SMC(IM,KM), &
& SLIMSK(IM), CANOPY(IM), PRECIP(IM), &
& RUNOFF(IM), SNOWMT(IM), ZSOIL(IM,KM), &
& SIGMAF(IM)
INTEGER SOILTYP(IM)
!
integer k, lond, i
real(kind=kind_phys) CNPY(IM), PRCP(IM), TSAT(IM), &
& INF(IM), INFMAX(IM), SMSOIL(IM,KM)
!
real(kind=kind_phys) cc, ctfil1, ctfil2, delt2, &
& drip, rffact, rhoh2o, &
!WRF & rzero, scanop, tdif, thsat, KSAT
& rzero, scanop, tdif, KSAT
!
LOGICAL FLAG(IM)
!c
PARAMETER (SCANOP=.5, RHOH2O=1000.)
PARAMETER (CTFIL1=.5, CTFIL2=1.-CTFIL1)
! PARAMETER (CTFIL1=1., CTFIL2=1.-CTFIL1)
PARAMETER (RFFACT=.15)
!
!##DG LATD = 44
LOND = 353
DELT2 = DELT * 2.
!
! PRECIPITATION RATE IS NEEDED IN UNIT OF KG M-2 S-1
!
DO I=1,IM
PRCP(I) = RHOH2O * (PRECIP(I)+SNOWMT(I)) / DELT
RUNOFF(I) = 0.
CNPY(I) = CANOPY(I)
ENDDO
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' BEFORE RUNOFF CAL, RHSMC =', RHSMC(1)
!##DG ENDIF
!
! UPDATE CANOPY WATER CONTENT
!
DO I=1,IM
IF(SLIMSK(I).EQ.1.) THEN
RHSCNPY(I) = RHSCNPY(I) + SIGMAF(I) * PRCP(I)
CANOPY(I) = CANOPY(I) + DELT * RHSCNPY(I)
CANOPY(I) = MAX(CANOPY(I),0._kind_phys)
PRCP(I) = PRCP(I) * (1. - SIGMAF(I))
IF(CANOPY(I).GT.SCANOP) THEN
DRIP = CANOPY(I) - SCANOP
CANOPY(I) = SCANOP
PRCP(I) = PRCP(I) + DRIP / DELT
ENDIF
!
! CALCULATE INFILTRATION RATE
!
INF(I) = PRCP(I)
TSAT(I) = THSAT(SOILTYP(I))
! DSAT = FUNCDF(TSAT(I),SOILTYP(I))
! KSAT = FUNCKT(TSAT(I),SOILTYP(I))
! INFMAX(I) = -DSAT * (TSAT(I) - SMC(I,1))
! & / (.5 * ZSOIL(I,1)) &
! & + KSAT
INFMAX(I) = (-ZSOIL(I,1)) * &
& ((TSAT(I) - SMC(I,1)) / DELT - RHSMC(I,1)) &
& * RHOH2O
INFMAX(I) = MAX(RFFACT*INFMAX(I),0._kind_phys)
! IF(SMC(I,1).GE.TSAT(I)) INFMAX(I) = KSAT
! IF(SMC(I,1).GE.TSAT(I)) INFMAX(I) = ZSOIL(I,1) * RHSMC(I,1)
IF(INF(I).GT.INFMAX(I)) THEN
RUNOFF(I) = INF(I) - INFMAX(I)
INF(I) = INFMAX(I)
ENDIF
INF(I) = INF(I) / RHOH2O
RHSMC(I,1) = RHSMC(I,1) - INF(I) / ZSOIL(I,1)
ENDIF
ENDDO
!!
!##DG IF(LAT.EQ.LATD) THEN
!##DG PRINT *, ' PRCP(I), INFMAX(I), RUNOFF =', PRCP(I),INFMAX(I),RUNOFF
!##DG PRINT *, ' SMSOIL =', SMC(1), SMC(2)
!##DG PRINT *, ' RHSMC =', RHSMC(1)
!##DG ENDIF
!
! WE CURRENTLY IGNORE THE EFFECT OF RAIN ON SEA ICE
!
DO I=1,IM
FLAG(I) = SLIMSK(I).EQ.1.
ENDDO
!!
!
! SOLVE THE TRI-DIAGONAL MATRIX
!
DO K = 1, KM
DO I=1,IM
IF(FLAG(I)) THEN
RHSMC(I,K) = RHSMC(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
! FORWARD ELIMINATION
DO I=1,IM
IF(FLAG(I)) THEN
CI(I,1) = -CI(I,1) / BI(I,1)
RHSMC(I,1) = RHSMC(I,1) / BI(I,1)
ENDIF
ENDDO
DO K = 2, KM
DO I=1,IM
IF(FLAG(I)) THEN
CC = 1. / (BI(I,K) + AI(I,K) * CI(I,K-1))
CI(I,K) = -CI(I,K) * CC
RHSMC(I,K)=(RHSMC(I,K)-AI(I,K)*RHSMC(I,K-1))*CC
ENDIF
ENDDO
ENDDO
! BACKWARD SUBSTITUTTION
DO I=1,IM
IF(FLAG(I)) THEN
CI(I,KM) = RHSMC(I,KM)
ENDIF
ENDDO
!!
DO K = KM-1, 1
DO I=1,IM
IF(FLAG(I)) THEN
CI(I,K) = CI(I,K) * CI(I,K+1) + RHSMC(I,K)
ENDIF
ENDDO
ENDDO
100 CONTINUE
!
! UPDATE SOIL MOISTURE
!
DO K = 1, KM
DO I=1,IM
IF(FLAG(I)) THEN
SMSOIL(I,K) = SMC(I,K) + CI(I,K)
SMSOIL(I,K) = MAX(SMSOIL(I,K),0._kind_phys)
TDIF = MAX(SMSOIL(I,K) - TSAT(I),0._kind_phys)
RUNOFF(I) = RUNOFF(I) - &
& RHOH2O * TDIF * ZSOIL(I,K) / DELT
SMSOIL(I,K) = SMSOIL(I,K) - TDIF
ENDIF
ENDDO
ENDDO
DO K = 1, KM
DO I=1,IM
IF(FLAG(I)) THEN
SMC(I,K) = CTFIL1 * SMSOIL(I,K) + CTFIL2 * SMC(I,K)
ENDIF
ENDDO
ENDDO
! IF(FLAG(I)) THEN
! CANOPY(I) = CTFIL1 * CANOPY(I) + CTFIL2 * CNPY(I)
! ENDIF
! I = 1
! PRINT *, ' SMC'
! PRINT 6000, SMC(1), SMC(2)
!6000 FORMAT(2(F8.5,','))
RETURN
END SUBROUTINE PROGT2
FUNCTION KTSOIL(THETA,KTYPE) 4,2
!
USE MODULE_GFS_MACHINE , ONLY : kind_phys
USE module_progtm , ONLY : TSAT, DFKT
implicit none
integer ktype,kw
real(kind=kind_phys) ktsoil, theta, w
!
W = (THETA / TSAT(KTYPE)) * 20. + 1.
KW = W
KW = MIN(KW,21)
KW = MAX(KW,1)
KTSOIL = DFKT(KW,KTYPE) &
& + (W - KW) * (DFKT(KW+1,KTYPE) - DFKT(KW,KTYPE))
RETURN
END FUNCTION KTSOIL
FUNCTION FUNCDF(THETA,KTYPE) 3,2
!
USE MODULE_GFS_MACHINE , ONLY : kind_phys
USE module_progtm , ONLY : TSAT, DFK
implicit none
integer ktype,kw
real(kind=kind_phys) funcdf,theta,w
!
W = (THETA / TSAT(KTYPE)) * 20. + 1.
KW = W
KW = MIN(KW,21)
KW = MAX(KW,1)
FUNCDF = DFK(KW,KTYPE) &
& + (W - KW) * (DFK(KW+1,KTYPE) - DFK(KW,KTYPE))
RETURN
END FUNCTION FUNCDF
FUNCTION FUNCKT(THETA,KTYPE) 4,2
!
USE MODULE_GFS_MACHINE , ONLY : kind_phys
USE module_progtm , ONLY : TSAT, KTK
implicit none
integer ktype,kw
real(kind=kind_phys) funckt,theta,w
!
W = (THETA / TSAT(KTYPE)) * 20. + 1.
KW = W
KW = MIN(KW,21)
KW = MAX(KW,1)
FUNCKT = KTK(KW,KTYPE) &
& + (W - KW) * (KTK(KW+1,KTYPE) - KTK(KW,KTYPE))
RETURN
END FUNCTION FUNCKT
FUNCTION THSAT(KTYPE) 1,2
!
USE MODULE_GFS_MACHINE , ONLY : kind_phys
USE module_progtm , ONLY : TSAT
implicit none
integer ktype
real(kind=kind_phys) thsat
!
THSAT = TSAT(KTYPE)
RETURN
END FUNCTION THSAT
FUNCTION TWLT(KTYPE),1
USE MODULE_GFS_MACHINE , ONLY : kind_phys
! USE module_progtm
implicit none
integer ktype
real(kind=kind_phys) twlt
!
TWLT = .1
RETURN
END FUNCTION TWLT
END MODULE module_sf_gfs