| blob | fa9853a2637604dabf82a0bfe06be6f0068f8c09 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_sf_mynn
5 !-------------------------------------------------------------------
6 !Modifications implemented by Joseph Olson NOAA/GSL
7 !The following overviews the current state of this scheme::
8 !
9 ! BOTH LAND AND WATER:
10 !1) Calculation of stability parameter (z/L) taken from Li et al. (2010 BLM)
11 ! for first iteration of first time step; afterwards, exact calculation
12 ! using a brute force iterative method described in Olson et al. (2021 NOAA
13 ! Tech memorandum). This method replaces the iterative technique used in
14 ! module_sf_sfclayrev.F (Jimenez et al. 2013) with mods. Either technique
15 ! gives about the same result. The former technique is retained in this
16 ! module (commented out) for potential subsequent benchmarking.
17 !2) Fixed isflux=0 option to turn off scalar fluxes, but keep momentum
18 ! fluxes for idealized studies (credit: Anna Fitch).
19 !3) Kinematic viscosity varies with temperature according to Andreas (1989).
20 !4) Uses the blended Monin-Obukhov flux-profile relationships COARE (Fairall
21 ! et al 2003) for the unstable regime (a blended mix of Dyer-Hicks 1974 and
22 ! Grachev et al (2000). Uses Cheng and Brutsaert (2005) for stable conditions.
23 !5) The following overviews the namelist variables that control the
24 ! aerodynamic roughness lengths (over water) and the thermal and moisture
25 ! roughness lengths (defaults are recommended):
26 !
27 ! LAND only:
28 ! "iz0tlnd" namelist option is used to select the following options:
29 ! (default) =0: Zilitinkevich (1995); Czil now set to 0.085
30 ! =1: Czil_new (modified according to Chen & Zhang 2008)
31 ! =2: Modified Yang et al (2002, 2008) - generalized for all landuse
32 ! =3: constant zt = z0/7.4 (original form; Garratt 1992)
33 !
34 ! WATER only:
35 ! "isftcflx" namelist option is used to select the following options:
36 ! (default) =0: z0, zt, and zq from the COARE algorithm. Set COARE_OPT (below) to
37 ! 3.0 (Fairall et al. 2003, default)
38 ! 3.5 (Edson et al 2013) - now with bug fix (Edson et al. 2014, JPO)
39 ! =1: z0 from Davis et al (2008), zt & zq from COARE 3.0/3.5
40 ! =2: z0 from Davis et al (2008), zt & zq from Garratt (1992)
41 ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5
42 !
43 ! SNOW/ICE only:
44 ! Andreas (2002) snow/ice parameterization for thermal and
45 ! moisture roughness is used over all gridpoints with snow deeper than
46 ! 0.1 m. This algorithm calculates a z0 for snow (Andreas et al. 2005, BLM),
47 ! which is only used as part of the thermal and moisture roughness
48 ! length calculation, not to directly impact the surface winds.
49 !
50 ! Misc:
51 !1) Added a more elaborate diagnostic for u10 & V10 for high vertical resolution
52 ! model configurations but for most model configurations with depth of
53 ! the lowest half-model level near 10 m, a neutral-log diagnostic is used.
54 !
55 !2) Option to activate stochastic parameter perturbations (SPP), which
56 ! perturb z0, zt, and zq, along with many other parameters in the MYNN-
57 ! EDMF scheme.
58 !
59 !NOTE: This code was primarily tested in combination with the RUC LSM.
60 ! Performance with the Noah (or other) LSM is relatively unknown.
61 !-------------------------------------------------------------------
62 !For WRF
63 USE module_model_constants, only: &
64 &p1000mb, ep_2
66 !-------------------------------------------------------------------
67 IMPLICIT NONE
68 !-------------------------------------------------------------------
69 !For non-WRF
70 ! REAL , PARAMETER :: g = 9.81
71 ! REAL , PARAMETER :: r_d = 287.
72 ! REAL , PARAMETER :: cp = 7.*r_d/2.
73 ! REAL , PARAMETER :: r_v = 461.6
74 ! REAL , PARAMETER :: cpv = 4.*r_v
75 ! REAL , PARAMETER :: rcp = r_d/cp
76 ! REAL , PARAMETER :: XLV = 2.5E6
77 ! REAL , PARAMETER :: XLF = 3.50E5
78 ! REAL , PARAMETER :: p1000mb = 100000.
79 ! REAL , PARAMETER :: EP_2 = r_d/r_v
81 REAL, PARAMETER :: ep_3=1.-ep_2
82 REAL, PARAMETER :: wmin=0.1 ! Minimum wind speed
83 REAL, PARAMETER :: VCONVC=1.25
84 REAL, PARAMETER :: SNOWZ0=0.011
85 REAL, PARAMETER :: COARE_OPT=3.0 ! 3.0 or 3.5
86 INTEGER, PARAMETER :: psi_opt = 0 ! 0 = stable: Cheng and Brustaert
87 ! unstable: blended COARE
88 ! 1 = GFS
89 REAL, PARAMETER :: qvmin = 1e-9
90 !For debugging purposes:
91 LOGICAL, PARAMETER :: debug_code = .false.
93 REAL, DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab, &
94 psih_stab,psih_unstab
96 CONTAINS
98 !-------------------------------------------------------------------
99 SUBROUTINE mynn_sf_init_driver(allowed_to_read)
101 LOGICAL, INTENT(in) :: allowed_to_read
103 !Fill the PSIM and PSIH tables. This code was leveraged from
104 !module_sf_sfclayrev.F, leveraging the work from Pedro Jimenez.
105 !This subroutine returns a blended form from Dyer and Hicks (1974)
106 !and Grachev et al (2000) for unstable conditions and the form
107 !from Cheng and Brutsaert (2005) for stable conditions.
109 CALL psi_init(psi_opt)
111 END SUBROUTINE mynn_sf_init_driver
113 !-------------------------------------------------------------------
114 SUBROUTINE SFCLAY_mynn( &
115 U3D,V3D,T3D,QV3D,P3D,dz8w,th3d,rho3d, &
116 CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM, &
117 ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
118 XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
119 U10,V10,TH2,T2,Q2,SNOWH, &
120 GZ1OZ0,WSPD,BR,ISFFLX,DX, &
121 SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
122 KARMAN,itimestep,ch,qcg, &
123 spp_pbl,pattern_spp_pbl, &
124 ids,ide, jds,jde, kds,kde, &
125 ims,ime, jms,jme, kms,kme, &
126 its,ite, jts,jte, kts,kte, &
127 ustm,ck,cka,cd,cda,isftcflx,iz0tlnd )
128 !-------------------------------------------------------------------
129 IMPLICIT NONE
130 !-------------------------------------------------------------------
131 !-- U3D 3D u-velocity interpolated to theta points (m/s)
132 !-- V3D 3D v-velocity interpolated to theta points (m/s)
133 !-- T3D 3D temperature (K)
134 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
135 !-- P3D 3D pressure (Pa)
136 !-- RHO3D 3D density (kg/m3)
137 !-- dz8w 3D dz between full levels (m)
138 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
139 !-- G acceleration due to gravity (m/s^2)
140 !-- ROVCP R/CP
141 !-- R gas constant for dry air (J/kg/K)
142 !-- XLV latent heat of vaporization for water (J/kg)
143 !-- PSFCPA surface pressure (Pa)
144 !-- ZNT roughness length (m)
145 !-- UST u* in similarity theory (m/s)
146 !-- USTM u* in similarity theory (m/s) w* added to WSPD. This is
147 ! used to couple with TKE scheme but not in MYNN.
148 ! (as of now, USTM = UST in this version)
149 !-- PBLH PBL height from previous time (m)
150 !-- MAVAIL surface moisture availability (between 0 and 1)
151 !-- ZOL z/L height over Monin-Obukhov length
152 !-- MOL T* (similarity theory) (K)
153 !-- RMOL Reciprocal of M-O length (/m)
154 !-- REGIME flag indicating PBL regime (stable, unstable, etc.)
155 !-- PSIM similarity stability function for momentum
156 !-- PSIH similarity stability function for heat
157 !-- XLAND land mask (1 for land, 2 for water)
158 !-- HFX upward heat flux at the surface (W/m^2)
159 !-- QFX upward moisture flux at the surface (kg/m^2/s)
160 !-- LH net upward latent heat flux at surface (W/m^2)
161 !-- TSK surface temperature (K)
162 !-- FLHC exchange coefficient for heat (W/m^2/K)
163 !-- FLQC exchange coefficient for moisture (kg/m^2/s)
164 !-- CHS heat/moisture exchange coefficient for LSM (m/s)
165 !-- QGH lowest-level saturated mixing ratio
166 !-- QSFC qv (specific humidity) at the surface
167 !-- QSFCMR qv (mixing ratio) at the surface
168 !-- U10 diagnostic 10m u wind
169 !-- V10 diagnostic 10m v wind
170 !-- TH2 diagnostic 2m theta (K)
171 !-- T2 diagnostic 2m temperature (K)
172 !-- Q2 diagnostic 2m mixing ratio (kg/kg)
173 !-- SNOWH Snow height (m)
174 !-- GZ1OZ0 log((z1+ZNT)/ZNT) where ZNT is roughness length
175 !-- WSPD wind speed at lowest model level (m/s)
176 !-- BR bulk Richardson number in surface layer
177 !-- ISFFLX isfflx=1 for surface heat and moisture fluxes
178 !-- DX horizontal grid size (m)
179 !-- SVP1 constant for saturation vapor pressure (=0.6112 kPa)
180 !-- SVP2 constant for saturation vapor pressure (=17.67 dimensionless)
181 !-- SVP3 constant for saturation vapor pressure (=29.65 K)
182 !-- SVPT0 constant for saturation vapor pressure (=273.15 K)
183 !-- EP1 constant for virtual temperature (Rv/Rd - 1) (dimensionless)
184 !-- EP2 constant for spec. hum. calc (Rd/Rv = 0.622) (dimensionless)
185 !-- EP3 constant for spec. hum. calc (1 - Rd/Rv = 0.378 ) (dimensionless)
186 !-- KARMAN Von Karman constant
187 !-- ck enthalpy exchange coeff at 10 meters
188 !-- cd momentum exchange coeff at 10 meters
189 !-- cka enthalpy exchange coeff at the lowest model level
190 !-- cda momentum exchange coeff at the lowest model level
191 !-- isftcflx =0: z0, zt, and zq from COARE3.0/3.5 (Fairall et al 2003/Edson et al 2013)
192 ! (water =1: z0 from Davis et al (2008), zt & zq from COARE3.0/3.5
193 ! only) =2: z0 from Davis et al (2008), zt & zq from Garratt (1992)
194 ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5
195 !-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.085,
196 ! (land =1: Czil_new (modified according to Chen & Zhang 2008)
197 ! only) =2: Modified Yang et al (2002, 2008) - generalized for all landuse
198 ! =3: constant zt = z0/7.4 (Garratt 1992)
199 !
200 !-- ids start index for i in domain
201 !-- ide end index for i in domain
202 !-- jds start index for j in domain
203 !-- jde end index for j in domain
204 !-- kds start index for k in domain
205 !-- kde end index for k in domain
206 !-- ims start index for i in memory
207 !-- ime end index for i in memory
208 !-- jms start index for j in memory
209 !-- jme end index for j in memory
210 !-- kms start index for k in memory
211 !-- kme end index for k in memory
212 !-- its start index for i in tile
213 !-- ite end index for i in tile
214 !-- jts start index for j in tile
215 !-- jte end index for j in tile
216 !-- kts start index for k in tile
217 !-- kte end index for k in tile
218 !=================================================================
219 ! SCALARS
220 !===================================
221 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
222 ims,ime, jms,jme, kms,kme, &
223 its,ite, jts,jte, kts,kte
224 INTEGER, INTENT(IN) :: itimestep
225 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
226 REAL, INTENT(IN) :: EP1,EP2,KARMAN
227 REAL, INTENT(IN) :: CP,G,ROVCP,R,XLV,DX
228 !NAMELIST OPTIONS:
229 INTEGER, INTENT(IN) :: ISFFLX
230 INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND
231 INTEGER, OPTIONAL, INTENT(IN) :: spp_pbl
233 !===================================
234 ! 3D VARIABLES
235 !===================================
236 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
237 INTENT(IN ) :: dz8w, &
238 QV3D, &
239 P3D, &
240 T3D, &
241 U3D,V3D, &
242 RHO3D,th3d
244 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
245 INTENT(IN) :: pattern_spp_pbl
246 !===================================
247 ! 2D VARIABLES
248 !===================================
249 REAL, DIMENSION( ims:ime, jms:jme ) , &
250 INTENT(IN ) :: MAVAIL, &
251 PBLH, &
252 XLAND, &
253 TSK, &
254 QCG, &
255 PSFCPA ,&
256 SNOWH
259 REAL, DIMENSION( ims:ime, jms:jme ) , &
260 INTENT(OUT ) :: U10,V10, &
261 TH2,T2,Q2
263 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , &
264 INTENT(OUT) :: ck,cka,cd,cda,ustm
265 !
266 REAL, DIMENSION( ims:ime, jms:jme ) , &
267 INTENT(INOUT) :: REGIME, &
268 HFX, &
269 QFX, &
270 LH, &
271 MOL,RMOL, &
272 QSFC, QGH, &
273 ZNT, &
274 ZOL, &
275 UST, &
276 CPM, &
277 CHS2, &
278 CQS2, &
279 CHS, &
280 CH, &
281 FLHC,FLQC, &
282 GZ1OZ0,WSPD,BR, &
283 PSIM,PSIH
285 !ADDITIONAL OUTPUT
286 REAL, DIMENSION( ims:ime, jms:jme ) :: wstar,qstar
287 !===================================
288 ! 1D LOCAL ARRAYS
289 !===================================
290 REAL, DIMENSION( its:ite ) :: U1D, &
291 V1D, &
292 U1D2,V1D2, & !level2 winds
293 QV1D, &
294 P1D, &
295 T1D, &
296 RHO1D, &
297 dz8w1d, & !level 1 height
298 dz2w1d !level 2 height
301 REAL, DIMENSION( its:ite ) :: rstoch1D
303 INTEGER :: I,J,K,itf,jtf,ktf
304 !-----------------------------------------------------------
306 itf=MIN0(ite,ide-1)
307 jtf=MIN0(jte,jde-1)
308 ktf=MIN0(kte,kde-1)
310 DO J=jts,jte
311 DO i=its,ite
312 dz8w1d(I) = dz8w(i,kts,j)
313 dz2w1d(I) = dz8w(i,kts+1,j)
314 U1D(i) =U3D(i,kts,j)
315 V1D(i) =V3D(i,kts,j)
316 !2nd model level winds - for diags with high-res grids
317 U1D2(i) =U3D(i,kts+1,j)
318 V1D2(i) =V3D(i,kts+1,j)
319 QV1D(i)=QV3D(i,kts,j)
320 P1D(i) =P3D(i,kts,j)
321 T1D(i) =T3D(i,kts,j)
322 RHO1D(i)=RHO3D(i,kts,j)
323 if (spp_pbl==1) then
324 rstoch1D(i)=pattern_spp_pbl(i,kts,j)
325 else
326 rstoch1D(i)=0.0
327 endif
328 ENDDO
330 IF (itimestep==1) THEN
331 DO i=its,ite
332 UST(i,j)=MAX(0.04*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i)),0.001)
333 MOL(i,j)=0. ! Tstar
334 QSFC(i,j)=QV3D(i,kts,j)/(1.+QV3D(i,kts,j))
335 qstar(i,j)=0.0
336 ENDDO
337 ENDIF
339 CALL SFCLAY1D_mynn( &
340 J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,rho1d, &
341 U1D2,V1D2,dz2w1d, &
342 CP,G,ROVCP,R,XLV,PSFCPA(ims,j),CHS(ims,j),CHS2(ims,j),&
343 CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
344 ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
345 MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
346 XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
347 U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), &
348 Q2(ims,j),FLHC(ims,j),FLQC(ims,j),SNOWH(ims,j), &
349 QGH(ims,j),QSFC(ims,j),LH(ims,j), &
350 GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, &
351 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, &
352 ch(ims,j),qcg(ims,j), &
353 itimestep, &
354 !JOE-begin additional output
355 wstar(ims,j),qstar(ims,j), &
356 !JOE-end
357 spp_pbl,rstoch1D, &
358 ids,ide, jds,jde, kds,kde, &
359 ims,ime, jms,jme, kms,kme, &
360 its,ite, jts,jte, kts,kte, &
361 isftcflx, iz0tlnd, &
362 USTM(ims,j),CK(ims,j),CKA(ims,j), &
363 CD(ims,j),CDA(ims,j) &
364 )
365 ENDDO
367 END SUBROUTINE SFCLAY_MYNN
369 !-------------------------------------------------------------------
370 SUBROUTINE SFCLAY1D_mynn( &
371 J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,rho1d, &
372 U1D2,V1D2,dz2w1d, &
373 CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM, &
374 PBLH,RMOL,ZNT,UST,MAVAIL,ZOL,MOL,REGIME, &
375 PSIM,PSIH,XLAND,HFX,QFX,TSK, &
376 U10,V10,TH2,T2,Q2,FLHC,FLQC,SNOWH,QGH, &
377 QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
378 SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
379 KARMAN,ch,qcg, &
380 itimestep, &
381 !JOE-additional output
382 wstar,qstar, &
383 !JOE-end
384 spp_pbl,rstoch1D, &
385 ids,ide, jds,jde, kds,kde, &
386 ims,ime, jms,jme, kms,kme, &
387 its,ite, jts,jte, kts,kte, &
388 isftcflx, iz0tlnd, &
389 ustm,ck,cka,cd,cda &
390 )
392 !-------------------------------------------------------------------
393 IMPLICIT NONE
394 !-------------------------------------------------------------------
395 ! SCALARS
396 !-----------------------------
397 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
398 ims,ime, jms,jme, kms,kme, &
399 its,ite, jts,jte, kts,kte, &
400 J, itimestep
402 REAL, PARAMETER :: XKA=2.4E-5 !molecular diffusivity
403 REAL, PARAMETER :: PRT=1. !prandlt number
404 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0,EP1,EP2
405 REAL, INTENT(IN) :: KARMAN,CP,G,ROVCP,R,XLV,DX
407 !-----------------------------
408 ! NAMELIST OPTIONS
409 !-----------------------------
410 INTEGER, INTENT(IN) :: ISFFLX
411 INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
412 INTEGER, INTENT(IN) :: spp_pbl
414 !-----------------------------
415 ! 1D ARRAYS
416 !-----------------------------
417 REAL, DIMENSION( ims:ime ), INTENT(IN) :: MAVAIL, &
418 PBLH, &
419 XLAND, &
420 TSK, &
421 PSFCPA, &
422 QCG, &
423 SNOWH
425 REAL, DIMENSION( its:ite ), INTENT(IN) :: U1D,V1D, &
426 U1D2,V1D2, &
427 QV1D,P1D, &
428 T1D, &
429 dz8w1d,dz2w1d, &
430 RHO1D
432 REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: REGIME, &
433 HFX,QFX,LH, &
434 MOL,RMOL, &
435 QGH,QSFC, &
436 ZNT, &
437 ZOL, &
438 UST, &
439 CPM, &
440 CHS2,CQS2, &
441 CHS,CH, &
442 FLHC,FLQC, &
443 GZ1OZ0, &
444 WSPD, &
445 BR, &
446 PSIM,PSIH
447 REAL, DIMENSION( its:ite ), INTENT(IN) :: rstoch1D
449 ! DIAGNOSTIC OUTPUT
450 REAL, DIMENSION( ims:ime ), INTENT(OUT) :: U10,V10, &
451 TH2,T2,Q2
453 REAL, OPTIONAL, DIMENSION( ims:ime ) , &
454 INTENT(OUT) :: ck,cka,cd,cda,ustm
455 !--------------------------------------------
456 !JOE-additinal output
457 REAL, DIMENSION( ims:ime ) :: wstar,qstar
458 !JOE-end
459 !----------------------------------------------------------------
460 ! LOCAL VARS
461 !----------------------------------------------------------------
462 REAL, DIMENSION(its:ite) :: &
463 ZA, & !Height of lowest 1/2 sigma level(m)
464 ZA2, & !Height of 2nd lowest 1/2 sigma level(m)
465 THV1D, & !Theta-v at lowest 1/2 sigma (K)
466 TH1D, & !Theta at lowest 1/2 sigma (K)
467 TC1D, & !T at lowest 1/2 sigma (Celsius)
468 TV1D, & !Tv at lowest 1/2 sigma (K)
469 QVSH, & !qv at lowest 1/2 sigma (spec humidity)
470 PSIH2, & !M-O stability functions at z=2 m
471 PSIM10, & !M-O stability functions at z=10 m
472 PSIH10, & !M-O stability functions at z=10 m
473 WSPDI, &
474 z_t,z_q, & !thermal & moisture roughness lengths
475 ZNTstoch, &
476 GOVRTH, & !g/theta
477 THGB, & !theta at ground
478 THVGB, & !theta-v at ground
479 PSFC, & !press at surface (Pa/1000)
480 QSFCMR, & !qv at surface (mixing ratio, kg/kg)
481 GZ2OZ0, & !LOG((2.0+ZNT(I))/ZNT(I))
482 GZ10OZ0, & !LOG((10.+ZNT(I))/ZNT(I))
483 GZ2OZt, & !LOG((2.0+z_t(i))/z_t(i))
484 GZ10OZt, & !LOG((10.+z_t(i))/z_t(i))
485 GZ1OZt, & !LOG((ZA(I)+z_t(i))/z_t(i))
486 zratio !z0/zt
488 INTEGER :: N,I,K,L,yesno
490 REAL :: PL,THCON,TVCON,E1
491 REAL :: DTHVDZ,DTHVM,VCONV,ZOL2,ZOL10,ZOLZA,ZOLZ0,ZOLZT
492 REAL :: DTG,PSIX,DTTHX,DTHDZ,PSIX10,PSIT,PSIT2, &
493 PSIQ,PSIQ2,PSIQ10,DZDT
494 REAL :: FLUXC,VSGD
495 REAL :: restar,VISC,DQG,OLDUST,OLDTST
497 !-------------------------------------------------------------------
499 DO I=its,ite
500 ! CONVERT GROUND & LOWEST LAYER TEMPERATURE TO POTENTIAL TEMPERATURE:
501 ! PSFC cmb
502 PSFC(I)=PSFCPA(I)/1000.
503 THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP !(K)
504 ! PL cmb
505 PL=P1D(I)/1000.
506 THCON=(100./PL)**ROVCP
507 TH1D(I)=T1D(I)*THCON !(Theta, K)
508 TC1D(I)=T1D(I)-273.15 !(T, Celsius)
510 ! CONVERT TO VIRTUAL TEMPERATURE
511 QVSH(I)=QV1D(I)/(1.+QV1D(I)) !CONVERT TO SPEC HUM (kg/kg)
512 TVCON=(1.+EP1*QVSH(I))
513 THV1D(I)=TH1D(I)*TVCON !(K)
514 TV1D(I)=T1D(I)*TVCON !(K)
516 !RHO1D(I)=PSFCPA(I)/(R*TV1D(I)) !now using value calculated in sfc driver
517 ZA(I)=0.5*dz8w1d(I) !height of first half-sigma level
518 ZA2(I)=dz8w1d(I) + 0.5*dz2w1d(I) !height of 2nd half-sigma level
519 GOVRTH(I)=G/TH1D(I)
520 ENDDO
522 DO I=its,ite
523 IF (TSK(I) .LT. 273.15) THEN
524 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
525 E1=SVP1*EXP(4648*(1./273.15 - 1./TSK(I)) - &
526 & 11.64*LOG(273.15/TSK(I)) + 0.02265*(273.15 - TSK(I)))
527 ELSE
528 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
529 E1=SVP1*EXP(SVP2*(TSK(I)-SVPT0)/(TSK(I)-SVP3))
530 ENDIF
531 !FOR LAND POINTS, QSFC can come from LSM, ONLY RECOMPUTE OVER WATER
532 IF (xland(i).gt.1.5 .or. QSFC(i).le.0.0) THEN !WATER
533 QSFC(I)=EP2*E1/(PSFC(I)-ep_3*E1) !specific humidity
534 QSFCMR(I)=EP2*E1/(PSFC(I)-E1) !mixing ratio
535 ELSE !LAND
536 QSFCMR(I)=QSFC(I)/(1.-QSFC(I))
537 ENDIF
539 ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
540 ! Q2SAT = QGH IN LSM
541 IF (TSK(I) .LT. 273.15) THEN
542 !SATURATION VAPOR PRESSURE WRT ICE
543 E1=SVP1*EXP(4648*(1./273.15 - 1./T1D(I)) - &
544 & 11.64*LOG(273.15/T1D(I)) + 0.02265*(273.15 - T1D(I)))
545 ELSE
546 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
547 E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
548 ENDIF
549 PL=P1D(I)/1000.
550 !QGH(I)=EP2*E1/(PL-ep_3*E1) !specific humidity
551 QGH(I)=EP2*E1/(PL-E1) !mixing ratio
552 CPM(I)=CP*(1.+0.84*QV1D(I))
553 ENDDO
555 DO I=its,ite
556 WSPD(I)=SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I))
558 !TGS:THVGB(I)=THGB(I)*(1.+EP1*QSFC(I)*MAVAIL(I))
559 THVGB(I)=THGB(I)*(1.+EP1*QSFC(I))
561 DTHDZ=(TH1D(I)-THGB(I))
562 DTHVDZ=(THV1D(I)-THVGB(I))
564 !--------------------------------------------------------
565 ! Calculate the convective velocity scale (WSTAR) and
566 ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS)
567 ! and Mahrt and Sun (1995, MWR), respectively
568 !-------------------------------------------------------
569 ! Use Beljaars over land and water
570 fluxc = max(hfx(i)/RHO1D(i)/cp &
571 & + ep1*THVGB(I)*qfx(i)/RHO1D(i),0.)
572 !WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**.33
573 IF (xland(i).gt.1.5 .or. QSFC(i).le.0.0) THEN !WATER
574 WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**.33
575 ELSE !LAND
576 !increase height scale, assuming that the non-local transoport
577 !from the mass-flux (plume) mixing exceedsd the PBLH.
578 WSTAR(I) = vconvc*(g/TSK(i)*MIN(1.5*pblh(i),4000.)*fluxc)**.33
579 ENDIF
580 !--------------------------------------------------------
581 ! Mahrt and Sun low-res correction
582 ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s)
583 !--------------------------------------------------------
584 VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
585 WSPD(I)=SQRT(WSPD(I)*WSPD(I)+WSTAR(I)*WSTAR(I)+vsgd*vsgd)
586 WSPD(I)=MAX(WSPD(I),wmin)
588 !--------------------------------------------------------
589 ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER,
590 ! ACCORDING TO AKB(1976), EQ(12).
591 !--------------------------------------------------------
592 BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
593 IF (ITIMESTEP == 1) THEN
594 !SET LIMITS ACCORDING TO Li et al. (2010) Boundary-Layer Meteorol (p.158)
595 BR(I)=MAX(BR(I),-2.0)
596 BR(I)=MIN(BR(I),2.0)
597 ELSE
598 BR(I)=MAX(BR(I),-4.0)
599 BR(I)=MIN(BR(I), 4.0)
600 ENDIF
602 ! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2 (STABLE)
603 !if (itimestep .GT. 1) THEN
604 ! IF(MOL(I).LT.0.)BR(I)=MIN(BR(I),0.0)
605 !ENDIF
607 ENDDO
609 1006 format(A,F7.3,A,f9.4,A,f9.5,A,f9.4)
610 1007 format(A,F2.0,A,f6.2,A,f7.3,A,f7.2)
612 !--------------------------------------------------------------------
613 !--------------------------------------------------------------------
614 !--- BEGIN I-LOOP
615 !--------------------------------------------------------------------
616 !--------------------------------------------------------------------
618 DO I=its,ite
620 !COMPUTE KINEMATIC VISCOSITY (m2/s) Andreas (1989) CRREL Rep. 89-11
621 !valid between -173 and 277 degrees C.
622 VISC=1.326e-5*(1. + 6.542e-3*TC1D(I) + 8.301e-6*TC1D(I)*TC1D(I) &
623 - 4.84e-9*TC1D(I)*TC1D(I)*TC1D(I))
625 IF ((XLAND(I)-1.5).GE.0) THEN
626 !--------------------------------------
627 ! WATER
628 !--------------------------------------
629 ! CALCULATE z0 (znt)
630 !--------------------------------------
631 IF ( PRESENT(ISFTCFLX) ) THEN
632 IF ( ISFTCFLX .EQ. 0 ) THEN
633 IF (COARE_OPT .EQ. 3.0) THEN
634 !COARE 3.0 (MISLEADING SUBROUTINE NAME)
635 CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc,ZA(I))
636 ELSE
637 !COARE 3.5
638 CALL edson_etal_2013(ZNT(i),UST(i),WSPD(i),visc,ZA(I))
639 ENDIF
640 ELSEIF ( ISFTCFLX .EQ. 1 .OR. ISFTCFLX .EQ. 2 ) THEN
641 CALL davis_etal_2008(ZNT(i),UST(i))
642 ELSEIF ( ISFTCFLX .EQ. 3 ) THEN
643 CALL Taylor_Yelland_2001(ZNT(i),UST(i),WSPD(i))
644 ELSEIF ( ISFTCFLX .EQ. 4 ) THEN
645 IF (COARE_OPT .EQ. 3.0) THEN
646 !COARE 3.0 (MISLEADING SUBROUTINE NAME)
647 CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc,ZA(I))
648 ELSE
649 !COARE 3.5
650 CALL edson_etal_2013(ZNT(i),UST(i),WSPD(i),visc,ZA(I))
651 ENDIF
652 ENDIF
653 ELSE
654 !DEFAULT TO COARE 3.0/3.5
655 IF (COARE_OPT .EQ. 3.0) THEN
656 !COARE 3.0
657 CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc,ZA(I))
658 ELSE
659 !COARE 3.5
660 CALL edson_etal_2013(ZNT(i),UST(i),WSPD(i),visc,ZA(I))
661 ENDIF
662 ENDIF
664 ! add stochastic perturbaction of ZNT
665 if (spp_pbl==1) then
666 ZNTstoch(I) = MAX(ZNT(I) + ZNT(I)*1.0*rstoch1D(i), 1e-6)
667 else
668 ZNTstoch(I) = ZNT(I)
669 endif
671 !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING NEW ZNT
672 ! AHW: Garrattt formula: Calculate roughness Reynolds number
673 ! Kinematic viscosity of air (linear approx to
674 ! temp dependence at sea level)
675 restar=MAX(ust(i)*ZNTstoch(i)/visc, 0.1)
677 !--------------------------------------
678 !CALCULATE z_t and z_q
679 !--------------------------------------
680 IF ( PRESENT(ISFTCFLX) ) THEN
681 IF ( ISFTCFLX .EQ. 0 ) THEN
682 IF (COARE_OPT .EQ. 3.0) THEN
683 CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
684 ELSE
685 !presumably, this will be published soon, but hasn't yet
686 CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
687 ENDIF
688 ELSEIF ( ISFTCFLX .EQ. 1 ) THEN
689 IF (COARE_OPT .EQ. 3.0) THEN
690 CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
691 ELSE
692 CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
693 ENDIF
694 ELSEIF ( ISFTCFLX .EQ. 2 ) THEN
695 CALL garratt_1992(z_t(i),z_q(i),ZNTstoch(i),restar,XLAND(I))
696 ELSEIF ( ISFTCFLX .EQ. 3 ) THEN
697 IF (COARE_OPT .EQ. 3.0) THEN
698 CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
699 ELSE
700 CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
701 ENDIF
702 ENDIF
703 ELSE
704 !DEFAULT TO COARE 3.0/3.5
705 IF (COARE_OPT .EQ. 3.0) THEN
706 CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
707 ELSE
708 CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl)
709 ENDIF
710 ENDIF
712 ELSE
714 ! add stochastic perturbaction of ZNT
715 if (spp_pbl==1) then
716 ZNTstoch(I) = MAX(ZNT(I) + ZNT(I)*1.0*rstoch1D(i), 1e-6)
717 else
718 ZNTstoch(I) = ZNT(I)
719 endif
721 !--------------------------------------
722 ! LAND
723 !--------------------------------------
724 !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT
725 restar=MAX(ust(i)*ZNTstoch(i)/visc, 0.1)
727 !--------------------------------------
728 !GET z_t and z_q
729 !--------------------------------------
730 !CHECK FOR SNOW/ICE POINTS OVER LAND
731 IF ( SNOWH(i) .GE. 0.1) THEN
732 CALL Andreas_2002(ZNTSTOCH(i),visc,ust(i),z_t(i),z_q(i))
733 ELSE
734 IF ( PRESENT(IZ0TLND) ) THEN
735 IF ( IZ0TLND .LE. 1 ) THEN
736 CALL zilitinkevich_1995(ZNTSTOCH(i),z_t(i),z_q(i),restar,&
737 UST(I),KARMAN,XLAND(I),IZ0TLND,spp_pbl,rstoch1D(i))
738 ELSEIF ( IZ0TLND .EQ. 2 ) THEN
739 CALL Yang_2008(ZNTSTOCH(i),z_t(i),z_q(i),UST(i),MOL(I),&
740 qstar(I),restar,visc,XLAND(I))
741 ELSEIF ( IZ0TLND .EQ. 3 ) THEN
742 !Original MYNN in WRF-ARW used this form:
743 CALL garratt_1992(z_t(i),z_q(i),ZNTSTOCH(i),restar,XLAND(I))
744 ENDIF
745 ELSE
746 !DEFAULT TO ZILITINKEVICH
747 CALL zilitinkevich_1995(ZNTSTOCH(i),z_t(i),z_q(i),restar,&
748 UST(I),KARMAN,XLAND(I),0,spp_pbl,rstoch1D(i))
749 ENDIF
750 ENDIF
752 ENDIF
753 zratio(i)=ZNTstoch(I)/z_t(I) !needed for Li et al.
755 GZ1OZ0(I)= LOG((ZA(I)+ZNTstoch(I))/ZNTstoch(I))
756 GZ1OZt(I)= LOG((ZA(I)+ZNTstoch(I))/z_t(i))
757 GZ2OZ0(I)= LOG((2.0+ZNTstoch(I))/ZNTstoch(I))
758 GZ2OZt(I)= LOG((2.0+ZNTstoch(I))/z_t(i))
759 GZ10OZ0(I)=LOG((10.+ZNTstoch(I))/ZNTstoch(I))
760 GZ10OZt(I)=LOG((10.+ZNTstoch(I))/z_t(i))
762 !--------------------------------------------------------------------
763 !--- DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATE STABILITY CLASS:
764 !
765 ! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.).
766 !
767 ! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
768 !
769 ! 1. BR .GE. 0.2;
770 ! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
771 !
772 ! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
773 ! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
774 ! (REGIME=2),
775 !
776 ! 3. BR .EQ. 0.0
777 ! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
778 !
779 ! 4. BR .LT. 0.0
780 ! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
781 !
782 !--------------------------------------------------------------------
783 IF (BR(I) .GT. 0.0) THEN
784 IF (BR(I) .GT. 0.2) THEN
785 !---CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
786 REGIME(I)=1.
787 ELSE
788 !---CLASS 2; DAMPED MECHANICAL TURBULENCE:
789 REGIME(I)=2.
790 ENDIF
792 !COMPUTE z/L first guess:
793 IF (itimestep .LE. 1) THEN
794 CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNTstoch(I),zratio(I))
795 ELSE
796 ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST(I)*UST(I),0.0001))
797 ZOL(I)=MAX(ZOL(I),0.0)
798 ZOL(I)=MIN(ZOL(I),20.)
799 ENDIF
801 !Use Pedros iterative function to find z/L
802 !zol(I)=zolri(br(I),ZA(I),ZNTstoch(I),z_t(I),ZOL(I),psi_opt)
803 !Use brute-force method
804 zol(I)=zolrib(br(I),ZA(I),ZNTstoch(I),z_t(I),GZ1OZ0(I),GZ1OZt(I),ZOL(I),psi_opt)
805 ZOL(I)=MAX(ZOL(I),0.0)
806 ZOL(I)=MIN(ZOL(I),20.)
808 zolzt = zol(I)*z_t(I)/ZA(I) ! zt/L
809 zolz0 = zol(I)*ZNTstoch(I)/ZA(I) ! z0/L
810 zolza = zol(I)*(za(I)+ZNTstoch(I))/za(I) ! (z+z0/L
811 zol10 = zol(I)*(10.+ZNTstoch(I))/za(I) ! (10+z0)/L
812 zol2 = zol(I)*(2.+ZNTstoch(I))/za(I) ! (2+z0)/L
814 !COMPUTE PSIM and PSIH
815 IF ((XLAND(I)-1.5).GE.0) THEN
816 ! WATER
817 !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I))
818 !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
819 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
820 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),z_t(I),ZNTstoch(I),ZA(I))
821 !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0)
822 ! or use tables
823 psim(I)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt)
824 psih(I)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt)
825 psim10(I)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt)
826 psih10(I)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt)
827 psih2(I)=psih_stable(zol2,psi_opt)-psih_stable(zolz0,psi_opt)
828 ELSE
829 ! LAND
830 !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
831 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
832 !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I))
833 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),z_t(I),ZNTstoch(I),ZA(I))
834 !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0)
835 ! or use tables
836 psim(I)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt)
837 psih(I)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt)
838 psim10(I)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt)
839 psih10(I)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt)
840 psih2(I)=psih_stable(zol2,psi_opt)-psih_stable(zolz0,psi_opt)
841 ENDIF
843 !PSIM10(I)=10./ZA(I)*PSIM(I)
844 !PSIH10(I)=10./ZA(I)*PSIH(I)
845 !PSIM2(I)=2./ZA(I)*PSIM(I)
846 !PSIH2(I)=2./ZA(I)*PSIH(I)
848 ! 1.0 over Monin-Obukhov length
849 RMOL(I)= ZOL(I)/ZA(I)
851 ELSEIF(BR(I) .EQ. 0.) THEN
852 !=========================================================
853 !-----CLASS 3; FORCED CONVECTION/NEUTRAL:
854 !=========================================================
855 REGIME(I)=3.
857 PSIM(I)=0.0
858 PSIH(I)=PSIM(I)
859 PSIM10(I)=0.
860 PSIH10(I)=0.
861 PSIH2(I)=0.
863 ZOL(I)=0.
864 !IF (UST(I) .LT. 0.01) THEN
865 ! ZOL(I)=BR(I)*GZ1OZ0(I)
866 !ELSE
867 ! ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(MAX(UST(I)*UST(I),0.001))
868 !ENDIF
869 RMOL(I) = ZOL(I)/ZA(I)
871 ELSEIF(BR(I) .LT. 0.)THEN
872 !==========================================================
873 !-----CLASS 4; FREE CONVECTION:
874 !==========================================================
875 REGIME(I)=4.
877 !COMPUTE z/L first guess:
878 IF (itimestep .LE. 1) THEN
879 CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNTstoch(I),zratio(I))
880 ELSE
881 ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST(I)*UST(I),0.001))
882 ZOL(I)=MAX(ZOL(I),-20.0)
883 ZOL(I)=MIN(ZOL(I),0.0)
884 ENDIF
886 !Use Pedros iterative function to find z/L
887 !zol(I)=zolri(br(I),ZA(I),ZNTstoch(I),z_t(I),ZOL(I),psi_opt)
888 !Use alternative method
889 zol(I)=zolrib(br(I),ZA(I),ZNTstoch(I),z_t(I),GZ1OZ0(I),GZ1OZt(I),ZOL(I),psi_opt)
890 ZOL(I)=MAX(ZOL(I),-20.0)
891 ZOL(I)=MIN(ZOL(I),0.0)
893 zolzt = zol(I)*z_t(I)/ZA(I) ! zt/L
894 zolz0 = zol(I)*ZNTstoch(I)/ZA(I) ! z0/L
895 zolza = zol(I)*(za(I)+ZNTstoch(I))/za(I) ! (z+z0/L
896 zol10 = zol(I)*(10.+ZNTstoch(I))/za(I) ! (10+z0)/L
897 zol2 = zol(I)*(2.+ZNTstoch(I))/za(I) ! (2+z0)/L
899 !COMPUTE PSIM and PSIH
900 IF ((XLAND(I)-1.5).GE.0) THEN
901 ! WATER
902 !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I))
903 !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNTstoch(I), ZA(I))
904 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
905 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),z_t(I),ZNTstoch(I),ZA(I))
906 ! use tables
907 psim(I)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt)
908 psih(I)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt)
909 psim10(I)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt)
910 psih10(I)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt)
911 psih2(I)=psih_unstable(zol2,psi_opt)-psih_unstable(zolz0,psi_opt)
912 ELSE
913 ! LAND
914 !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNTstoch(I), ZA(I))
915 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
916 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),z_t(I),ZNTstoch(I),ZA(I))
917 ! use tables
918 psim(I)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt)
919 psih(I)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt)
920 psim10(I)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt)
921 psih10(I)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt)
922 psih2(I)=psih_unstable(zol2,psi_opt)-psih_unstable(zolz0,psi_opt)
923 ENDIF
925 !PSIM10(I)=10./ZA(I)*PSIM(I)
926 !PSIH2(I)=2./ZA(I)*PSIH(I)
928 !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND
929 !---HIGH ROUGHNESS. THIS PREVENTS DENOMINATOR IN FLUXES
930 !---FROM GETTING TOO SMALL
931 PSIH(I)=MIN(PSIH(I),0.9*GZ1OZt(I))
932 PSIM(I)=MIN(PSIM(I),0.9*GZ1OZ0(I))
933 PSIH2(I)=MIN(PSIH2(I),0.9*GZ2OZt(I))
934 PSIM10(I)=MIN(PSIM10(I),0.9*GZ10OZ0(I))
935 PSIH10(I)=MIN(PSIH10(I),0.9*GZ10OZt(I))
937 RMOL(I) = ZOL(I)/ZA(I)
939 ENDIF
941 !------------------------------------------------------------
942 !-----COMPUTE THE FRICTIONAL VELOCITY:
943 !------------------------------------------------------------
944 ! ZA(1982) EQS(2.60),(2.61).
945 PSIX=GZ1OZ0(I)-PSIM(I)
946 PSIX10=GZ10OZ0(I)-PSIM10(I)
947 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
948 OLDUST = UST(I)
949 UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
950 !NON-AVERAGED: UST(I)=KARMAN*WSPD(I)/PSIX
952 ! Compute u* without vconv for use in HFX calc when isftcflx > 0
953 WSPDI(I)=MAX(SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)), wmin)
954 IF ( PRESENT(USTM) ) THEN
955 USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
956 ENDIF
958 IF ((XLAND(I)-1.5).LT.0.) THEN !LAND
959 UST(I)=MAX(UST(I),0.005) !Further relaxing this limit - no need to go lower
960 !Keep ustm = ust over land.
961 IF ( PRESENT(USTM) ) USTM(I)=UST(I)
962 ENDIF
964 !------------------------------------------------------------
965 !-----COMPUTE THE THERMAL AND MOISTURE RESISTANCE (PSIQ AND PSIT):
966 !------------------------------------------------------------
967 ! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
968 ! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
969 GZ1OZt(I)= LOG((ZA(I)+ZNTstoch(I))/z_t(i))
970 GZ2OZt(I)= LOG((2.0+ZNTstoch(I))/z_t(i))
972 PSIT =MAX(GZ1OZt(I)-PSIH(I) ,1.)
973 PSIT2=MAX(GZ2OZt(I)-PSIH2(I),1.)
975 PSIQ=MAX(LOG((ZA(I)+ZNTstoch(I))/z_q(I))-PSIH(I) ,1.0)
976 PSIQ2=MAX(LOG((2.0+ZNTstoch(I))/z_q(I))-PSIH2(I) ,1.0)
977 PSIQ10=MAX(LOG((10.0+ZNTstoch(I))/z_q(I))-PSIH10(I) ,1.0)
978 !----------------------------------------------------
979 !COMPUTE THE TEMPERATURE SCALE (or FRICTION TEMPERATURE, T*)
980 !----------------------------------------------------
981 DTG=THV1D(I)-THVGB(I)
982 OLDTST=MOL(I)
983 MOL(I)=KARMAN*DTG/PSIT/PRT
984 !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I))
985 !t_star(I) = MOL(I)
986 !----------------------------------------------------
987 !COMPUTE THE MOISTURE SCALE (or q*)
988 DQG=(QVSH(i)-qsfc(i))*1000. !(kg/kg -> g/kg)
989 qstar(I)=KARMAN*DQG/PSIQ/PRT
991 !IF () THEN
992 ! write(*,1001)"REGIME:",REGIME(I)," z/L:",ZOL(I)," U*:",UST(I)," Tstar:",MOL(I)
993 ! write(*,1002)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I)," DTHV:",THV1D(I)-THVGB(I)
994 ! write(*,1003)"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",ZOL(I)/ZA(I)," DTH:",TH1D(I)-THGB(I)
995 ! write(*,1004)"Z0/Zt:",zratio(I)," Z0:",ZNTstoch(I)," Zt:",z_t(I)," za:",za(I)
996 ! write(*,1005)"Re:",restar," MAVAIL:",MAVAIL(I)," QSFC(I):",QSFC(I)," QVSH(I):",QVSH(I)
997 ! print*,"VISC=",VISC," Z0:",ZNTstoch(I)," T1D(i):",T1D(i)
998 ! write(*,*)"============================================="
999 !ENDIF
1001 ENDDO ! end i-loop
1003 1000 format(A,F6.1, A,f6.1, A,f5.1, A,f7.1)
1004 1001 format(A,F2.0, A,f10.4,A,f5.3, A,f11.5)
1005 1002 format(A,f7.2, A,f7.2, A,f7.2, A,f10.3)
1006 1003 format(A,f7.2, A,f7.2, A,f10.3,A,f10.3)
1007 1004 format(A,f11.3,A,f9.7, A,f9.7, A,f6.2, A,f10.3)
1008 1005 format(A,f9.2,A,f6.4,A,f7.4,A,f7.4)
1010 !----------------------------------------------------------
1011 ! COMPUTE SURFACE HEAT AND MOISTURE FLUXES
1012 !----------------------------------------------------------
1013 DO I=its,ite
1015 !For computing the diagnostics and fluxes (below), whether the fluxes
1016 !are turned off or on, we need the following:
1017 PSIX=GZ1OZ0(I)-PSIM(I)
1018 PSIX10=GZ10OZ0(I)-PSIM10(I)
1020 PSIT =MAX(GZ1OZt(I)-PSIH(I), 1.0)
1021 PSIT2=MAX(GZ2OZt(I)-PSIH2(I),1.0)
1023 PSIQ=MAX(LOG((ZA(I)+z_q(i))/z_q(I))-PSIH(I) ,1.0)
1024 PSIQ2=MAX(LOG((2.0+z_q(i))/z_q(I))-PSIH2(I) ,1.0)
1025 PSIQ10=MAX(LOG((10.0+z_q(i))/z_q(I))-PSIH10(I) ,1.0)
1027 IF (ISFFLX .LT. 1) THEN
1029 QFX(i) = 0.
1030 HFX(i) = 0.
1031 FLHC(I) = 0.
1032 FLQC(I) = 0.
1033 LH(I) = 0.
1034 CHS(I) = 0.
1035 CH(I) = 0.
1036 CHS2(i) = 0.
1037 CQS2(i) = 0.
1038 IF(PRESENT(ck) .and. PRESENT(cd) .and. &
1039 &PRESENT(cka) .and. PRESENT(cda)) THEN
1040 Ck(I) = 0.
1041 Cd(I) = 0.
1042 Cka(I)= 0.
1043 Cda(I)= 0.
1044 ENDIF
1045 ELSE
1047 !------------------------------------------
1048 ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC)
1049 ! AND MOISTURE (FLQC)
1050 !------------------------------------------
1051 FLQC(I)=RHO1D(I)*MAVAIL(I)*UST(I)*KARMAN/PSIQ
1052 FLHC(I)=RHO1D(I)*CPM(I)*UST(I)*KARMAN/PSIT
1054 !----------------------------------
1055 ! COMPUTE SURFACE MOISTURE FLUX:
1056 !----------------------------------
1057 QFX(I)=FLQC(I)*(QSFCMR(I)-QV1D(I))
1058 !JOE: QFX(I)=MAX(QFX(I),0.) !originally did not allow neg QFX
1059 QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX, like MYJ
1060 LH(I)=XLV*QFX(I)
1062 !----------------------------------
1063 ! COMPUTE SURFACE HEAT FLUX:
1064 !----------------------------------
1065 IF(XLAND(I)-1.5.GT.0.)THEN !WATER
1066 HFX(I)=FLHC(I)*(THGB(I)-TH1D(I))
1067 IF ( PRESENT(ISFTCFLX) ) THEN
1068 IF ( ISFTCFLX.NE.0 ) THEN
1069 ! AHW: add dissipative heating term
1070 HFX(I)=HFX(I)+RHO1D(I)*USTM(I)*USTM(I)*WSPDI(I)
1071 ENDIF
1072 ENDIF
1073 ELSEIF(XLAND(I)-1.5.LT.0.)THEN !LAND
1074 HFX(I)=FLHC(I)*(THGB(I)-TH1D(I))
1075 HFX(I)=MAX(HFX(I),-250.)
1076 ENDIF
1078 !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
1079 ! /XKA+ZA(I)/ZL)-PSIH(I))
1081 CHS(I)=UST(I)*KARMAN/PSIT
1083 ! The exchange coefficient for cloud water is assumed to be the
1084 ! same as that for heat. CH is multiplied by WSPD.
1086 !ch(i)=chs(i)
1087 ch(i)=flhc(i)/( cpm(i)*RHO1D(i) )
1089 !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY
1090 CQS2(I)=UST(I)*KARMAN/PSIQ2
1091 CHS2(I)=UST(I)*KARMAN/PSIT2
1093 IF(PRESENT(ck) .and. PRESENT(cd) .and. &
1094 &PRESENT(cka) .and. PRESENT(cda)) THEN
1095 Ck(I)=(karman/psix10)*(karman/psiq10)
1096 Cd(I)=(karman/psix10)*(karman/psix10)
1097 Cka(I)=(karman/psix)*(karman/psiq)
1098 Cda(I)=(karman/psix)*(karman/psix)
1099 ENDIF
1101 ENDIF !end ISFFLX option
1103 !-----------------------------------------------------
1104 !COMPUTE DIAGNOSTICS
1105 !-----------------------------------------------------
1106 !COMPUTE 10 M WNDS
1107 !-----------------------------------------------------
1108 ! If the lowest model level is close to 10-m, use it
1109 ! instead of the flux-based diagnostic formula.
1110 if (ZA(i) .le. 7.0) then
1111 ! high vertical resolution
1112 if(ZA2(i) .gt. 7.0 .and. ZA2(i) .lt. 13.0) then
1113 !use 2nd model level
1114 U10(I)=U1D2(I)
1115 V10(I)=V1D2(I)
1116 else
1117 U10(I)=U1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I))
1118 V10(I)=V1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I))
1119 endif
1120 elseif(ZA(i) .gt. 7.0 .and. ZA(i) .lt. 13.0) then
1121 !moderate vertical resolution
1122 !U10(I)=U1D(I)*PSIX10/PSIX
1123 !V10(I)=V1D(I)*PSIX10/PSIX
1124 !use neutral-log:
1125 U10(I)=U1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I))
1126 V10(I)=V1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I))
1127 else
1128 ! very coarse vertical resolution
1129 U10(I)=U1D(I)*PSIX10/PSIX
1130 V10(I)=V1D(I)*PSIX10/PSIX
1131 endif
1133 !-----------------------------------------------------
1134 !COMPUTE 2m T, TH, AND Q
1135 !THESE WILL BE OVERWRITTEN FOR LAND POINTS IN THE LSM
1136 !-----------------------------------------------------
1137 DTG=TH1D(I)-THGB(I)
1138 TH2(I)=THGB(I)+DTG*PSIT2/PSIT
1139 !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY
1140 !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL.
1141 IF ((TH1D(I)>THGB(I) .AND. (TH2(I)<THGB(I) .OR. TH2(I)>TH1D(I))) .OR. &
1142 (TH1D(I)<THGB(I) .AND. (TH2(I)>THGB(I) .OR. TH2(I)<TH1D(I)))) THEN
1143 TH2(I)=THGB(I) + 2.*(TH1D(I)-THGB(I))/ZA(I)
1144 ENDIF
1145 T2(I)=TH2(I)*(PSFC(I)/100.)**ROVCP
1147 Q2(I)=QSFCMR(I)+(QV1D(I)-QSFCMR(I))*PSIQ2/PSIQ
1148 Q2(I)= MAX(Q2(I), MIN(QSFCMR(I), QV1D(I)))
1149 Q2(I)= MIN(Q2(I), 1.05*QV1D(I))
1151 IF ( debug_code ) THEN
1152 yesno = 0
1153 IF (HFX(I) > 1200. .OR. HFX(I) < -700.)THEN
1154 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1155 I,J, "HFX: ",HFX(I)
1156 yesno = 1
1157 ENDIF
1158 IF (LH(I) > 1200. .OR. LH(I) < -700.)THEN
1159 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1160 I,J, "LH: ",LH(I)
1161 yesno = 1
1162 ENDIF
1163 IF (UST(I) < 0.0 .OR. UST(I) > 4.0 )THEN
1164 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1165 I,J, "UST: ",UST(I)
1166 yesno = 1
1167 ENDIF
1168 IF (WSTAR(I)<0.0 .OR. WSTAR(I) > 6.0)THEN
1169 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1170 I,J, "WSTAR: ",WSTAR(I)
1171 yesno = 1
1172 ENDIF
1173 IF (RHO1D(I)<0.0 .OR. RHO1D(I) > 1.6 )THEN
1174 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1175 I,J, "rho: ",RHO1D(I)
1176 yesno = 1
1177 ENDIF
1178 IF (QSFC(I)*1000. <0.0 .OR. QSFC(I)*1000. >40.)THEN
1179 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1180 I,J, "QSFC: ",QSFC(I)
1181 yesno = 1
1182 ENDIF
1183 IF (PBLH(I)<0. .OR. PBLH(I)>6000.)THEN
1184 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
1185 I,J, "PBLH: ",PBLH(I)
1186 yesno = 1
1187 ENDIF
1189 IF (yesno == 1) THEN
1190 print*," OTHER INFO:"
1191 write(*,1001)"REGIME:",REGIME(I)," z/L:",ZOL(I)," U*:",UST(I),&
1192 " Tstar:",MOL(I)
1193 write(*,1002)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),&
1194 " DTHV:",THV1D(I)-THVGB(I)
1195 write(*,1003)"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",&
1196 ZOL(I)/ZA(I)," DTH:",TH1D(I)-THGB(I)
1197 write(*,*)" Z0:",ZNTstoch(I)," Zt:",z_t(I)," za:",za(I)
1198 write(*,1005)"Re:",restar," MAVAIL:",MAVAIL(I)," QSFC(I):",&
1199 QSFC(I)," QVSH(I):",QVSH(I)
1200 print*,"PSIX=",PSIX," Z0:",ZNTstoch(I)," T1D(i):",T1D(i)
1201 write(*,*)"============================================="
1202 ENDIF
1203 ENDIF
1205 ENDDO !end i-loop
1207 END SUBROUTINE SFCLAY1D_mynn
1208 !-------------------------------------------------------------------
1209 SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,&
1210 & landsea,IZ0TLND2,spp_pbl,rstoch)
1212 ! This subroutine returns the thermal and moisture roughness lengths
1213 ! from Zilitinkevich (1995) and Zilitinkevich et al. (2001) over
1214 ! land and water, respectively.
1215 !
1216 ! MODS:
1217 ! 20120705 : added IZ0TLND option. Note: This option was designed
1218 ! to work with the Noah LSM and may be specific for that
1219 ! LSM only. Tests with RUC LSM showed no improvements.
1221 IMPLICIT NONE
1222 REAL, INTENT(IN) :: Z_0,restar,ustar,KARMAN,landsea
1223 INTEGER, OPTIONAL, INTENT(IN):: IZ0TLND2
1224 REAL, INTENT(OUT) :: Zt,Zq
1225 REAL :: CZIL !=0.100 in Chen et al. (1997)
1226 !=0.075 in Zilitinkevich (1995)
1227 !=0.500 in Lemone et al. (2008)
1228 INTEGER, INTENT(IN) :: spp_pbl
1229 REAL, INTENT(IN) :: rstoch
1232 IF (landsea-1.5 .GT. 0) THEN !WATER
1234 !THIS IS BASED ON Zilitinkevich, Grachev, and Fairall (2001;
1235 !Their equations 15 and 16).
1236 IF (restar .LT. 0.1) THEN
1237 Zt = Z_0*EXP(KARMAN*2.0)
1238 Zt = MIN( Zt, 6.0e-5)
1239 Zt = MAX( Zt, 2.0e-9)
1240 Zq = Z_0*EXP(KARMAN*3.0)
1241 Zq = MIN( Zq, 6.0e-5)
1242 Zq = MAX( Zq, 2.0e-9)
1243 ELSE
1244 Zt = Z_0*EXP(-KARMAN*(4.0*SQRT(restar)-3.2))
1245 Zt = MIN( Zt, 6.0e-5)
1246 Zt = MAX( Zt, 2.0e-9)
1247 Zq = Z_0*EXP(-KARMAN*(4.0*SQRT(restar)-4.2))
1248 Zq = MIN( Zt, 6.0e-5)
1249 Zq = MAX( Zt, 2.0e-9)
1250 ENDIF
1252 ELSE !LAND
1254 !Option to modify CZIL according to Chen & Zhang, 2009
1255 IF ( IZ0TLND2 .EQ. 1 ) THEN
1256 CZIL = 10.0 ** ( -0.40 * ( Z_0 / 0.07 ) )
1257 ELSE
1258 CZIL = 0.085 !0.075 !0.10
1259 END IF
1261 Zt = Z_0*EXP(-KARMAN*CZIL*SQRT(restar))
1262 Zt = MIN( Zt, 0.75*Z_0)
1264 Zq = Z_0*EXP(-KARMAN*CZIL*SQRT(restar))
1265 Zq = MIN( Zq, 0.75*Z_0)
1267 ! stochastically perturb thermal and moisture roughness length.
1268 ! currently set to half the amplitude:
1269 if (spp_pbl==1) then
1270 Zt = Zt + Zt * 0.5 * rstoch
1271 Zt = MAX(Zt, 0.0001)
1272 Zq = Zt
1273 endif
1275 ENDIF
1277 return
1279 END SUBROUTINE zilitinkevich_1995
1280 !--------------------------------------------------------------------
1281 SUBROUTINE davis_etal_2008(Z_0,ustar)
1283 !a.k.a. : Donelan et al. (2004)
1284 !This formulation for roughness length was designed to match
1285 !the labratory experiments of Donelan et al. (2004).
1286 !This is an update version from Davis et al. 2008, which
1287 !corrects a small-bias in Z_0 (AHW real-time 2012).
1289 IMPLICIT NONE
1290 REAL, INTENT(IN) :: ustar
1291 REAL, INTENT(OUT) :: Z_0
1292 REAL :: ZW, ZN1, ZN2
1293 REAL, PARAMETER :: G=9.81, OZO=1.59E-5
1295 !OLD FORM: Z_0 = 10.*EXP(-10./(ustar**(1./3.)))
1296 !NEW FORM:
1298 ZW = MIN((ustar/1.06)**(0.3),1.0)
1299 ZN1 = 0.011*ustar*ustar/G + OZO
1300 ZN2 = 10.*exp(-9.5*ustar**(-.3333)) + &
1301 0.11*1.5E-5/AMAX1(ustar,0.01)
1302 Z_0 = (1.0-ZW) * ZN1 + ZW * ZN2
1304 Z_0 = MAX( Z_0, 1.27e-7) !These max/mins were suggested by
1305 Z_0 = MIN( Z_0, 2.85e-3) !Davis et al. (2008)
1307 return
1309 END SUBROUTINE davis_etal_2008
1310 !--------------------------------------------------------------------
1311 SUBROUTINE Taylor_Yelland_2001(Z_0,ustar,wsp10)
1313 !This formulation for roughness length was designed account for
1314 !wave steepness.
1316 IMPLICIT NONE
1317 REAL, INTENT(IN) :: ustar,wsp10
1318 REAL, INTENT(OUT) :: Z_0
1319 REAL, parameter :: g=9.81, pi=3.14159265
1320 REAL :: hs, Tp, Lp
1322 !hs is the significant wave height
1323 hs = 0.0248*(wsp10**2.)
1324 !Tp dominant wave period
1325 Tp = 0.729*MAX(wsp10,0.1)
1326 !Lp is the wavelength of the dominant wave
1327 Lp = g*Tp**2/(2*pi)
1329 Z_0 = 1200.*hs*(hs/Lp)**4.5
1330 Z_0 = MAX( Z_0, 1.27e-7) !These max/mins were suggested by
1331 Z_0 = MIN( Z_0, 2.85e-3) !Davis et al. (2008)
1333 return
1335 END SUBROUTINE Taylor_Yelland_2001
1336 !--------------------------------------------------------------------
1337 SUBROUTINE charnock_1955(Z_0,ustar,wsp10,visc,zu)
1339 !This version of Charnock's relation employs a varying
1340 !Charnock parameter, similar to COARE3.0 [Fairall et al. (2003)].
1341 !The Charnock parameter CZC is varied from .011 to .018
1342 !between 10-m wsp = 10 and 18.
1344 IMPLICIT NONE
1345 REAL, INTENT(IN) :: ustar, visc, wsp10, zu
1346 REAL, INTENT(OUT) :: Z_0
1347 REAL, PARAMETER :: G=9.81, CZO2=0.011
1348 REAL :: CZC !variable charnock "constant"
1349 REAL :: wsp10m ! logarithmically calculated 10 m
1351 wsp10m = wsp10*log(10./1e-4)/log(zu/1e-4)
1352 CZC = CZO2 + 0.007*MIN(MAX((wsp10m-10.)/8., 0.), 1.0)
1354 Z_0 = CZC*ustar*ustar/G + (0.11*visc/MAX(ustar,0.05))
1355 Z_0 = MAX( Z_0, 1.27e-7) !These max/mins were suggested by
1356 Z_0 = MIN( Z_0, 2.85e-3) !Davis et al. (2008)
1358 return
1360 END SUBROUTINE charnock_1955
1361 !--------------------------------------------------------------------
1362 SUBROUTINE edson_etal_2013(Z_0,ustar,wsp10,visc,zu)
1364 !This version of Charnock's relation employs a varying
1365 !Charnock parameter, taken from COARE 3.5 [Edson et al. (2001, JPO)].
1366 !The Charnock parameter CZC is varied from about .005 to .028
1367 !between 10-m wind speeds of 6 and 19 m/s.
1368 !11 Nov 2021: Note that this was finally fixed according to the
1369 ! Edson et al (2014) corrigendum, where "m" was corrected.
1371 IMPLICIT NONE
1372 REAL, INTENT(IN) :: ustar, visc, wsp10, zu
1373 REAL, INTENT(OUT) :: Z_0
1374 REAL, PARAMETER :: G=9.81
1375 REAL, PARAMETER :: m=0.0017, b=-0.005
1376 REAL :: CZC ! variable charnock "constant"
1377 REAL :: wsp10m ! logarithmically calculated 10 m
1379 wsp10m = wsp10*log(10/1e-4)/log(zu/1e-4)
1380 wsp10m = MIN(19., wsp10m)
1381 CZC = m*wsp10m + b
1382 CZC = MAX(CZC, 0.0)
1384 Z_0 = CZC*ustar*ustar/G + (0.11*visc/MAX(ustar,0.07))
1385 Z_0 = MAX( Z_0, 1.27e-7) !These max/mins were suggested by
1386 Z_0 = MIN( Z_0, 2.85e-3) !Davis et al. (2008)
1388 return
1390 END SUBROUTINE edson_etal_2013
1391 !--------------------------------------------------------------------
1392 SUBROUTINE garratt_1992(Zt,Zq,Z_0,Ren,landsea)
1394 !This formulation for the thermal and moisture roughness lengths
1395 !(Zt and Zq) relates them to Z0 via the roughness Reynolds number (Ren).
1396 !This formula comes from Fairall et al. (2003). It is modified from
1397 !the original Garratt-Brutsaert model to better fit the COARE/HEXMAX
1398 !data. The formula for land uses a constant ratio (Z_0/7.4) taken
1399 !from Garratt (1992).
1401 IMPLICIT NONE
1402 REAL, INTENT(IN) :: Ren, Z_0,landsea
1403 REAL, INTENT(OUT) :: Zt,Zq
1404 REAL :: Rq
1405 REAL, PARAMETER :: e=2.71828183
1407 IF (landsea-1.5 .GT. 0) THEN !WATER
1409 Zt = Z_0*EXP(2.0 - (2.48*(Ren**0.25)))
1410 Zq = Z_0*EXP(2.0 - (2.28*(Ren**0.25)))
1412 Zq = MIN( Zq, 5.5e-5)
1413 Zq = MAX( Zq, 2.0e-9)
1414 Zt = MIN( Zt, 5.5e-5)
1415 Zt = MAX( Zt, 2.0e-9) !same lower limit as ECMWF
1416 ELSE !LAND
1417 Zq = Z_0/(e**2.) !taken from Garratt (1980,1992)
1418 Zt = Zq
1419 ENDIF
1421 return
1423 END SUBROUTINE garratt_1992
1424 !--------------------------------------------------------------------
1425 SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl)
1427 !This formulation for thermal and moisture roughness length (Zt and Zq)
1428 !as a function of the roughness Reynolds number (Ren) comes from the
1429 !COARE3.0 formulation, empirically derived from COARE and HEXMAX data
1430 ![Fairall et al. (2003)]. Edson et al. (2004; JGR) suspected that this
1431 !relationship overestimated the scalar roughness lengths for low Reynolds
1432 !number flows, so an optional smooth flow relationship, taken from Garratt
1433 !(1992, p. 102), is available for flows with Ren < 2.
1434 !
1435 !This is for use over water only.
1437 IMPLICIT NONE
1438 REAL, INTENT(IN) :: Ren,ustar,visc,rstoch
1439 INTEGER, INTENT(IN):: spp_pbl
1440 REAL, INTENT(OUT) :: Zt,Zq
1442 IF (Ren .le. 2.) then
1444 Zt = (5.5e-5)*(Ren**(-0.60))
1445 Zq = Zt
1446 !FOR SMOOTH SEAS, CAN USE GARRATT
1447 !Zq = 0.2*visc/MAX(ustar,0.1)
1448 !Zq = 0.3*visc/MAX(ustar,0.1)
1450 ELSE
1452 !FOR ROUGH SEAS, USE COARE
1453 Zt = (5.5e-5)*(Ren**(-0.60))
1454 Zq = Zt
1456 ENDIF
1458 if (spp_pbl==1) then
1459 Zt = Zt + Zt * 0.5 * rstoch
1460 Zq = Zt
1461 endif
1463 Zt = MIN(Zt,1.0e-4)
1464 Zt = MAX(Zt,2.0e-9)
1466 Zq = MIN(Zt,1.0e-4)
1467 Zq = MAX(Zt,2.0e-9)
1469 return
1471 END SUBROUTINE fairall_etal_2003
1472 !--------------------------------------------------------------------
1473 SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl)
1475 !This formulation for thermal and moisture roughness length (Zt and Zq)
1476 !as a function of the roughness Reynolds number (Ren) comes from the
1477 !COARE 3.5/4.0 formulation, empirically derived from COARE and HEXMAX data
1478 ![Fairall et al. (2014? coming soon, not yet published as of July 2014)].
1479 !This is for use over water only.
1481 IMPLICIT NONE
1482 REAL, INTENT(IN) :: Ren,ustar,visc,rstoch
1483 INTEGER, INTENT(IN):: spp_pbl
1484 REAL, INTENT(OUT) :: Zt,Zq
1486 !Zt = (5.5e-5)*(Ren**(-0.60))
1487 Zt = MIN(1.6E-4, 5.8E-5/(Ren**0.72))
1488 Zq = Zt
1490 IF (spp_pbl ==1) THEN
1491 Zt = MAX(Zt + Zt*0.5*rstoch,2.0e-9)
1492 Zq = MAX(Zt + Zt*0.5*rstoch,2.0e-9)
1493 ELSE
1494 Zt = MAX(Zt,2.0e-9)
1495 Zq = MAX(Zt,2.0e-9)
1496 ENDIF
1498 return
1500 END SUBROUTINE fairall_etal_2014
1501 !--------------------------------------------------------------------
1502 SUBROUTINE Yang_2008(Z_0,Zt,Zq,ustar,tstar,qst,Ren,visc,landsea)
1504 !This is a modified version of Yang et al (2002 QJRMS, 2008 JAMC)
1505 !and Chen et al (2010, J of Hydromet). Although it was originally
1506 !designed for arid regions with bare soil, it is modified
1507 !here to perform over a broader spectrum of vegetation.
1508 !
1509 !The original formulation relates the thermal roughness length (Zt)
1510 !to u* and T*:
1511 !
1512 ! Zt = ht * EXP(-beta*(ustar**0.5)*(ABS(tstar)**0.25))
1513 !
1514 !where ht = Renc*visc/ustar and the critical Reynolds number
1515 !(Renc) = 70. Beta was originally = 10 (2002 paper) but was revised
1516 !to 7.2 (in 2008 paper). Their form typically varies the
1517 !ratio Z0/Zt by a few orders of magnitude (1-1E4).
1518 !
1519 !This modified form uses beta = 1.5 and a variable Renc (function of Z_0),
1520 !so zt generally varies similarly to the Zilitinkevich form (with Czil ~ 0.1)
1521 !for very small or negative surface heat fluxes but can become close to the
1522 !Zilitinkevich with Czil = 0.2 for very large HFX (large negative T*).
1523 !Also, the exponent (0.25) on tstar was changed to 1.0, since we found
1524 !Zt was reduced too much for low-moderate positive heat fluxes.
1525 !
1526 !This should only be used over land!
1528 IMPLICIT NONE
1529 REAL, INTENT(IN) :: Z_0, Ren, ustar, tstar, qst, visc, landsea
1530 REAL :: ht, &! roughness height at critical Reynolds number
1531 tstar2, &! bounded T*, forced to be non-positive
1532 qstar2, &! bounded q*, forced to be non-positive
1533 Z_02, &! bounded Z_0 for variable Renc2 calc
1534 Renc2 ! variable Renc, function of Z_0
1535 REAL, INTENT(OUT) :: Zt,Zq
1536 REAL, PARAMETER :: Renc=300., & !old constant Renc
1537 beta=1.5, & !important for diurnal variation
1538 m=170., & !slope for Renc2 function
1539 b=691. !y-intercept for Renc2 function
1541 Z_02 = MIN(Z_0,0.5)
1542 Z_02 = MAX(Z_02,0.04)
1543 Renc2= b + m*log(Z_02)
1544 ht = Renc2*visc/MAX(ustar,0.01)
1545 tstar2 = MIN(tstar, 0.0)
1546 qstar2 = MIN(qst,0.0)
1548 Zt = ht * EXP(-beta*(ustar**0.5)*(ABS(tstar2)**1.0))
1549 Zq = ht * EXP(-beta*(ustar**0.5)*(ABS(qstar2)**1.0))
1550 !Zq = Zt
1552 Zt = MIN(Zt, Z_0/2.0)
1553 Zq = MIN(Zq, Z_0/2.0)
1555 return
1557 END SUBROUTINE Yang_2008
1558 !--------------------------------------------------------------------
1559 SUBROUTINE Andreas_2002(Z_0,bvisc,ustar,Zt,Zq)
1561 ! This is taken from Andreas (2002; J. of Hydromet) and
1562 ! Andreas et al. (2005; BLM).
1563 !
1564 ! This should only be used over snow/ice!
1566 IMPLICIT NONE
1567 REAL, INTENT(IN) :: Z_0, bvisc, ustar
1568 REAL, INTENT(OUT) :: Zt, Zq
1569 REAL :: Ren2, zntsno
1571 REAL, PARAMETER :: bt0_s=1.25, bt0_t=0.149, bt0_r=0.317, &
1572 bt1_s=0.0, bt1_t=-0.55, bt1_r=-0.565, &
1573 bt2_s=0.0, bt2_t=0.0, bt2_r=-0.183
1575 REAL, PARAMETER :: bq0_s=1.61, bq0_t=0.351, bq0_r=0.396, &
1576 bq1_s=0.0, bq1_t=-0.628, bq1_r=-0.512, &
1577 bq2_s=0.0, bq2_t=0.0, bq2_r=-0.180
1579 !Calculate zo for snow (Andreas et al. 2005, BLM)
1580 zntsno = 0.135*bvisc/ustar + &
1581 (0.035*(ustar*ustar)/9.8) * &
1582 (5.*exp(-1.*(((ustar - 0.18)/0.1)*((ustar - 0.18)/0.1))) + 1.)
1583 Ren2 = ustar*zntsno/bvisc
1585 ! Make sure that Re is not outside of the range of validity
1586 ! for using their equations
1587 IF (Ren2 .gt. 1000.) Ren2 = 1000.
1589 IF (Ren2 .le. 0.135) then
1591 Zt = zntsno*EXP(bt0_s + bt1_s*LOG(Ren2) + bt2_s*LOG(Ren2)**2)
1592 Zq = zntsno*EXP(bq0_s + bq1_s*LOG(Ren2) + bq2_s*LOG(Ren2)**2)
1594 ELSE IF (Ren2 .gt. 0.135 .AND. Ren2 .lt. 2.5) then
1596 Zt = zntsno*EXP(bt0_t + bt1_t*LOG(Ren2) + bt2_t*LOG(Ren2)**2)
1597 Zq = zntsno*EXP(bq0_t + bq1_t*LOG(Ren2) + bq2_t*LOG(Ren2)**2)
1599 ELSE
1601 Zt = zntsno*EXP(bt0_r + bt1_r*LOG(Ren2) + bt2_r*LOG(Ren2)**2)
1602 Zq = zntsno*EXP(bq0_r + bq1_r*LOG(Ren2) + bq2_r*LOG(Ren2)**2)
1604 ENDIF
1606 return
1608 END SUBROUTINE Andreas_2002
1609 !--------------------------------------------------------------------
1610 SUBROUTINE PSI_Hogstrom_1996(psi_m, psi_h, zL, Zt, Z_0, Za)
1612 ! This subroutine returns the stability functions based off
1613 ! of Hogstrom (1996).
1615 IMPLICIT NONE
1616 REAL, INTENT(IN) :: zL, Zt, Z_0, Za
1617 REAL, INTENT(OUT) :: psi_m, psi_h
1618 REAL :: x, x0, y, y0, zmL, zhL
1620 zmL = Z_0*zL/Za
1621 zhL = Zt*zL/Za
1623 IF (zL .gt. 0.) THEN !STABLE (not well tested - seem large)
1625 psi_m = -5.3*(zL - zmL)
1626 psi_h = -8.0*(zL - zhL)
1628 ELSE !UNSTABLE
1630 x = (1.-19.0*zL)**0.25
1631 x0= (1.-19.0*zmL)**0.25
1632 y = (1.-11.6*zL)**0.5
1633 y0= (1.-11.6*zhL)**0.5
1635 psi_m = 2.*LOG((1.+x)/(1.+x0)) + &
1636 &LOG((1.+x**2.)/(1.+x0**2.)) - &
1637 &2.0*ATAN(x) + 2.0*ATAN(x0)
1638 psi_h = 2.*LOG((1.+y)/(1.+y0))
1640 ENDIF
1642 return
1644 END SUBROUTINE PSI_Hogstrom_1996
1645 !--------------------------------------------------------------------
1646 SUBROUTINE PSI_DyerHicks(psi_m, psi_h, zL, Zt, Z_0, Za)
1648 ! This subroutine returns the stability functions based off
1649 ! of Hogstrom (1996), but with different constants compatible
1650 ! with Dyer and Hicks (1970/74?). This formulation is used for
1651 ! testing/development by Nakanishi (personal communication).
1653 IMPLICIT NONE
1654 REAL, INTENT(IN) :: zL, Zt, Z_0, Za
1655 REAL, INTENT(OUT) :: psi_m, psi_h
1656 REAL :: x, x0, y, y0, zmL, zhL
1658 zmL = Z_0*zL/Za !Zo/L
1659 zhL = Zt*zL/Za !Zt/L
1661 IF (zL .gt. 0.) THEN !STABLE
1663 psi_m = -5.0*(zL - zmL)
1664 psi_h = -5.0*(zL - zhL)
1666 ELSE !UNSTABLE
1668 x = (1.-16.*zL)**0.25
1669 x0= (1.-16.*zmL)**0.25
1671 y = (1.-16.*zL)**0.5
1672 y0= (1.-16.*zhL)**0.5
1674 psi_m = 2.*LOG((1.+x)/(1.+x0)) + &
1675 &LOG((1.+x**2.)/(1.+x0**2.)) - &
1676 &2.0*ATAN(x) + 2.0*ATAN(x0)
1677 psi_h = 2.*LOG((1.+y)/(1.+y0))
1679 ENDIF
1681 return
1683 END SUBROUTINE PSI_DyerHicks
1684 !--------------------------------------------------------------------
1685 SUBROUTINE PSI_Beljaars_Holtslag_1991(psi_m, psi_h, zL)
1687 ! This subroutine returns the stability functions based off
1688 ! of Beljaar and Holtslag 1991, which is an extension of Holtslag
1689 ! and Debruin 1989.
1691 IMPLICIT NONE
1692 REAL, INTENT(IN) :: zL
1693 REAL, INTENT(OUT) :: psi_m, psi_h
1694 REAL, PARAMETER :: a=1., b=0.666, c=5., d=0.35
1696 IF (zL .lt. 0.) THEN !UNSTABLE
1698 WRITE(*,*)"WARNING: Universal stability functions from"
1699 WRITE(*,*)" Beljaars and Holtslag (1991) should only"
1700 WRITE(*,*)" be used in the stable regime!"
1701 psi_m = 0.
1702 psi_h = 0.
1704 ELSE !STABLE
1706 psi_m = -(a*zL + b*(zL -(c/d))*exp(-d*zL) + (b*c/d))
1707 psi_h = -((1.+.666*a*zL)**1.5 + &
1708 b*(zL - (c/d))*exp(-d*zL) + (b*c/d) -1.)
1710 ENDIF
1712 return
1714 END SUBROUTINE PSI_Beljaars_Holtslag_1991
1715 !--------------------------------------------------------------------
1716 SUBROUTINE PSI_Zilitinkevich_Esau_2007(psi_m, psi_h, zL)
1718 ! This subroutine returns the stability functions come from
1719 ! Zilitinkevich and Esau (2007, BM), which are formulatioed from the
1720 ! "generalized similarity theory" and tuned to the LES DATABASE64
1721 ! to determine their dependence on z/L.
1723 IMPLICIT NONE
1724 REAL, INTENT(IN) :: zL
1725 REAL, INTENT(OUT) :: psi_m, psi_h
1726 REAL, PARAMETER :: Cm=3.0, Ct=2.5
1728 IF (zL .lt. 0.) THEN !UNSTABLE
1730 WRITE(*,*)"WARNING: Universal stability function from"
1731 WRITE(*,*)" Zilitinkevich and Esau (2007) should only"
1732 WRITE(*,*)" be used in the stable regime!"
1733 psi_m = 0.
1734 psi_h = 0.
1736 ELSE !STABLE
1738 psi_m = -Cm*(zL**(5./6.))
1739 psi_h = -Ct*(zL**(4./5.))
1741 ENDIF
1743 return
1745 END SUBROUTINE PSI_Zilitinkevich_Esau_2007
1746 !--------------------------------------------------------------------
1747 SUBROUTINE PSI_Businger_1971(psi_m, psi_h, zL)
1749 ! This subroutine returns the flux-profile relationships
1750 ! of Businger el al. 1971.
1752 IMPLICIT NONE
1753 REAL, INTENT(IN) :: zL
1754 REAL, INTENT(OUT) :: psi_m, psi_h
1755 REAL :: x, y
1756 REAL, PARAMETER :: Pi180 = 3.14159265/180.
1758 IF (zL .lt. 0.) THEN !UNSTABLE
1760 x = (1. - 15.0*zL)**0.25
1761 y = (1. - 9.0*zL)**0.5
1763 psi_m = LOG(((1.+x)/2.)**2.) + &
1764 &LOG((1.+x**2.)/2.) - &
1765 &2.0*ATAN(x) + Pi180*90.
1766 psi_h = 2.*LOG((1.+y)/2.)
1768 ELSE !STABLE
1770 psi_m = -4.7*zL
1771 psi_h = -(4.7/0.74)*zL
1773 ENDIF
1775 return
1777 END SUBROUTINE PSI_Businger_1971
1778 !--------------------------------------------------------------------
1779 SUBROUTINE PSI_Suselj_Sood_2010(psi_m, psi_h, zL)
1781 !This subroutine returns flux-profile relatioships based off
1782 !of Lobocki (1993), which is derived from the MY-level 2 model.
1783 !Suselj and Sood (2010) applied the surface layer length scales
1784 !from Nakanishi (2001) to get this new relationship. These functions
1785 !are more agressive (larger magnitude) than most formulations. They
1786 !showed improvement over water, but untested over land.
1788 IMPLICIT NONE
1789 REAL, INTENT(IN) :: zL
1790 REAL, INTENT(OUT) :: psi_m, psi_h
1791 REAL, PARAMETER :: Rfc=0.19, Ric=0.183, PHIT=0.8
1793 IF (zL .gt. 0.) THEN !STABLE
1795 psi_m = -(zL/Rfc + 1.1223*EXP(1.-1.6666/zL))
1796 !psi_h = -zL*Ric/((Rfc**2.)*PHIT) + 8.209*(zL**1.1091)
1797 !THEIR EQ FOR PSI_H CRASHES THE MODEL AND DOES NOT MATCH
1798 !THEIR FIG 1. THIS EQ (BELOW) MATCHES THEIR FIG 1 BETTER:
1799 psi_h = -(zL*Ric/((Rfc**2.)*5.) + 7.09*(zL**1.1091))
1801 ELSE !UNSTABLE
1803 psi_m = 0.9904*LOG(1. - 14.264*zL)
1804 psi_h = 1.0103*LOG(1. - 16.3066*zL)
1806 ENDIF
1808 return
1810 END SUBROUTINE PSI_Suselj_Sood_2010
1811 !--------------------------------------------------------------------
1812 SUBROUTINE PSI_CB2005(psim1,psih1,zL,z0L)
1814 ! This subroutine returns the stability functions based off
1815 ! of Cheng and Brutseart (2005, BLM), for use in stable conditions only.
1816 ! The returned values are the combination of psi((za+zo)/L) - psi(z0/L)
1818 IMPLICIT NONE
1819 REAL, INTENT(IN) :: zL,z0L
1820 REAL, INTENT(OUT) :: psim1,psih1
1822 psim1 = -6.1*LOG(zL + (1.+ zL **2.5)**0.4) &
1823 -6.1*LOG(z0L + (1.+ z0L**2.5)**0.4)
1824 psih1 = -5.5*LOG(zL + (1.+ zL **1.1)**0.90909090909) &
1825 -5.5*LOG(z0L + (1.+ z0L**1.1)**0.90909090909)
1827 return
1829 END SUBROUTINE PSI_CB2005
1830 !--------------------------------------------------------------------
1831 SUBROUTINE Li_etal_2010(zL, Rib, zaz0, z0zt)
1833 !This subroutine returns a more robust z/L that best matches
1834 !the z/L from Hogstrom (1996) for unstable conditions and Beljaars
1835 !and Holtslag (1991) for stable conditions.
1837 IMPLICIT NONE
1838 REAL, INTENT(OUT) :: zL
1839 REAL, INTENT(IN) :: Rib, zaz0, z0zt
1840 REAL :: alfa, beta, zaz02, z0zt2
1841 REAL, PARAMETER :: au11=0.045, bu11=0.003, bu12=0.0059, &
1842 &bu21=-0.0828, bu22=0.8845, bu31=0.1739, &
1843 &bu32=-0.9213, bu33=-0.1057
1844 REAL, PARAMETER :: aw11=0.5738, aw12=-0.4399, aw21=-4.901,&
1845 &aw22=52.50, bw11=-0.0539, bw12=1.540, &
1846 &bw21=-0.669, bw22=-3.282
1847 REAL, PARAMETER :: as11=0.7529, as21=14.94, bs11=0.1569,&
1848 &bs21=-0.3091, bs22=-1.303
1850 !set limits according to Li et al (2010), p 157.
1851 zaz02=zaz0
1852 IF (zaz0 .lt. 100.0) zaz02=100.
1853 IF (zaz0 .gt. 100000.0) zaz02=100000.
1855 !set more limits according to Li et al (2010)
1856 z0zt2=z0zt
1857 IF (z0zt .lt. 0.5) z0zt2=0.5
1858 IF (z0zt .gt. 100.0) z0zt2=100.
1860 alfa = LOG(zaz02)
1861 beta = LOG(z0zt2)
1863 IF (Rib .le. 0.0) THEN
1864 zL = au11*alfa*Rib**2 + ( &
1865 & (bu11*beta + bu12)*alfa**2 + &
1866 & (bu21*beta + bu22)*alfa + &
1867 & (bu31*beta**2 + bu32*beta + bu33))*Rib
1868 !if(zL .LT. -15 .OR. zl .GT. 0.)print*,"VIOLATION Rib<0:",zL
1869 zL = MAX(zL,-15.) !LIMITS SET ACCORDING TO Li et al (2010)
1870 zL = MIN(zL,0.) !Figure 1.
1871 ELSEIF (Rib .gt. 0.0 .AND. Rib .le. 0.2) THEN
1872 zL = ((aw11*beta + aw12)*alfa + &
1873 & (aw21*beta + aw22))*Rib**2 + &
1874 & ((bw11*beta + bw12)*alfa + &
1875 & (bw21*beta + bw22))*Rib
1876 !if(zL .LT. 0 .OR. zl .GT. 4)print*,"VIOLATION 0<Rib<0.2:",zL
1877 zL = MIN(zL,4.) !LIMITS APPROX SET ACCORDING TO Li et al (2010)
1878 zL = MAX(zL,0.) !THEIR FIGURE 1B.
1879 ELSE
1880 zL = (as11*alfa + as21)*Rib + bs11*alfa + &
1881 & bs21*beta + bs22
1882 !if(zL .LE. 1 .OR. zl .GT. 23)print*,"VIOLATION Rib>0.2:",zL
1883 zL = MIN(zL,20.) !LIMITS ACCORDING TO Li et al (2010), THIER
1884 !FIGUE 1C.
1885 zL = MAX(zL,1.)
1886 ENDIF
1888 return
1890 END SUBROUTINE Li_etal_2010
1891 !-------------------------------------------------------------------
1892 REAL function zolri(ri,za,z0,zt,zol1,psi_opt)
1894 ! This iterative algorithm is a two-point secant method taken from the revised
1895 ! surface layer scheme in WRF-ARW, written by Pedro Jimenez and Jimy Dudhia and
1896 ! summarized in Jimenez et al. (2012, MWR). This function was adapted
1897 ! to input the thermal roughness length, zt, (as well as z0) and use initial
1898 ! estimate of z/L.
1900 IMPLICIT NONE
1901 REAL, INTENT(IN) :: ri,za,z0,zt,zol1
1902 INTEGER, INTENT(IN) :: psi_opt
1903 REAL :: x1,x2,fx1,fx2
1904 INTEGER :: n
1905 INTEGER, PARAMETER :: nmax = 20
1907 if (ri.lt.0.)then
1908 x1=zol1 - 0.02 !-5.
1909 x2=0.
1910 else
1911 x1=0.
1912 x2=zol1 + 0.02 !5.
1913 endif
1915 n=0
1916 fx1=zolri2(x1,ri,za,z0,zt,psi_opt)
1917 fx2=zolri2(x2,ri,za,z0,zt,psi_opt)
1919 Do While (abs(x1 - x2) > 0.01 .and. n < nmax)
1920 if(abs(fx2) .lt. abs(fx1))then
1921 x1=x1-fx1/(fx2-fx1)*(x2-x1)
1922 fx1=zolri2(x1,ri,za,z0,zt,psi_opt)
1923 zolri=x1
1924 else
1925 x2=x2-fx2/(fx2-fx1)*(x2-x1)
1926 fx2=zolri2(x2,ri,za,z0,zt,psi_opt)
1927 zolri=x2
1928 endif
1929 n=n+1
1930 enddo
1932 if (n==nmax .and. abs(x1 - x2) >= 0.01) then
1933 !if convergence fails, use approximate values:
1934 CALL Li_etal_2010(zolri, ri, za/z0, z0/zt)
1935 !print*,"FAILED, n=",n," Ri=",ri," zt=",zt
1936 else
1937 !print*,"SUCCESS,n=",n," Ri=",ri," z/L=",zolri
1938 endif
1940 return
1941 end function
1942 !-------------------------------------------------------------------
1943 REAL function zolri2(zol2,ri2,za,z0,zt,psi_opt)
1945 ! INPUT: =================================
1946 ! zol2 - estimated z/L
1947 ! ri2 - calculated bulk Richardson number
1948 ! za - 1/2 depth of first model layer
1949 ! z0 - aerodynamic roughness length
1950 ! zt - thermal roughness length
1951 ! OUTPUT: ================================
1952 ! zolri2 - delta Ri
1954 IMPLICIT NONE
1955 INTEGER, INTENT(IN) :: psi_opt
1956 REAL, INTENT(IN) :: ri2,za,z0,zt
1957 REAL, INTENT(INOUT) :: zol2
1958 REAL :: zol20,zol3,psim1,psih1,psix2,psit2,zolt
1960 if(zol2*ri2 .lt. 0.) THEN
1961 !print*,"WRONG QUADRANTS: z/L=",zol2," ri=",ri2
1962 zol2=0.
1963 endif
1965 zol20=zol2*z0/za ! z0/L
1966 zol3=zol2+zol20 ! (z+z0)/L
1967 zolt=zol2*zt/za ! zt/L
1969 if (ri2.lt.0) then
1970 psit2=MAX(log((za+z0)/zt)-(psih_unstable(zol3,psi_opt)-psih_unstable(zolt,psi_opt)), 1.0)
1971 psix2=MAX(log((za+z0)/z0)-(psim_unstable(zol3,psi_opt)-psim_unstable(zol20,psi_opt)),1.0)
1972 else
1973 psit2=MAX(log((za+z0)/zt)-(psih_stable(zol3,psi_opt)-psih_stable(zolt,psi_opt)), 1.0)
1974 psix2=MAX(log((za+z0)/z0)-(psim_stable(zol3,psi_opt)-psim_stable(zol20,psi_opt)),1.0)
1975 endif
1977 zolri2=zol2*psit2/psix2**2 - ri2
1978 !print*," target ri=",ri2," est ri=",zol2*psit2/psix2**2
1980 return
1981 end function
1982 !====================================================================
1984 REAL function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
1986 ! This iterative algorithm to compute z/L from bulk-Ri
1988 IMPLICIT NONE
1989 REAL, INTENT(IN) :: ri,za,z0,zt,logz0,logzt
1990 INTEGER, INTENT(IN) :: psi_opt
1991 REAL, INTENT(INOUT) :: zol1
1992 REAL :: zol20,zol3,zolt,zolold
1993 INTEGER :: n
1994 INTEGER, PARAMETER :: nmax = 20
1995 !REAL, DIMENSION(nmax):: zLhux
1996 REAL :: psit2,psix2
1998 if(zol1*ri .lt. 0.) THEN
1999 !print*,"WRONG QUADRANTS: z/L=",zol1," ri=",ri
2000 zol1=0.
2001 endif
2003 if (ri .lt. 0.) then
2004 zolold=-99999.
2005 zolrib=-66666.
2006 else
2007 zolold=99999.
2008 zolrib=66666.
2009 endif
2010 n=1
2012 DO While (abs(zolold - zolrib) > 0.01 .and. n < nmax)
2014 if(n==1)then
2015 zolold=zol1
2016 else
2017 zolold=zolrib
2018 endif
2019 zol20=zolold*z0/za ! z0/L
2020 zol3=zolold+zol20 ! (z+z0)/L
2021 zolt=zolold*zt/za ! zt/L
2023 if (ri.lt.0) then
2024 psit2=MAX(logzt-(psih_unstable(zol3,psi_opt)-psih_unstable(zolt,psi_opt)), 1.0)
2025 psix2=MAX(logz0-(psim_unstable(zol3,psi_opt)-psim_unstable(zol20,psi_opt)), 1.0)
2026 else
2027 psit2=MAX(logzt-(psih_stable(zol3,psi_opt)-psih_stable(zolt,psi_opt)), 1.0)
2028 psix2=MAX(logz0-(psim_stable(zol3,psi_opt)-psim_stable(zol20,psi_opt)), 1.0)
2029 endif
2031 zolrib=ri*psix2**2/psit2
2032 !zLhux(n)=zolrib
2033 n=n+1
2034 enddo
2036 if (n==nmax .and. abs(zolold - zolrib) > 0.01 ) then
2037 !print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri
2038 !if convergence fails, use approximate values:
2039 CALL Li_etal_2010(zolrib, ri, za/z0, z0/zt)
2040 !zLhux(n)=zolri
2041 !print*,"FAILED, n=",n," Ri=",ri," zt=",zt
2042 !print*,"z/L=",zLhux(1:nmax)
2043 else
2044 !print*,"SUCCESS,n=",n," Ri=",ri," z/L=",zolrib
2045 endif
2047 return
2048 end function
2049 !====================================================================
2051 SUBROUTINE psi_init(psi_opt)
2053 ! Define tables from -10 <= z/L <= 10
2055 INTEGER :: N,psi_opt
2056 REAL :: zolf
2058 if (psi_opt == 0) then
2059 DO N=0,1000
2060 ! stable function tables
2061 zolf = float(n)*0.01
2062 psim_stab(n)=psim_stable_full(zolf)
2063 psih_stab(n)=psih_stable_full(zolf)
2065 ! unstable function tables
2066 zolf = -float(n)*0.01
2067 psim_unstab(n)=psim_unstable_full(zolf)
2068 psih_unstab(n)=psih_unstable_full(zolf)
2069 ENDDO
2070 else
2071 DO N=0,1000
2072 ! stable function tables
2073 zolf = float(n)*0.01
2074 psim_stab(n)=psim_stable_full_gfs(zolf)
2075 psih_stab(n)=psih_stable_full_gfs(zolf)
2077 ! unstable function tables
2078 zolf = -float(n)*0.01
2079 psim_unstab(n)=psim_unstable_full_gfs(zolf)
2080 psih_unstab(n)=psih_unstable_full_gfs(zolf)
2081 ENDDO
2082 endif
2084 END SUBROUTINE psi_init
2085 ! ==================================================================
2086 ! ... Full equations for the integrated similarity functions ...
2087 ! ==================================================================
2088 REAL function psim_stable_full(zolf)
2089 REAL :: zolf
2091 psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
2093 return
2094 end function
2096 REAL function psih_stable_full(zolf)
2097 REAL :: zolf
2099 psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
2101 return
2102 end function
2104 REAL function psim_unstable_full(zolf)
2105 REAL :: zolf,x,ym,psimc,psimk
2107 x=(1.-16.*zolf)**.25
2108 psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*ATAN(1.)
2110 ym=(1.-10.*zolf)**0.33
2111 psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
2113 psim_unstable_full=(psimk+zolf**2*(psimc))/(1+zolf**2.)
2115 return
2116 end function
2118 REAL function psih_unstable_full(zolf)
2119 REAL :: zolf,y,yh,psihc,psihk
2121 y=(1.-16.*zolf)**.5
2122 psihk=2.*log((1+y)/2.)
2124 yh=(1.-34.*zolf)**0.33
2125 psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
2127 psih_unstable_full=(psihk+zolf**2*(psihc))/(1+zolf**2.)
2129 return
2130 end function
2132 ! ==================================================================
2133 ! ... integrated similarity functions from GFS...
2134 !
2135 REAL function psim_stable_full_gfs(zolf)
2136 REAL :: zolf
2137 REAL, PARAMETER :: alpha4 = 20.
2138 REAL :: aa
2140 aa = sqrt(1. + alpha4 * zolf)
2141 psim_stable_full_gfs = -1.*aa + log(aa + 1.)
2143 return
2144 end function
2146 REAL function psih_stable_full_gfs(zolf)
2147 REAL :: zolf
2148 REAL, PARAMETER :: alpha4 = 20.
2149 REAL :: bb
2151 bb = sqrt(1. + alpha4 * zolf)
2152 psih_stable_full_gfs = -1.*bb + log(bb + 1.)
2154 return
2155 end function
2157 REAL function psim_unstable_full_gfs(zolf)
2158 REAL :: zolf
2159 REAL :: hl1,tem1
2160 REAL, PARAMETER :: a0=-3.975, a1=12.32, &
2161 b1=-7.755, b2=6.041
2163 if (zolf .ge. -0.5) then
2164 hl1 = zolf
2165 psim_unstable_full_gfs = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1)
2166 else
2167 hl1 = -zolf
2168 tem1 = 1.0 / sqrt(hl1)
2169 psim_unstable_full_gfs = log(hl1) + 2. * sqrt(tem1) - .8776
2170 end if
2172 return
2173 end function
2175 REAL function psih_unstable_full_gfs(zolf)
2176 REAL :: zolf
2177 REAL :: hl1,tem1
2178 REAL, PARAMETER :: a0p=-7.941, a1p=24.75, &
2179 b1p=-8.705, b2p=7.899
2181 if (zolf .ge. -0.5) then
2182 hl1 = zolf
2183 psih_unstable_full_gfs = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1)
2184 else
2185 hl1 = -zolf
2186 tem1 = 1.0 / sqrt(hl1)
2187 psih_unstable_full_gfs = log(hl1) + .5 * tem1 + 1.386
2188 end if
2190 return
2191 end function
2193 !=================================================================
2194 ! These functions use the look-up table functions when |z/L| <= 10
2195 ! but default to the full equations when |z/L| > 10.
2196 !=================================================================
2197 REAL function psim_stable(zolf,psi_opt)
2198 integer :: nzol,psi_opt
2199 real :: rzol,zolf
2201 nzol = int(zolf*100.)
2202 rzol = zolf*100. - nzol
2203 if(nzol+1 .le. 1000)then
2204 psim_stable = psim_stab(nzol) + rzol*(psim_stab(nzol+1)-psim_stab(nzol))
2205 else
2206 if (psi_opt == 0) then
2207 psim_stable = psim_stable_full(zolf)
2208 else
2209 psim_stable = psim_stable_full_gfs(zolf)
2210 endif
2211 endif
2213 return
2214 end function
2216 REAL function psih_stable(zolf,psi_opt)
2217 integer :: nzol,psi_opt
2218 real :: rzol,zolf
2220 nzol = int(zolf*100.)
2221 rzol = zolf*100. - nzol
2222 if(nzol+1 .le. 1000)then
2223 psih_stable = psih_stab(nzol) + rzol*(psih_stab(nzol+1)-psih_stab(nzol))
2224 else
2225 if (psi_opt == 0) then
2226 psih_stable = psih_stable_full(zolf)
2227 else
2228 psih_stable = psih_stable_full_gfs(zolf)
2229 endif
2230 endif
2232 return
2233 end function
2235 REAL function psim_unstable(zolf,psi_opt)
2236 integer :: nzol,psi_opt
2237 real :: rzol,zolf
2239 nzol = int(-zolf*100.)
2240 rzol = -zolf*100. - nzol
2241 if(nzol+1 .le. 1000)then
2242 psim_unstable = psim_unstab(nzol) + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol))
2243 else
2244 if (psi_opt == 0) then
2245 psim_unstable = psim_unstable_full(zolf)
2246 else
2247 psim_unstable = psim_unstable_full_gfs(zolf)
2248 endif
2249 endif
2251 return
2252 end function
2254 REAL function psih_unstable(zolf,psi_opt)
2255 integer :: nzol,psi_opt
2256 real :: rzol,zolf
2258 nzol = int(-zolf*100.)
2259 rzol = -zolf*100. - nzol
2260 if(nzol+1 .le. 1000)then
2261 psih_unstable = psih_unstab(nzol) + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol))
2262 else
2263 if (psi_opt == 0) then
2264 psih_unstable = psih_unstable_full(zolf)
2265 else
2266 psih_unstable = psih_unstable_full_gfs(zolf)
2267 endif
2268 endif
2270 return
2271 end function
2272 !========================================================================
2274 END MODULE module_sf_mynn