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