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