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