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