| blob | 8411b5b78b03664b57686c506e83e7983c03d63b |
1 !!WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_sf_gfs
6 CONTAINS
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 !-------------------------------------------------------------------
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
80 REAL, INTENT(IN) :: &
81 CP, &
82 EP1, &
83 EP2, &
84 KARMAN, &
85 R, &
86 ROVCP, &
87 XLV
89 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
90 P3D, &
91 QV3D, &
92 T3D, &
93 U3D, &
94 V3D
96 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
97 TSK, &
98 PSFC, &
99 XLAND
101 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
102 UST, &
103 ZNT
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
126 !--------------------------- LOCAL VARS ------------------------------
128 REAL :: ESAT
130 REAL (kind=kind_phys) :: &
131 RHOX
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
174 INTEGER :: &
175 I, &
176 IM, &
177 J, &
178 K, &
179 KM
182 if(itimestep.eq.0) then
183 CALL GFUNCPHYS
184 endif
186 IM=ITE-ITS+1
187 KM=KTE-KTS+1
189 DO J=jts,jte
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
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./(R*T1(I)*(1.+EP1*Q1(I)))
214 Q1(I)=Q1(I)/(1.+Q1(I))
215 ENDDO
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 )
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.8*QV3D(i,kts,j))
241 esat = fpvs(t1(i))
242 QGH(I,J)=ep2*esat/(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
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
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
279 ENDDO
282 END SUBROUTINE SF_GFS
285 !-------------------------------------------------------------------
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 !
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 !
381 PARAMETER (CHARNOCK=.014,CA=.4)!C CA IS THE VON KARMAN CONSTANT
382 PARAMETER (G=grav,sigma=sbc)
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)
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/5*30.,65.,4*100.,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 !
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
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
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
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 !
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
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
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
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
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
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,ETPFAC*HVAP
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)*BFACT*EP(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 & -HFUS*SNOWMT(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
1432 5555 CONTINUE
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 !
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
1500 !
1501 !##DG IF(LAT.EQ.LATD) THEN
1502 !C RBAL = -SLWD-SIGMA*TSKIN**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
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)
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
1767 END MODULE module_sf_gfs