Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_gfs.F @ 1544

Last change on this file since 1544 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 61.5 KB

Line
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	UST, &
103	ZNT
104
105	REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
106	BR, &
107	CHS, &
108	CHS2, &
109	CPM, &
110	CQS2, &
111	FLHC, &
112	FLQC, &
113	GZ1OZ0, &
114	HFX, &
115	LH, &
116	PSIM, &
117	PSIH, &
118	QFX, &
119	QGH, &
120	QSFC, &
121	U10, &
122	V10, &
123	WSPD
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