module_sf_myjsfc.F
References to this file elsewhere.
1 !----------------------------------------------------------------------
2 !
3 MODULE MODULE_SF_MYJSFC
4 !
5 !----------------------------------------------------------------------
6 !
7 USE MODULE_MODEL_CONSTANTS
8 USE MODULE_DM, ONLY : WRF_DM_MAXVAL
9 !
10 !----------------------------------------------------------------------
11 !
12 ! REFERENCES: Janjic (2002), NCEP Office Note 437
13 !
14 ! ABSTRACT:
15 ! MYJSFC GENERATES THE SURFACE EXCHANGE COEFFICIENTS FOR VERTICAL
16 ! TURBULENT EXCHANGE BASED UPON MONIN_OBUKHOV THEORY WITH
17 ! VARIOUS REFINEMENTS.
18 !
19 !----------------------------------------------------------------------
20 !
21 INTEGER :: ITRMX=5 ! Iteration count for sfc layer computations
22 !
23 REAL,PARAMETER :: VKARMAN=0.4
24 REAL,PARAMETER :: CAPA=R_D/CP,ELOCP=2.72E6/CP,RCAP=1./CAPA
25 REAL,PARAMETER :: GOCP02=G/CP*2.,GOCP10=G/CP*10.
26 ! REAL,PARAMETER :: EPSU2=1.E-4,EPSUST=0.07,EPSZT=1.E-5
27 ! ECMWF sets lower limit on ln(Z0h)=-20 --> EPSZT=2.e-9
28 REAL,PARAMETER :: EPSU2=1.E-6,EPSUST=1.e-9,EPSZT=1.E-28
29 REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86
30 REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC
31 REAL,PARAMETER :: BETA=1./273.,CZIL=0.1,EXCML=0.0001,EXCMS=0.0001 &
32 !old REAL,PARAMETER :: BETA=1./273.,CZIL=0.1,EXCML=0.001,EXCMS=0.001 &
33 & ,GLKBR=10.,GLKBS=30.,PI=3.1415926 &
34 & ,QVISC=2.1E-5,RIC=0.505,SMALL=0.35 &
35 & ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5 &
36 & ,USTC=0.7,USTR=0.225,VISC=1.5E-5,FH=1.01 &
37 & ,WWST=1.2,ZTFC=1.,TOPOFAC=9.0e-6
38 !
39 REAL,PARAMETER :: BTG=BETA*G,CZIV=SMALL*GLKBS &
40 & ,GRRS=GLKBR/GLKBS &
41 & ,RTVISC=1./TVISC,RVISC=1./VISC &
42 & ,ZQRZT=SQSC/SQPR &
43 & ,FZQ1=RTVISC*QVISC*ZQRZT &
44 & ,FZQ2=RTVISC*QVISC*ZQRZT &
45 & ,FZT1=RVISC *TVISC*SQPR &
46 & ,FZT2=CZIV*GRRS*TVISC*SQPR &
47 & ,FZU1=CZIV*VISC &
48 & ,PIHF=0.5*PI &
49 & ,RQVISC=1./QVISC &
50 & ,USTFC=0.018/G &
51 & ,WWST2=WWST*WWST &
52 & ,ZILFC=-CZIL*VKARMAN*SQVISC
53 !
54 !----------------------------------------------------------------------
55 INTEGER, PARAMETER :: KZTM=10001,KZTM2=KZTM-2
56 !
57 REAL :: DZETA1,DZETA2,FH01,FH02,ZTMAX1,ZTMAX2,ZTMIN1,ZTMIN2
58 !
59 REAL,DIMENSION(KZTM) :: PSIH1,PSIH2,PSIM1,PSIM2
60 !
61 !----------------------------------------------------------------------
62 !
63 CONTAINS
64 !
65 !----------------------------------------------------------------------
66 SUBROUTINE MYJSFC(ITIMESTEP,HT,DZ &
67 & ,PMID,PINT,TH,T,QV,QC,U,V,Q2 &
68 & ,TSK,QSFC,THZ0,QZ0,UZ0,VZ0 &
69 & ,LOWLYR,XLAND &
70 & ,USTAR,ZNT,Z0BASE,PBLH,MAVAIL,RMOL &
71 & ,AKHS,AKMS &
72 & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC &
73 & ,QGH,CPM,CT &
74 & ,U10,V10,T02,TH02,TSHLTR,TH10,Q02,QSHLTR,Q10,PSHLTR &
75 & ,P1000mb &
76 & ,IDS,IDE,JDS,JDE,KDS,KDE &
77 & ,IMS,IME,JMS,JME,KMS,KME &
78 & ,ITS,ITE,JTS,JTE,KTS,KTE)
79 !----------------------------------------------------------------------
80 !
81 IMPLICIT NONE
82 !
83 !----------------------------------------------------------------------
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 !
88 INTEGER,INTENT(IN) :: ITIMESTEP
89 !
90 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
91 !
92 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,MAVAIL,TSK &
93 & ,XLAND,Z0BASE
94 !
95 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ &
96 & ,PMID,PINT &
97 & ,Q2,QC,QV &
98 & ,T,TH &
99 & ,U,V
100 !
101 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: FLX_LH,HFX,PSHLTR &
102 & ,QFX,Q10,QSHLTR &
103 & ,TH10,TSHLTR,T02 &
104 & ,U10,V10,TH02,Q02
105 !
106 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS &
107 & ,PBLH,QSFC
108 !
109 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: QZ0,RMOL,THZ0 &
110 & ,USTAR,UZ0,VZ0 &
111 & ,ZNT
112 !
113 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CHS,CHS2,CQS2 &
114 & ,CPM,CT,FLHC,FLQC &
115 & ,QGH
116 !
117 REAL,INTENT(IN) :: P1000mb
118 !----------------------------------------------------------------------
119 !***
120 !*** LOCAL VARIABLES
121 !***
122 INTEGER :: I,J,K,KFLIP,LMH,LPBL,NTSD
123 !
124 REAL :: A,APESFC,B,BTGX,CWMLOW &
125 & ,DQDT,DTDIF,DTDT,DUDT,DVDT &
126 & ,FIS &
127 & ,P02P,P10P,PLOW,PSFC,PTOP,QLOW,QS02,QS10 &
128 & ,RAPA,RAPA02,RAPA10,RATIOMX,RDZ,SEAMASK,SM &
129 & ,T02P,T10P,TEM,TH02P,TH10P,THLOW,THELOW,THM &
130 & ,TLOW,TZ0,ULOW,VLOW,ZSL
131 !
132 REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,THK,TK,UK,VK
133 !
134 REAL,DIMENSION(KTS:KTE-1) :: EL,ELM
135 !
136 REAL,DIMENSION(KTS:KTE+1) :: ZHK
137 !
138 REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
139 !
140 REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
141 !
142 !----------------------------------------------------------------------
143 !**********************************************************************
144 !----------------------------------------------------------------------
145 !
146 !*** MAKE PREPARATIONS
147 !
148 !----------------------------------------------------------------------
149 DO J=JTS,JTE
150 DO K=KTS,KTE+1
151 DO I=ITS,ITE
152 ZINT(I,K,J)=0.
153 ENDDO
154 ENDDO
155 ENDDO
156 !
157 DO J=JTS,JTE
158 DO I=ITS,ITE
159 ZINT(I,KTE+1,J)=HT(I,J) ! Z at bottom of lowest sigma layer
160 PBLH(I,J)=-1.
161 !
162 !!!!!!!!!
163 !!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
164 !!!!!!!!!
165 !!!!!! ZINT(I,KTE+1,J)=1.E-4 ! Z of bottom of lowest eta layer
166 !!!!!! ZHK(KTE+1)=1.E-4 ! Z of bottom of lowest eta layer
167 !
168 ENDDO
169 ENDDO
170 !
171 DO J=JTS,JTE
172 DO K=KTE,KTS,-1
173 KFLIP=KTE+1-K
174 DO I=ITS,ITE
175 ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
176 ENDDO
177 ENDDO
178 ENDDO
179 !
180 NTSD=ITIMESTEP
181 !
182 #if ( NMM_CORE == 1 )
183 if(NTSD+1.eq.1) then
184 #else
185 IF(NTSD==1)THEN
186 #endif
187 !tgs IF(NTSD==1)THEN
188 DO J=JTS,JTE
189 DO I=ITS,ITE
190 USTAR(I,J)=0.1
191 FIS=HT(I,J)*G
192 SM=XLAND(I,J)-1.
193 !!! Z0(I,J)=SM*Z0SEA+(1.-SM)*(Z0(I,J)*Z0MAX+FIS*FCM+Z0LAND)
194 !!! ZNT(I,J)=SM*Z0SEA+(1.-SM)*(ZNT(I,J)*Z0MAX+FIS*FCM+Z0LAND)
195 ENDDO
196 ENDDO
197 ENDIF
198 !
199 !!!! IF(NTSD==1)THEN
200 DO J=JTS,JTE
201 DO I=ITS,ITE
202 CT(I,J)=0.
203 ENDDO
204 ENDDO
205 !!!! ENDIF
206 !
207 !----------------------------------------------------------------------
208 setup_integration: DO J=JTS,JTE
209 !----------------------------------------------------------------------
210 !
211 DO I=ITS,ITE
212 !
213 !*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
214 !
215 LMH=KTE-LOWLYR(I,J)+1
216 !
217 PTOP=PINT(I,KTE+1,J) ! KTE+1=KME
218 PSFC=PINT(I,LOWLYR(I,J),J)
219 ! Define THSK here (for first timestep mostly)
220 ! THSK(I,J)=TSK(I,J)/(PSFC*1.E-5)**CAPA
221 THSK(I,J)=TSK(I,J)/(PSFC/P1000mb)**CAPA
222 !
223 !*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
224 !
225 SEAMASK=XLAND(I,J)-1.
226 !
227 !*** FILL 1-D VERTICAL ARRAYS
228 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
229 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
230 !
231 DO K=KTE,KTS,-1
232 KFLIP=KTE+1-K
233 THK(K)=TH(I,KFLIP,J)
234 TK(K)=T(I,KFLIP,J)
235 RATIOMX=QV(I,KFLIP,J)
236 QK(K)=RATIOMX/(1.+RATIOMX)
237 PK(K)=PMID(I,KFLIP,J)
238 CWMK(K)=QC(I,KFLIP,J)
239 THEK(K)=(CWMK(K)*(-ELOCP/TK(K))+1.)*THK(K)
240 Q2K(K)=2.*Q2(I,KFLIP,J)
241 !
242 !
243 !*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
244 !
245 ZHK(K)=ZINT(I,K,J)
246 !
247 ENDDO
248 ZHK(KTE+1)=HT(I,J) ! Z at bottom of lowest sigma layer
249 !
250 DO K=KTE,KTS,-1
251 KFLIP=KTE+1-K
252 UK(K)=U(I,KFLIP,J)
253 VK(K)=V(I,KFLIP,J)
254 ENDDO
255 !
256 !*** FIND THE HEIGHT OF THE PBL
257 !
258 LPBL=LMH
259 DO K=LMH-1,1,-1
260 IF(Q2K(K)<=EPSQ2*FH) THEN
261 LPBL=K
262 GO TO 110
263 ENDIF
264 ENDDO
265 !
266 LPBL=1
267 !
268 !-----------------------------------------------------------------------
269 !--------------THE HEIGHT OF THE PBL------------------------------------
270 !-----------------------------------------------------------------------
271 !
272 110 PBLH(I,J)=ZHK(LPBL)-ZHK(LMH+1)
273 !
274 !----------------------------------------------------------------------
275 !***
276 !*** FIND THE SURFACE EXCHANGE COEFFICIENTS
277 !***
278 !----------------------------------------------------------------------
279 PLOW=PK(LMH)
280 TLOW=TK(LMH)
281 THLOW=THK(LMH)
282 THELOW=THEK(LMH)
283 QLOW=QK(LMH)
284 CWMLOW=CWMK(LMH)
285 ULOW=UK(LMH)
286 VLOW=VK(LMH)
287 ZSL=(ZHK(LMH)-ZHK(LMH+1))*0.5
288 ! APESFC=(PSFC*1.E-5)**CAPA
289 APESFC=(PSFC/P1000mb)**CAPA
290 !tgs - in ARW THZ0 is not initialized when MYJSFC is called first time
291 #if ( NMM_CORE == 1 )
292 if(NTSD+1.eq.1) then
293 #else
294 IF(NTSD==1)THEN
295 #endif
296 ! if(itimestep.le.1) then
297 TZ0=TSK(I,J)
298 else
299 TZ0=THZ0(I,J)*APESFC
300 endif
301 !
302 CALL SFCDIF(NTSD,SEAMASK,THSK(I,J),QSFC(I,J),PSFC &
303 & ,UZ0(I,J),VZ0(I,J),TZ0,THZ0(I,J),QZ0(I,J) &
304 & ,USTAR(I,J),ZNT(I,J),Z0BASE(I,J),CT(I,J),RMOL(I,J) &
305 & ,AKMS(I,J),AKHS(I,J),PBLH(I,J),MAVAIL(I,J) &
306 & ,CHS(I,J),CHS2(I,J),CQS2(I,J) &
307 & ,HFX(I,J),QFX(I,J),FLX_LH(I,J) &
308 & ,FLHC(I,J),FLQC(I,J),QGH(I,J),CPM(I,J) &
309 & ,ULOW,VLOW,TLOW,THLOW,THELOW,QLOW,CWMLOW &
310 & ,ZSL,PLOW &
311 & ,U10(I,J),V10(I,J),TSHLTR(I,J),TH10(I,J) &
312 & ,QSHLTR(I,J),Q10(I,J),PSHLTR(I,J) &
313 & ,IDS,IDE,JDS,JDE,KDS,KDE &
314 & ,IMS,IME,JMS,JME,KMS,KME &
315 & ,ITS,ITE,JTS,JTE,KTS,KTE,i,j,ZHK(LMH+1))
316 !
317 !*** REMOVE SUPERATURATION AT 2M AND 10M
318 !
319 RAPA=APESFC
320 TH02P=TSHLTR(I,J)
321 TH10P=TH10(I,J)
322 TH02(I,J)=TSHLTR(I,J)
323 !
324 RAPA02=RAPA-GOCP02/TH02P
325 RAPA10=RAPA-GOCP10/TH10P
326 !
327 T02P=TH02P*RAPA02
328 T10P=TH10P*RAPA10
329 ! 1 may 06 tgs T02(I,J) = T02P
330 T02(I,J) = TH02(I,J)*APESFC
331 !
332 ! P02P=(RAPA02**RCAP)*1.E5
333 ! P10P=(RAPA10**RCAP)*1.E5
334 P02P=(RAPA02**RCAP)*P1000mb
335 P10P=(RAPA10**RCAP)*P1000mb
336 !
337 QS02=PQ0/P02P*EXP(A2*(T02P-A3)/(T02P-A4))
338 QS10=PQ0/P10P*EXP(A2*(T10P-A3)/(T10P-A4))
339 !
340 IF(QSHLTR(I,J)>QS02)QSHLTR(I,J)=QS02
341 IF(Q10 (I,J)>QS10)Q10 (I,J)=QS10
342 Q02(I,J)=QSHLTR(I,J)/(1.-QSHLTR(I,J))
343 !----------------------------------------------------------------------
344 !
345 ENDDO
346 !
347 !----------------------------------------------------------------------
348 ENDDO setup_integration
349 !----------------------------------------------------------------------
350
351 END SUBROUTINE MYJSFC
352 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
353 !----------------------------------------------------------------------
354 SUBROUTINE SFCDIF(NTSD,SEAMASK,THS,QS,PSFC &
355 & ,UZ0,VZ0,TZ0,THZ0,QZ0 &
356 & ,USTAR,Z0,Z0BASE,CT,RLMO,AKMS,AKHS,PBLH,WETM &
357 & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC,QGH,CPM &
358 & ,ULOW,VLOW,TLOW,THLOW,THELOW,QLOW,CWMLOW &
359 & ,ZSL,PLOW &
360 & ,U10,V10,TH02,TH10,Q02,Q10,PSHLTR &
361 & ,IDS,IDE,JDS,JDE,KDS,KDE &
362 & ,IMS,IME,JMS,JME,KMS,KME &
363 & ,ITS,ITE,JTS,JTE,KTS,KTE,i,j,ZSFC)
364 ! ****************************************************************
365 ! * *
366 ! * SURFACE LAYER *
367 ! * *
368 ! ****************************************************************
369 !----------------------------------------------------------------------
370 !
371 IMPLICIT NONE
372 !
373 !----------------------------------------------------------------------
374 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
375 & ,IMS,IME,JMS,JME,KMS,KME &
376 & ,ITS,ITE,JTS,JTE,KTS,KTE,i,j
377 !
378 INTEGER,INTENT(IN) :: NTSD
379 !
380 REAL,INTENT(IN) :: CWMLOW,PBLH,PLOW,QLOW,PSFC,SEAMASK,ZSFC &
381 & ,THELOW,THLOW,THS,TLOW,TZ0,ULOW,VLOW,WETM,ZSL &
382 & ,Z0BASE
383 !
384 REAL,INTENT(OUT) :: CHS,CHS2,CPM,CQS2,CT,FLHC,FLQC,FLX_LH,HFX &
385 & ,PSHLTR,Q02,Q10,QFX,QGH,RLMO,TH02,TH10,U10,V10
386 !
387 REAL,INTENT(INOUT) :: AKHS,AKMS,QZ0,THZ0,USTAR,UZ0,VZ0,Z0,QS
388 !----------------------------------------------------------------------
389 !***
390 !*** LOCAL VARIABLES
391 !***
392 INTEGER :: ITR,K
393 !
394 REAL :: A,B,BTGH,BTGX,CXCHL,CXCHS,DTHV,DU2,ELFC,FCT &
395 & ,HLFLX,HSFLX,HV,PSH02,PSH10,PSHZ,PSHZL,PSM10,PSMZ,PSMZL &
396 & ,RDZ,RDZT,RIB,RLMA,RLMN,RLMP &
397 & ,RLOGT,RLOGU,RWGH,RZ,RZST,RZSU,SIMH,SIMM,TEM,THM &
398 & ,UMFLX,USTARK,VMFLX,WGHT,WGHTT,WGHTQ,WSTAR2 &
399 & ,X,XLT,XLT4,XLU,XLU4,XT,XT4,XU,XU4,ZETALT,ZETALU &
400 & ,ZETAT,ZETAU,ZQ,ZSLT,ZSLU,ZT,ZU,TOPOTERM,ZZIL
401 !
402 !*** DIAGNOSTICS
403 !
404 REAL :: AKHS02,AKHS10,AKMS02,AKMS10,EKMS10,QSAT10,QSAT2 &
405 & ,RLNT02,RLNT10,RLNU10,SIMH02,SIMH10,SIMM10,T02,T10 &
406 & ,TERM1,RLOW,U10E,V10E,WSTAR,XLT02,XLT024,XLT10 &
407 & ,XLT104,XLU10,XLU104,XU10,XU104,ZT02,ZT10,ZTAT02,ZTAT10 &
408 & ,ZTAU,ZTAU10,ZU10,ZUUZ
409 !----------------------------------------------------------------------
410 !**********************************************************************
411 !----------------------------------------------------------------------
412 RDZ=1./ZSL
413 CXCHL=EXCML*RDZ
414 CXCHS=EXCMS*RDZ
415 !
416 BTGX=G/THLOW
417 ELFC=VKARMAN*BTGX
418 !
419 IF(PBLH>1000.)THEN
420 BTGH=BTGX*PBLH
421 ELSE
422 BTGH=BTGX*1000.
423 ENDIF
424 !
425 !----------------------------------------------------------------------
426 !
427 !*** SEA POINTS
428 !
429 !----------------------------------------------------------------------
430 !
431 IF(SEAMASK>0.5)THEN
432 !
433 !----------------------------------------------------------------------
434 DO ITR=1,ITRMX
435 !----------------------------------------------------------------------
436 Z0=MAX(USTFC*USTAR*USTAR,1.59E-5)
437 !
438 !*** VISCOUS SUBLAYER, JANJIC MWR 1994
439 !
440 !----------------------------------------------------------------------
441 IF(USTAR<USTC)THEN
442 !----------------------------------------------------------------------
443 !
444 IF(USTAR<USTR)THEN
445 !
446 #if ( NMM_CORE == 1 )
447 if(NTSD+1.eq.1) then
448 #else
449 IF(NTSD==1)THEN
450 #endif
451 !tgs IF(NTSD==1)THEN
452 AKMS=CXCHS
453 AKHS=CXCHS
454 QS=QLOW
455 ENDIF
456 !
457 ZU=FZU1*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR
458 WGHT=AKMS*ZU*RVISC
459 RWGH=WGHT/(WGHT+1.)
460 UZ0=(ULOW*RWGH+UZ0)*0.5
461 VZ0=(VLOW*RWGH+VZ0)*0.5
462 !
463 ZT=FZT1*ZU
464 ZQ=FZQ1*ZT
465 WGHTT=AKHS*ZT*RTVISC
466 WGHTQ=AKHS*ZQ*RQVISC
467 !
468 #if ( NMM_CORE == 1 )
469 if(NTSD+1>1) then
470 #else
471 IF(NTSD>1)THEN
472 #endif
473 !tgs IF(NTSD>1)THEN
474 THZ0=((WGHTT*THLOW+THS)/(WGHTT+1.)+THZ0)*0.5
475 QZ0=((WGHTQ*QLOW+QS)/(WGHTQ+1.)+QZ0)*0.5
476 ELSE
477 THZ0=(WGHTT*THLOW+THS)/(WGHTT+1.)
478 QZ0=(WGHTQ*QLOW+QS)/(WGHTQ+1.)
479 ENDIF
480 !
481 ENDIF
482 !
483 IF(USTAR>=USTR.AND.USTAR<USTC)THEN
484 ZU=Z0
485 UZ0=0.
486 VZ0=0.
487 !
488 ZT=FZT2*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR
489 ZQ=FZQ2*ZT
490 WGHTT=AKHS*ZT*RTVISC
491 WGHTQ=AKHS*ZQ*RQVISC
492 !
493 #if ( NMM_CORE == 1 )
494 if(NTSD+1>1) then
495 #else
496 IF(NTSD>1)THEN
497 #endif
498
499 !tgs IF(NTSD>1)THEN
500 THZ0=((WGHTT*THLOW+THS)/(WGHTT+1.)+THZ0)*0.5
501 QZ0=((WGHTQ*QLOW+QS)/(WGHTQ+1.)+QZ0)*0.5
502 ELSE
503 THZ0=(WGHTT*THLOW+THS)/(WGHTT+1.)
504 QZ0=(WGHTQ*QLOW+QS)/(WGHTQ+1.)
505 ENDIF
506 !
507 ENDIF
508 !----------------------------------------------------------------------
509 ELSE
510 !----------------------------------------------------------------------
511 ZU=Z0
512 UZ0=0.
513 VZ0=0.
514 !
515 ZT=Z0
516 THZ0=THS
517 !
518 ZQ=Z0
519 QZ0=QS
520 !----------------------------------------------------------------------
521 ENDIF
522 !----------------------------------------------------------------------
523 TEM=(TLOW+TZ0)*0.5
524 THM=(THELOW+THZ0)*0.5
525 !
526 A=THM*P608
527 B=(ELOCP/TEM-1.-P608)*THM
528 !
529 DTHV=((THELOW-THZ0)*((QLOW+QZ0+CWMLOW)*(0.5*P608)+1.) &
530 & +(QLOW-QZ0+CWMLOW)*A+CWMLOW*B)
531 !
532 DU2=MAX((ULOW-UZ0)**2+(VLOW-VZ0)**2,EPSU2)
533 RIB=BTGX*DTHV*ZSL/DU2
534 !----------------------------------------------------------------------
535 ! IF(RIB>=RIC)THEN
536 !----------------------------------------------------------------------
537 ! AKMS=MAX( VISC*RDZ,CXCHS)
538 ! AKHS=MAX(TVISC*RDZ,CXCHS)
539 !----------------------------------------------------------------------
540 ! ELSE ! turbulent branch
541 !----------------------------------------------------------------------
542 ZSLU=ZSL+ZU
543 ZSLT=ZSL+ZT
544 !
545 RZSU=ZSLU/ZU
546 RZST=ZSLT/ZT
547 !
548 RLOGU=LOG(RZSU)
549 RLOGT=LOG(RZST)
550 !
551 !----------------------------------------------------------------------
552 !*** 1./MONIN-OBUKHOV LENGTH
553 !----------------------------------------------------------------------
554 !
555 RLMO=ELFC*AKHS*DTHV/USTAR**3
556 !
557 ZETALU=ZSLU*RLMO
558 ZETALT=ZSLT*RLMO
559 ZETAU=ZU*RLMO
560 ZETAT=ZT*RLMO
561 !
562 ZETALU=MIN(MAX(ZETALU,ZTMIN1),ZTMAX1)
563 ZETALT=MIN(MAX(ZETALT,ZTMIN1),ZTMAX1)
564 ZETAU=MIN(MAX(ZETAU,ZTMIN1/RZSU),ZTMAX1/RZSU)
565 ZETAT=MIN(MAX(ZETAT,ZTMIN1/RZST),ZTMAX1/RZST)
566 !
567 !----------------------------------------------------------------------
568 !*** WATER FUNCTIONS
569 !----------------------------------------------------------------------
570 !
571 RZ=(ZETAU-ZTMIN1)/DZETA1
572 K=INT(RZ)
573 RDZT=RZ-REAL(K)
574 K=MIN(K,KZTM2)
575 K=MAX(K,0)
576 PSMZ=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1)
577 !
578 RZ=(ZETALU-ZTMIN1)/DZETA1
579 K=INT(RZ)
580 RDZT=RZ-REAL(K)
581 K=MIN(K,KZTM2)
582 K=MAX(K,0)
583 PSMZL=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1)
584 !
585 SIMM=PSMZL-PSMZ+RLOGU
586 !
587 RZ=(ZETAT-ZTMIN1)/DZETA1
588 K=INT(RZ)
589 RDZT=RZ-REAL(K)
590 K=MIN(K,KZTM2)
591 K=MAX(K,0)
592 PSHZ=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
593 !
594 RZ=(ZETALT-ZTMIN1)/DZETA1
595 K=INT(RZ)
596 RDZT=RZ-REAL(K)
597 K=MIN(K,KZTM2)
598 K=MAX(K,0)
599 PSHZL=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
600 !
601 SIMH=(PSHZL-PSHZ+RLOGT)*FH01
602 !----------------------------------------------------------------------
603 USTARK=USTAR*VKARMAN
604 AKMS=MAX(USTARK/SIMM,CXCHS)
605 AKHS=MAX(USTARK/SIMH,CXCHS)
606 !
607 !----------------------------------------------------------------------
608 !*** BELJAARS CORRECTION FOR USTAR
609 !----------------------------------------------------------------------
610 !
611 IF(DTHV<=0.)THEN !zj
612 WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.) !zj
613 ELSE !zj
614 WSTAR2=0. !zj
615 ENDIF !zj
616 USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
617 !
618 !----------------------------------------------------------------------
619 ! ENDIF ! End of turbulent branch
620 !----------------------------------------------------------------------
621 !
622 ENDDO ! End of the iteration loop over sea points
623 !
624 !----------------------------------------------------------------------
625 !
626 !*** LAND POINTS
627 !
628 !----------------------------------------------------------------------
629 !
630 ELSE
631 !
632 !----------------------------------------------------------------------
633 IF(NTSD==1)THEN
634 QS=QLOW
635 ENDIF
636 !
637 ZU=Z0
638 UZ0=0.
639 VZ0=0.
640 !
641 ZT=ZU*ZTFC
642 THZ0=THS
643 !
644 ZQ=ZT
645 QZ0=QS
646 !----------------------------------------------------------------------
647 TEM=(TLOW+TZ0)*0.5
648 THM=(THELOW+THZ0)*0.5
649 !
650 A=THM*P608
651 B=(ELOCP/TEM-1.-P608)*THM
652 !
653 DTHV=((THELOW-THZ0)*((QLOW+QZ0+CWMLOW)*(0.5*P608)+1.) &
654 & +(QLOW-QZ0+CWMLOW)*A+CWMLOW*B)
655 !
656 DU2=MAX((ULOW-UZ0)**2+(VLOW-VZ0)**2,EPSU2)
657 RIB=BTGX*DTHV*ZSL/DU2
658 !----------------------------------------------------------------------
659 ! IF(RIB>=RIC)THEN
660 ! AKMS=MAX( VISC*RDZ,CXCHL)
661 ! AKHS=MAX(TVISC*RDZ,CXCHL)
662 !----------------------------------------------------------------------
663 ! ELSE ! Turbulent branch
664 !----------------------------------------------------------------------
665 ZSLU=ZSL+ZU
666 !
667 RZSU=ZSLU/ZU
668 !
669 RLOGU=LOG(RZSU)
670
671 ZSLT=ZSL+ZU ! u,v and t are at the same level
672 !----------------------------------------------------------------------
673 !
674 !mp Topo modification of ZILFC term
675 !
676 TOPOTERM=TOPOFAC*ZSFC**2.
677 TOPOTERM=MAX(TOPOTERM,3.0)
678 !
679 IF(DTHV>0.)THEN
680 ZZIL=ZILFC*TOPOTERM
681 ELSE
682 ZZIL=ZILFC
683 ENDIF
684 !
685 !----------------------------------------------------------------------
686 !
687 land_point_iteration: DO ITR=1,ITRMX
688 !
689 !----------------------------------------------------------------------
690 !*** ZILITINKEVITCH FIX FOR ZT
691 !----------------------------------------------------------------------
692 !
693 ! oldform ZT=MAX(EXP(ZZIL*SQRT(USTAR*ZU))*ZU,EPSZT)
694 ZT=MAX(EXP(ZZIL*SQRT(USTAR*Z0BASE))*Z0BASE,EPSZT)
695 !
696 RZST=ZSLT/ZT
697 RLOGT=LOG(RZST)
698 !
699 !----------------------------------------------------------------------
700 !*** 1./MONIN-OBUKHOV LENGTH-SCALE
701 !----------------------------------------------------------------------
702 !
703 RLMO=ELFC*AKHS*DTHV/USTAR**3
704 ZETALU=ZSLU*RLMO
705 ZETALT=ZSLT*RLMO
706 ZETAU=ZU*RLMO
707 ZETAT=ZT*RLMO
708 !
709 ZETALU=MIN(MAX(ZETALU,ZTMIN2),ZTMAX2)
710 ZETALT=MIN(MAX(ZETALT,ZTMIN2),ZTMAX2)
711 ZETAU=MIN(MAX(ZETAU,ZTMIN2/RZSU),ZTMAX2/RZSU)
712 ZETAT=MIN(MAX(ZETAT,ZTMIN2/RZST),ZTMAX2/RZST)
713 !
714 !----------------------------------------------------------------------
715 !*** LAND FUNCTIONS
716 !----------------------------------------------------------------------
717 !
718 RZ=(ZETAU-ZTMIN2)/DZETA2
719 K=INT(RZ)
720 RDZT=RZ-REAL(K)
721 K=MIN(K,KZTM2)
722 K=MAX(K,0)
723 PSMZ=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
724 !
725 RZ=(ZETALU-ZTMIN2)/DZETA2
726 K=INT(RZ)
727 RDZT=RZ-REAL(K)
728 K=MIN(K,KZTM2)
729 K=MAX(K,0)
730 PSMZL=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
731 !
732 SIMM=PSMZL-PSMZ+RLOGU
733 !
734 RZ=(ZETAT-ZTMIN2)/DZETA2
735 K=INT(RZ)
736 RDZT=RZ-REAL(K)
737 K=MIN(K,KZTM2)
738 K=MAX(K,0)
739 PSHZ=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
740 !
741 RZ=(ZETALT-ZTMIN2)/DZETA2
742 K=INT(RZ)
743 RDZT=RZ-REAL(K)
744 K=MIN(K,KZTM2)
745 K=MAX(K,0)
746 PSHZL=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
747 !
748 SIMH=(PSHZL-PSHZ+RLOGT)*FH02
749 !----------------------------------------------------------------------
750 USTARK=USTAR*VKARMAN
751 AKMS=MAX(USTARK/SIMM,CXCHL)
752 AKHS=MAX(USTARK/SIMH,CXCHL)
753 !
754 !----------------------------------------------------------------------
755 !*** BELJAARS CORRECTION FOR USTAR
756 !----------------------------------------------------------------------
757 !
758 IF(DTHV<=0.)THEN !zj
759 WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.) !zj
760 ELSE !zj
761 WSTAR2=0. !zj
762 ENDIF !zj
763 !
764 USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
765 !
766 !----------------------------------------------------------------------
767 ENDDO land_point_iteration
768 !----------------------------------------------------------------------
769 !
770 ! ENDIF ! End of turbulant branch over land
771 !
772 !----------------------------------------------------------------------
773 !
774 ENDIF ! End of land/sea branch
775 !
776 !----------------------------------------------------------------------
777 !*** COUNTERGRADIENT FIX
778 !----------------------------------------------------------------------
779 !
780 ! HV=-AKHS*DTHV
781 ! IF(HV>0.)THEN
782 ! FCT=-10.*(BTGX)**(-1./3.)
783 ! CT=FCT*(HV/(PBLH*PBLH))**(2./3.)
784 ! ELSE
785 CT=0.
786 ! ENDIF
787 !
788 !----------------------------------------------------------------------
789 !----------------------------------------------------------------------
790 !*** THE FOLLOWING DIAGNOSTIC BLOCK PRODUCES 2-m and 10-m VALUES
791 !*** FOR TEMPERATURE, MOISTURE, AND WINDS. IT IS DONE HERE SINCE
792 !*** THE VARIOUS QUANTITIES NEEDED FOR THE COMPUTATION ARE LOST
793 !*** UPON EXIT FROM THE ROTUINE.
794 !----------------------------------------------------------------------
795 !----------------------------------------------------------------------
796 !
797 WSTAR=SQRT(WSTAR2)/WWST
798 !
799 UMFLX=AKMS*(ULOW -UZ0 )
800 VMFLX=AKMS*(VLOW -VZ0 )
801 HSFLX=AKHS*(THLOW-THZ0)
802 HLFLX=AKHS*(QLOW -QZ0 )
803 !----------------------------------------------------------------------
804 ! IF(RIB>=RIC)THEN
805 !----------------------------------------------------------------------
806 ! IF(SEAMASK>0.5)THEN
807 ! AKMS10=MAX( VISC/10.,CXCHS)
808 ! AKHS02=MAX(TVISC/02.,CXCHS)
809 ! AKHS10=MAX(TVISC/10.,CXCHS)
810 ! ELSE
811 ! AKMS10=MAX( VISC/10.,CXCHL)
812 ! AKHS02=MAX(TVISC/02.,CXCHL)
813 ! AKHS10=MAX(TVISC/10.,CXCHL)
814 ! ENDIF
815 !----------------------------------------------------------------------
816 ! ELSE
817 !----------------------------------------------------------------------
818 ZU10=ZU+10.
819 ZT02=ZT+02.
820 ZT10=ZT+10.
821 !
822 RLNU10=LOG(ZU10/ZU)
823 RLNT02=LOG(ZT02/ZT)
824 RLNT10=LOG(ZT10/ZT)
825 !
826 ZTAU10=ZU10*RLMO
827 ZTAT02=ZT02*RLMO
828 ZTAT10=ZT10*RLMO
829 !
830 !----------------------------------------------------------------------
831 !*** SEA
832 !----------------------------------------------------------------------
833 !
834 IF(SEAMASK>0.5)THEN
835 !
836 !----------------------------------------------------------------------
837 ZTAU10=MIN(MAX(ZTAU10,ZTMIN1),ZTMAX1)
838 ZTAT02=MIN(MAX(ZTAT02,ZTMIN1),ZTMAX1)
839 ZTAT10=MIN(MAX(ZTAT10,ZTMIN1),ZTMAX1)
840 !----------------------------------------------------------------------
841 RZ=(ZTAU10-ZTMIN1)/DZETA1
842 K=INT(RZ)
843 RDZT=RZ-REAL(K)
844 K=MIN(K,KZTM2)
845 K=MAX(K,0)
846 PSM10=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1)
847 !
848 SIMM10=PSM10-PSMZ+RLNU10
849 !
850 RZ=(ZTAT02-ZTMIN1)/DZETA1
851 K=INT(RZ)
852 RDZT=RZ-REAL(K)
853 K=MIN(K,KZTM2)
854 K=MAX(K,0)
855 PSH02=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
856 !
857 SIMH02=(PSH02-PSHZ+RLNT02)*FH01
858 !
859 RZ=(ZTAT10-ZTMIN1)/DZETA1
860 K=INT(RZ)
861 RDZT=RZ-REAL(K)
862 K=MIN(K,KZTM2)
863 K=MAX(K,0)
864 PSH10=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
865 !
866 SIMH10=(PSH10-PSHZ+RLNT10)*FH01
867 !
868 AKMS10=MAX(USTARK/SIMM10,CXCHS)
869 AKHS02=MAX(USTARK/SIMH02,CXCHS)
870 AKHS10=MAX(USTARK/SIMH10,CXCHS)
871 !
872 !----------------------------------------------------------------------
873 !*** LAND
874 !----------------------------------------------------------------------
875 !
876 ELSE
877 !
878 !----------------------------------------------------------------------
879 ZTAU10=MIN(MAX(ZTAU10,ZTMIN2),ZTMAX2)
880 ZTAT02=MIN(MAX(ZTAT02,ZTMIN2),ZTMAX2)
881 ZTAT10=MIN(MAX(ZTAT10,ZTMIN2),ZTMAX2)
882 !----------------------------------------------------------------------
883 RZ=(ZTAU10-ZTMIN2)/DZETA2
884 K=INT(RZ)
885 RDZT=RZ-REAL(K)
886 K=MIN(K,KZTM2)
887 K=MAX(K,0)
888 PSM10=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
889 !
890 SIMM10=PSM10-PSMZ+RLNU10
891 !
892 RZ=(ZTAT02-ZTMIN2)/DZETA2
893 K=INT(RZ)
894 RDZT=RZ-REAL(K)
895 K=MIN(K,KZTM2)
896 K=MAX(K,0)
897 PSH02=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
898 !
899 SIMH02=(PSH02-PSHZ+RLNT02)*FH02
900 !
901 RZ=(ZTAT10-ZTMIN2)/DZETA2
902 K=INT(RZ)
903 RDZT=RZ-REAL(K)
904 K=MIN(K,KZTM2)
905 K=MAX(K,0)
906 PSH10=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
907 !
908 SIMH10=(PSH10-PSHZ+RLNT10)*FH02
909 !
910 AKMS10=MAX(USTARK/SIMM10,CXCHL)
911 AKHS02=MAX(USTARK/SIMH02,CXCHL)
912 AKHS10=MAX(USTARK/SIMH10,CXCHL)
913 !----------------------------------------------------------------------
914 ENDIF
915 !----------------------------------------------------------------------
916 ! ENDIF
917 !----------------------------------------------------------------------
918 U10 =UMFLX/AKMS10+UZ0
919 V10 =VMFLX/AKMS10+VZ0
920 TH02=HSFLX/AKHS02+THZ0
921 !
922 !*** BE CERTAIN THAT THE 2-M THETA AND 10-M THETA ARE BRACKETED BY
923 !*** THE VALUES OF THZ0 AND THLOW.
924 !
925 IF(THLOW>THZ0.AND.(TH02<THZ0.OR.TH02>THLOW).OR. &
926 THLOW<THZ0.AND.(TH02>THZ0.OR.TH02<THLOW))THEN
927 TH02=THZ0+2.*RDZ*(THLOW-THZ0)
928 ENDIF
929 !
930 TH10=HSFLX/AKHS10+THZ0
931 !
932 IF(THLOW>THZ0.AND.(TH10<THZ0.OR.TH10>THLOW).OR. &
933 THLOW<THZ0.AND.(TH10>THZ0.OR.TH10<THLOW))THEN
934 TH10=THZ0+10.*RDZ*(THLOW-THZ0)
935 ENDIF
936 !
937 Q02 =HLFLX/AKHS02+QZ0
938 Q10 =HLFLX/AKHS10+QZ0
939 TERM1=-0.068283/TLOW
940 PSHLTR=PSFC*EXP(TERM1)
941 !
942 !----------------------------------------------------------------------
943 !*** COMPUTE "EQUIVALENT" Z0 TO APPROXIMATE LOCAL SHELTER READINGS.
944 !----------------------------------------------------------------------
945 !
946 U10E=U10
947 V10E=V10
948 !
949 IF(SEAMASK<0.5)THEN
950
951 !1st ZUUZ=MIN(0.5*ZU,0.1)
952 !1st ZU=MAX(0.1*ZU,ZUUZ)
953 !tst ZUUZ=amin1(ZU*0.50,0.3)
954 !tst ZU=amax1(ZU*0.3,ZUUZ)
955
956 ZUUZ=AMIN1(ZU*0.50,0.18)
957 ZU=AMAX1(ZU*0.35,ZUUZ)
958 !
959 ZU10=ZU+10.
960 RZSU=ZU10/ZU
961 RLNU10=LOG(RZSU)
962
963 ZETAU=ZU*RLMO
964 ZTAU10=ZU10*RLMO
965
966 ZTAU10=MIN(MAX(ZTAU10,ZTMIN2),ZTMAX2)
967 ZETAU=MIN(MAX(ZETAU,ZTMIN2/RZSU),ZTMAX2/RZSU)
968
969 RZ=(ZTAU10-ZTMIN2)/DZETA2
970 K=INT(RZ)
971 RDZT=RZ-REAL(K)
972 K=MIN(K,KZTM2)
973 K=MAX(K,0)
974 PSM10=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
975 SIMM10=PSM10-PSMZ+RLNU10
976 EKMS10=MAX(USTARK/SIMM10,CXCHL)
977
978 U10E=UMFLX/EKMS10+UZ0
979 V10E=VMFLX/EKMS10+VZ0
980
981 ENDIF
982 !
983 U10=U10E
984 V10=V10E
985 !
986 !----------------------------------------------------------------------
987 !*** SET OTHER WRF DRIVER ARRAYS
988 !----------------------------------------------------------------------
989 !
990 RLOW=PLOW/(R_D*TLOW)
991 CHS=AKHS
992 CHS2=AKHS02
993 CQS2=AKHS02
994 HFX=-RLOW*CP*HSFLX
995 QFX=-RLOW*HLFLX*WETM
996 FLX_LH=XLV*QFX
997 FLHC=RLOW*CP*AKHS
998 FLQC=RLOW*AKHS*WETM
999 !!! QGH=PQ0/PSHLTR*EXP(A2S*(TSK-A3S)/(TSK-A4S))
1000 QGH=((1.-SEAMASK)*PQ0+SEAMASK*PQ0SEA) &
1001 & /PLOW*EXP(A2S*(TLOW-A3S)/(TLOW-A4S))
1002 QGH=QGH/(1.-QGH) !Convert to mixing ratio
1003 CPM=CP*(1.+0.8*QLOW)
1004 !
1005 !*** DO NOT COMPUTE QS OVER LAND POINTS HERE SINCE IT IS
1006 !*** A PROGNOSTIC VARIABLE THERE. IT IS OKAY TO USE IT
1007 !*** AS A DIAGNOSTIC OVER WATER SINCE IT WILL CAUSE NO
1008 !*** INTERFERENCE BEFORE BEING RECOMPUTED IN MYJPBL.
1009 !
1010 IF(SEAMASK>0.5)THEN
1011 QS=QLOW+QFX/(RLOW*AKHS)
1012 QS=QS/(1.-QS)
1013 ENDIF
1014 !----------------------------------------------------------------------
1015 !
1016 END SUBROUTINE SFCDIF
1017 !
1018 !----------------------------------------------------------------------
1019 SUBROUTINE MYJSFCINIT(LOWLYR,USTAR,Z0 &
1020 & ,SEAMASK,XICE,IVGTYP,RESTART &
1021 & ,ALLOWED_TO_READ &
1022 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1023 & ,IMS,IME,JMS,JME,KMS,KME &
1024 & ,ITS,ITE,JTS,JTE,KTS,KTE)
1025 !----------------------------------------------------------------------
1026 IMPLICIT NONE
1027 !----------------------------------------------------------------------
1028 LOGICAL,INTENT(IN) :: RESTART,ALLOWED_TO_READ
1029 !
1030 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1031 & ,IMS,IME,JMS,JME,KMS,KME &
1032 & ,ITS,ITE,JTS,JTE,KTS,KTE
1033 !
1034 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IVGTYP
1035 !
1036 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: LOWLYR
1037 !
1038 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: SEAMASK,XICE
1039 !
1040 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: USTAR,Z0
1041 !
1042 REAL,DIMENSION(0:30) :: VZ0TBL
1043 REAL,DIMENSION(0:30) :: VZ0TBL_24
1044 !
1045 INTEGER :: I,IDUM,IRECV,J,JDUM,K,ITF,JTF,KTF,MAXGBL_IVGTYP &
1046 &, MAXLOC_IVGTYP,MPI_COMM_COMP
1047 !
1048 ! INTEGER :: MPI_INTEGER,MPI_MAX
1049 !
1050 REAL :: SM,X,ZETA1,ZETA2,ZRNG1,ZRNG2
1051 !
1052 REAL :: PIHF=3.1415926/2.,EPS=1.E-6
1053 !----------------------------------------------------------------------
1054 VZ0TBL= &
1055 & (/0., &
1056 & 2.653,0.826,0.563,1.089,0.854,0.856,0.035,0.238,0.065,0.076 &
1057 & ,0.011,0.035,0.011,0.000,0.000,0.000,0.000,0.000,0.000,0.000 &
1058 & ,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/)
1059
1060 VZ0TBL_24= (/0., &
1061 & 1.00, 0.07, 0.07, 0.07, 0.07, 0.15, &
1062 & 0.08, 0.03, 0.05, 0.86, 0.80, 0.85, &
1063 & 2.65, 1.09, 0.80, 0.001, 0.04, 0.05, &
1064 & 0.01, 0.04, 0.06, 0.05, 0.03, 0.001, &
1065 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 /)
1066
1067 !----------------------------------------------------------------------
1068 !
1069 JTF=MIN0(JTE,JDE-1)
1070 KTF=MIN0(KTE,KDE-1)
1071 ITF=MIN0(ITE,IDE-1)
1072 !
1073 !
1074 !*** FOR NOW, ASSUME SIGMA MODE FOR LOWEST MODEL LAYER
1075 !
1076 DO J=JTS,JTF
1077 DO I=ITS,ITF
1078 LOWLYR(I,J)=1
1079 ! USTAR(I,J)=EPSUST
1080 ENDDO
1081 ENDDO
1082 !----------------------------------------------------------------------
1083 #if (NMM_CORE == 1)
1084 !
1085 IF(.NOT.RESTART)THEN
1086 ! CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
1087 ! MAXLOC_IVGTYP=MAXVAL(IVGTYP)
1088 ! CALL MPI_ALLREDUCE(MAXLOC_IVGTYP,MAXGBL_IVGTYP,1,MPI_INTEGER &
1089 ! &, MPI_MAX,MPI_COMM_COMP,IRECV)
1090 MAXGBL_IVGTYP=MAXVAL(IVGTYP)
1091 CALL WRF_DM_MAXVAL(MAXGBL_IVGTYP,IDUM,JDUM)
1092 !
1093 IF (MAXGBL_IVGTYP<13) THEN
1094 DO J=JTS,JTE
1095 DO I=ITS,ITE
1096 SM=SEAMASK(I,J)-1.
1097 IF(SM+XICE(I,J)<0.5)THEN
1098 Z0(I,J)=VZ0TBL(IVGTYP(I,J))
1099 ENDIF
1100 ENDDO
1101 ENDDO
1102 !
1103 ELSE
1104 !
1105 DO J=JTS,JTE
1106 DO I=ITS,ITE
1107 SM=SEAMASK(I,J)-1.
1108 IF(SM+XICE(I,J)<0.5)THEN
1109 Z0(I,J)=VZ0TBL_24(IVGTYP(I,J))
1110 ENDIF
1111 ENDDO
1112 ENDDO
1113 !
1114 ENDIF ! Vegtype check
1115 !
1116 ENDIF ! Restart check
1117
1118 #endif
1119 !----------------------------------------------------------------------
1120 IF(.NOT.RESTART)THEN
1121 DO J=JTS,JTE
1122 DO I=ITS,ITF
1123 USTAR(I,J)=0.1
1124 ENDDO
1125 ENDDO
1126 ENDIF
1127 !----------------------------------------------------------------------
1128 !
1129 !*** COMPUTE SURFACE LAYER INTEGRAL FUNCTIONS
1130 !
1131 !----------------------------------------------------------------------
1132 FH01=1.
1133 FH02=1.
1134 !
1135 ! ZTMIN1=-10.0
1136 ! ZTMAX1=2.0
1137 ! ZTMIN2=-10.0
1138 ! ZTMAX2=2.0
1139 ZTMIN1=-5.0
1140 ZTMAX1=1.0
1141 ZTMIN2=-5.0
1142 ZTMAX2=1.0
1143 !
1144 ZRNG1=ZTMAX1-ZTMIN1
1145 ZRNG2=ZTMAX2-ZTMIN2
1146 !
1147 DZETA1=ZRNG1/(KZTM-1)
1148 DZETA2=ZRNG2/(KZTM-1)
1149 !
1150 !----------------------------------------------------------------------
1151 !*** FUNCTION DEFINITION LOOP
1152 !----------------------------------------------------------------------
1153 !
1154 ZETA1=ZTMIN1
1155 ZETA2=ZTMIN2
1156 !
1157 DO K=1,KZTM
1158 !
1159 !----------------------------------------------------------------------
1160 !*** UNSTABLE RANGE
1161 !----------------------------------------------------------------------
1162 !
1163 IF(ZETA1<0.)THEN
1164 !
1165 !----------------------------------------------------------------------
1166 !*** PAULSON 1970 FUNCTIONS
1167 !----------------------------------------------------------------------
1168 X=SQRT(SQRT(1.-16.*ZETA1))
1169 !
1170 PSIM1(K)=-2.*LOG((X+1.)/2.)-LOG((X*X+1.)/2.)+2.*ATAN(X)-PIHF
1171 PSIH1(K)=-2.*LOG((X*X+1.)/2.)
1172 !
1173 !----------------------------------------------------------------------
1174 !*** STABLE RANGE
1175 !----------------------------------------------------------------------
1176 !
1177 ELSE
1178 !
1179 !----------------------------------------------------------------------
1180 !*** PAULSON 1970 FUNCTIONS
1181 !----------------------------------------------------------------------
1182 !
1183 ! PSIM1(K)=5.*ZETA1
1184 ! PSIH1(K)=5.*ZETA1
1185 !----------------------------------------------------------------------
1186 !*** HOLTSLAG AND DE BRUIN 1988
1187 !----------------------------------------------------------------------
1188 !
1189 PSIM1(K)=0.7*ZETA1+0.75*ZETA1*(6.-0.35*ZETA1)*EXP(-0.35*ZETA1)
1190 PSIH1(K)=0.7*ZETA1+0.75*ZETA1*(6.-0.35*ZETA1)*EXP(-0.35*ZETA1)
1191 !----------------------------------------------------------------------
1192 !
1193 ENDIF
1194 !
1195 !----------------------------------------------------------------------
1196 !*** UNSTABLE RANGE
1197 !----------------------------------------------------------------------
1198 !
1199 IF(ZETA2<0.)THEN
1200 !
1201 !----------------------------------------------------------------------
1202 !*** PAULSON 1970 FUNCTIONS
1203 !----------------------------------------------------------------------
1204 !
1205 X=SQRT(SQRT(1.-16.*ZETA2))
1206 !
1207 PSIM2(K)=-2.*LOG((X+1.)/2.)-LOG((X*X+1.)/2.)+2.*ATAN(X)-PIHF
1208 PSIH2(K)=-2.*LOG((X*X+1.)/2.)
1209 !----------------------------------------------------------------------
1210 !*** STABLE RANGE
1211 !----------------------------------------------------------------------
1212 !
1213 ELSE
1214 !
1215 !----------------------------------------------------------------------
1216 !*** PAULSON 1970 FUNCTIONS
1217 !----------------------------------------------------------------------
1218 !
1219 ! PSIM2(K)=5.*ZETA2
1220 ! PSIH2(K)=5.*ZETA2
1221 !
1222 !----------------------------------------------------------------------
1223 !*** HOLTSLAG AND DE BRUIN 1988
1224 !----------------------------------------------------------------------
1225 !
1226 PSIM2(K)=0.7*ZETA2+0.75*ZETA2*(6.-0.35*ZETA2)*EXP(-0.35*ZETA2)
1227 PSIH2(K)=0.7*ZETA2+0.75*ZETA2*(6.-0.35*ZETA2)*EXP(-0.35*ZETA2)
1228 !----------------------------------------------------------------------
1229 !
1230 ENDIF
1231 !
1232 !----------------------------------------------------------------------
1233 IF(K==KZTM)THEN
1234 ZTMAX1=ZETA1
1235 ZTMAX2=ZETA2
1236 ENDIF
1237 !
1238 ZETA1=ZETA1+DZETA1
1239 ZETA2=ZETA2+DZETA2
1240 !----------------------------------------------------------------------
1241 ENDDO
1242 !----------------------------------------------------------------------
1243 ZTMAX1=ZTMAX1-EPS
1244 ZTMAX2=ZTMAX2-EPS
1245 !----------------------------------------------------------------------
1246 !
1247 END SUBROUTINE MYJSFCINIT
1248 !
1249 !----------------------------------------------------------------------
1250 !
1251 END MODULE MODULE_SF_MYJSFC
1252 !
1253 !----------------------------------------------------------------------