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