Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_gfs.F @ 3584

Last change on this file since 3584 was 2759, checked in by aslmd, 3 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 61.5 KB

Rev	Line
[2759]	1	!!WRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_sf_gfs
	4
	5
	6	CONTAINS
	7
	8	!-------------------------------------------------------------------
	9	SUBROUTINE SF_GFS(U3D,V3D,T3D,QV3D,P3D, &
	10	CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
	11	ZNT,UST,PSIM,PSIH, &
	12	XLAND,HFX,QFX,LH,TSK,FLHC,FLQC, &
	13	QGH,QSFC,U10,V10, &
	14	GZ1OZ0,WSPD,BR,ISFFLX, &
	15	EP1,EP2,KARMAN,itimestep, &
	16	ids,ide, jds,jde, kds,kde, &
	17	ims,ime, jms,jme, kms,kme, &
	18	its,ite, jts,jte, kts,kte )
	19	!-------------------------------------------------------------------
	20	USE MODULE_GFS_MACHINE, ONLY : kind_phys
	21	USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
	22	!-------------------------------------------------------------------
	23	IMPLICIT NONE
	24	!-------------------------------------------------------------------
	25	!-- U3D 3D u-velocity interpolated to theta points (m/s)
	26	!-- V3D 3D v-velocity interpolated to theta points (m/s)
	27	!-- T3D temperature (K)
	28	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	29	!-- P3D 3D pressure (Pa)
	30	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
	31	!-- ROVCP R/CP
	32	!-- R gas constant for dry air (J/kg/K)
	33	!-- XLV latent heat of vaporization for water (J/kg)
	34	!-- PSFC surface pressure (Pa)
	35	!-- ZNT roughness length (m)
	36	!-- UST u* in similarity theory (m/s)
	37	!-- PSIM similarity stability function for momentum
	38	!-- PSIH similarity stability function for heat
	39	!-- XLAND land mask (1 for land, 2 for water)
	40	!-- HFX upward heat flux at the surface (W/m^2)
	41	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	42	!-- LH net upward latent heat flux at surface (W/m^2)
	43	!-- TSK surface temperature (K)
	44	!-- FLHC exchange coefficient for heat (m/s)
	45	!-- FLQC exchange coefficient for moisture (m/s)
	46	!-- QGH lowest-level saturated mixing ratio
	47	!-- U10 diagnostic 10m u wind
	48	!-- V10 diagnostic 10m v wind
	49	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
	50	!-- WSPD wind speed at lowest model level (m/s)
	51	!-- BR bulk Richardson number in surface layer
	52	!-- ISFFLX isfflx=1 for surface heat and moisture fluxes
	53	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
	54	!-- KARMAN Von Karman constant
	55	!-- ids start index for i in domain
	56	!-- ide end index for i in domain
	57	!-- jds start index for j in domain
	58	!-- jde end index for j in domain
	59	!-- kds start index for k in domain
	60	!-- kde end index for k in domain
	61	!-- ims start index for i in memory
	62	!-- ime end index for i in memory
	63	!-- jms start index for j in memory
	64	!-- jme end index for j in memory
	65	!-- kms start index for k in memory
	66	!-- kme end index for k in memory
	67	!-- its start index for i in tile
	68	!-- ite end index for i in tile
	69	!-- jts start index for j in tile
	70	!-- jte end index for j in tile
	71	!-- kts start index for k in tile
	72	!-- kte end index for k in tile
	73	!-------------------------------------------------------------------
	74
	75	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
	76	ims,ime, jms,jme, kms,kme, &
	77	its,ite, jts,jte, kts,kte, &
	78	ISFFLX,itimestep
	79
	80	REAL, INTENT(IN) :: &
	81	CP, &
	82	EP1, &
	83	EP2, &
	84	KARMAN, &
	85	R, &
	86	ROVCP, &
	87	XLV
	88
	89	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
	90	P3D, &
	91	QV3D, &
	92	T3D, &
	93	U3D, &
	94	V3D
	95
	96	REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
	97	TSK, &
	98	PSFC, &
	99	XLAND
	100
	101	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
	102	BR, &
	103	CHS, &
	104	CHS2, &
	105	CPM, &
	106	CQS2, &
	107	FLHC, &
	108	FLQC, &
	109	GZ1OZ0, &
	110	HFX, &
	111	LH, &
	112	PSIM, &
	113	PSIH, &
	114	QFX, &
	115	QGH, &
	116	QSFC, &
	117	UST, &
	118	ZNT, &
	119	WSPD
	120
	121	REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
	122	U10, &
	123	V10
	124
	125
	126	!--------------------------- LOCAL VARS ------------------------------
	127
	128	REAL :: ESAT
	129
	130	REAL (kind=kind_phys) :: &
	131	RHOX
	132
	133	REAL (kind=kind_phys), DIMENSION(its:ite) :: &
	134	CH, &
	135	CM, &
	136	DDVEL, &
	137	DRAIN, &
	138	EP, &
	139	EVAP, &
	140	FH, &
	141	FH2, &
	142	FM, &
	143	HFLX, &
	144	PH, &
	145	PM, &
	146	PRSL1, &
	147	PRSLKI, &
	148	PS, &
	149	Q1, &
	150	Q2M, &
	151	QSS, &
	152	QSURF, &
	153	RB, &
	154	RCL, &
	155	RHO1, &
	156	SLIMSK, &
	157	STRESS, &
	158	T1, &
	159	T2M, &
	160	THGB, &
	161	THX, &
	162	TSKIN, &
	163	SHELEG, &
	164	U1, &
	165	U10M, &
	166	USTAR, &
	167	V1, &
	168	V10M, &
	169	WIND, &
	170	Z0RL, &
	171	Z1
	172
	173
	174	INTEGER :: &
	175	I, &
	176	IM, &
	177	J, &
	178	K, &
	179	KM
	180
	181
	182	if(itimestep.eq.0) then
	183	CALL GFUNCPHYS
	184	endif
	185
	186	IM=ITE-ITS+1
	187	KM=KTE-KTS+1
	188
	189	DO J=jts,jte
	190
	191	DO i=its,ite
	192	DDVEL(I)=0.
	193	RCL(i)=1.
	194	PRSL1(i)=P3D(i,kts,j)*.001
	195	PS(i)=PSFC(i,j)*.001
	196	Q1(I) = QV3D(i,kts,j)
	197	! QSURF(I)=QSFC(I,J)
	198	QSURF(I)=0.
	199	SHELEG(I)=0.
	200	SLIMSK(i)=ABS(XLAND(i,j)-2.)
	201	TSKIN(i)=TSK(i,j)
	202	T1(I) = T3D(i,kts,j)
	203	U1(I) = U3D(i,kts,j)
	204	USTAR(I) = UST(i,j)
	205	V1(I) = V3D(i,kts,j)
	206	Z0RL(I) = ZNT(i,j)*100.
	207	ENDDO
	208
	209	DO i=its,ite
	210	PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
	211	THGB(I)=TSKIN(i)(100./PS(I))*ROVCP
	212	THX(I)=T1(i)(100./PRSL1(I))*ROVCP
	213	RHO1(I)=PRSL1(I)1000./(RT1(I)(1.+EP1Q1(I)))
	214	Q1(I)=Q1(I)/(1.+Q1(I))
	215	ENDDO
	216
	217
	218	CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
	219	SHELEG,TSKIN,QSURF, &
	220	!WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, &
	221	!WRF SLRAD,SNOWMT,DELT, &
	222	Z0RL, &
	223	!WRF TG3,GFLUX,F10M, &
	224	U10M,V10M,T2M,Q2M, &
	225	!WRF ZSOIL, &
	226	CM,CH,RB, &
	227	!WRF RHSCNPY,RHSMC,AIM,BIM,CIM, &
	228	RCL,PRSL1,PRSLKI,SLIMSK, &
	229	DRAIN,EVAP,HFLX,STRESS,EP, &
	230	FM,FH,USTAR,WIND,DDVEL, &
	231	PM,PH,FH2,QSS,Z1 )
	232
	233
	234	DO i=its,ite
	235	U10(i,j)=U10M(i)
	236	V10(i,j)=V10M(i)
	237	BR(i,j)=RB(i)
	238	CHS(I,J)=CH(I)*WIND(I)
	239	CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
	240	CPM(I,J)=CP(1.+0.8QV3D(i,kts,j))
	241	esat = fpvs(t1(i))
	242	QGH(I,J)=ep2esat/(1000.ps(i)-esat)
	243	QSFC(I,J)=qss(i)
	244	PSIH(i,j)=PH(i)
	245	PSIM(i,j)=PM(i)
	246	UST(i,j)=ustar(i)
	247	WSPD(i,j)=WIND(i)
	248	ZNT(i,j)=Z0RL(i)*.01
	249	ENDDO
	250
	251	DO i=its,ite
	252	FLHC(i,j)=CPM(I,J)RHO1(I)CHS(I,J)
	253	FLQC(i,j)=RHO1(I)*CHS(I,J)
	254	GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
	255	CQS2(i,j)=CHS2(I,J)
	256	ENDDO
	257
	258	IF (ISFFLX.EQ.0) THEN
	259	DO i=its,ite
	260	HFX(i,j)=0.
	261	LH(i,j)=0.
	262	QFX(i,j)=0.
	263	ENDDO
	264	ELSE
	265	DO i=its,ite
	266	IF(XLAND(I,J)-1.5.GT.0.)THEN
	267	HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
	268	ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
	269	HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
	270	HFX(I,J)=AMAX1(HFX(I,J),-250.)
	271	ENDIF
	272	QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
	273	QFX(I,J)=AMAX1(QFX(I,J),0.)
	274	LH(I,J)=XLV*QFX(I,J)
	275	ENDDO
	276	ENDIF
	277
	278
	279	ENDDO
	280
	281
	282	END SUBROUTINE SF_GFS
	283
	284
	285	!-------------------------------------------------------------------
	286
	287	SUBROUTINE PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
	288	& SHELEG,TSKIN,QSURF, &
	289	!WRF & SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY, &
	290	!WRF & DLWFLX,SLRAD,SNOWMT,DELT, &
	291	& Z0RL, &
	292	!WRF & TG3,GFLUX,F10M, &
	293	& U10M,V10M,T2M,Q2M, &
	294	!WRF & ZSOIL, &
	295	& CM, CH, RB, &
	296	!WRF & RHSCNPY,RHSMC,AIM,BIM,CIM, &
	297	& RCL,PRSL1,PRSLKI,SLIMSK, &
	298	& DRAIN,EVAP,HFLX,STRESS,EP, &
	299	& FM,FH,USTAR,WIND,DDVEL, &
	300	& PM,PH,FH2,QSS,Z1 )
	301	!
	302
	303	USE MODULE_GFS_MACHINE, ONLY : kind_phys
	304	USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
	305	USE MODULE_GFS_PHYSCONS, grav => con_g, SBC => con_sbc, HVAP => con_HVAP &
	306	&, CP => con_CP, HFUS => con_HFUS, JCAL => con_JCAL &
	307	&, EPS => con_eps, EPSM1 => con_epsm1, t0c => con_t0c &
	308	&, RVRDM1 => con_FVirt, RD => con_RD
	309	implicit none
	310	!
	311	! include 'constant.h'
	312	!
	313	integer IM, km
	314	!
	315	real(kind=kind_phys), parameter :: cpinv=1.0/cp, HVAPI=1.0/HVAP
	316	real(kind=kind_phys) DELT
	317	INTEGER SOILTYP(IM), VEGTYPE(IM)
	318	real(kind=kind_phys) PS(IM), U1(IM), V1(IM), &
	319	& T1(IM), Q1(IM), SHELEG(IM), &
	320	& TSKIN(IM), QSURF(IM), SMC(IM,KM), &
	321	& STC(IM,KM), DM(IM), SIGMAF(IM), &
	322	& CANOPY(IM), DLWFLX(IM), SLRAD(IM), &
	323	& SNOWMT(IM), Z0RL(IM), TG3(IM), &
	324	& GFLUX(IM), F10M(IM), U10M(IM), &
	325	& V10M(IM), T2M(IM), Q2M(IM), &
	326	& ZSOIL(IM,KM), CM(IM), CH(IM), &
	327	& RB(IM), RHSCNPY(IM), RHSMC(IM,KM), &
	328	& AIM(IM,KM), BIM(IM,KM), CIM(IM,KM), &
	329	& RCL(IM), PRSL1(IM), PRSLKI(IM), &
	330	& SLIMSK(IM), DRAIN(IM), EVAP(IM), &
	331	& HFLX(IM), RNET(IM), EP(IM), &
	332	& FM(IM), FH(IM), USTAR(IM), &
	333	& WIND(IM), DDVEL(IM), STRESS(IM)
	334	!
	335	! Locals
	336	!
	337	integer k,i
	338	!
	339	real(kind=kind_phys) CANFAC(IM), &
	340	& DDZ(IM), DDZ2(IM), DELTA(IM), &
	341	& DEW(IM), DF1(IM), DFT0(IM), &
	342	& DFT2(IM), DFT1(IM), &
	343	& DMDZ(IM), DMDZ2(IM), DTDZ1(IM), &
	344	& DTDZ2(IM), DTV(IM), EC(IM), &
	345	& EDIR(IM), ETPFAC(IM), &
	346	& FACTSNW(IM), FH2(IM), FM10(IM), &
	347	& FX(IM), GX(IM), &
	348	& HCPCT(IM), HL1(IM), HL12(IM), &
	349	& HLINF(IM), PARTLND(IM), PH(IM), &
	350	& PH2(IM), PM(IM), PM10(IM), &
	351	& PSURF(IM), Q0(IM), QS1(IM), &
	352	& QSS(IM), RAT(IM), RCAP(IM), &
	353	& RCH(IM), RHO(IM), RS(IM), &
	354	& RSMALL(IM), SLWD(IM), SMCZ(IM), &
	355	& SNET(IM), SNOEVP(IM), SNOWD(IM), &
	356	& T1O(IM), T2MO(IM), TERM1(IM), &
	357	& TERM2(IM), THETA1(IM), THV1(IM), &
	358	& TREF(IM), TSURF(IM), TV1(IM), &
	359	& TVS(IM), TSURFO(IM), TWILT(IM), &
	360	& XX(IM), XRCL(IM), YY(IM), &
	361	& Z0(IM), Z0MAX(IM), Z1(IM), &
	362	& ZTMAX(IM), ZZ(IM), PS1(IM)
	363	!
	364	real(kind=kind_phys) a0, a0p, a1, a1p, aa, aa0, &
	365	& aa1, adtv, alpha, arnu, b1, b1p, &
	366	& b2, b2p, bb, bb0, bb1, bb2, &
	367	& bfact, ca, cc, cc1, cc2, cfactr, &
	368	& ch2o, charnock, cice, convrad, cq, csoil, &
	369	& ctfil1,ctfil2, delt2, df2, dfsnow, &
	370	& elocp, eth, ff, FMS, &
	371	!WRF & fhs, funcdf, funckt,g, hl0, hl0inf, &
	372	& fhs, g, hl0, hl0inf, &
	373	& hl110, hlt, hltinf,OLINF, rcq, rcs, &
	374	& rct, restar, rhoh2o,rnu, RSI, &
	375	& rss, scanop, sig2k, sigma, smcdry, &
	376	& t12, t14, tflx, tgice, topt, &
	377	& val, vis, zbot, snomin, tem
	378	!
	379	!
	380
	381	PARAMETER (CHARNOCK=.014,CA=.4)!C CA IS THE VON KARMAN CONSTANT
	382	PARAMETER (G=grav,sigma=sbc)
	383
	384	PARAMETER (ALPHA=5.,A0=-3.975,A1=12.32,B1=-7.755,B2=6.041)
	385	PARAMETER (A0P=-7.941,A1P=24.75,B1P=-8.705,B2P=7.899,VIS=1.4E-5)
	386	PARAMETER (AA1=-1.076,BB1=.7045,CC1=-.05808)
	387	PARAMETER (BB2=-.1954,CC2=.009999)
	388	PARAMETER (ELOCP=HVAP/CP,DFSNOW=.31,CH2O=4.2E6,CSOIL=1.26E6)
	389	PARAMETER (SCANOP=.5,CFACTR=.5,ZBOT=-3.,TGICE=271.2)
	390	PARAMETER (CICE=1880.*917.,topt=298.)
	391	PARAMETER (RHOH2O=1000.,CONVRAD=JCAL*1.E4/60.)
	392	PARAMETER (CTFIL1=.5,CTFIL2=1.-CTFIL1)
	393	PARAMETER (RNU=1.51E-5,ARNU=.135*RNU)
	394	parameter (snomin=1.0e-9)
	395	!
	396	LOGICAL FLAG(IM), FLAGSNW(IM)
	397	!WRF real(kind=kind_phys) KT1(IM), KT2(IM), KTSOIL, &
	398	real(kind=kind_phys) KT1(IM), KT2(IM), &
	399	& ET(IM,KM), &
	400	& STSOIL(IM,KM), AI(IM,KM), BI(IM,KM), &
	401	& CI(IM,KM), RHSTC(IM,KM)
	402	real(kind=kind_phys) rsmax(13), rgl(13), rsmin(13), hs(13), &
	403	& smmax(9), smdry(9), smref(9), smwlt(9)
	404
	405	!
	406	! the 13 vegetation types are:
	407	!
	408	! 1 ... broadleave-evergreen trees (tropical forest)
	409	! 2 ... broadleave-deciduous trees
	410	! 3 ... broadleave and needle leave trees (mixed forest)
	411	! 4 ... needleleave-evergreen trees
	412	! 5 ... needleleave-deciduous trees (larch)
	413	! 6 ... broadleave trees with groundcover (savanna)
	414	! 7 ... groundcover only (perenial)
	415	! 8 ... broadleave shrubs with perenial groundcover
	416	! 9 ... broadleave shrubs with bare soil
	417	! 10 ... dwarf trees and shrubs with ground cover (trunda)
	418	! 11 ... bare soil
	419	! 12 ... cultivations (use parameters from type 7)
	420	! 13 ... glacial
	421	!
	422	data rsmax/13*5000./
	423	data rsmin/150.,100.,125.,150.,100.,70.,40., &
	424	& 300.,400.,150.,999.,40.,999./
	425	data rgl/530.,65.,4100.,999.,100.,999./
	426	data hs/41.69,54.53,51.93,47.35,47.35,54.53,36.35, &
	427	& 3*42.00,999.,36.35,999./
	428	data smmax/.421,.464,.468,.434,.406,.465,.404,.439,.421/
	429	data smdry/.07,.14,.22,.08,.18,.16,.12,.10,.07/
	430	data smref/.283,.387,.412,.312,.338,.382,.315,.329,.283/
	431	data smwlt/.029,.119,.139,.047,.010,.103,.069,.066,.029/
	432	!
	433	!!! save rsmax, rsmin, rgl, hs, smmax, smdry, smref, smwlt
	434	!
	435
	436	!WRF DELT2 = DELT * 2.
	437	!
	438	! ESTIMATE SIGMA ** K AT 2 M
	439	!
	440	SIG2K = 1. - 4. * G * 2. / (CP * 280.)
	441	!
	442	! INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIE
	443	! PSURF IS IN PASCALS
	444	! WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
	445	! RHO IS DENSITY, QS1 IS SAT. HUM. AT LEVEL1 AND QSS IS SAT. HUM. AT
	446	! SURFACE
	447	! CONVERT SLRAD TO THE CIVILIZED UNIT FROM LANGLEY MINUTE-1 K-4
	448	! SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM
	449	!
	450	!!
	451	! qs1 = fpvs(t1)
	452	! qss = fpvs(tskin)
	453	DO I=1,IM
	454	XRCL(I) = SQRT(RCL(I))
	455	PSURF(I) = 1000. * PS(I)
	456	PS1(I) = 1000. * PRSL1(I)
	457	! SLWD(I) = SLRAD(I) * CONVRAD
	458	!WRF SLWD(I) = SLRAD(I)
	459	!
	460	! DLWFLX has been given a negative sign for downward longwave
	461	! snet is the net shortwave flux
	462	!
	463	!WRF SNET(I) = -SLWD(I) - DLWFLX(I)
	464	WIND(I) = XRCL(I) * SQRT(U1(I) * U1(I) + V1(I) * V1(I)) &
	465	& + MAX(0.0_kind_phys, MIN(DDVEL(I), 30.0_kind_phys))
	466	WIND(I) = MAX(WIND(I),1._kind_phys)
	467	Q0(I) = MAX(Q1(I),1.E-8_kind_phys)
	468	TSURF(I) = TSKIN(I)
	469	THETA1(I) = T1(I) * PRSLKI(I)
	470	TV1(I) = T1(I) * (1. + RVRDM1 * Q0(I))
	471	THV1(I) = THETA1(I) * (1. + RVRDM1 * Q0(I))
	472	TVS(I) = TSURF(I) * (1. + RVRDM1 * Q0(I))
	473	RHO(I) = PS1(I) / (RD * TV1(I))
	474	!jfe QS1(I) = 1000. * FPVS(T1(I))
	475	qs1(i) = fpvs(t1(i))
	476	QS1(I) = EPS * QS1(I) / (PS1(I) + EPSM1 * QS1(I))
	477	QS1(I) = MAX(QS1(I), 1.E-8_kind_phys)
	478	Q0(I) = min(QS1(I),Q0(I))
	479	!jfe QSS(I) = 1000. * FPVS(TSURF(I))
	480	qss(i) = fpvs(tskin(i))
	481	QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
	482	! RS = PLANTR
	483	RS(I) = 0.
	484	!WRF if(VEGTYPE(I).gt.0.) RS(I) = rsmin(VEGTYPE(I))
	485	Z0(I) = .01 * Z0RL(i)
	486	!WRF CANOPY(I)= MAX(CANOPY(I),0._kind_phys)
	487	DM(I) = 1.
	488	!WRF
	489	GOTO 1111
	490	!WRF
	491	FACTSNW(I) = 10.
	492	IF(SLIMSK(I).EQ.2.) FACTSNW(I) = 3.
	493	!
	494	! SNOW DEPTH IN WATER EQUIVALENT IS CONVERTED FROM MM TO M UNIT
	495	!
	496	SNOWD(I) = SHELEG(I) / 1000.
	497	FLAGSNW(I) = .FALSE.
	498	!
	499	! WHEN SNOW DEPTH IS LESS THAN 1 MM, A PATCHY SNOW IS ASSUMED AND
	500	! SOIL IS ALLOWED TO INTERACT WITH THE ATMOSPHERE.
	501	! WE SHOULD EVENTUALLY MOVE TO A LINEAR COMBINATION OF SOIL AND
	502	! SNOW UNDER THE CONDITION OF PATCHY SNOW.
	503	!
	504	IF(SNOWD(I).GT..001.OR.SLIMSK(I).EQ.2.) RS(I) = 0.
	505	IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE.
	506	!##DG IF(LAT.EQ.LATD) THEN
	507	!##DG PRINT *, ' WIND,TV1,TVS,Q1,QS1,SNOW,SLIMSK=',
	508	!##DG& WIND,TV1,TVS,Q1,QS1,SNOWD,SLIMSK
	509	!##DG PRINT *, ' SNET, SLWD =', SNET, SLWD(I)
	510	!##DG ENDIF
	511	IF(SLIMSK(I).EQ.0.) THEN
	512	ZSOIL(I,1) = 0.
	513	ELSEIF(SLIMSK(I).EQ.1.) THEN
	514	ZSOIL(I,1) = -.10
	515	ELSE
	516	ZSOIL(I,1) = -3. / KM
	517	ENDIF
	518	!WRF
	519	1111 CONTINUE
	520	!WRF
	521	ENDDO
	522
	523	!!
	524	!WRF
	525	GOTO 2222
	526	!WRF
	527	DO K = 2, KM
	528	DO I=1,IM
	529	IF(SLIMSK(I).EQ.0.) THEN
	530	ZSOIL(I,K) = 0.
	531	ELSEIF(SLIMSK(I).EQ.1.) THEN
	532	ZSOIL(I,K) = ZSOIL(I,K-1) &
	533	& + (-2. - ZSOIL(I,1)) / (KM - 1)
	534	ELSE
	535	ZSOIL(I,K) = - 3. * FLOAT(K) / FLOAT(KM)
	536	ENDIF
	537	ENDDO
	538	ENDDO
	539	!WRF
	540	2222 CONTINUE
	541	!WRF
	542	!!
	543	DO I=1,IM
	544	Z1(I) = -RD * TV1(I) * LOG(PS1(I)/PSURF(I)) / G
	545	DRAIN(I) = 0.
	546	ENDDO
	547
	548	!!
	549	DO K = 1, KM
	550	DO I=1,IM
	551	ET(I,K) = 0.
	552	RHSMC(I,K) = 0.
	553	AIM(I,K) = 0.
	554	BIM(I,K) = 1.
	555	CIM(I,K) = 0.
	556	STSOIL(I,K) = STC(I,K)
	557	ENDDO
	558	ENDDO
	559
	560	DO I=1,IM
	561	EDIR(I) = 0.
	562	EC(I) = 0.
	563	EVAP(I) = 0.
	564	EP(I) = 0.
	565	SNOWMT(I) = 0.
	566	GFLUX(I) = 0.
	567	RHSCNPY(I) = 0.
	568	FX(I) = 0.
	569	ETPFAC(I) = 0.
	570	CANFAC(I) = 0.
	571	ENDDO
	572	!
	573	! COMPUTE STABILITY DEPENDENT EXCHANGE COEFFICIENTS
	574	!
	575	! THIS PORTION OF THE CODE IS PRESENTLY SUPPRESSED
	576	!
	577	DO I=1,IM
	578	IF(SLIMSK(I).EQ.0.) THEN
	579	USTAR(I) = SQRT(G * Z0(I) / CHARNOCK)
	580	ENDIF
	581	!
	582	! COMPUTE STABILITY INDICES (RB AND HLINF)
	583	!
	584
	585	Z0MAX(I) = MIN(Z0(I),0.1 * Z1(I))
	586	ZTMAX(I) = Z0MAX(I)
	587	IF(SLIMSK(I).EQ.0.) THEN
	588	RESTAR = USTAR(I) * Z0MAX(I) / VIS
	589	RESTAR = MAX(RESTAR,.000001_kind_phys)
	590	! RESTAR = ALOG(RESTAR)
	591	! RESTAR = MIN(RESTAR,5.)
	592	! RESTAR = MAX(RESTAR,-5.)
	593	! RAT(I) = AA1 + BB1 * RESTAR + CC1 * RESTAR ** 2
	594	! RAT(I) = RAT(I) / (1. + BB2 * RESTAR
	595	! & + CC2 * RESTAR ** 2)
	596	! Rat taken from Zeng, Zhao and Dickinson 1997
	597	RAT(I) = 2.67 * restar ** .25 - 2.57
	598	RAT(I) = min(RAT(I),7._kind_phys)
	599	ZTMAX(I) = Z0MAX(I) * EXP(-RAT(I))
	600	ENDIF
	601	ENDDO
	602
	603	!##DG IF(LAT.EQ.LATD) THEN
	604	!##DG PRINT *, ' z0max, ztmax, restar, RAT(I) =',
	605	!##DG & z0max, ztmax, restar, RAT(I)
	606	!##DG ENDIF
	607	DO I = 1, IM
	608	DTV(I) = THV1(I) - TVS(I)
	609	ADTV = ABS(DTV(I))
	610	ADTV = MAX(ADTV,.001_kind_phys)
	611	DTV(I) = SIGN(1._kind_phys,DTV(I)) * ADTV
	612	RB(I) = G * DTV(I) * Z1(I) / (.5 * (THV1(I) + TVS(I)) &
	613	& * WIND(I) * WIND(I))
	614	RB(I) = MAX(RB(I),-5000._kind_phys)
	615	! FM(I) = LOG((Z0MAX(I)+Z1(I)) / Z0MAX(I))
	616	! FH(I) = LOG((ZTMAX(I)+Z1(I)) / ZTMAX(I))
	617	FM(I) = LOG((Z1(I)) / Z0MAX(I))
	618	FH(I) = LOG((Z1(I)) / ZTMAX(I))
	619	HLINF(I) = RB(I) * FM(I) * FM(I) / FH(I)
	620	FM10(I) = LOG((Z0MAX(I)+10.) / Z0MAX(I))
	621	FH2(I) = LOG((ZTMAX(I)+2.) / ZTMAX(I))
	622	ENDDO
	623	!##DG IF(LAT.EQ.LATD) THEN
	624	!##DG PRINT *, ' DTV, RB(I), FM(I), FH(I), HLINF =',
	625	!##DG & dtv, rb, FM(I), FH(I), hlinf
	626	!##DG ENDIF
	627	!
	628	! STABLE CASE
	629	!
	630	DO I = 1, IM
	631	IF(DTV(I).GE.0.) THEN
	632	HL1(I) = HLINF(I)
	633	ENDIF
	634	IF(DTV(I).GE.0..AND.HLINF(I).GT..25) THEN
	635	HL0INF = Z0MAX(I) * HLINF(I) / Z1(I)
	636	HLTINF = ZTMAX(I) * HLINF(I) / Z1(I)
	637	AA = SQRT(1. + 4. * ALPHA * HLINF(I))
	638	AA0 = SQRT(1. + 4. * ALPHA * HL0INF)
	639	BB = AA
	640	BB0 = SQRT(1. + 4. * ALPHA * HLTINF)
	641	PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
	642	PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
	643	FMS = FM(I) - PM(I)
	644	FHS = FH(I) - PH(I)
	645	HL1(I) = FMS * FMS * RB(I) / FHS
	646	ENDIF
	647	ENDDO
	648	!
	649	! SECOND ITERATION
	650	!
	651	DO I = 1, IM
	652	IF(DTV(I).GE.0.) THEN
	653	HL0 = Z0MAX(I) * HL1(I) / Z1(I)
	654	HLT = ZTMAX(I) * HL1(I) / Z1(I)
	655	AA = SQRT(1. + 4. * ALPHA * HL1(I))
	656	AA0 = SQRT(1. + 4. * ALPHA * HL0)
	657	BB = AA
	658	BB0 = SQRT(1. + 4. * ALPHA * HLT)
	659	PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
	660	PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
	661	HL110 = HL1(I) * 10. / Z1(I)
	662	AA = SQRT(1. + 4. * ALPHA * HL110)
	663	PM10(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
	664	HL12(I) = HL1(I) * 2. / Z1(I)
	665	! AA = SQRT(1. + 4. * ALPHA * HL12(I))
	666	BB = SQRT(1. + 4. * ALPHA * HL12(I))
	667	PH2(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
	668	ENDIF
	669	ENDDO
	670	!!
	671	!##DG IF(LAT.EQ.LATD) THEN
	672	!##DG PRINT *, ' HL1(I), PM, PH =',
	673	!##DG & HL1(I), pm, ph
	674	!##DG ENDIF
	675	!
	676	! UNSTABLE CASE
	677	!
	678	!
	679	! CHECK FOR UNPHYSICAL OBUKHOV LENGTH
	680	!
	681	DO I=1,IM
	682	IF(DTV(I).LT.0.) THEN
	683	OLINF = Z1(I) / HLINF(I)
	684	IF(ABS(OLINF).LE.50. * Z0MAX(I)) THEN
	685	HLINF(I) = -Z1(I) / (50. * Z0MAX(I))
	686	ENDIF
	687	ENDIF
	688	ENDDO
	689	!
	690	! GET PM AND PH
	691	!
	692	DO I = 1, IM
	693	IF(DTV(I).LT.0..AND.HLINF(I).GE.-.5) THEN
	694	HL1(I) = HLINF(I)
	695	PM(I) = (A0 + A1 * HL1(I)) * HL1(I) &
	696	& / (1. + B1 * HL1(I) + B2 * HL1(I) * HL1(I))
	697	PH(I) = (A0P + A1P * HL1(I)) * HL1(I) &
	698	& / (1. + B1P * HL1(I) + B2P * HL1(I) * HL1(I))
	699	HL110 = HL1(I) * 10. / Z1(I)
	700	PM10(I) = (A0 + A1 * HL110) * HL110 &
	701	& / (1. + B1 * HL110 + B2 * HL110 * HL110)
	702	HL12(I) = HL1(I) * 2. / Z1(I)
	703	PH2(I) = (A0P + A1P * HL12(I)) * HL12(I) &
	704	& / (1. + B1P * HL12(I) + B2P * HL12(I) * HL12(I))
	705	ENDIF
	706	IF(DTV(I).LT.0.AND.HLINF(I).LT.-.5) THEN
	707	HL1(I) = -HLINF(I)
	708	PM(I) = LOG(HL1(I)) + 2. * HL1(I) ** (-.25) - .8776
	709	PH(I) = LOG(HL1(I)) + .5 * HL1(I) ** (-.5) + 1.386
	710	HL110 = HL1(I) * 10. / Z1(I)
	711	PM10(I) = LOG(HL110) + 2. * HL110 ** (-.25) - .8776
	712	HL12(I) = HL1(I) * 2. / Z1(I)
	713	PH2(I) = LOG(HL12(I)) + .5 * HL12(I) ** (-.5) + 1.386
	714	ENDIF
	715	ENDDO
	716	!
	717	! FINISH THE EXCHANGE COEFFICIENT COMPUTATION TO PROVIDE FM AND FH
	718	!
	719	DO I = 1, IM
	720
	721	FM(I) = FM(I) - PM(I)
	722	FH(I) = FH(I) - PH(I)
	723	FM10(I) = FM10(I) - PM10(I)
	724	FH2(I) = FH2(I) - PH2(I)
	725	CM(I) = CA * CA / (FM(I) * FM(I))
	726	CH(I) = CA * CA / (FM(I) * FH(I))
	727	CQ = CH(I)
	728	STRESS(I) = CM(I) * WIND(I) * WIND(I)
	729	USTAR(I) = SQRT(STRESS(I))
	730	! USTAR(I) = SQRT(CM(I) * WIND(I) * WIND(I))
	731	ENDDO
	732	!##DG IF(LAT.EQ.LATD) THEN
	733	!##DG PRINT *, ' FM, FH, CM, CH(I), USTAR =',
	734	!##DG & FM, FH, CM, ch, USTAR
	735	!##DG ENDIF
	736	!
	737	! UPDATE Z0 OVER OCEAN
	738	!
	739	DO I = 1, IM
	740	IF(SLIMSK(I).EQ.0.) THEN
	741	Z0(I) = (CHARNOCK / G) * USTAR(I) ** 2
	742	! NEW IMPLEMENTATION OF Z0
	743	! CC = USTAR(I) * Z0 / RNU
	744	! PP = CC / (1. + CC)
	745	! FF = G * ARNU / (CHARNOCK * USTAR(I) ** 3)
	746	! Z0 = ARNU / (USTAR(I) * FF ** PP)
	747	Z0(I) = MIN(Z0(I),.1_kind_phys)
	748	Z0(I) = MAX(Z0(I),1.E-7_kind_phys)
	749	Z0RL(I) = 100. * Z0(I)
	750	ENDIF
	751	ENDDO
	752
	753	GOTO 5555
	754	!
	755	! RCP = RHO CP CH V
	756	!
	757	DO I = 1, IM
	758	RCH(I) = RHO(I) * CP * CH(I) * WIND(I)
	759	ENDDO
	760
	761
	762	!
	763	! SENSIBLE AND LATENT HEAT FLUX OVER OPEN WATER
	764	!
	765	DO I = 1, IM
	766	IF(SLIMSK(I).EQ.0.) THEN
	767	EVAP(I) = ELOCP * RCH(I) * (QSS(I) - Q1(I))
	768	DM(I) = 1.
	769	QSURF(I) = QSS(I)
	770	ENDIF
	771	ENDDO
	772
	773	!
	774	! COMPUTE SOIL/SNOW/ICE HEAT FLUX IN PREPARATION FOR SURFACE ENERGY
	775	! BALANCE CALCULATION
	776	!
	777	DO I = 1, IM
	778	GFLUX(I) = 0.
	779	IF(SLIMSK(I).EQ.1.) THEN
	780	SMCZ(I) = .5 * (SMC(I,1) + .20)
	781	DFT0(I) = KTSOIL(SMCZ(I),SOILTYP(I))
	782	ELSEIF(SLIMSK(I).EQ.2.) THEN
	783	! DF FOR ICE IS TAKEN FROM MAYKUT AND UNTERSTEINER
	784	! DF IS IN SI UNIT OF W K-1 M-1
	785	DFT0(I) = 2.2
	786	ENDIF
	787	ENDDO
	788	!!
	789	DO I=1,IM
	790	IF(SLIMSK(I).NE.0.) THEN
	791	! IF(SNOWD(I).GT..001) THEN
	792	IF(FLAGSNW(I)) THEN
	793	!
	794	! WHEN SNOW COVERED, GROUND HEAT FLUX COMES FROM SNOW
	795	!
	796	TFLX = MIN(T1(I), TSURF(I))
	797	GFLUX(I) = -DFSNOW * (TFLX - STSOIL(I,1)) &
	798	& / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
	799	ELSE
	800	GFLUX(I) = DFT0(I) * (STSOIL(I,1) - TSURF(I)) &
	801	& / (-.5 * ZSOIL(I,1))
	802	ENDIF
	803	GFLUX(I) = MAX(GFLUX(I),-200._kind_phys)
	804	GFLUX(I) = MIN(GFLUX(I),+200._kind_phys)
	805	ENDIF
	806	ENDDO
	807	DO I = 1, IM
	808	FLAG(I) = SLIMSK(I).NE.0.
	809	PARTLND(I) = 1.
	810	IF(SNOWD(I).GT.0..AND.SNOWD(I).LE..001) THEN
	811	PARTLND(I) = 1. - SNOWD(I) / .001
	812	ENDIF
	813	ENDDO
	814	DO I = 1, IM
	815	SNOEVP(I) = 0.
	816	if(SNOWD(I).gt..001) PARTLND(I) = 0.
	817	ENDDO
	818	!
	819	! COMPUTE POTENTIAL EVAPORATION FOR LAND AND SEA ICE
	820	!
	821	DO I = 1, IM
	822	IF(FLAG(I)) THEN
	823	T12 = T1(I) * T1(I)
	824	T14 = T12 * T12
	825	!
	826	! RCAP = FNET - SIGMA T**4 + GFLX - RHO CP CH V (T1-THETA1)
	827	!
	828	RCAP(I) = -SLWD(I) - SIGMA * T14 + GFLUX(I) &
	829	& - RCH(I) * (T1(I) - THETA1(I))
	830	!
	831	! RSMALL = 4 SIGMA T**3 / RCH(I) + 1
	832	!
	833	RSMALL(I) = 4. * SIGMA * T1(I) * T12 / RCH(I) + 1.
	834	!
	835	! DELTA = L / CP * DQS/DT
	836	!
	837	DELTA(I) = ELOCP * EPS * HVAP * QS1(I) / (RD * T12)
	838	!
	839	! POTENTIAL EVAPOTRANSPIRATION ( WATTS / M**2 ) AND
	840	! POTENTIAL EVAPORATION
	841	!
	842	TERM1(I) = ELOCP * RSMALL(I) * RCH(I)*(QS1(I)-Q0(I))
	843	TERM2(I) = RCAP(I) * DELTA(I)
	844	EP(I) = (ELOCP * RSMALL(I) * RCH(I) * (QS1(I) - Q0(I)) &
	845	& + RCAP(I) * DELTA(I))
	846	EP(I) = EP(I) / (RSMALL(I) + DELTA(I))
	847	ENDIF
	848	ENDDO
	849	!
	850	! ACTUAL EVAPORATION OVER LAND IN THREE PARTS : EDIR, ET, AND EC
	851	!
	852	! DIRECT EVAPORATION FROM SOIL, THE UNIT GOES FROM M S-1 TO KG M-2 S-1
	853	!
	854	DO I = 1, IM
	855	FLAG(I) = SLIMSK(I).EQ.1..AND.EP(I).GT.0.
	856	ENDDO
	857	DO I = 1, IM
	858	IF(FLAG(I)) THEN
	859	DF1(I) = FUNCDF(SMC(I,1),SOILTYP(I))
	860	KT1(I) = FUNCKT(SMC(I,1),SOILTYP(I))
	861	endif
	862	if(FLAG(I).and.STC(I,1).lt.t0c) then
	863	DF1(I) = 0.
	864	KT1(I) = 0.
	865	endif
	866	IF(FLAG(I)) THEN
	867	! TREF = .75 * THSAT(SOILTYP(I))
	868	TREF(I) = smref(SOILTYP(I))
	869	! TWILT = TWLT(SOILTYP(I))
	870	TWILT(I) = smwlt(SOILTYP(I))
	871	smcdry = smdry(SOILTYP(I))
	872	! FX(I) = -2. * DF1(I) * (SMC(I,1) - .23) / ZSOIL(I,1)
	873	! & - KT1(I)
	874	FX(I) = -2. * DF1(I) * (SMC(I,1) - smcdry) / ZSOIL(I,1) &
	875	& - KT1(I)
	876	FX(I) = MIN(FX(I), EP(I)/HVAP)
	877	FX(I) = MAX(FX(I),0._kind_phys)
	878	!
	879	! SIGMAF IS THE FRACTION OF AREA COVERED BY VEGETATION
	880	!
	881	EDIR(I) = FX(I) * (1. - SIGMAF(I)) * PARTLND(I)
	882	ENDIF
	883	ENDDO
	884	!
	885	! calculate stomatal resistance
	886	!
	887	DO I = 1, IM
	888	if(FLAG(I)) then
	889	!
	890	! resistance due to PAR. We use net solar flux as proxy at the present time
	891	!
	892	ff = .55 * 2. * SNET(I) / rgl(VEGTYPE(I))
	893	rcs = (ff + RS(I)/rsmax(VEGTYPE(I))) / (1. + ff)
	894	rcs = max(rcs,.0001_kind_phys)
	895	rct = 1.
	896	rcq = 1.
	897	!
	898	! resistance due to thermal effect
	899	!
	900	! rct = 1. - .0016 * (topt - theta1) ** 2
	901	! rct = max(rct,.0001)
	902	!
	903	! resistance due to humidity
	904	!
	905	! rcq = 1. / (1. + hs(VEGTYPE(I)) * (QS1(I) - Q0(I)))
	906	! rcq = max(rcq,.0001)
	907	!
	908	! compute resistance without the effect of soil moisture
	909	!
	910	RS(I) = RS(I) / (rcs * rct * rcq)
	911	endif
	912	ENDDO
	913	!
	914	! TRANSPIRATION FROM ALL LEVELS OF THE SOIL
	915	!
	916	DO I = 1, IM
	917	IF(FLAG(I)) THEN
	918	CANFAC(I) = (CANOPY(I) / SCANOP) ** CFACTR
	919	endif
	920	IF(FLAG(I)) THEN
	921	ETPFAC(I) = SIGMAF(I) &
	922	& * (1. - CANFAC(I)) / HVAP
	923	GX(I) = (SMC(I,1) - TWILT(I)) / (TREF(I) - TWILT(I))
	924	GX(I) = MAX(GX(I),0._kind_phys)
	925	GX(I) = MIN(GX(I),1._kind_phys)
	926	!
	927	! resistance due to soil moisture deficit
	928	!
	929	rss = GX(I) * (ZSOIL(I,1) / ZSOIL(I,km))
	930	rss = max(rss,.0001_kind_phys)
	931	RSI = RS(I) / rss
	932	!
	933	! transpiration a la Monteith
	934	!
	935	eth = (TERM1(I) + TERM2(I)) / &
	936	& (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
	937	ET(I,1) = ETPFAC(I) * eth &
	938	& * PARTLND(I)
	939	ENDIF
	940	ENDDO
	941	!!
	942	DO K = 2, KM
	943	DO I=1,IM
	944	IF(FLAG(I)) THEN
	945	GX(I) = (SMC(I,K) - TWILT(I)) / (TREF(I) - TWILT(I))
	946	GX(I) = MAX(GX(I),0._kind_phys)
	947	GX(I) = MIN(GX(I),1._kind_phys)
	948	!
	949	! resistance due to soil moisture deficit
	950	!
	951	rss = GX(I) * ((ZSOIL(I,k) - ZSOIL(I,k-1))/ZSOIL(I,km))
	952	rss = max(rss,1.e-6_kind_phys)
	953	RSI = RS(I) / rss
	954	!
	955	! transpiration a la Monteith
	956	!
	957	eth = (TERM1(I) + TERM2(I)) / &
	958	& (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
	959	ET(I,K) = eth &
	960	& * ETPFAC(I) * PARTLND(I)
	961	ENDIF
	962	ENDDO
	963	ENDDO
	964	!!
	965	400 CONTINUE
	966	!
	967	! CANOPY RE-EVAPORATION
	968	!
	969	DO I=1,IM
	970	IF(FLAG(I)) THEN
	971	EC(I) = SIGMAF(I) * CANFAC(I) * EP(I) / HVAP
	972	EC(I) = EC(I) * PARTLND(I)
	973	EC(I) = min(EC(I),CANOPY(I)/delt)
	974	ENDIF
	975	ENDDO
	976	!
	977	! SUM UP TOTAL EVAPORATION
	978	!
	979	DO I = 1, IM
	980	IF(FLAG(I)) THEN
	981	EVAP(I) = EDIR(I) + EC(I)
	982	ENDIF
	983	ENDDO
	984	!!
	985	DO K = 1, KM
	986	DO I=1,IM
	987	IF(FLAG(I)) THEN
	988	EVAP(I) = EVAP(I) + ET(I,K)
	989	ENDIF
	990	ENDDO
	991	ENDDO
	992	!!
	993	!
	994	! RETURN EVAP UNIT FROM KG M-2 S-1 TO WATTS M-2
	995	!
	996	DO I=1,IM
	997	IF(FLAG(I)) THEN
	998	EVAP(I) = MIN(EVAP(I)*HVAP,EP(I))
	999	ENDIF
	1000	ENDDO
	1001	!##DG IF(LAT.EQ.LATD) THEN
	1002	!##DG PRINT , 'FX(I), SIGMAF, EDIR(I), ETPFAC=', FX(I)HVAP,SIGMAF,
	1003	!##DG& EDIR(I)HVAP,ETPFACHVAP
	1004	!##DG PRINT , ' ET =', (ET(K)HVAP,K=1,KM)
	1005	!##DG PRINT , ' CANFAC(I), EC(I), EVAP', CANFAC(I),EC(I)HVAP,EVAP
	1006	!##DG ENDIF
	1007	!
	1008	! EVAPORATION OVER BARE SEA ICE
	1009	!
	1010	DO I = 1, IM
	1011	! IF(SLIMSK(I).EQ.2.AND.SNOWD(I).LE..001) THEN
	1012	IF(SLIMSK(I).EQ.2.) THEN
	1013	EVAP(I) = PARTLND(I) * EP(I)
	1014	ENDIF
	1015	ENDDO
	1016	!
	1017	! TREAT DOWNWARD MOISTURE FLUX SITUATION
	1018	! (EVAP WAS PRESET TO ZERO SO NO UPDATE NEEDED)
	1019	! DEW IS CONVERTED FROM KG M-2 TO M TO CONFORM TO PRECIP UNIT
	1020	!
	1021	DO I = 1, IM
	1022	FLAG(I) = SLIMSK(I).NE.0..AND.EP(I).LE.0.
	1023	DEW(I) = 0.
	1024	ENDDO
	1025	DO I = 1, IM
	1026	IF(FLAG(I)) THEN
	1027	DEW(I) = -EP(I) * DELT / (HVAP * RHOH2O)
	1028	EVAP(I) = EP(I)
	1029	DEW(I) = DEW(I) * PARTLND(I)
	1030	EVAP(I) = EVAP(I) * PARTLND(I)
	1031	DM(I) = 1.
	1032	ENDIF
	1033	ENDDO
	1034	!
	1035	! SNOW COVERED LAND AND SEA ICE
	1036	!
	1037	DO I = 1, IM
	1038	FLAG(I) = SLIMSK(I).NE.0..AND.SNOWD(I).GT.0.
	1039	ENDDO
	1040	!
	1041	! CHANGE OF SNOW DEPTH DUE TO EVAPORATION OR SUBLIMATION
	1042	!
	1043	! CONVERT EVAP FROM KG M-2 S-1 TO M S-1 TO DETERMINE THE REDUCTION OF S
	1044	!
	1045	DO I = 1, IM
	1046	IF(FLAG(I)) THEN
	1047	BFACT = SNOWD(I) / (DELT * EP(I) / (HVAP * RHOH2O))
	1048	BFACT = MIN(BFACT,1._kind_phys)
	1049	!
	1050	! THE EVAPORATION OF SNOW
	1051	!
	1052	IF(EP(I).LE.0.) BFACT = 1.
	1053	IF(SNOWD(I).LE..001) THEN
	1054	! EVAP = (SNOWD(I)/.001)BFACTEP(I) + EVAP
	1055	! SNOEVP(I) = bfact * EP(I) * (1. - PARTLND(I))
	1056	! EVAP = EVAP + SNOEVP(I)
	1057	SNOEVP(I) = bfact * EP(I)
	1058	! EVAP = EVAP + SNOEVP(I) * (1. - PARTLND(I))
	1059	EVAP(I)=EVAP(I)+SNOEVP(I)*(1.-PARTLND(I))
	1060	ELSE
	1061	! EVAP(I) = BFACT * EP(I)
	1062	SNOEVP(I) = bfact * EP(I)
	1063	EVAP(I) = SNOEVP(I)
	1064	ENDIF
	1065	TSURF(I) = T1(I) + &
	1066	& (RCAP(I) - GFLUX(I) - DFSNOW * (T1(I) - STSOIL(I,1)) &
	1067	& /(FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys)) &
	1068	! & + THETA1 - T1 &
	1069	! & - BFACT * EP(I)) / (RSMALL(I) * RCH(I) &
	1070	& - SNOEVP(I)) / (RSMALL(I) * RCH(I) &
	1071	& + DFSNOW / (FACTSNW(I)* MAX(SNOWD(I),.001_kind_phys)))
	1072	! SNOWD(I) = SNOWD(I) - BFACT * EP(I) * DELT / (RHOH2O * HVAP)
	1073	SNOWD(I) = SNOWD(I) - SNOEVP(I) * delt / (rhoh2o * hvap)
	1074	SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
	1075	ENDIF
	1076	ENDDO
	1077	!
	1078	! SNOW MELT (M)
	1079	!
	1080	500 CONTINUE
	1081	DO I = 1, IM
	1082	FLAG(I) = SLIMSK(I).NE.0. &
	1083	& .AND.SNOWD(I).GT..0
	1084	ENDDO
	1085	DO I = 1, IM
	1086	IF(FLAG(I).AND.TSURF(I).GT.T0C) THEN
	1087	SNOWMT(I) = RCH(I) * RSMALL(I) * DELT &
	1088	& * (TSURF(I) - T0C) / (RHOH2O * HFUS)
	1089	SNOWMT(I) = min(SNOWMT(I),SNOWD(I))
	1090	SNOWD(I) = SNOWD(I) - SNOWMT(I)
	1091	SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
	1092	TSURF(I) = MAX(T0C,TSURF(I) &
	1093	& -HFUSSNOWMT(I)RHOH2O/(RCH(I)RSMALL(I)DELT))
	1094	ENDIF
	1095	ENDDO
	1096	!
	1097	! We need to re-evaluate evaporation because of snow melt
	1098	! the skin temperature is now bounded to 0 deg C
	1099	!
	1100	! qss = fpvs(tsurf)
	1101	DO I = 1, IM
	1102	! IF (SNOWD(I) .GT. 0.0) THEN
	1103	IF (SNOWD(I) .GT. snomin) THEN
	1104	!jfe QSS(I) = 1000. * FPVS(TSURF(I))
	1105	qss(i) = fpvs(tsurf(i))
	1106	QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
	1107	EVAP(I) = elocp * RCH(I) * (QSS(I) - Q0(I))
	1108	ENDIF
	1109	ENDDO
	1110	!
	1111	! PREPARE TENDENCY TERMS FOR THE SOIL MOISTURE FIELD WITHOUT PRECIPITAT
	1112	! THE UNIT OF MOISTURE FLUX NEEDS TO BECOME M S-1 FOR SOIL MOISTURE
	1113	! HENCE THE FACTOR OF RHOH2O
	1114	!
	1115	DO I = 1, IM
	1116	FLAG(I) = SLIMSK(I).EQ.1.
	1117	if(FLAG(I)) then
	1118	DF1(I) = FUNCDF(SMCZ(I),SOILTYP(I))
	1119	KT1(I) = FUNCKT(SMCZ(I),SOILTYP(I))
	1120	endif
	1121	if(FLAG(I).and.STC(I,1).lt.t0c) then
	1122	DF1(I) = 0.
	1123	KT1(I) = 0.
	1124	endif
	1125	IF(FLAG(I)) THEN
	1126	RHSCNPY(I) = -EC(I) + SIGMAF(I) * RHOH2O * DEW(I) / DELT
	1127	SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
	1128	DMDZ(I) = (SMC(I,1) - SMC(I,2)) / (-.5 * ZSOIL(I,2))
	1129	RHSMC(I,1) = (DF1(I) * DMDZ(I) + KT1(I) &
	1130	& + (EDIR(I) + ET(I,1))) / (ZSOIL(I,1) * RHOH2O)
	1131	RHSMC(I,1) = RHSMC(I,1) - (1. - SIGMAF(I)) * DEW(I) / &
	1132	& ( ZSOIL(I,1) * delt)
	1133	DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
	1134	!
	1135	! AIM, BIM, AND CIM ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
	1136	! IMPLICIT UPDATE OF THE SOIL MOISTURE
	1137	!
	1138	AIM(I,1) = 0.
	1139	BIM(I,1) = DF1(I) * DDZ(I) / (-ZSOIL(I,1) * RHOH2O)
	1140	CIM(I,1) = -BIM(I,1)
	1141	ENDIF
	1142	ENDDO
	1143	!!
	1144	DO K = 2, KM
	1145	IF(K.LT.KM) THEN
	1146	DO I=1,IM
	1147	IF(FLAG(I)) THEN
	1148	DF2 = FUNCDF(SMCZ(I),SOILTYP(I))
	1149	KT2(I) = FUNCKT(SMCZ(I),SOILTYP(I))
	1150	ENDIF
	1151	IF(FLAG(I).and.STC(I,k).lt.t0c) THEN
	1152	df2 = 0.
	1153	KT2(I) = 0.
	1154	ENDIF
	1155	IF(FLAG(I)) THEN
	1156	DMDZ2(I) = (SMC(I,K) - SMC(I,K+1)) &
	1157	& / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
	1158	SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
	1159	RHSMC(I,K) = (DF2 * DMDZ2(I) + KT2(I) &
	1160	& - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K)) &
	1161	& / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
	1162	DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
	1163	CIM(I,K) = -DF2 * DDZ2(I) &
	1164	& / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
	1165	ENDIF
	1166	ENDDO
	1167	ELSE
	1168	DO I = 1, IM
	1169	IF(FLAG(I)) THEN
	1170	KT2(I) = FUNCKT(SMC(I,K),SOILTYP(I))
	1171	ENDIF
	1172	if(FLAG(I).and.STC(I,k).lt.t0c) KT2(I) = 0.
	1173	IF(FLAG(I)) THEN
	1174	RHSMC(I,K) = (KT2(I) &
	1175	& - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K)) &
	1176	& / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
	1177	DRAIN(I) = KT2(I)
	1178	CIM(I,K) = 0.
	1179	ENDIF
	1180	ENDDO
	1181	ENDIF
	1182	DO I = 1, IM
	1183	IF(FLAG(I)) THEN
	1184	AIM(I,K) = -DF1(I) * DDZ(I) &
	1185	& / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
	1186	BIM(I,K) = -(AIM(I,K) + CIM(I,K))
	1187	DF1(I) = DF2
	1188	KT1(I) = KT2(I)
	1189	DMDZ(I) = DMDZ2(I)
	1190	DDZ(I) = DDZ2(I)
	1191	ENDIF
	1192	ENDDO
	1193	ENDDO
	1194	!!
	1195	600 CONTINUE
	1196	!
	1197	! UPDATE SOIL TEMPERATURE AND SEA ICE TEMPERATURE
	1198	!
	1199	DO I=1,IM
	1200	FLAG(I) = SLIMSK(I).NE.0.
	1201	ENDDO
	1202	!
	1203	! SURFACE TEMPERATURE IS PART OF THE UPDATE WHEN SNOW IS ABSENT
	1204	!
	1205	DO I=1,IM
	1206	! IF(FLAG(I).AND.SNOWD(I).LE..001) THEN
	1207	IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
	1208	YY(I) = T1(I) + &
	1209	! & (RCAP(I)-GFLUX(I) + THETA1 - T1(I) &
	1210	& (RCAP(I)-GFLUX(I) &
	1211	& - EVAP(I)) / (RSMALL(I) * RCH(I))
	1212	ZZ(I) = 1. + DFT0(I) / (-.5 * ZSOIL(I,1) * RCH(I) * RSMALL(I))
	1213	XX(I) = DFT0(I) * (STSOIL(I,1) - YY(I)) / &
	1214	& (.5 * ZSOIL(I,1) * ZZ(I))
	1215	ENDIF
	1216	! IF(FLAG(I).AND.SNOWD(I).GT..001) THEN
	1217	IF(FLAG(I).AND.FLAGSNW(I)) THEN
	1218	YY(I) = STSOIL(I,1)
	1219	!
	1220	! HEAT FLUX FROM SNOW IS EXPLICIT IN TIME
	1221	!
	1222	ZZ(I) = 1.
	1223	XX(I) = DFSNOW * (STSOIL(I,1) - TSURF(I)) &
	1224	& / (-FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
	1225	ENDIF
	1226	ENDDO
	1227	!
	1228	! COMPUTE THE FORCING AND THE IMPLICIT MATRIX ELEMENTS FOR UPDATE
	1229	!
	1230	! CH2O IS THE HEAT CAPACITY OF WATER AND CSOIL IS THE HEAT CAPACITY OF
	1231	!
	1232	DO I = 1, IM
	1233	IF(FLAG(I)) THEN
	1234	SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
	1235	DTDZ1(I) = (STSOIL(I,1) - STSOIL(I,2)) / (-.5 * ZSOIL(I,2))
	1236	IF(SLIMSK(I).EQ.1.) THEN
	1237	DFT1(I) = KTSOIL(SMCZ(I),SOILTYP(I))
	1238	HCPCT(I) = SMC(I,1) * CH2O + (1. - SMC(I,1)) * CSOIL
	1239	ELSE
	1240	DFT1(I) = DFT0(I)
	1241	HCPCT(I) = CICE
	1242	ENDIF
	1243	DFT2(I) = DFT1(I)
	1244	DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
	1245	!
	1246	! AI, BI, AND CI ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
	1247	! IMPLICIT UPDATE OF THE SOIL TEMPERATURE
	1248	!
	1249	AI(I,1) = 0.
	1250	BI(I,1) = DFT1(I) * DDZ(I) / (-ZSOIL(I,1) * HCPCT(I))
	1251	CI(I,1) = -BI(I,1)
	1252	BI(I,1) = BI(I,1) &
	1253	& + DFT0(I) / (.5 * ZSOIL(I,1) *2 HCPCT(I) * ZZ(I))
	1254	! SS = DFT0(I) * (STSOIL(I,1) - YY(I)) &
	1255	! & / (.5 * ZSOIL(I,1) * ZZ(I))
	1256	! RHSTC(1) = (DFT1(I) * DTDZ1(I) - SS)
	1257	RHSTC(I,1) = (DFT1(I) * DTDZ1(I) - XX(I)) &
	1258	& / (ZSOIL(I,1) * HCPCT(I))
	1259	ENDIF
	1260	ENDDO
	1261	!!
	1262	DO K = 2, KM
	1263	DO I=1,IM
	1264	IF(SLIMSK(I).EQ.1.) THEN
	1265	HCPCT(I) = SMC(I,K) * CH2O + (1. - SMC(I,K)) * CSOIL
	1266	ELSEIF(SLIMSK(I).EQ.2.) THEN
	1267	HCPCT(I) = CICE
	1268	ENDIF
	1269	ENDDO
	1270	IF(K.LT.KM) THEN
	1271	DO I = 1, IM
	1272	IF(FLAG(I)) THEN
	1273	DTDZ2(I) = (STSOIL(I,K) - STSOIL(I,K+1)) &
	1274	& / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
	1275	SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
	1276	IF(SLIMSK(I).EQ.1.) THEN
	1277	DFT2(I) = KTSOIL(SMCZ(I),SOILTYP(I))
	1278	ENDIF
	1279	DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
	1280	CI(I,K) = -DFT2(I) * DDZ2(I) &
	1281	& / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
	1282	ENDIF
	1283	ENDDO
	1284	ELSE
	1285	!
	1286	! AT THE BOTTOM, CLIMATOLOGY IS ASSUMED AT 2M DEPTH FOR LAND AND
	1287	! FREEZING TEMPERATURE IS ASSUMED FOR SEA ICE AT Z(KM)
	1288	DO I = 1, IM
	1289	IF(SLIMSK(I).EQ.1.) THEN
	1290	DTDZ2(I) = (STSOIL(I,K) - TG3(I)) &
	1291	& / (.5 * (ZSOIL(I,K-1) + ZSOIL(I,K)) - ZBOT)
	1292	DFT2(I) = KTSOIL(SMC(I,K),SOILTYP(I))
	1293	CI(I,K) = 0.
	1294	ENDIF
	1295	IF(SLIMSK(I).EQ.2.) THEN
	1296	DTDZ2(I) = (STSOIL(I,K) - TGICE) &
	1297	& / (.5 * ZSOIL(I,K-1) - .5 * ZSOIL(I,K))
	1298	DFT2(I) = DFT1(I)
	1299	CI(I,K) = 0.
	1300	ENDIF
	1301	ENDDO
	1302	ENDIF
	1303	DO I = 1, IM
	1304	IF(FLAG(I)) THEN
	1305	RHSTC(I,K) = (DFT2(I) * DTDZ2(I) - DFT1(I) * DTDZ1(I)) &
	1306	& / ((ZSOIL(I,K) - ZSOIL(I,K-1)) * HCPCT(I))
	1307	AI(I,K) = -DFT1(I) * DDZ(I) &
	1308	& / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
	1309	BI(I,K) = -(AI(I,K) + CI(I,K))
	1310	DFT1(I) = DFT2(I)
	1311	DTDZ1(I) = DTDZ2(I)
	1312	DDZ(I) = DDZ2(I)
	1313	ENDIF
	1314	ENDDO
	1315	ENDDO
	1316	!!
	1317	700 CONTINUE
	1318	!
	1319	! SOLVE THE TRI-DIAGONAL MATRIX
	1320	!
	1321	DO K = 1, KM
	1322	DO I=1,IM
	1323	IF(FLAG(I)) THEN
	1324	RHSTC(I,K) = RHSTC(I,K) * DELT2
	1325	AI(I,K) = AI(I,K) * DELT2
	1326	BI(I,K) = 1. + BI(I,K) * DELT2
	1327	CI(I,K) = CI(I,K) * DELT2
	1328	ENDIF
	1329	ENDDO
	1330	ENDDO
	1331	! FORWARD ELIMINATION
	1332	DO I=1,IM
	1333	IF(FLAG(I)) THEN
	1334	CI(I,1) = -CI(I,1) / BI(I,1)
	1335	RHSTC(I,1) = RHSTC(I,1) / BI(I,1)
	1336	ENDIF
	1337	ENDDO
	1338	!!
	1339	DO K = 2, KM
	1340	DO I=1,IM
	1341	IF(FLAG(I)) THEN
	1342	CC = 1. / (BI(I,K) + AI(I,K) * CI(I,K-1))
	1343	CI(I,K) = -CI(I,K) * CC
	1344	RHSTC(I,K) = (RHSTC(I,K) - AI(I,K) * RHSTC(I,K-1)) * CC
	1345	ENDIF
	1346	ENDDO
	1347	ENDDO
	1348	!!
	1349	! BACKWARD SUBSTITUTTION
	1350	DO I=1,IM
	1351	IF(FLAG(I)) THEN
	1352	CI(I,KM) = RHSTC(I,KM)
	1353	ENDIF
	1354	ENDDO
	1355	!!
	1356	DO K = KM-1, 1
	1357	DO I=1,IM
	1358	IF(FLAG(I)) THEN
	1359	CI(I,K) = CI(I,K) * CI(I,K+1) + RHSTC(I,K)
	1360	ENDIF
	1361	ENDDO
	1362	ENDDO
	1363	!
	1364	! UPDATE SOIL AND ICE TEMPERATURE
	1365	!
	1366	DO K = 1, KM
	1367	DO I=1,IM
	1368	IF(FLAG(I)) THEN
	1369	STSOIL(I,K) = STSOIL(I,K) + CI(I,K)
	1370	ENDIF
	1371	ENDDO
	1372	ENDDO
	1373	!
	1374	! UPDATE SURFACE TEMPERATURE FOR SNOW FREE SURFACES
	1375	!
	1376	DO I=1,IM
	1377	! IF(SLIMSK(I).NE.0..AND.SNOWD(I).LE..001) THEN
	1378	IF(SLIMSK(I).NE.0..AND..NOT.FLAGSNW(I)) THEN
	1379	TSURF(I) = (YY(I) + (ZZ(I) - 1.) * STSOIL(I,1)) / ZZ(I)
	1380	ENDIF
	1381	! IF(SLIMSK(I).EQ.2..AND.SNOWD(I).LE..001) THEN
	1382	IF(SLIMSK(I).EQ.2..AND..NOT.FLAGSNW(I)) THEN
	1383	TSURF(I) = MIN(TSURF(I),T0C)
	1384	ENDIF
	1385	ENDDO
	1386	!!
	1387	DO K = 1, KM
	1388	DO I=1,IM
	1389	IF(SLIMSK(I).EQ.2) THEN
	1390	STSOIL(I,K) = MIN(STSOIL(I,K),T0C)
	1391	ENDIF
	1392	ENDDO
	1393	ENDDO
	1394	!
	1395	! TIME FILTER FOR SOIL AND SKIN TEMPERATURE
	1396	!
	1397	DO I=1,IM
	1398	IF(SLIMSK(I).NE.0.) THEN
	1399	TSKIN(I) = CTFIL1 * TSURF(I) + CTFIL2 * TSKIN(I)
	1400	ENDIF
	1401	ENDDO
	1402	DO K = 1, KM
	1403	DO I=1,IM
	1404	IF(SLIMSK(I).NE.0.) THEN
	1405	STC(I,K) = CTFIL1 * STSOIL(I,K) + CTFIL2 * STC(I,K)
	1406	ENDIF
	1407	ENDDO
	1408	ENDDO
	1409	!
	1410	! GFLUX CALCULATION
	1411	!
	1412	DO I=1,IM
	1413	FLAG(I) = SLIMSK(I).NE.0. &
	1414	! & .AND.SNOWD(I).GT..001 &
	1415	& .AND.FLAGSNW(I)
	1416	ENDDO
	1417	DO I = 1, IM
	1418	IF(FLAG(I)) THEN
	1419	GFLUX(I) = -DFSNOW * (TSKIN(I) - STC(I,1)) &
	1420	& / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
	1421	ENDIF
	1422	ENDDO
	1423	DO I = 1, IM
	1424	! IF(SLIMSK(I).NE.0..AND.SNOWD(I).LE..001) THEN
	1425	IF( SLIMSK(I).NE.0..AND..NOT.FLAGSNW(I)) THEN
	1426	GFLUX(I) = DFT0(I) * (STC(I,1) - TSKIN(I)) &
	1427	& / (-.5 * ZSOIL(I,1))
	1428	ENDIF
	1429	ENDDO
	1430
	1431
	1432	5555 CONTINUE
	1433
	1434	!
	1435	! CALCULATE SENSIBLE HEAT FLUX
	1436	!
	1437	!WRF DO I = 1, IM
	1438	!WRF HFLX(I) = RCH(I) * (TSKIN(I) - THETA1(I))
	1439	!WRF ENDDO
	1440	!
	1441	! THE REST OF THE OUTPUT
	1442	!
	1443	!WRF DO I = 1, IM
	1444	!WRF QSURF(I) = Q1(I) + EVAP(I) / (ELOCP * RCH(I))
	1445	!WRF DM(I) = 1.
	1446	!
	1447	! CONVERT SNOW DEPTH BACK TO MM OF WATER EQUIVALENT
	1448	!
	1449	!WRF SHELEG(I) = SNOWD(I) * 1000.
	1450	!WRF ENDDO
	1451	!
	1452
	1453	DO I = 1, IM
	1454	F10M(I) = FM10(I) / FM(I)
	1455	F10M(I) = min(F10M(I),1._kind_phys)
	1456	U10M(I) = F10M(I) * XRCL(I) * U1(I)
	1457	V10M(I) = F10M(I) * XRCL(I) * V1(I)
	1458	!WRF T2M(I) = TSKIN(I) * (1. - FH2(I) / FH(I)) &
	1459	!WRF & + THETA1(I) * FH2(I) / FH(I)
	1460	!WRF T2M(I) = T2M(I) * SIG2K
	1461	! Q2M(I) = QSURF(I) * (1. - FH2(I) / FH(I)) &
	1462	! & + Q1(I) * FH2(I) / FH(I)
	1463	! T2M(I) = T1
	1464	! Q2M(I) = Q1
	1465	!WRF IF(EVAP(I).GE.0.) THEN
	1466	!
	1467	! IN CASE OF EVAPORATION, USE THE INFERRED QSURF TO DEDUCE Q2M
	1468	!
	1469	!WRF Q2M(I) = QSURF(I) * (1. - FH2(I) / FH(I)) &
	1470	!WRF & + Q1(I) * FH2(I) / FH(I)
	1471	!WRF ELSE
	1472	!
	1473	! FOR DEW FORMATION SITUATION, USE SATURATED Q AT TSKIN
	1474	!
	1475	!jfe QSS(I) = 1000. * FPVS(TSKIN(I))
	1476	!WRF qss(I) = fpvs(tskin(I))
	1477	!WRF QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
	1478	!WRF Q2M(I) = QSS(I) * (1. - FH2(I) / FH(I)) &
	1479	!WRF & + Q1(I) * FH2(I) / FH(I)
	1480	!WRF ENDIF
	1481	!jfe QSS(I) = 1000. * FPVS(T2M(I))
	1482	!WRF QSS(I) = fpvs(t2m(I))
	1483	! QSS(I) = 1000. * T2MO(I)
	1484	!WRF QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
	1485	!WRF Q2M(I) = MIN(Q2M(I),QSS(I))
	1486	ENDDO
	1487	!!
	1488	! DO I = 1, IM
	1489	! RNET(I) = -SLWD(I) - SIGMA * TSKIN(I) **4
	1490	! ENDDO
	1491	!!
	1492	!
	1493	!WRF do i=1,im
	1494	!WRF tem = 1.0 / rho(i)
	1495	!WRF hflx(i) = hflx(i) * tem * cpinv
	1496	!WRF evap(i) = evap(i) * tem * hvapi
	1497	!WRF enddo
	1498
	1499
	1500	!
	1501	!##DG IF(LAT.EQ.LATD) THEN
	1502	!C RBAL = -SLWD-SIGMATSKIN*4+GFLUX
	1503	!C & -EVAP - HFLX
	1504	!##DG PRINT 6000,HFLX,EVAP,GFLUX,
	1505	!##DG& STC(1), STC(2),TSKIN,RNET,SLWD
	1506	!##DG PRINT *, ' T1 =', T1
	1507	6000 FORMAT(8(F8.2,','))
	1508	!C PRINT *, ' EP, ETP,T2M(I) =', EP, ETP,T2M(I)
	1509	!C PRINT *, ' FH, FH2 =', FH, FH2
	1510	!C PRINT *, ' PH, PH2 =', PH, PH2
	1511	!C PRINT *, ' CH, RCH =', CH, RCH
	1512	!C PRINT *, ' TERM1, TERM2 =', TERM1, TERM2
	1513	!C PRINT *, ' RS(I), PLANTR =', RS(I), PLANTR
	1514	!##DG ENDIF
	1515
	1516	RETURN
	1517	END SUBROUTINE PROGTM
	1518	!
	1519	! PROGT2 IS THE SECOND PART OF THE SOIL MODEL THAT IS EXECUTED
	1520	! AFTER PRECIPITATION FOR THE TIME STEP HAS BEEN CALCULATED
	1521	!
	1522	!FPP$ NOCONCUR R
	1523	!FPP$ EXPAND(FUNCDF,FUNCKT,THSAT)
	1524	SUBROUTINE PROGT2(IM,KM,RHSCNPY, &
	1525	& RHSMC,AI,BI,CI,SMC,SLIMSK, &
	1526	& CANOPY,PRECIP,RUNOFF,SNOWMT, &
	1527	& ZSOIL,SOILTYP,SIGMAF,DELT,me)
	1528	!c
	1529	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1530	implicit none
	1531	integer km, IM, me
	1532	real(kind=kind_phys) delt
	1533	real(kind=kind_phys) RHSCNPY(IM), RHSMC(IM,KM), AI(IM,KM), &
	1534	& BI(IM,KM), CI(IM,KM), SMC(IM,KM), &
	1535	& SLIMSK(IM), CANOPY(IM), PRECIP(IM), &
	1536	& RUNOFF(IM), SNOWMT(IM), ZSOIL(IM,KM), &
	1537	& SIGMAF(IM)
	1538	INTEGER SOILTYP(IM)
	1539	!
	1540	integer k, lond, i
	1541	real(kind=kind_phys) CNPY(IM), PRCP(IM), TSAT(IM), &
	1542	& INF(IM), INFMAX(IM), SMSOIL(IM,KM)
	1543	!
	1544	real(kind=kind_phys) cc, ctfil1, ctfil2, delt2, &
	1545	& drip, rffact, rhoh2o, &
	1546	!WRF & rzero, scanop, tdif, thsat, KSAT
	1547	& rzero, scanop, tdif, KSAT
	1548	!
	1549	LOGICAL FLAG(IM)
	1550	!c
	1551	PARAMETER (SCANOP=.5, RHOH2O=1000.)
	1552	PARAMETER (CTFIL1=.5, CTFIL2=1.-CTFIL1)
	1553	! PARAMETER (CTFIL1=1., CTFIL2=1.-CTFIL1)
	1554	PARAMETER (RFFACT=.15)
	1555	!
	1556	!##DG LATD = 44
	1557	LOND = 353
	1558	DELT2 = DELT * 2.
	1559	!
	1560	! PRECIPITATION RATE IS NEEDED IN UNIT OF KG M-2 S-1
	1561	!
	1562	DO I=1,IM
	1563	PRCP(I) = RHOH2O * (PRECIP(I)+SNOWMT(I)) / DELT
	1564	RUNOFF(I) = 0.
	1565	CNPY(I) = CANOPY(I)
	1566	ENDDO
	1567	!##DG IF(LAT.EQ.LATD) THEN
	1568	!##DG PRINT *, ' BEFORE RUNOFF CAL, RHSMC =', RHSMC(1)
	1569	!##DG ENDIF
	1570	!
	1571	! UPDATE CANOPY WATER CONTENT
	1572	!
	1573	DO I=1,IM
	1574	IF(SLIMSK(I).EQ.1.) THEN
	1575	RHSCNPY(I) = RHSCNPY(I) + SIGMAF(I) * PRCP(I)
	1576	CANOPY(I) = CANOPY(I) + DELT * RHSCNPY(I)
	1577	CANOPY(I) = MAX(CANOPY(I),0._kind_phys)
	1578	PRCP(I) = PRCP(I) * (1. - SIGMAF(I))
	1579	IF(CANOPY(I).GT.SCANOP) THEN
	1580	DRIP = CANOPY(I) - SCANOP
	1581	CANOPY(I) = SCANOP
	1582	PRCP(I) = PRCP(I) + DRIP / DELT
	1583	ENDIF
	1584	!
	1585	! CALCULATE INFILTRATION RATE
	1586	!
	1587	INF(I) = PRCP(I)
	1588	TSAT(I) = THSAT(SOILTYP(I))
	1589	! DSAT = FUNCDF(TSAT(I),SOILTYP(I))
	1590	! KSAT = FUNCKT(TSAT(I),SOILTYP(I))
	1591	! INFMAX(I) = -DSAT * (TSAT(I) - SMC(I,1))
	1592	! & / (.5 * ZSOIL(I,1)) &
	1593	! & + KSAT
	1594	INFMAX(I) = (-ZSOIL(I,1)) * &
	1595	& ((TSAT(I) - SMC(I,1)) / DELT - RHSMC(I,1)) &
	1596	& * RHOH2O
	1597	INFMAX(I) = MAX(RFFACT*INFMAX(I),0._kind_phys)
	1598	! IF(SMC(I,1).GE.TSAT(I)) INFMAX(I) = KSAT
	1599	! IF(SMC(I,1).GE.TSAT(I)) INFMAX(I) = ZSOIL(I,1) * RHSMC(I,1)
	1600	IF(INF(I).GT.INFMAX(I)) THEN
	1601	RUNOFF(I) = INF(I) - INFMAX(I)
	1602	INF(I) = INFMAX(I)
	1603	ENDIF
	1604	INF(I) = INF(I) / RHOH2O
	1605	RHSMC(I,1) = RHSMC(I,1) - INF(I) / ZSOIL(I,1)
	1606	ENDIF
	1607	ENDDO
	1608	!!
	1609	!##DG IF(LAT.EQ.LATD) THEN
	1610	!##DG PRINT *, ' PRCP(I), INFMAX(I), RUNOFF =', PRCP(I),INFMAX(I),RUNOFF
	1611	!##DG PRINT *, ' SMSOIL =', SMC(1), SMC(2)
	1612	!##DG PRINT *, ' RHSMC =', RHSMC(1)
	1613	!##DG ENDIF
	1614	!
	1615	! WE CURRENTLY IGNORE THE EFFECT OF RAIN ON SEA ICE
	1616	!
	1617	DO I=1,IM
	1618	FLAG(I) = SLIMSK(I).EQ.1.
	1619	ENDDO
	1620	!!
	1621	!
	1622	! SOLVE THE TRI-DIAGONAL MATRIX
	1623	!
	1624	DO K = 1, KM
	1625	DO I=1,IM
	1626	IF(FLAG(I)) THEN
	1627	RHSMC(I,K) = RHSMC(I,K) * DELT2
	1628	AI(I,K) = AI(I,K) * DELT2
	1629	BI(I,K) = 1. + BI(I,K) * DELT2
	1630	CI(I,K) = CI(I,K) * DELT2
	1631	ENDIF
	1632	ENDDO
	1633	ENDDO
	1634	! FORWARD ELIMINATION
	1635	DO I=1,IM
	1636	IF(FLAG(I)) THEN
	1637	CI(I,1) = -CI(I,1) / BI(I,1)
	1638	RHSMC(I,1) = RHSMC(I,1) / BI(I,1)
	1639	ENDIF
	1640	ENDDO
	1641	DO K = 2, KM
	1642	DO I=1,IM
	1643	IF(FLAG(I)) THEN
	1644	CC = 1. / (BI(I,K) + AI(I,K) * CI(I,K-1))
	1645	CI(I,K) = -CI(I,K) * CC
	1646	RHSMC(I,K)=(RHSMC(I,K)-AI(I,K)RHSMC(I,K-1))CC
	1647	ENDIF
	1648	ENDDO
	1649	ENDDO
	1650	! BACKWARD SUBSTITUTTION
	1651	DO I=1,IM
	1652	IF(FLAG(I)) THEN
	1653	CI(I,KM) = RHSMC(I,KM)
	1654	ENDIF
	1655	ENDDO
	1656	!!
	1657	DO K = KM-1, 1
	1658	DO I=1,IM
	1659	IF(FLAG(I)) THEN
	1660	CI(I,K) = CI(I,K) * CI(I,K+1) + RHSMC(I,K)
	1661	ENDIF
	1662	ENDDO
	1663	ENDDO
	1664	100 CONTINUE
	1665	!
	1666	! UPDATE SOIL MOISTURE
	1667	!
	1668	DO K = 1, KM
	1669	DO I=1,IM
	1670	IF(FLAG(I)) THEN
	1671	SMSOIL(I,K) = SMC(I,K) + CI(I,K)
	1672	SMSOIL(I,K) = MAX(SMSOIL(I,K),0._kind_phys)
	1673	TDIF = MAX(SMSOIL(I,K) - TSAT(I),0._kind_phys)
	1674	RUNOFF(I) = RUNOFF(I) - &
	1675	& RHOH2O * TDIF * ZSOIL(I,K) / DELT
	1676	SMSOIL(I,K) = SMSOIL(I,K) - TDIF
	1677	ENDIF
	1678	ENDDO
	1679	ENDDO
	1680	DO K = 1, KM
	1681	DO I=1,IM
	1682	IF(FLAG(I)) THEN
	1683	SMC(I,K) = CTFIL1 * SMSOIL(I,K) + CTFIL2 * SMC(I,K)
	1684	ENDIF
	1685	ENDDO
	1686	ENDDO
	1687	! IF(FLAG(I)) THEN
	1688	! CANOPY(I) = CTFIL1 * CANOPY(I) + CTFIL2 * CNPY(I)
	1689	! ENDIF
	1690	! I = 1
	1691	! PRINT *, ' SMC'
	1692	! PRINT 6000, SMC(1), SMC(2)
	1693	!6000 FORMAT(2(F8.5,','))
	1694	RETURN
	1695	END SUBROUTINE PROGT2
	1696	FUNCTION KTSOIL(THETA,KTYPE)
	1697	!
	1698	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1699	USE module_progtm , ONLY : TSAT, DFKT
	1700	implicit none
	1701	integer ktype,kw
	1702	real(kind=kind_phys) ktsoil, theta, w
	1703	!
	1704	W = (THETA / TSAT(KTYPE)) * 20. + 1.
	1705	KW = W
	1706	KW = MIN(KW,21)
	1707	KW = MAX(KW,1)
	1708	KTSOIL = DFKT(KW,KTYPE) &
	1709	& + (W - KW) * (DFKT(KW+1,KTYPE) - DFKT(KW,KTYPE))
	1710	RETURN
	1711	END FUNCTION KTSOIL
	1712	FUNCTION FUNCDF(THETA,KTYPE)
	1713	!
	1714	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1715	USE module_progtm , ONLY : TSAT, DFK
	1716	implicit none
	1717	integer ktype,kw
	1718	real(kind=kind_phys) funcdf,theta,w
	1719	!
	1720	W = (THETA / TSAT(KTYPE)) * 20. + 1.
	1721	KW = W
	1722	KW = MIN(KW,21)
	1723	KW = MAX(KW,1)
	1724	FUNCDF = DFK(KW,KTYPE) &
	1725	& + (W - KW) * (DFK(KW+1,KTYPE) - DFK(KW,KTYPE))
	1726	RETURN
	1727	END FUNCTION FUNCDF
	1728	FUNCTION FUNCKT(THETA,KTYPE)
	1729	!
	1730	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1731	USE module_progtm , ONLY : TSAT, KTK
	1732	implicit none
	1733	integer ktype,kw
	1734	real(kind=kind_phys) funckt,theta,w
	1735	!
	1736	W = (THETA / TSAT(KTYPE)) * 20. + 1.
	1737	KW = W
	1738	KW = MIN(KW,21)
	1739	KW = MAX(KW,1)
	1740	FUNCKT = KTK(KW,KTYPE) &
	1741	& + (W - KW) * (KTK(KW+1,KTYPE) - KTK(KW,KTYPE))
	1742	RETURN
	1743	END FUNCTION FUNCKT
	1744	FUNCTION THSAT(KTYPE)
	1745	!
	1746	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1747	USE module_progtm , ONLY : TSAT
	1748	implicit none
	1749	integer ktype
	1750	real(kind=kind_phys) thsat
	1751	!
	1752	THSAT = TSAT(KTYPE)
	1753	RETURN
	1754	END FUNCTION THSAT
	1755	FUNCTION TWLT(KTYPE)
	1756
	1757	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1758	! USE module_progtm
	1759	implicit none
	1760	integer ktype
	1761	real(kind=kind_phys) twlt
	1762	!
	1763	TWLT = .1
	1764	RETURN
	1765	END FUNCTION TWLT
	1766
	1767	END MODULE module_sf_gfs

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format