module_sf_sfclay.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_sf_sfclay
4
5 REAL , PARAMETER :: VCONVC=1.
6 REAL , PARAMETER :: CZO=0.0185
7 REAL , PARAMETER :: OZO=1.59E-5
8
9 REAL, DIMENSION(0:1000 ),SAVE :: PSIMTB,PSIHTB
10
11 CONTAINS
12
13 !-------------------------------------------------------------------
14 SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w, &
15 CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
16 ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
17 XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
18 U10,V10,TH2,T2,Q2, &
19 GZ1OZ0,WSPD,BR,ISFFLX,DX, &
20 SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
21 KARMAN,EOMEG,STBOLT, &
22 P1000mb, &
23 ids,ide, jds,jde, kds,kde, &
24 ims,ime, jms,jme, kms,kme, &
25 its,ite, jts,jte, kts,kte, &
26 uratx,vratx,tratx )
27 !-------------------------------------------------------------------
28 IMPLICIT NONE
29 !-------------------------------------------------------------------
30 !-- U3D 3D u-velocity interpolated to theta points (m/s)
31 !-- V3D 3D v-velocity interpolated to theta points (m/s)
32 !-- T3D temperature (K)
33 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
34 !-- P3D 3D pressure (Pa)
35 !-- dz8w dz between full levels (m)
36 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
37 !-- G acceleration due to gravity (m/s^2)
38 !-- ROVCP R/CP
39 !-- R gas constant for dry air (J/kg/K)
40 !-- XLV latent heat of vaporization for water (J/kg)
41 !-- PSFC surface pressure (Pa)
42 !-- ZNT roughness length (m)
43 !-- UST u* in similarity theory (m/s)
44 !-- PBLH PBL height from previous time (m)
45 !-- MAVAIL surface moisture availability (between 0 and 1)
46 !-- ZOL z/L height over Monin-Obukhov length
47 !-- MOL T* (similarity theory) (K)
48 !-- REGIME flag indicating PBL regime (stable, unstable, etc.)
49 !-- PSIM similarity stability function for momentum
50 !-- PSIH similarity stability function for heat
51 !-- XLAND land mask (1 for land, 2 for water)
52 !-- HFX upward heat flux at the surface (W/m^2)
53 !-- QFX upward moisture flux at the surface (kg/m^2/s)
54 !-- LH net upward latent heat flux at surface (W/m^2)
55 !-- TSK surface temperature (K)
56 !-- FLHC exchange coefficient for heat (W/m^2/K)
57 !-- FLQC exchange coefficient for moisture (kg/m^2/s)
58 !-- CHS heat/moisture exchange coefficient for LSM (m/s)
59 !-- QGH lowest-level saturated mixing ratio
60 !-- QSFC ground saturated mixing ratio
61 !-- uratx ratio of surface U to U10
62 !-- vratx ratio of surface V to V10
63 !-- tratx ratio of surface T to TH2
64 !-- U10 diagnostic 10m u wind
65 !-- V10 diagnostic 10m v wind
66 !-- TH2 diagnostic 2m theta (K)
67 !-- T2 diagnostic 2m temperature (K)
68 !-- Q2 diagnostic 2m mixing ratio (kg/kg)
69 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
70 !-- WSPD wind speed at lowest model level (m/s)
71 !-- BR bulk Richardson number in surface layer
72 !-- ISFFLX isfflx=1 for surface heat and moisture fluxes
73 !-- DX horizontal grid size (m)
74 !-- SVP1 constant for saturation vapor pressure (kPa)
75 !-- SVP2 constant for saturation vapor pressure (dimensionless)
76 !-- SVP3 constant for saturation vapor pressure (K)
77 !-- SVPT0 constant for saturation vapor pressure (K)
78 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
79 !-- EP2 constant for specific humidity calculation
80 ! (R_d/R_v) (dimensionless)
81 !-- KARMAN Von Karman constant
82 !-- EOMEG angular velocity of earth's rotation (rad/s)
83 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
84 !-- ids start index for i in domain
85 !-- ide end index for i in domain
86 !-- jds start index for j in domain
87 !-- jde end index for j in domain
88 !-- kds start index for k in domain
89 !-- kde end index for k in domain
90 !-- ims start index for i in memory
91 !-- ime end index for i in memory
92 !-- jms start index for j in memory
93 !-- jme end index for j in memory
94 !-- kms start index for k in memory
95 !-- kme end index for k in memory
96 !-- its start index for i in tile
97 !-- ite end index for i in tile
98 !-- jts start index for j in tile
99 !-- jte end index for j in tile
100 !-- kts start index for k in tile
101 !-- kte end index for k in tile
102 !-------------------------------------------------------------------
103 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
104 ims,ime, jms,jme, kms,kme, &
105 its,ite, jts,jte, kts,kte
106 !
107 INTEGER, INTENT(IN ) :: ISFFLX
108 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
109 REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
110 REAL, INTENT(IN ) :: P1000mb
111 !
112 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
113 INTENT(IN ) :: dz8w
114
115 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
116 INTENT(IN ) :: QV3D, &
117 P3D, &
118 T3D
119
120 REAL, DIMENSION( ims:ime, jms:jme ) , &
121 INTENT(IN ) :: MAVAIL, &
122 PBLH, &
123 XLAND, &
124 TSK
125 REAL, DIMENSION( ims:ime, jms:jme ) , &
126 INTENT(OUT ) :: U10, &
127 V10, &
128 TH2, &
129 T2, &
130 Q2, &
131 QSFC
132
133 !
134 REAL, DIMENSION( ims:ime, jms:jme ) , &
135 INTENT(INOUT) :: REGIME, &
136 HFX, &
137 QFX, &
138 LH, &
139 MOL,RMOL
140 !m the following 5 are change to memory size
141 !
142 REAL, DIMENSION( ims:ime, jms:jme ) , &
143 INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
144 PSIM,PSIH
145
146 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
147 INTENT(IN ) :: U3D, &
148 V3D
149
150 REAL, DIMENSION( ims:ime, jms:jme ) , &
151 INTENT(IN ) :: PSFC
152
153 REAL, DIMENSION( ims:ime, jms:jme ) , &
154 INTENT(INOUT) :: ZNT, &
155 ZOL, &
156 UST, &
157 CPM, &
158 CHS2, &
159 CQS2, &
160 CHS
161
162 REAL, DIMENSION( ims:ime, jms:jme ) , &
163 INTENT(INOUT) :: FLHC,FLQC
164
165 REAL, DIMENSION( ims:ime, jms:jme ) , &
166 INTENT(INOUT) :: &
167 QGH
168
169
170
171 REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
172
173 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , &
174 INTENT(OUT) :: uratx,vratx,tratx
175 ! LOCAL VARS
176
177 REAL, DIMENSION( its:ite ) :: U1D, &
178 V1D, &
179 QV1D, &
180 P1D, &
181 T1D
182
183 REAL, DIMENSION( its:ite ) :: dz8w1d
184
185 INTEGER :: I,J
186
187 DO J=jts,jte
188 DO i=its,ite
189 dz8w1d(I) = dz8w(i,1,j)
190 ENDDO
191
192 DO i=its,ite
193 U1D(i) =U3D(i,1,j)
194 V1D(i) =V3D(i,1,j)
195 QV1D(i)=QV3D(i,1,j)
196 P1D(i) =P3D(i,1,j)
197 T1D(i) =T3D(i,1,j)
198 ENDDO
199
200 CALL SFCLAY1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, &
201 CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
202 CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
203 ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
204 MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
205 XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
206 U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), &
207 Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j), &
208 QSFC(ims,j),LH(ims,j), &
209 GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, &
210 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
211 P1000mb, &
212 ids,ide, jds,jde, kds,kde, &
213 ims,ime, jms,jme, kms,kme, &
214 its,ite, jts,jte, kts,kte, &
215 uratx(ims,j),vratx(ims,j),tratx(ims,j) )
216 ENDDO
217
218
219 END SUBROUTINE SFCLAY
220
221
222 !-------------------------------------------------------------------
223 SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d, &
224 CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
225 ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
226 XLAND,HFX,QFX,TSK, &
227 U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH, &
228 QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
229 SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
230 KARMAN,EOMEG,STBOLT, &
231 P1000mb, &
232 ids,ide, jds,jde, kds,kde, &
233 ims,ime, jms,jme, kms,kme, &
234 its,ite, jts,jte, kts,kte, &
235 uratx,vratx,tratx )
236 !-------------------------------------------------------------------
237 IMPLICIT NONE
238 !-------------------------------------------------------------------
239 REAL, PARAMETER :: XKA=2.4E-5
240 REAL, PARAMETER :: PRT=1.
241
242 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
243 ims,ime, jms,jme, kms,kme, &
244 its,ite, jts,jte, kts,kte, &
245 J
246 !
247 INTEGER, INTENT(IN ) :: ISFFLX
248 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
249 REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
250 REAL, INTENT(IN ) :: P1000mb
251
252 !
253 REAL, DIMENSION( ims:ime ) , &
254 INTENT(IN ) :: MAVAIL, &
255 PBLH, &
256 XLAND, &
257 TSK
258 !
259 REAL, DIMENSION( ims:ime ) , &
260 INTENT(IN ) :: PSFCPA
261
262 REAL, DIMENSION( ims:ime ) , &
263 INTENT(INOUT) :: REGIME, &
264 HFX, &
265 QFX, &
266 MOL,RMOL
267 !m the following 5 are changed to memory size---
268 !
269 REAL, DIMENSION( ims:ime ) , &
270 INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
271 PSIM,PSIH
272
273 REAL, DIMENSION( ims:ime ) , &
274 INTENT(INOUT) :: ZNT, &
275 ZOL, &
276 UST, &
277 CPM, &
278 CHS2, &
279 CQS2, &
280 CHS
281
282 REAL, DIMENSION( ims:ime ) , &
283 INTENT(INOUT) :: FLHC,FLQC
284
285 REAL, DIMENSION( ims:ime ) , &
286 INTENT(INOUT) :: &
287 QGH
288
289 REAL, DIMENSION( ims:ime ) , &
290 INTENT(OUT) :: U10,V10, &
291 TH2,T2,Q2,QSFC,LH
292
293
294 REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
295
296 ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
297 REAL, DIMENSION( its:ite ), INTENT(IN ) :: dz8w1d
298
299 REAL, DIMENSION( its:ite ), INTENT(IN ) :: UX, &
300 VX, &
301 QV1D, &
302 P1D, &
303 T1D
304
305 REAL, OPTIONAL, DIMENSION( ims:ime ) , &
306 INTENT(OUT) :: uratx,vratx,tratx
307
308 ! LOCAL VARS
309
310 REAL, DIMENSION( its:ite ) :: ZA, &
311 THVX,ZQKL, &
312 ZQKLP1, &
313 THX,QX, &
314 PSIH2, &
315 PSIM2, &
316 PSIH10, &
317 PSIM10, &
318 GZ2OZ0, &
319 GZ10OZ0
320 !
321 REAL, DIMENSION( its:ite ) :: &
322 RHOX,GOVRTH, &
323 TGDSA
324 !
325 REAL, DIMENSION( its:ite) :: SCR3,SCR4
326 REAL, DIMENSION( its:ite ) :: THGB, PSFC
327 !
328 INTEGER :: KL
329
330 INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
331
332 REAL :: PL,THCON,TVCON,E1
333 REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
334 REAL :: DTG,PSIX,USTM,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2
335 REAL :: FLUXC,VSGD
336 !-------------------------------------------------------------------
337 KL=kte
338
339 DO i=its,ite
340 ! PSFC cmb
341 PSFC(I)=PSFCPA(I)/1000.
342 ENDDO
343 !
344 !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
345 !
346 DO 5 I=its,ite
347 TGDSA(I)=TSK(I)
348 ! PSFC cmb
349 ! THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP
350 THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
351 5 CONTINUE
352 !
353 !-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
354 ! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
355 !
356 ! *** NOTE ***
357 ! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
358 ! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
359 ! TENDENCIES.
360 !
361 10 CONTINUE
362
363 ! DO 24 I=its,ite
364 ! UX(I)=U1D(I)
365 ! VX(I)=V1D(I)
366 ! 24 CONTINUE
367
368 26 CONTINUE
369
370 !.....SCR3(I,K) STORE TEMPERATURE,
371 ! SCR4(I,K) STORE VIRTUAL TEMPERATURE.
372
373 DO 30 I=its,ite
374 ! PL cmb
375 PL=P1D(I)/1000.
376 SCR3(I)=T1D(I)
377 ! THCON=(100./PL)**ROVCP
378 THCON=(P1000mb*0.001/PL)**ROVCP
379 THX(I)=SCR3(I)*THCON
380 SCR4(I)=SCR3(I)
381 THVX(I)=THX(I)
382 QX(I)=0.
383 30 CONTINUE
384 !
385 DO I=its,ite
386 QGH(I)=0.
387 FLHC(I)=0.
388 FLQC(I)=0.
389 CPM(I)=CP
390 ENDDO
391 !
392 ! IF(IDRY.EQ.1)GOTO 80
393 DO 50 I=its,ite
394 QX(I)=QV1D(I)
395 TVCON=(1.+EP1*QX(I))
396 THVX(I)=THX(I)*TVCON
397 SCR4(I)=SCR3(I)*TVCON
398 50 CONTINUE
399 !
400 DO 60 I=its,ite
401 E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
402 QSFC(I)=EP2*E1/(PSFC(I)-E1)
403 ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
404 ! Q2SAT = QGH IN LSM
405 E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
406 QGH(I)=EP2*E1/(PSFC(I)-E1)
407 CPM(I)=CP*(1.+0.8*QX(I))
408 60 CONTINUE
409 80 CONTINUE
410
411 !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
412 ! LEVEL, AND THE LAYER THICKNESSES.
413
414 DO 90 I=its,ite
415 ZQKLP1(I)=0.
416 RHOX(I)=PSFC(I)*1000./(R*SCR4(I))
417 90 CONTINUE
418 !
419 DO 110 I=its,ite
420 ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
421 110 CONTINUE
422 !
423 DO 120 I=its,ite
424 ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
425 120 CONTINUE
426 !
427 DO 160 I=its,ite
428 GOVRTH(I)=G/THX(I)
429 160 CONTINUE
430
431 !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
432 ! AKB(1976), EQ(12).
433
434 DO 260 I=its,ite
435 GZ1OZ0(I)=ALOG(ZA(I)/ZNT(I))
436 GZ2OZ0(I)=ALOG(2./ZNT(I))
437 GZ10OZ0(I)=ALOG(10./ZNT(I))
438 IF((XLAND(I)-1.5).GE.0)THEN
439 ZL=ZNT(I)
440 ELSE
441 ZL=0.01
442 ENDIF
443 WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
444
445 TSKV=THGB(I)*(1.+EP1*QSFC(I)*MAVAIL(I))
446 DTHVDZ=(THVX(I)-TSKV)
447 ! Convective velocity scale Vc and subgrid-scale velocity Vsg
448 ! following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
449 ! ... HONG Aug. 2001
450 !
451 ! VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
452 ! Use Beljaars over land, old MM5 (Wyngaard) formula over water
453 if (xland(i).lt.1.5) then
454 fluxc = max(hfx(i)/rhox(i)/cp &
455 + ep1*tskv*qfx(i)/rhox(i),0.)
456 VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
457 else
458 IF(-DTHVDZ.GE.0)THEN
459 DTHVM=-DTHVDZ
460 ELSE
461 DTHVM=0.
462 ENDIF
463 VCONV = 2.*SQRT(DTHVM)
464 endif
465 ! Mahrt and Sun low-res correction
466 VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
467 WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
468 WSPD(I)=AMAX1(WSPD(I),0.1)
469 BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
470 ! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
471 IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
472 !jdf
473 RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
474 !jdf
475
476 260 CONTINUE
477
478 !
479 !-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
480 !
481 !
482 ! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
483 ! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
484 !
485 ! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
486 !
487 ! 1. BR .GE. 0.2;
488 ! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
489 !
490 ! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
491 ! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
492 ! (REGIME=2),
493 !
494 ! 3. BR .EQ. 0.0
495 ! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
496 !
497 ! 4. BR .LT. 0.0
498 ! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
499 !
500 !CCCCC
501
502 DO 320 I=its,ite
503 !CCCCC
504 !CC REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
505 !CC IF(BR(I).LT.0..AND.HOL(I,J).GT.1.5)GOTO 310
506 IF(BR(I).LT.0.)GOTO 310
507 !
508 !-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
509 !
510 IF(BR(I).LT.0.2)GOTO 270
511 REGIME(I)=1.
512 PSIM(I)=-10.*GZ1OZ0(I)
513 ! LOWER LIMIT ON PSI IN STABLE CONDITIONS
514 PSIM(I)=AMAX1(PSIM(I),-10.)
515 PSIH(I)=PSIM(I)
516 PSIM10(I)=10./ZA(I)*PSIM(I)
517 PSIM10(I)=AMAX1(PSIM10(I),-10.)
518 PSIH10(I)=PSIM10(I)
519 PSIM2(I)=2./ZA(I)*PSIM(I)
520 PSIM2(I)=AMAX1(PSIM2(I),-10.)
521 PSIH2(I)=PSIM2(I)
522
523 ! 1.0 over Monin-Obukhov length
524 IF(UST(I).LT.0.01)THEN
525 RMOL(I)=BR(I)*GZ1OZ0(I) !ZA/L
526 ELSE
527 RMOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I)) !ZA/L
528 ENDIF
529 RMOL(I)=AMIN1(RMOL(I),9.999) ! ZA/L
530 RMOL(I) = RMOL(I)/ZA(I) !1.0/L
531
532 GOTO 320
533 !
534 !-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
535 !
536 270 IF(BR(I).EQ.0.0)GOTO 280
537 REGIME(I)=2.
538 PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))
539 ! LOWER LIMIT ON PSI IN STABLE CONDITIONS
540 PSIM(I)=AMAX1(PSIM(I),-10.)
541 !.....AKB(1976), EQ(16).
542 PSIH(I)=PSIM(I)
543 PSIM10(I)=10./ZA(I)*PSIM(I)
544 PSIM10(I)=AMAX1(PSIM10(I),-10.)
545 PSIH10(I)=PSIM10(I)
546 PSIM2(I)=2./ZA(I)*PSIM(I)
547 PSIM2(I)=AMAX1(PSIM2(I),-10.)
548 PSIH2(I)=PSIM2(I)
549
550 ! Linear form: PSIM = -0.5*ZA/L; e.g, see eqn 16 of
551 ! Blackadar, Modeling the nocturnal boundary layer, Preprints,
552 ! Third Symposium on Atmospheric Turbulence Diffusion and Air Quality,
553 ! Raleigh, NC, 1976
554 ZOL(I) = BR(I)*GZ1OZ0(I)/(1.00001-5.0*BR(I))
555
556 if ( ZOL(I) .GT. 0.5 ) then ! linear form ok
557 ! Holtslag and de Bruin, J. App. Meteor 27, 689-704, 1988;
558 ! see also, Launiainen, Boundary-Layer Meteor 76,165-179, 1995
559 ! Eqn (8) of Launiainen, 1995
560 ZOL(I) = ( 1.89*GZ1OZ0(I) + 44.2 ) * BR(I)*BR(I) &
561 + ( 1.18*GZ1OZ0(I) - 1.37 ) * BR(I)
562 ZOL(I)=AMIN1(ZOL(I),9.999)
563 end if
564
565 ! 1.0 over Monin-Obukhov length
566 RMOL(I)= ZOL(I)/ZA(I)
567
568 GOTO 320
569 !
570 !-----CLASS 3; FORCED CONVECTION:
571 !
572 280 REGIME(I)=3.
573 PSIM(I)=0.0
574 PSIH(I)=PSIM(I)
575 PSIM10(I)=0.
576 PSIH10(I)=PSIM10(I)
577 PSIM2(I)=0.
578 PSIH2(I)=PSIM2(I)
579
580
581 IF(UST(I).LT.0.01)THEN
582 ZOL(I)=BR(I)*GZ1OZ0(I)
583 ELSE
584 ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
585 ENDIF
586
587 RMOL(I) = ZOL(I)/ZA(I)
588
589 GOTO 320
590 !
591 !-----CLASS 4; FREE CONVECTION:
592 !
593 310 CONTINUE
594 REGIME(I)=4.
595 IF(UST(I).LT.0.01)THEN
596 ZOL(I)=BR(I)*GZ1OZ0(I)
597 ELSE
598 ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
599 ENDIF
600 ZOL10=10./ZA(I)*ZOL(I)
601 ZOL2=2./ZA(I)*ZOL(I)
602 ZOL(I)=AMIN1(ZOL(I),0.)
603 ZOL(I)=AMAX1(ZOL(I),-9.9999)
604 ZOL10=AMIN1(ZOL10,0.)
605 ZOL10=AMAX1(ZOL10,-9.9999)
606 ZOL2=AMIN1(ZOL2,0.)
607 ZOL2=AMAX1(ZOL2,-9.9999)
608 NZOL=INT(-ZOL(I)*100.)
609 RZOL=-ZOL(I)*100.-NZOL
610 NZOL10=INT(-ZOL10*100.)
611 RZOL10=-ZOL10*100.-NZOL10
612 NZOL2=INT(-ZOL2*100.)
613 RZOL2=-ZOL2*100.-NZOL2
614 PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
615 PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
616 PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))
617 PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
618 PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))
619 PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))
620
621 !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
622 !--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
623 ! PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
624 ! PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
625 PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
626 PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
627 PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
628 PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
629
630 RMOL(I) = ZOL(I)/ZA(I)
631
632 320 CONTINUE
633 !
634 !-----COMPUTE THE FRICTIONAL VELOCITY:
635 ! ZA(1982) EQS(2.60),(2.61).
636 !
637 DO 330 I=its,ite
638 DTG=THX(I)-THGB(I)
639 PSIX=GZ1OZ0(I)-PSIM(I)
640 PSIX10=GZ10OZ0(I)-PSIM10(I)
641 ! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
642 ! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
643 PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
644
645 IF((XLAND(I)-1.5).GE.0)THEN
646 ZL=ZNT(I)
647 ELSE
648 ZL=0.01
649 ENDIF
650 PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)
651 PSIT2=GZ2OZ0(I)-PSIH2(I)
652 PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-PSIH2(I)
653 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
654 UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
655 U10(I)=UX(I)*PSIX10/PSIX
656 V10(I)=VX(I)*PSIX10/PSIX
657 TH2(I)=THGB(I)+DTG*PSIT2/PSIT
658 Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
659 ! T2(I) = TH2(I)*(PSFC(I)/100.)**ROVCP
660 T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP
661 ! LATER Q2 WILL BE OVERWRITTEN FOR LAND POINTS IN SURFCE
662 ! QA2(I,J) = Q2(I)
663 ! UA10(I,J) = U10(I)
664 ! VA10(I,J) = V10(I)
665 ! write(*,1002)UST(I),KARMAN*WSPD(I),PSIX,KARMAN*WSPD(I)/PSIX
666 !
667 IF(PRESENT(uratx) .and. PRESENT(vratx) .and. PRESENT(tratx))THEN
668 IF(ABS(U10(I)) .GT. 1.E-10) THEN
669 uratx(I) = UX(I)/U10(I)
670 ELSE
671 uratx(I) = 1.2
672 END IF
673 IF(ABS(V10(I)) .GT. 1.E-10) THEN
674 vratx(I) = VX(I)/V10(I)
675 ELSE
676 vratx(I) = 1.2
677 END IF
678 tratx(I) = THX(I)/TH2(I)
679 ENDIF
680
681 USTM=AMAX1(UST(I),0.1)
682 IF((XLAND(I)-1.5).GE.0)THEN
683 UST(I)=UST(I)
684 ELSE
685 UST(I)=USTM
686 ENDIF
687 ! write(*,1002)UST(I),USTM,I,J
688 1002 format(f15.12,2x,f15.12,2x,f15.12,2x,f15.12,2x,f15.12)
689 MOL(I)=KARMAN*DTG/PSIT/PRT
690 330 CONTINUE
691 !
692 335 CONTINUE
693
694 !-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:
695
696 DO i=its,ite
697 QFX(i)=0.
698 HFX(i)=0.
699 ENDDO
700
701 IF (ISFFLX.EQ.0) GOTO 410
702
703 !-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
704
705 DO 360 I=its,ite
706 IF((XLAND(I)-1.5).GE.0)THEN
707 ZNT(I)=CZO*UST(I)*UST(I)/G+OZO
708 ENDIF
709 IF((XLAND(I)-1.5).GE.0)THEN
710 ZL=ZNT(I)
711 ELSE
712 ZL=0.01
713 ENDIF
714 FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/( &
715 ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
716 DTTHX=ABS(THX(I)-THGB(I))
717 IF(DTTHX.GT.1.E-5)THEN
718 FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))
719 ! write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
720 1001 format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
721 ELSE
722 FLHC(I)=0.
723 ENDIF
724 360 CONTINUE
725
726 !
727 !-----COMPUTE SURFACE MOIST FLUX:
728 !
729 ! IF(IDRY.EQ.1)GOTO 390
730 !
731 DO 370 I=its,ite
732 QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
733 QFX(I)=AMAX1(QFX(I),0.)
734 LH(I)=XLV*QFX(I)
735 370 CONTINUE
736
737 !-----COMPUTE SURFACE HEAT FLUX:
738 !
739 390 CONTINUE
740 DO 400 I=its,ite
741 IF(XLAND(I)-1.5.GT.0.)THEN
742 HFX(I)=FLHC(I)*(THGB(I)-THX(I))
743 ELSEIF(XLAND(I)-1.5.LT.0.)THEN
744 HFX(I)=FLHC(I)*(THGB(I)-THX(I))
745 HFX(I)=AMAX1(HFX(I),-250.)
746 ENDIF
747 400 CONTINUE
748
749 DO I=its,ite
750 IF((XLAND(I)-1.5).GE.0)THEN
751 ZL=ZNT(I)
752 ELSE
753 ZL=0.01
754 ENDIF
755 CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
756 /XKA+ZA(I)/ZL)-PSIH(I))
757 ! GZ2OZ0(I)=ALOG(2./ZNT(I))
758 ! PSIM2(I)=-10.*GZ2OZ0(I)
759 ! PSIM2(I)=AMAX1(PSIM2(I),-10.)
760 ! PSIH2(I)=PSIM2(I)
761 CQS2(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*2.0 &
762 /XKA+2.0/ZL)-PSIH2(I))
763 CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I))
764 ENDDO
765
766 410 CONTINUE
767 !jdf
768 ! DO I=its,ite
769 ! IF(UST(I).GE.0.1) THEN
770 ! RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I))
771 ! ELSE
772 ! RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1)
773 ! ENDIF
774 ! ENDDO
775 !jdf
776
777 !
778 END SUBROUTINE SFCLAY1D
779
780 !====================================================================
781 SUBROUTINE sfclayinit( allowed_to_read )
782
783 LOGICAL , INTENT(IN) :: allowed_to_read
784 INTEGER :: N
785 REAL :: ZOLN,X,Y
786
787 DO N=0,1000
788 ZOLN=-FLOAT(N)*0.01
789 X=(1-16.*ZOLN)**0.25
790 PSIMTB(N)=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))- &
791 2.*ATAN(X)+2.*ATAN(1.)
792 Y=(1-16*ZOLN)**0.5
793 PSIHTB(N)=2*ALOG(0.5*(1+Y))
794 ENDDO
795
796 END SUBROUTINE sfclayinit
797
798 !-------------------------------------------------------------------
799
800 END MODULE module_sf_sfclay