| blob | 3033569a400a563d89474cf1e70ad6824d0b76bc |
1 !-- XLV latent heat of vaporization for water (J/kg)
2 !
3 MODULE module_sf_gfdl
5 !real, dimension(-100:2000,-100:2000), save :: z00
7 CONTAINS
9 !-------------------------------------------------------------------
10 SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, &
11 CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, &
12 DT, SMOIS,num_soil_layers,ISLTYP,ZNT,UST,PSIM,PSIH, &
13 XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, & ! gopal's doing for Ocean coupling
14 QGH,QSFC,U10,V10, &
15 GZ1OZ0,WSPD,BR,ISFFLX, &
16 EP1,EP2,KARMAN,NTSFLG,SFENTH, &
17 ids,ide, jds,jde, kds,kde, &
18 ims,ime, jms,jme, kms,kme, &
19 its,ite, jts,jte, kts,kte )
20 !-------------------------------------------------------------------
21 USE MODULE_GFS_MACHINE, ONLY : kind_phys
22 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
23 USE MODULE_GFS_PHYSCONS, grav => con_g
24 !-------------------------------------------------------------------
25 IMPLICIT NONE
26 !-------------------------------------------------------------------
27 !-- U3D 3D u-velocity interpolated to theta points (m/s)
28 !-- V3D 3D v-velocity interpolated to theta points (m/s)
29 !-- T3D temperature (K)
30 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
31 !-- P3D 3D pressure (Pa)
32 !-- DT time step (second)
33 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
34 !-- ROVCP R/CP
35 !-- R gas constant for dry air (J/kg/K)
36 !-- XLV latent heat of vaporization for water (J/kg)
37 !-- PSFC surface pressure (Pa)
38 !-- ZNT roughness length (m)
39 !-- MAVAIL surface moisture availability (between 0 and 1)
40 !-- UST u* in similarity theory (m/s)
41 !-- PSIM similarity stability function for momentum
42 !-- PSIH similarity stability function for heat
43 !-- XLAND land mask (1 for land, 2 for water)
44 !-- HFX upward heat flux at the surface (W/m^2)
45 !-- QFX upward moisture flux at the surface (kg/m^2/s)
46 !-- TAUX RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
47 !-- TAUY RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
48 !-- LH net upward latent heat flux at surface (W/m^2)
49 !-- GSW downward short wave flux at ground surface (W/m^2)
50 !-- GLW downward long wave flux at ground surface (W/m^2)
51 !-- TSK surface temperature (K)
52 !-- FLHC exchange coefficient for heat (m/s)
53 !-- FLQC exchange coefficient for moisture (m/s)
54 !-- QGH lowest-level saturated mixing ratio
55 !-- U10 diagnostic 10m u wind
56 !-- V10 diagnostic 10m v wind
57 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
58 !-- WSPD wind speed at lowest model level (m/s)
59 !-- BR bulk Richardson number in surface layer
60 !-- ISFFLX isfflx=1 for surface heat and moisture fluxes
61 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
62 !-- KARMAN Von Karman constant
63 !-- SFENTH enthalpy flux factor 0 zot via charnock ..>0 zot enhanced>15m/s
64 !-- ids start index for i in domain
65 !-- ide end index for i in domain
66 !-- jds start index for j in domain
67 !-- jde end index for j in domain
68 !-- kds start index for k in domain
69 !-- kde end index for k in domain
70 !-- ims start index for i in memory
71 !-- ime end index for i in memory
72 !-- jms start index for j in memory
73 !-- jme end index for j in memory
74 !-- kms start index for k in memory
75 !-- kme end index for k in memory
76 !-- its start index for i in tile
77 !-- ite end index for i in tile
78 !-- jts start index for j in tile
79 !-- jte end index for j in tile
80 !-- kts start index for k in tile
81 !-- kte end index for k in tile
82 !-------------------------------------------------------------------
84 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
85 ims,ime, jms,jme, kms,kme, &
86 its,ite, jts,jte, kts,kte, &
87 ISFFLX,NUM_SOIL_LAYERS,NTSFLG
89 REAL, INTENT(IN) :: &
90 CP, &
91 EP1, &
92 EP2, &
93 KARMAN, &
94 R, &
95 ROVCP, &
96 DT, &
97 SFENTH, &
98 XLV
100 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
101 P3D, &
102 QV3D, &
103 T3D, &
104 U3D, &
105 V3D
106 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP
107 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS
109 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
110 PSFC, &
111 GLW, &
112 GSW, &
113 XLAND
114 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
115 TSK, &
116 BR, &
117 CHS, &
118 CHS2, &
119 CPM, &
120 CQS2, &
121 FLHC, &
122 FLQC, &
123 GZ1OZ0, &
124 HFX, &
125 LH, &
126 PSIM, &
127 PSIH, &
128 QFX, &
129 QGH, &
130 QSFC, &
131 UST, &
132 ZNT, &
133 WSPD, &
134 TAUX, & ! gopal's doing for Ocean coupling
135 TAUY
137 REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
138 U10, &
139 V10
142 !--------------------------- LOCAL VARS ------------------------------
144 REAL :: ESAT, &
145 cpcgs, &
146 smc, &
147 smcdry, &
148 smcmax
150 REAL (kind=kind_phys) :: &
151 RHOX
152 REAL, DIMENSION(1:30) :: MAXSMC, &
153 DRYSMC
154 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
155 CH, &
156 CM, &
157 DDVEL, &
158 DRAIN, &
159 EP, &
160 EVAP, &
161 FH, &
162 FH2, &
163 FM, &
164 HFLX, &
165 PH, &
166 PM, &
167 PRSL1, &
168 PRSLKI, &
169 PS, &
170 Q1, &
171 Q2M, &
172 QSS, &
173 QSURF, &
174 RB, &
175 RCL, &
176 RHO1, &
177 SLIMSK, &
178 STRESS, &
179 T1, &
180 T2M, &
181 THGB, &
182 THX, &
183 TSKIN, &
184 SHELEG, &
185 U1, &
186 U10M, &
187 USTAR, &
188 V1, &
189 V10M, &
190 WIND, &
191 Z0RL, &
192 Z1
194 REAL, DIMENSION(kms:kme, ims:ime) :: &
195 rpc, &
196 tpc, &
197 upc, &
198 vpc
200 REAL, DIMENSION(ims:ime) :: &
201 pspc, &
202 pkmax, &
203 tstrc, &
204 zoc, &
205 wetc, &
206 slwdc, &
207 rib, &
208 zkmax, &
209 tkmax, &
210 fxmx, &
211 fxmy, &
212 cdm, &
213 fxh, &
214 fxe, &
215 xxfh, &
216 xxfh2, &
217 wind10, &
218 tjloc
220 INTEGER :: &
221 I, &
222 II, &
223 IGPVS, &
224 IM, &
225 J, &
226 K, &
227 KM
230 DATA MAXSMC/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
231 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
232 0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
233 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
234 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
236 DATA DRYSMC/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
237 0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
238 0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
239 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
240 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
241 DATA IGPVS/0/
242 save igpvs
245 if(igpvs.eq.0) then
246 ! call readzo(glat,glon,6,ims,ime,jms,jme,its,ite,jts,jte,z00)
247 endif
248 igpvs=1
250 IM=ITE-ITS+1
251 KM=KTE-KTS+1
253 WRITE(0,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB 2010 UPGRADS',NTSFLG
255 DO J=jts,jte
257 DO i=its,ite
258 DDVEL(I)=0.
259 RCL(i)=1.
260 PRSL1(i)=P3D(i,kts,j)*.001
261 wetc(i)=1.0
262 if(xland(i,j).lt.1.99) then
263 smc=smois(i,1,j)
264 smcdry=drysmc(isltyp(i,j))
265 smcmax=maxsmc(isltyp(i,j))
266 wetc(i)=(smc-smcdry)/(smcmax-smcdry)
267 wetc(i)=amin1(1.,amax1(wetc(i),0.))
268 endif
269 ! convert from Pa to cgs...
270 pspc(i)=PSFC(i,j)*10.
271 pkmax(i)=P3D(i,kts,j)*10.
272 PS(i)=PSFC(i,j)*.001
273 Q1(I) = QV3D(i,kts,j)
274 rpc(kts,i)=QV3D(i,kts,j)
275 ! QSURF(I)=QSFC(I,J)
276 QSURF(I)=0.
277 SHELEG(I)=0.
278 SLIMSK(i)=ABS(XLAND(i,j)-2.)
279 TSKIN(i)=TSK(i,j)
280 tstrc(i)=TSK(i,j)
281 T1(I) = T3D(i,kts,j)
282 tpc(kts,i)=T3D(i,kts,j)
283 U1(I) = U3D(i,kts,j)
284 upc(kts,i)=U3D(i,kts,j) * 100.
285 USTAR(I) = UST(i,j)
286 V1(I) = V3D(i,kts,j)
287 vpc(kts,i)=v3D(i,kts,j) * 100.
288 Z0RL(I) = ZNT(i,j)*100.
289 zoc(i)=ZNT(i,j)*100.
290 if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i)
291 ! Z0RL(I) = z00(i,j)*100.
292 ! slwdc... GFDL downward net flux in units of cal/(cm**2/min)
293 ! also divide by 10**4 to convert from /m**2 to /cm**2
294 slwdc(i)=gsw(i,j)+glw(i,j)
295 slwdc(i)=0.239*60.*slwdc(i)*1.e-4
296 tjloc(i)=float(j)
298 ENDDO
300 DO i=its,ite
301 PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
302 THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
303 THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
304 RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
305 Q1(I)=Q1(I)/(1.+Q1(I))
306 ENDDO
308 ! if(j==2)then
309 ! write(0,*)'--------------------------------------------'
310 ! write(0,*) 'u, v, t, r, pkmax, pspc,wetc, tjloc,zoc,tstr'
311 ! write(0,*)'--------------------------------------------'
312 ! endif
314 ! do i = its,ite
315 ! WRITE(0,1010)i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), &
316 ! pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i)
317 ! enddo
319 CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,tstrc, &
320 pspc,pkmax,wetc,slwdc,tjloc, &
321 upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, &
322 ids,ide, jds,jde, kds,kde, &
323 ims,ime, jms,jme, kms,kme, &
324 its,ite, jts,jte, kts,kte )
326 ! if(j==2)then
327 ! write(0,*)'--------------------------------------------'
328 ! write(0,*) 'fxh, fxe, fxmx, fxmy, cdm, xxfh zoc,tstrc'
329 ! write(0,*)'--------------------------------------------'
330 ! endif
332 ! do i = its,ite
333 ! WRITE(0,1010)i,j,fxh(i),fxe(i),fxmx(i),fxmy(i),cdm(i),rib(i),xxfh(i),zoc(i),tstrc(i)
334 ! enddo
336 1010 format(2I4,9F11.6)
340 !GFDL CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
341 !GFDL SHELEG,TSKIN,QSURF, &
342 !WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, &
343 !WRF SLRAD,SNOWMT,DELT, &
344 !GFDL Z0RL, &
345 !WRF TG3,GFLUX,F10M, &
346 !GFDL U10M,V10M,T2M,Q2M, &
347 !WRF ZSOIL, &
348 !GFDL CM,CH,RB, &
349 !WRF RHSCNPY,RHSMC,AIM,BIM,CIM, &
350 !GFDL RCL,PRSL1,PRSLKI,SLIMSK, &
351 !GFDL DRAIN,EVAP,HFLX,STRESS,EP, &
352 !GFDL FM,FH,USTAR,WIND,DDVEL, &
353 !GFDL PM,PH,FH2,QSS,Z1 )
356 DO i=its,ite
357 ! update skin temp only when using GFDL slab...
359 IF(NTSFLG==1) then
360 tsk(i,j) = tstrc(i) ! gopal's doing
361 !bugfix 4
362 ! bob's doing patch tsk with neigboring values... are grid boundaries
363 if(j.eq.jde) then
364 tsk(i,j) = tsk(i,j-1)
365 endif
367 if(j.eq.jds) then
368 tsk(i,j) = tsk(i,j+1)
369 endif
371 if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
372 if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
373 endif
376 znt(i,j)= 0.01*abs(zoc(i))
377 wspd(i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
378 wspd(i,j) = amax1(wspd(i,j) ,100.)/100.
379 u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100.
380 v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100.
381 ! br =0.0001*zfull(i,kmax)*dthv/
382 ! & (gmks*theta(i,kmax)*wspd *wspd )
383 ! zkmax = rgas*tpc(kmax,i)*qqlog(kmax)*og
384 zkmax(i) = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
385 !------------------------------------------------------------------------
387 gz1oz0(i,j)=alog(zkmax(i)/znt(i,j))
388 ustar (i)= 0.01*sqrt(cdm(i)* &
389 (upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)))
390 ! convert from g/(cm*cm*sec) to kg/(m*m*sec)
391 qfx (i,j)=-10.*fxe(i) ! BOB: qfx (i,1)=-10.*fxe(i)
392 ! cpcgs = 1.00464e7
393 ! convert from ergs/gram/K to J/kg/K cpmks=1004
394 ! hfx (i,1)=-0.001*cpcgs*fxh(i)
395 hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i)
396 taux (i,j)= fxmx(i)/10. ! gopal's doing for Ocean coupling
397 tauy (i,j)= fxmy(i)/10. ! gopal's doing for Ocean coupling
398 fm(i) = karman/sqrt(cdm(i))
399 fh(i) = karman*xxfh(i)
400 PSIM(i,j)=GZ1OZ0(i,j)-FM(i)
401 PSIH(i,j)=GZ1OZ0(i,j)-FH(i)
402 fh2(i) = karman*xxfh2(i)
403 ch(i) = karman*karman/(fm(i) * fh(i))
404 cm(i) = cdm(i)
406 U10(i,j)=U10M(i)
407 V10(i,j)=V10M(i)
408 BR(i,j)=rib(i)
409 CHS(I,J)=CH(I)*wspd (i,j)
410 CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
411 CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
412 esat = fpvs(t1(i))
413 QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
414 esat = fpvs(tskin(i))
415 qss(i) = ep2*esat/(1000.*ps(i)-esat)
416 QSFC(I,J)=qss(i)
417 ! PSIH(i,j)=PH(i)
418 ! PSIM(i,j)=PM(i)
419 UST(i,j)=ustar(i)
420 ! wspd (i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
421 ! wspd (i,j) = amax1(wspd (i,j) ,100.)/100.
422 ! WSPD(i,j)=WIND(i)
423 ! ZNT(i,j)=Z0RL(i)*.01
424 ENDDO
426 ! write(0,*)'fm,fh,cm,ch(125)', fm(125),fh(125),cm(125),ch(125)
428 DO i=its,ite
429 FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
430 FLQC(i,j)=RHO1(I)*CHS(I,J)
431 ! GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
432 CQS2(i,j)=CHS2(I,J)
433 ENDDO
435 IF (ISFFLX.EQ.0) THEN
436 DO i=its,ite
437 HFX(i,j)=0.
438 LH(i,j)=0.
439 QFX(i,j)=0.
440 ENDDO
441 ELSE
442 DO i=its,ite
443 IF(XLAND(I,J)-1.5.GT.0.)THEN
444 ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
445 ! cpcgs = 1.00464e7
446 ! convert from ergs/gram/K to J/kg/K cpmks=1004
447 ! hfx (i,j)=-0.001*cpcgs*fxh(i)
448 hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i)
449 ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
450 ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
451 ! cpcgs = 1.00464e7
452 ! convert from ergs/gram/K to J/kg/K cpmks=1004
453 ! hfx (i,j)=-0.001*cpcgs*fxh(i)
454 hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,j)= -10.*CP*fxh(i)
455 HFX(I,J)=AMAX1(HFX(I,J),-250.)
456 ENDIF
457 ! QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
458 ! convert from g/(cm*cm*sec) to kg/(m*m*sec)
459 qfx(i,j)=-10.*fxe(i)
460 QFX(I,J)=AMAX1(QFX(I,J),0.)
461 LH(I,J)=XLV*QFX(I,J)
462 ENDDO
463 ENDIF
464 ! if(j.eq.2) write(0,*) 'u3d ,ustar,cdm at end of gfdlsfcmod'
465 ! write(0,*) j,(u3d(ii,1,j),ii=70,90)
466 ! write(0,*) j,(ustar(ii),ii=70,90)
467 ! write(0,*) j,(cdm(ii),ii=70,90)
468 if(j.eq.jds.or.j.eq.jde) then
470 write(0,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde
471 write(0,*) "TSFC", (TSK(i,j),i=its,ite)
472 endif
474 ENDDO ! FOR THE J LOOP I PRESUME
475 ! if(100.ge.its.and.100.le.ite.and.100.ge.jts.and.100.le.jte) then
476 ! write(0,*) 'output vars of sf_gfdl at i,j=100'
477 ! write(0,*) 'TSK',TSK(100,100)
478 ! write(0,*) 'PSFC',PSFC(100,100)
479 ! write(0,*) 'GLW',GLW(100,100)
480 ! write(0,*) 'GSW',GSW(100,100)
481 ! write(0,*) 'XLAND',XLAND(100,100)
482 ! write(0,*) 'BR',BR(100,100)
483 ! write(0,*) 'CHS',CHS(100,100)
484 ! write(0,*) 'CHS2',CHS2(100,100)
485 ! write(0,*) 'CPM',CPM(100,100)
486 ! write(0,*) 'FLHC',FLHC(100,100)
487 ! write(0,*) 'FLQC',FLQC(100,100)
488 ! write(0,*) 'GZ1OZ0',GZ1OZ0(100,100)
489 ! write(0,*) 'HFX',HFX(100,100)
490 ! write(0,*) 'LH',LH(100,100)
491 ! write(0,*) 'PSIM',PSIM(100,100)
492 ! write(0,*) 'PSIH',PSIH(100,100)
493 ! write(0,*) 'QFX',QFX(100,100)
494 ! write(0,*) 'QGH',QGH(100,100)
495 ! write(0,*) 'QSFC',QSFC(100,100)
496 ! write(0,*) 'UST',UST(100,100)
497 ! write(0,*) 'ZNT',ZNT(100,100)
498 ! write(0,*) 'wet',wet(100)
499 ! write(0,*) 'smois',smois(100,1,100)
500 ! write(0,*) 'WSPD',WSPD(100,100)
501 ! write(0,*) 'U10',U10(100,100)
502 ! write(0,*) 'V10',V10(100,100)
503 ! endif
506 END SUBROUTINE SF_GFDL
508 !-------------------------------------------------------------------
509 SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,tstrc, &
510 pspc,pkmax,wetc,slwdc,tjloc, &
511 upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, &
512 ids,ide, jds,jde, kds,kde, &
513 ims,ime, jms,jme, kms,kme, &
514 its,ite, jts,jte, kts,kte )
516 !------------------------------------------------------------------------
517 !
518 ! MFLUX2 computes surface fluxes of momentum, heat,and moisture
519 ! using monin-obukhov. the roughness length "z0" is prescribed
520 ! over land and over ocean "z0" is computed using charnocks formula.
521 ! the universal functions (from similarity theory approach) are
522 ! those of hicks. This is Bob's doing.
523 !
524 !------------------------------------------------------------------------
526 IMPLICIT NONE
527 ! use allocate_mod
528 ! use module_TLDATA , ONLY : tab,table,cp,g,rgas,og
530 ! include 'RESOLUTION.h'
531 ! include 'PARAMETERS.h'
532 ! include 'STDUNITS.h' stdout
534 ! include 'FLAGS.h'
535 ! include 'BKINFO.h' nstep
536 ! include 'CONSLEV.h'
537 ! include 'CONMLEV.h'
538 ! include 'ESTAB.h'
539 ! include 'FILEC.h'
540 ! include 'FILEPC.h'
541 ! include 'GDINFO.h' ngd
542 ! include 'LIMIT.h'
543 ! include 'QLOGS.h'
544 ! include 'TIME.h' dt(nnst)
545 ! include 'WINDD.h'
546 ! include 'ZLDATA.h' old MOBFLX?
548 !-----------------------------------------------------------------------
549 ! user interface variables
550 !-----------------------------------------------------------------------
551 integer,intent(in) :: ids,ide, jds,jde, kds,kde
552 integer,intent(in) :: ims,ime, jms,jme, kms,kme
553 integer,intent(in) :: its,ite, jts,jte, kts,kte
554 integer,intent(in) :: jfix,ntsflg
556 real, intent (out), dimension (ims :ime ) :: fxh
557 real, intent (out), dimension (ims :ime ) :: fxe
558 real, intent (out), dimension (ims :ime ) :: fxmx
559 real, intent (out), dimension (ims :ime ) :: fxmy
560 real, intent (out), dimension (ims :ime ) :: cdm
561 ! real, intent (out), dimension (ims :ime ) :: cdm2
562 real, intent (out), dimension (ims :ime ) :: rib
563 real, intent (out), dimension (ims :ime ) :: xxfh
564 real, intent (out), dimension (ims :ime ) :: xxfh2
565 real, intent (out), dimension (ims :ime ) :: wind10
567 real, intent ( inout), dimension (ims :ime ) :: zoc
568 real, intent ( inout), dimension (ims :ime ) :: tstrc
570 real, intent ( in) :: dt
571 real, intent ( in) :: sfenth
572 real, intent ( in), dimension (ims :ime ) :: pspc
573 real, intent ( in), dimension (ims :ime ) :: pkmax
574 real, intent ( in), dimension (ims :ime ) :: wetc
575 real, intent ( in), dimension (ims :ime ) :: slwdc
576 real, intent ( in), dimension (ims :ime ) :: tjloc
578 real, intent ( in), dimension (kms:kme, ims :ime ) :: upc
579 real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc
580 real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc
581 real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc
583 !-----------------------------------------------------------------------
584 ! internal variables
585 !-----------------------------------------------------------------------
587 integer, parameter :: icntx = 30
589 integer, dimension(1 :ime) :: ifz
590 integer, dimension(1 :ime) :: indx
591 integer, dimension(1 :ime) :: istb
592 integer, dimension(1 :ime) :: it
593 integer, dimension(1 :ime) :: iutb
595 real, dimension(1 :ime) :: aap
596 real, dimension(1 :ime) :: bq1
597 real, dimension(1 :ime) :: bq1p
598 real, dimension(1 :ime) :: delsrad
599 real, dimension(1 :ime) :: ecof
600 real, dimension(1 :ime) :: ecofp
601 real, dimension(1 :ime) :: estso
602 real, dimension(1 :ime) :: estsop
603 real, dimension(1 :ime) :: fmz1
604 real, dimension(1 :ime) :: fmz10
605 real, dimension(1 :ime) :: fmz2
606 real, dimension(1 :ime) :: fmzo1
607 real, dimension(1 :ime) :: foft
608 real, dimension(1 :ime) :: foftm
609 real, dimension(1 :ime) :: frac
610 real, dimension(1 :ime) :: land
611 real, dimension(1 :ime) :: pssp
612 real, dimension(1 :ime) :: qf
613 real, dimension(1 :ime) :: rdiff
614 real, dimension(1 :ime) :: rho
615 real, dimension(1 :ime) :: rkmaxp
616 real, dimension(1 :ime) :: rstso
617 real, dimension(1 :ime) :: rstsop
618 real, dimension(1 :ime) :: sf10
619 real, dimension(1 :ime) :: sf2
620 real, dimension(1 :ime) :: sfm
621 real, dimension(1 :ime) :: sfzo
622 real, dimension(1 :ime) :: sgzm
623 real, dimension(1 :ime) :: slwa
624 real, dimension(1 :ime) :: szeta
625 real, dimension(1 :ime) :: szetam
626 real, dimension(1 :ime) :: t1
627 real, dimension(1 :ime) :: t2
628 real, dimension(1 :ime) :: tab1
629 real, dimension(1 :ime) :: tab2
630 real, dimension(1 :ime) :: tempa1
631 real, dimension(1 :ime) :: tempa2
632 real, dimension(1 :ime) :: theta
633 real, dimension(1 :ime) :: thetap
634 real, dimension(1 :ime) :: tsg
635 real, dimension(1 :ime) :: tsm
636 real, dimension(1 :ime) :: tsp
637 real, dimension(1 :ime) :: tss
638 real, dimension(1 :ime) :: ucom
639 real, dimension(1 :ime) :: uf10
640 real, dimension(1 :ime) :: uf2
641 real, dimension(1 :ime) :: ufh
642 real, dimension(1 :ime) :: ufm
643 real, dimension(1 :ime) :: ufzo
644 real, dimension(1 :ime) :: ugzm
645 real, dimension(1 :ime) :: uzeta
646 real, dimension(1 :ime) :: uzetam
647 real, dimension(1 :ime) :: vcom
648 real, dimension(1 :ime) :: vrtkx
649 real, dimension(1 :ime) :: vrts
650 real, dimension(1 :ime) :: wind
651 real, dimension(1 :ime) :: windp
652 ! real, dimension(1 :ime) :: xxfh
653 real, dimension(1 :ime) :: xxfm
654 real, dimension(1 :ime) :: xxsh
655 real, dimension(1 :ime) :: z10
656 real, dimension(1 :ime) :: z2
657 real, dimension(1 :ime) :: zeta
658 real, dimension(1 :ime) :: zkmax
660 real, dimension(1 :ime) :: pss
661 real, dimension(1 :ime) :: tstar
662 real, dimension(1 :ime) :: ukmax
663 real, dimension(1 :ime) :: vkmax
664 real, dimension(1 :ime) :: tkmax
665 real, dimension(1 :ime) :: rkmax
666 real, dimension(1 :ime) :: zot
667 real, dimension(1 :ime) :: fhzo1
668 real, dimension(1 :ime) :: sfh
670 real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
671 real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
672 real :: szet2, zal2,ugz2
673 real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
674 real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
675 real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
676 real :: shfx,sigt4,reflect
677 real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
678 real :: wndm,ckg
679 real :: yz,y1,y2,y3,y4,windmks,znott,znotm
680 integer:: i,j,ii,iq,nnest,icnt,ngd,ip
681 data amask/ -98.0/
683 !-----------------------------------------------------------------------
684 ! internal variables
685 !-----------------------------------------------------------------------
687 real, dimension (223) :: tab
688 real, dimension (223) :: table
689 real, dimension (101) :: tab11
690 real, dimension (41) :: table4
691 real, dimension (42) :: tab3
692 real, dimension (54) :: table2
693 real, dimension (54) :: table3
694 real, dimension (74) :: table1
695 real, dimension (80) :: tab22
697 equivalence (tab(1),tab11(1))
698 equivalence (tab(102),tab22(1))
699 equivalence (tab(182),tab3(1))
700 equivalence (table(1),table1(1))
701 equivalence (table(75),table2(1))
702 equivalence (table(129),table3(1))
703 equivalence (table(183),table4(1))
705 !-----------------------------------------------------------------------
706 ! tables used to obtain the vapor pressures or saturated vapor
707 ! pressure
708 !-----------------------------------------------------------------------
710 data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
711 &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
712 &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
713 &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
714 &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
715 &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
716 &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
717 &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
719 data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
720 &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
721 &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
722 &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
723 &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
724 &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
725 &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
726 &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
727 &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
728 &190220.,199260./
730 data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
731 &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
732 &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
733 &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
734 &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
735 &1013250.,1049940.,1087740./
737 data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
738 & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
739 &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
740 &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
741 &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
742 &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
743 &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
744 &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
745 &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
746 &.7220e+01/
748 data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
749 &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
750 &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
751 &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
752 &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
753 &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
754 &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
755 &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
756 &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
757 &.7090e+03/
759 data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
760 &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
761 &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
762 &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
763 &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
764 &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
765 &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
766 &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
767 &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
768 &.9780e+04/
770 data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
771 &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
772 &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
773 &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
774 &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
775 &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
776 &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
777 !
778 ! spcify constants needed by MFLUX2
779 !
780 real,parameter :: cp = 1.00464e7
781 real,parameter :: g = 980.6
782 real,parameter :: rgas = 2.87e6
783 real,parameter :: og = 1./g
785 !
786 ! character*10 routine
787 ! routine = 'mflux2'
788 !
789 !------------------------------------------------------------------------
790 ! set water availability constant "ecof" and land mask "land".
791 ! limit minimum wind speed to 100 cm/s
792 !------------------------------------------------------------------------
793 j = IFIX(tjloc(its))
794 ! constants for 10 m winds (correction for knots
795 !
796 cor1 = .120
797 cor2 = 720.
798 yz= -0.0001344
799 y1= 3.015e-05
800 y2= 1.517e-06
801 y3= -3.567e-08
802 y4= 2.046e-10
804 do i = its,ite
805 z10(i) = 1000.
806 z2 (i) = 200.
807 pss(i) = pspc(i)
808 tstar(i) = tstrc(i)
809 ukmax(i) = upc(1,i)
810 vkmax(i) = vpc(1,i)
811 tkmax(i) = tpc(1,i)
812 rkmax(i) = rpc(1,i)
813 enddo
815 ! write(0,*)'z10,pss,tstar,u...rkmax(ite)', &
816 ! z10(ite), pss(ite),tstar(ite),ukmax(ite), &
817 ! vkmax(ite),tkmax(ite),rkmax(ite)
819 do i = its,ite
820 windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
821 wind (i) = amax1(windp(i),100.)
822 if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind(i)*wind(i)*og
823 if (zoc(i) .GT. 0.0) then
824 ecof(i) = wetc(i)
825 land(i) = 1.0
826 zot (i) = zoc(i)
827 else
828 ecof(i) = wetc(i)
829 land(i) = 0.0
830 #ifdef HWRF
831 zot (i) = zoc(i)
832 ! now use 2 regime fit for znot thermal
833 windmks=wind(i)*.01
834 znott=0.2375*exp(-0.5250*windmks) + 0.0025*exp(-0.0211*windmks)
835 znott=0.01*znott
836 ! go back to mon/obv et al for below 7m/s
837 if(windmks.le. 7.) &
838 znott = 0.01*zot (i)
839 ! end of kwon correction....
840 ! put in patch to increase znott from thermal to mom vaules 20-40m/s
841 if(windmks.ge.20.) then
842 ! znot4020=(.163e-2)-(.104e-4)
843 znott =1.636e-14 * windmks**6.8905
844 endif
845 ! back to cgs
846 zot (i) = 100.*znott
847 #else
848 ! zot (i) = zoc(i)
849 ! now use 2 regime fit for znot thermal
850 windmks=wind(i)*.01
851 znott=1.9551e-5 - 2.6338e-7 * windmks
852 if(windmks.le.10.) znott=0.0025542 * windmks **(-1.8023)
853 znott=amax1(1.e-6,znott)
854 ! go back to moon et al for below 7m/s
855 if(windmks.le. 7.) &
856 znott = (0.0185/9.8*(7.59e-8*wind(i)**2+ &
857 2.46e-4*wind(i))**2)
858 ! back to cgs
859 zot (i) = 100.*znott
860 #endif
861 ! in hwrf, thermal znot(zot) is passed as argument zoc
862 ! in hwrf, momentum znot is recalculated internally
863 zoc(i)=-(0.0185/9.8*(7.59e-8*wind(i)**2+ &
864 2.46e-4*wind(i))**2)*100.
865 if(wind(i).ge.1250.0) &
866 zoc(i)=-(.000739793 * wind(i) -0.58)/10
867 if(wind(i).ge.3000.) then
868 ! kwon znotm
869 znotm = yz +windmks*y1 +windmks**2*y2 +windmks**3*y3 +windmks**4*y4 !powell 2003
870 ! back to cgs
871 zoc(i) = 100.*znotm
872 endif
873 endif
875 #ifdef HWRF
876 !
877 ! in this version set thermal znot = momentum znot of powell 2003 above 40m/s
878 if(windmks.ge.40.) zot (i) = zoc(i)
879 #endif
881 !------------------------------------------------------------------------
882 ! where necessary modify zo values over ocean.
883 !------------------------------------------------------------------------
885 enddo
887 !------------------------------------------------------------------------
888 ! define constants:
889 ! a and c = constants used in evaluating universal function for
890 ! stable case
891 ! ca = karmen constant
892 ! cm01 = constant part of vertical integral of universal
893 ! function; stable case ( 0.5 < zeta < or = 10.0)
894 ! cm02 = constant part of vertical integral of universal
895 ! function; stable case ( zeta > 10.0)
896 !------------------------------------------------------------------------
898 en = 2.
899 c = .76
900 a = 5.
901 ca = .4
902 cmo1 = .5*a - 1.648
903 cmo2 = 17.193 + .5*a - 10.*c
904 boycon = .61
905 rovcp=rgas/cp
906 ! write(0,*)'rgas,cp,rovcp ', rgas,cp,rovcp
908 ! write(0,*)'--------------------------------------------------'
909 ! write(0,*)'pkmax, pspc, theta, zkmax, zoc'
910 ! write(0,*)'--------------------------------------------------'
912 do i = its,ite
913 ! theta(i) = tkmax(i)*rqc9
914 theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
915 vrtkx(i) = 1.0 + boycon*rkmax(i)
916 ! zkmax(i) = rgas*tkmax(i)*qqlog(kmax)*og
917 zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og
918 ! IF(I==78)write(0,*)I,JFIX,pkmax(i),pspc(i),theta(i),zkmax(i),zoc(i)
919 enddo
921 ! write(0,*)'pkmax,pspc ', pkmax,pspc
922 ! write(0,*)'theta, zkmax, zoc ', theta, zkmax, zoc
924 !------------------------------------------------------------------------
925 ! get saturation mixing ratios at surface
926 !------------------------------------------------------------------------
928 do i = its,ite
929 tsg (i) = tstar(i)
930 tab1 (i) = tstar(i) - 153.16
931 it (i) = IFIX(tab1(i))
932 tab2 (i) = tab1(i) - FLOAT(it(i))
933 t1 (i) = tab(it(i) + 1)
934 t2 (i) = table(it(i) + 1)
935 estso(i) = t1(i) + tab2(i)*t2(i)
936 psps1 = (pss(i) - estso(i))
937 if(psps1 .EQ. 0.0)then
938 psps1 = .1
939 endif
940 rstso(i) = 0.622*estso(i)/psps1
941 vrts (i) = 1. + boycon*ecof(i)*rstso(i)
942 enddo
944 !------------------------------------------------------------------------
945 ! check if consideration of virtual temperature changes stability.
946 ! if so, set "dthetav" to near neutral value (1.0e-4). also check
947 ! for very small lapse rates; if ABS(tempa1) <1.0e-4 then
948 ! tempa1=1.0e-4
949 !------------------------------------------------------------------------
951 do i = its,ite
952 tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
953 tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
954 if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
955 tab1(i) = ABS(tempa1(i))
956 if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
957 !------------------------------------------------------------------------
958 ! compute bulk richardson number "rib" at each point. if "rib"
959 ! exceeds 95% of critical richardson number "tab1" then "rib = tab1"
960 !------------------------------------------------------------------------
962 rib (i) = g*zkmax(i)*tempa1(i)/ &
963 (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
964 tab2(i) = ABS(zoc(i))
965 tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
966 if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
967 enddo
969 do i = its,ite
970 zeta(i) = ca*rib(i)/0.03
971 enddo
973 ! write(0,*)'rib,zeta,vrtkx,vrts(ite) ', rib(ite),zeta(ite), &
974 ! vrtkx(ite),vrts(ite)
975 !------------------------------------------------------------------------
976 ! begin looping through points on line, solving wegsteins iteration
977 ! for zeta at each point, and using hicks functions
978 !------------------------------------------------------------------------
980 !------------------------------------------------------------------------
981 ! set initial guess of zeta=non - dimensional height "szeta" for
982 ! stable points
983 !------------------------------------------------------------------------
985 rca = 1./ca
986 enrca = en*rca
987 ! turn off interfacial layer by zeroing out enrca
988 enrca = 0.0
989 zog = .0185*og
991 !------------------------------------------------------------------------
992 ! stable points
993 !------------------------------------------------------------------------
995 ip = 0
996 do i = its,ite
997 if (zeta(i) .GE. 0.0) then
998 ip = ip + 1
999 istb(ip) = i
1000 endif
1001 enddo
1003 if (ip .EQ. 0) go to 170
1004 do i = 1,ip
1005 szetam(i) = 1.0e+30
1006 sgzm(i) = 0.0e+00
1007 szeta(i) = zeta(istb(i))
1008 ifz(i) = 1
1009 enddo
1011 !------------------------------------------------------------------------
1012 ! begin wegstein iteration for "zeta" at stable points using
1013 ! hicks(1976)
1014 !------------------------------------------------------------------------
1016 do icnt = 1,icntx
1017 do i = 1,ip
1018 if (ifz(i) .EQ. 0) go to 80
1019 zal1g = ALOG(szeta(i))
1020 if (szeta(i) .LE. 0.5) then
1021 fmz1(i) = (zal1g + a*szeta(i))*rca
1022 else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
1023 rzeta1 = 1./szeta(i)
1024 fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
1025 0.5*rzeta1*rzeta1 + cmo1)*rca
1026 else if (szeta(i) .GT. 10.) then
1027 fmz1(i) = (c*szeta(i) + cmo2)*rca
1028 endif
1029 szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
1030 zal2g = ALOG(szetao)
1031 if (szetao .LE. 0.5) then
1032 fmzo1(i) = (zal2g + a*szetao)*rca
1033 sfzo (i) = 1. + a*szetao
1034 else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
1035 rzeta2 = 1./szetao
1036 fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
1037 0.5*rzeta2*rzeta2 + cmo1)*rca
1038 sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1039 else if (szetao .GT. 10.) then
1040 fmzo1(i) = (c*szetao + cmo2)*rca
1041 sfzo (i) = c*szetao
1042 endif
1045 ! compute heat & moisture parts of zot.. for calculation of sfh
1047 szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i)
1048 zal2gh = ALOG(szetho)
1049 if (szetho .LE. 0.5) then
1050 fhzo1(i) = (zal2gh + a*szetho)*rca
1051 sfzo (i) = 1. + a*szetho
1052 else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
1053 rzeta2 = 1./szetho
1054 fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
1055 0.5*rzeta2*rzeta2 + cmo1)*rca
1056 sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1057 else if (szetho .GT. 10.) then
1058 fhzo1(i) = (c*szetho + cmo2)*rca
1059 sfzo (i) = c*szetho
1060 endif
1062 !------------------------------------------------------------------------
1063 ! compute universal function at 10 meters for diagnostic purposes
1064 !------------------------------------------------------------------------
1066 !!!! if (ngd .EQ. nNEST) then
1067 szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i)
1068 zal10 = ALOG(szet10)
1069 if (szet10 .LE. 0.5) then
1070 fmz10(i) = (zal10 + a*szet10)*rca
1071 else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
1072 rzeta2 = 1./szet10
1073 fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
1074 0.5*rzeta2*rzeta2 + cmo1)*rca
1075 else if (szet10 .GT. 10.) then
1076 fmz10(i) = (c*szet10 + cmo2)*rca
1077 endif
1078 sf10(i) = fmz10(i) - fmzo1(i)
1079 ! compute 2m values for diagnostics in HWRF
1080 szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i)
1081 zal2 = ALOG(szet2 )
1082 if (szet2 .LE. 0.5) then
1083 fmz2 (i) = (zal2 + a*szet2 )*rca
1084 else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
1085 rzeta2 = 1./szet2
1086 fmz2 (i) = (8.*zal2 + 4.25*rzeta2 - &
1087 0.5*rzeta2*rzeta2 + cmo1)*rca
1088 else if (szet2 .GT. 2.) then
1089 fmz2 (i) = (c*szet2 + cmo2)*rca
1090 endif
1091 sf2 (i) = fmz2 (i) - fmzo1(i)
1093 !!!! endif
1094 sfm(i) = fmz1(i) - fmzo1(i)
1095 sfh(i) = fmz1(i) - fhzo1(i)
1096 sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
1097 (sfh(i) + enrca*sfzo(i))
1098 fmz = (sgz - szeta(i))/szeta(i)
1099 fmzo = ABS(fmz)
1100 if (fmzo .GE. 5.0e-5) then
1101 sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
1102 if(sq .EQ. 1) then
1103 write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (STABLE CASE)'
1104 write(0,*)'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
1105 endif
1106 szetam(i) = szeta(i)
1107 szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
1108 sgzm (i) = sgz
1109 else
1110 ifz(i) = 0
1111 endif
1112 80 continue
1113 enddo
1114 enddo
1116 do i = 1,ip
1117 if (ifz(i) .GE. 1) go to 110
1118 enddo
1120 go to 130
1122 110 continue
1123 write(6,120)
1124 120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
1125 ! call MPI_CLOSE(1,routine)
1127 !------------------------------------------------------------------------
1128 ! update "zo" for ocean points. "zo"cannot be updated within the
1129 ! wegsteins iteration as the scheme (for the near neutral case)
1130 ! can become unstable
1131 !------------------------------------------------------------------------
1133 130 continue
1134 do i = 1,ip
1135 szo = zoc(istb(i))
1136 if (szo .LT. 0.0) then
1137 wndm=wind(istb(i))*0.01
1138 if(wndm.lt.15.0) then
1139 ckg=0.0185*og
1140 else
1141 !! ckg=(0.000308*wndm+0.00925)*og
1142 !! ckg=(0.000616*wndm)*og
1143 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1144 endif
1146 szo = - ckg*wind(istb(i))*wind(istb(i))/ &
1147 (sfm(i)*sfm(i))
1148 cons_p000001 = .000001
1149 cons_7 = 7.
1150 vis = 1.4E-1
1152 ustar = sqrt( -szo / zog)
1153 restar = -ustar * szo / vis
1154 restar = max(restar,cons_p000001)
1155 ! Rat taken from Zeng, Zhao and Dickinson 1997
1156 rat = 2.67 * restar ** .25 - 2.57
1157 rat = min(rat ,cons_7) !constant
1158 rat=0.
1159 zot(istb(i)) = szo * exp(-rat)
1160 else
1161 zot(istb(i)) = zoc(istb(i))
1162 endif
1164 ! in hwrf thermal znot is loaded back into the zoc array for next step
1165 zoc(istb(i)) = szo
1166 enddo
1168 do i = 1,ip
1169 xxfm(istb(i)) = sfm(i)
1170 xxfh(istb(i)) = sfh(i)
1171 xxfh2(istb(i)) = sf2 (i)
1172 xxsh(istb(i)) = sfzo(i)
1173 enddo
1175 !------------------------------------------------------------------------
1176 ! obtain wind at 10 meters for diagnostic purposes
1177 !------------------------------------------------------------------------
1179 !!! if (ngd .EQ. nNEST) then
1180 do i = 1,ip
1181 wind10(istb(i)) = sf10(i)*wind(istb(i))/sfm(i)
1182 wind10(istb(i)) = wind10(istb(i)) * 1.944
1183 if(wind10(istb(i)) .GT. 6000.0) then
1184 wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
1185 - cor2
1186 endif
1187 ! the above correction done by GFDL in centi-kts!!!-change back
1188 wind10(istb(i)) = wind10(istb(i)) / 1.944
1189 enddo
1190 !!! endif
1191 !!! if (ngd .EQ. nNEST-1 .AND. llwe .EQ. 1 ) then
1192 !!! do i = 1,ip
1193 !!! wind10c(istb(i),j) = sf10(i)*wind(istb(i))/sfm(i)
1194 !!! wind10c(istb(i),j) = wind10c(istb(i),j) * 1.944
1195 !!! if(wind10c(istb(i),j) .GT. 6000.0) then
1196 !!! wind10c(istb(i),j)=wind10c(istb(i),j)+wind10c(istb(i),j)*cor1
1197 !!! * - cor2
1198 !!! endif
1199 !!! enddo
1200 !!! endif
1202 !------------------------------------------------------------------------
1203 ! unstable points
1204 !------------------------------------------------------------------------
1206 170 continue
1208 iq = 0
1209 do i = its,ite
1210 if (zeta(i) .LT. 0.0) then
1211 iq = iq + 1
1212 iutb(iq) = i
1213 endif
1214 enddo
1216 if (iq .EQ. 0) go to 290
1217 do i = 1,iq
1218 uzeta (i) = zeta(iutb(i))
1219 ifz (i) = 1
1220 uzetam(i) = 1.0e+30
1221 ugzm (i) = 0.0e+00
1222 enddo
1224 !------------------------------------------------------------------------
1225 ! begin wegstein iteration for "zeta" at unstable points using
1226 ! hicks functions
1227 !------------------------------------------------------------------------
1229 do icnt = 1,icntx
1230 do i = 1,iq
1231 if (ifz(i) .EQ. 0) go to 200
1232 ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i))))
1233 uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1234 ux11 = 1. - 16.*uzeta(i)
1235 ux12 = 1. - 16.*uzetao
1236 y = SQRT(ux11)
1237 yo = SQRT(ux12)
1238 ufzo(i) = 1./yo
1239 ux13 = (1. + y)/(1. + yo)
1240 ux21 = ALOG(ux13)
1241 ufh(i) = (ugzzo - 2.*ux21)*rca
1242 ! recompute scalers for ufm in terms of mom znot... zoc
1243 ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i))))
1244 uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1245 ux11 = 1. - 16.*uzeta(i)
1246 ux12 = 1. - 16.*uzetao
1247 y = SQRT(ux11)
1248 yo = SQRT(ux12)
1249 ux13 = (1. + y)/(1. + yo)
1250 ux21 = ALOG(ux13)
1251 ! ufzo(i) = 1./yo
1252 x = SQRT(y)
1253 xo = SQRT(yo)
1254 xnum = (x**2 + 1.)*((x + 1.)**2)
1255 xden = (xo**2 + 1.)*((xo + 1.)**2)
1256 xtan = ATAN(x) - ATAN(xo)
1257 ux3 = ALOG(xnum/xden)
1258 ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca
1259 !!!! if (ngd .EQ. nNEST) then
1261 !------------------------------------------------------------------------
1262 ! obtain ten meter winds for diagnostic purposes
1263 !------------------------------------------------------------------------
1265 ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i))))
1266 uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1267 uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1268 ux11 = 1. - 16.*uzet1o
1269 ux12 = 1. - 16.*uzetao
1270 y = SQRT(ux11)
1271 y10 = SQRT(ux12)
1272 ux13 = (1. + y)/(1. + y10)
1273 ux21 = ALOG(ux13)
1274 x = SQRT(y)
1275 x10 = SQRT(y10)
1276 xnum = (x**2 + 1.)*((x + 1.)**2)
1277 xden = (x10**2 + 1.)*((x10 + 1.)**2)
1278 xtan = ATAN(x) - ATAN(x10)
1279 ux3 = ALOG(xnum/xden)
1280 uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca
1282 ! obtain 2m values for diagnostics...
1285 ugz2 = ALOG(z2 (iutb(i))/ABS(zoc(iutb(i))))
1286 uzet1o = ABS(z2 (iutb(i)))/zkmax(iutb(i))*uzeta(i)
1287 uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1288 ux11 = 1. - 16.*uzet1o
1289 ux12 = 1. - 16.*uzetao
1290 y = SQRT(ux11)
1291 yo = SQRT(ux12)
1292 ux13 = (1. + y)/(1. + yo)
1293 ux21 = ALOG(ux13)
1294 uf2 (i) = (ugzzo - 2.*ux21)*rca
1297 !!! endif
1298 ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
1299 ux1 = (ugz - uzeta(i))/uzeta(i)
1300 ux2 = ABS(ux1)
1301 if (ux2 .GE. 5.0e-5) then
1302 uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
1303 uzetam(i) = uzeta(i)
1304 if(uq .EQ. 1) then
1305 write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (UNSTABLE CASE)'
1306 write(0,*)'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1307 endif
1308 uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
1309 ugzm (i) = ugz
1310 else
1311 ifz(i) = 0
1312 endif
1313 200 continue
1314 enddo
1315 enddo
1318 do i = 1,iq
1319 if (ifz(i) .GE. 1) go to 230
1320 enddo
1322 go to 250
1324 230 continue
1325 write(6,240)
1326 240 format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ')
1327 ! call MPI_CLOSE(1,routine)
1329 !------------------------------------------------------------------------
1330 ! gather unstable values
1331 !------------------------------------------------------------------------
1333 250 continue
1335 !------------------------------------------------------------------------
1336 ! update "zo" for ocean points. zo cannot be updated within the
1337 ! wegsteins iteration as the scheme (for the near neutral case)
1338 ! can become unstable.
1339 !------------------------------------------------------------------------
1341 do i = 1,iq
1342 uzo = zoc(iutb(i))
1343 if (zoc(iutb(i)) .LT. 0.0) then
1344 wndm=wind(iutb(i))*0.01
1345 if(wndm.lt.15.0) then
1346 ckg=0.0185*og
1347 else
1348 !! ckg=(0.000308*wndm+0.00925)*og <
1349 !! ckg=(0.000616*wndm)*og <
1350 ckg=(4*0.000308*wndm)*og
1351 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1352 endif
1353 uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
1354 cons_p000001 = .000001
1355 cons_7 = 7.
1356 vis = 1.4E-1
1358 ustar = sqrt( -uzo / zog)
1359 restar = -ustar * uzo / vis
1360 restar = max(restar,cons_p000001)
1361 ! Rat taken from Zeng, Zhao and Dickinson 1997
1362 rat = 2.67 * restar ** .25 - 2.57
1363 rat = min(rat ,cons_7) !constant
1364 rat=0.0
1365 zot(iutb(i)) = uzo * exp(-rat)
1366 else
1367 zot(iutb(i)) = zoc(iutb(i))
1368 endif
1369 ! in hwrf thermal znot is loaded back into the zoc array for next step
1370 zoc(iutb(i)) = uzo
1371 enddo
1373 !------------------------------------------------------------------------
1374 ! obtain wind at ten meters for diagnostic purposes
1375 !------------------------------------------------------------------------
1377 !!! if (ngd .EQ. nNEST) then
1378 do i = 1,iq
1379 wind10(iutb(i)) = uf10(i)*wind(iutb(i))/ufm(i)
1380 wind10(iutb(i)) = wind10(iutb(i)) * 1.944
1381 if(wind10(iutb(i)) .GT. 6000.0) then
1382 wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
1383 - cor2
1384 endif
1385 ! the above correction done by GFDL in centi-kts!!!-change back
1386 wind10(iutb(i)) = wind10(iutb(i)) / 1.944
1387 enddo
1388 !!! endif
1389 !!! if (ngd .EQ. nNEST-1) then
1390 !!! do i = 1,iq
1391 !!! wind10c(iutb(i),j) = uf10(i)*wind(iutb(i))/ufm(i)
1392 !!! wind10c(iutb(i),j) = wind10c(iutb(i),j) * 1.944
1393 !!! if(wind10c(iutb(i),j) .GT. 6000.0) then
1394 !!! wind10c(iutb(i),j)=wind10c(iutb(i),j)+wind10c(iutb(i),j)*cor1
1395 !!! * - cor2
1396 !!! endif
1397 !!! enddo
1398 !!! endif
1400 do i = 1,iq
1401 xxfm(iutb(i)) = ufm(i)
1402 xxfh(iutb(i)) = ufh(i)
1403 xxfh2(iutb(i)) = uf2 (i)
1404 xxsh(iutb(i)) = ufzo(i)
1405 enddo
1407 290 continue
1409 do i = its,ite
1410 ucom(i) = ukmax(i)
1411 vcom(i) = vkmax(i)
1412 if (windp(i) .EQ. 0.0) then
1413 windp(i) = 100.0
1414 ucom (i) = 100.0/SQRT(2.0)
1415 vcom (i) = 100.0/SQRT(2.0)
1416 endif
1417 rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - &
1418 tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i))))
1419 bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i)))
1420 enddo
1422 ! do land sfc temperature prediction if ntsflg=1
1423 ! ntsflg = 1 ! gopal's doing
1425 if (ntsflg .EQ. 0) go to 370
1426 alll = 600.
1427 xks = 0.01
1428 hcap = .5/2.39e-8
1429 pith = SQRT(4.*ATAN(1.0))
1430 alfus = alll/2.39e-8
1431 teps = 0.1
1432 ! slwdc... in units of cal/min ????
1433 ! slwa... in units of ergs/sec/cm*2
1434 ! 1 erg=2.39e-8 cal
1435 !------------------------------------------------------------------------
1436 ! pack land and sea ice points
1437 !------------------------------------------------------------------------
1439 ip = 0
1440 do i = its,ite
1441 if (land(i) .EQ. 1) then
1442 ip = ip + 1
1443 indx (ip) = i
1444 ! slwa is defined as positive down....
1445 slwa (ip) = slwdc(i)/(2.39e-8*60.)
1446 tss (ip) = tstar(i)
1447 thetap (ip) = theta(i)
1448 rkmaxp (ip) = rkmax(i)
1449 aap (ip) = 5.673e-5
1450 pssp (ip) = pss(i)
1451 ecofp (ip) = ecof(i)
1452 estsop (ip) = estso(i)
1453 rstsop (ip) = rstso(i)
1454 bq1p (ip) = bq1(i)
1455 bq1p (ip) = amax1(bq1p(ip),0.1e-3)
1456 delsrad(ip) = dt *pith/(hcap*SQRT(3600.*24.*xks))
1457 endif
1458 enddo
1460 !------------------------------------------------------------------------
1461 ! initialize variables for first pass of iteration
1462 !------------------------------------------------------------------------
1464 do i = 1,ip
1465 ifz (i) = 1
1466 tsm (i) = tss(i)
1467 rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1468 !!! if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. &
1469 !!! nstep .EQ. -99 .AND. ngd .EQ. 1) then
1471 !!! if (j .EQ. 1 .AND. i .EQ. 1) write(6,300)
1472 300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ')
1474 !! foftm(i) = thetap(i) + 1./(cp*bq1p(i))*(slwa(i) - aap(i)* &
1475 !! tsm(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i))
1476 !! else
1478 foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - &
1479 cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1480 rdiff(i))
1481 !! endif
1482 tsp(i) = foftm(i)
1483 enddo
1485 !------------------------------------------------------------------------
1486 ! do iteration to determine "tstar" at new time level
1487 !------------------------------------------------------------------------
1489 do icnt = 1,icntx
1490 do i = 1,ip
1491 if (ifz(i) .EQ. 0) go to 330
1492 tab1 (i) = tsp(i) - 153.16
1493 it (i) = IFIX(tab1(i))
1494 tab2 (i) = tab1(i) - FLOAT(it(i))
1495 t1 (i) = tab(it(i) + 1)
1496 t2 (i) = table(it(i) + 1)
1497 estsop(i) = t1(i) + tab2(i)*t2(i)
1498 psps2 = (pssp(i) - estsop(i))
1499 if(psps2 .EQ. 0.0)then
1500 psps2 = .1
1501 endif
1502 rstsop(i) = 0.622*estsop(i)/psps2
1503 rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1504 !!! if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. &
1505 !!! nstep .EQ. -99 .AND. ngd .EQ. 1) then
1506 !!! foft(i) = thetap(i) + (1./(cp*bq1p(i)))*(slwa(i) - aap(i)* &
1507 !!! tsp(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i))
1508 !!! else
1509 foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - &
1510 cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1511 rdiff(i))
1512 !!! endif
1513 !!! if (ngd .EQ. 1 .AND. j .EQ. 48 .AND. i .EQ. 19) then
1514 !!! reflect = slwa(i)
1515 !!! sigt4 = -aap(i)*tsp(i)**4
1516 !!! shfx = -cp*bq1p(i)*(tsp(i) - thetap(i))
1517 !!! alevp = ecofp(i)*alfus*bq1p(i)*rdiff(i)
1518 !!! delten = delsrad(i)
1519 !!! diffot = foft(i) - tss(i)
1520 !!! endif
1521 frac(i) = ABS((foft(i) - tsp(i))/tsp(i))
1523 !------------------------------------------------------------------------
1524 ! check for convergence of all points use wegstein iteration
1525 !------------------------------------------------------------------------
1527 if (frac(i) .GE. teps) then
1528 qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
1529 tsm (i) = tsp(i)
1530 tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
1531 foftm(i) = foft(i)
1532 else
1533 ifz(i) = 0
1534 endif
1535 330 continue
1536 enddo
1537 enddo
1539 !------------------------------------------------------------------------
1540 ! check for convergence of "t star" prediction
1541 !------------------------------------------------------------------------
1543 do i = 1,ip
1544 if (ifz(i) .EQ. 1) then
1545 write(6, 340) tsp(i), i, j
1546 340 format(2X, ' NON-CONVERGENCE OF T* PREDICTED (T*,I,J) = ', E14.8, &
1547 2I5)
1549 write(6,345) indx(i), j, tstar(indx(i)), tsp(i), ip
1550 345 format(2X, ' I, J, OLD T*, NEW T*, NPTS ', 2I5, 2E14.8, I5)
1552 ! write(6,350) reflect, sigt4, shfx, alevp, delten, diffot
1553 350 format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', &
1554 6E14.8)
1556 ! call MPI_CLOSE(1,routine)
1557 endif
1558 enddo
1560 do i = 1,ip
1561 ii = indx(i)
1562 tstrc(ii) = tsp (i)
1563 enddo
1565 !------------------------------------------------------------------------
1566 ! compute fluxes and momentum drag coef
1567 !------------------------------------------------------------------------
1569 370 continue
1570 do i = its,ite
1571 fxh(i) = bq1(i)*(theta(i) - tsg(i))
1572 fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i))
1573 if (fxe(i) .GT. 0.0) fxe(i) = 0.0
1574 fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ &
1575 windp(i)
1576 fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
1577 windp(i)
1578 cdm(i) = 1./(xxfm(i)*xxfm(i))
1579 ! print *, 'i,zot,zoc,cdm,cdm2,tsg,wind', &
1580 ! i, zot(i),zoc(i), cdm(i),cdm2(i), tsg(i),wind(i)
1581 enddo
1583 return
1584 end subroutine MFLUX2
1586 SUBROUTINE hwrfsfcinit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
1587 SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
1588 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
1589 TMN, &
1590 num_soil_layers, &
1591 allowed_to_read, &
1592 ids,ide, jds,jde, kds,kde, &
1593 ims,ime, jms,jme, kms,kme, &
1594 its,ite, jts,jte, kts,kte )
1596 IMPLICIT NONE
1598 ! Arguments
1599 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1600 ims,ime, jms,jme, kms,kme, &
1601 its,ite, jts,jte, kts,kte
1603 INTEGER, INTENT(IN) :: num_soil_layers
1605 REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
1607 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1608 INTENT(INOUT) :: SMOIS, &
1609 TSLB !STEMP
1611 REAL, DIMENSION( ims:ime, jms:jme ) , &
1612 INTENT(INOUT) :: SNOW, &
1613 SNOWC, &
1614 CANWAT, &
1615 SMSTAV, &
1616 SMSTOT, &
1617 SFCRUNOFF, &
1618 UDRUNOFF, &
1619 SFCEVP, &
1620 GRDFLX, &
1621 ACSNOW, &
1622 XICE, &
1623 VEGFRA, &
1624 TMN, &
1625 ACSNOM
1627 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1628 INTENT(INOUT) :: IVGTYP, &
1629 ISLTYP
1631 !
1633 INTEGER, INTENT(IN) :: isn
1634 LOGICAL, INTENT(IN) :: allowed_to_read
1635 ! Local
1636 INTEGER :: iseason
1637 INTEGER :: icm,jcm,itf,jtf
1638 INTEGER :: I,J,L
1641 itf=min0(ite,ide-1)
1642 jtf=min0(jte,jde-1)
1644 icm = ide/2
1645 jcm = jde/2
1647 iseason=isn
1649 DO J=jts,jtf
1650 DO I=its,itf
1651 ! SNOW(i,j)=0.
1652 SNOWC(i,j)=0.
1653 ! SMSTAV(i,j)=
1654 ! SMSTOT(i,j)=
1655 ! SFCRUNOFF(i,j)=
1656 ! UDRUNOFF(i,j)=
1657 ! GRDFLX(i,j)=
1658 ! ACSNOW(i,j)=
1659 ! ACSNOM(i,j)=
1660 ENDDO
1661 ENDDO
1663 END SUBROUTINE hwrfsfcinit
1665 END MODULE module_sf_gfdl