CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_sf_mynn.F90
1
3!WRF:MODEL_LAYER:PHYSICS
4!
8MODULE module_sf_mynn
9
10!-------------------------------------------------------------------
11!Modifications implemented by Joseph Olson NOAA/GSL
12!The following overviews the current state of this scheme::
13!
14! BOTH LAND AND WATER:
15!1) Calculation of stability parameter (z/L) taken from Li et al. (2010 BLM)
16! for first iteration of first time step; afterwards, exact calculation
17! using basically the same iterative technique in the module_sf_sfclayrev.F,
18! which leverages Pedro Jimenez's code, and is adapted for MYNN.
19!2) Fixed isflux=0 option to turn off scalar fluxes, but keep momentum
20! fluxes for idealized studies (credit: Anna Fitch).
21!3) Kinematic viscosity varies with temperature according to Andreas (1989).
22!4) Uses the blended Monin-Obukhov flux-profile relationships COARE (Fairall
23! et al 2003) for the unstable regime (a blended mix of Dyer-Hicks 1974 and
24! Grachev et al (2000). Uses Cheng and Brutsaert (2005) for stable conditions.
25!5) The following overviews the namelist variables that control the
26! aerodynamic roughness lengths (over water) and the thermal and moisture
27! roughness lengths (defaults are recommended):
28!
29! LAND only:
30! "iz0tlnd" namelist option is used to select the following momentum options:
31! (default) =0: Zilitinkevich (1995); Czil now set to 0.095
32! =1: Czil_new (modified according to Chen & Zhang 2008)
33! =2: Modified Yang et al (2002, 2008) - generalized for all landuse
34! =3: constant zt = z0/7.4 (original form; Garratt 1992)
35! =4: GFS - taken from sfc_diff.f, for comparison/testing
36!
37! WATER only:
38! "isftcflx" namelist option is used to select the following scalar options:
39! (default) =0: z0, zt, and zq from the COARE algorithm. Set COARE_OPT (below) to
40! 3.0 (Fairall et al. 2003, default)
41! 3.5 (Edson et al 2013)
42! =1: z0 from Davis et al (2008), zt & zq from COARE 3.0/3.5
43! =2: z0 from Davis et al (2008), zt & zq from Garratt (1992)
44! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5
45! =4: GFS - taken from sfc_diff.f, for comparison/testing
46!
47! SNOW/ICE only:
48! Andreas (2002) snow/ice parameterization for thermal and
49! moisture roughness is used over all gridpoints with snow deeper than
50! 0.1 m. This algorithm calculates a z0 for snow (Andreas et al. 2005, BLM),
51! which is only used as part of the thermal and moisture roughness
52! length calculation, not to directly impact the surface winds.
53!
54! Misc:
55!1) Added a more elaborate diagnostic for u10 & V10 for high vertical resolution
56! model configurations but for most model configurations with depth of
57! the lowest half-model level near 10 m, a neutral-log diagnostic is used.
58!
59!2) Option to activate stochastic parameter perturbations (SPP), which
60! perturb z0, zt, and zq, along with many other parameters in the MYNN-
61! EDMF scheme.
62!
63!NOTE: This code was primarily tested in combination with the RUC LSM.
64! Performance with the Noah (or other) LSM is relatively unknown.
65!-------------------------------------------------------------------
66!Include host model constants
67 use physcons, only : cp => con_cp, & !=7*Rd/2
68 & grav => con_g, & !=9.81
69 & rd => con_rd, & !=287.
70 & rv => con_rv, & !=461.6
71! & cpv => con_cvap, & !=4*Rv
72 & rovcp => con_rocp, & !=Rd/cp
73 & xlv => con_hvap, & !2.5e6
74 & xlf => con_hfus, & !3.5e5
75 & ep1 => con_fvirt, & !Rv/Rd - 1
76 & ep2 => con_eps !Rd/Rv
77
78!use kind_phys for real-types
79 use machine , only : kind_phys
80
81!-------------------------------------------------------------------
82 IMPLICIT NONE
83!-------------------------------------------------------------------
84!Drive and/or define more constant:
85 real(kind_phys), parameter :: ep3 = 1.-ep2
86 real(kind_phys), parameter :: g_inv = 1.0/grav
87 real(kind_phys), parameter :: rvovrd = rv/rd
88 real(kind_phys), parameter :: wmin = 0.1 ! Minimum wind speed
89 real(kind_phys), parameter :: karman = 0.4
90 real(kind_phys), parameter :: svp1 = 0.6112
91 real(kind_phys), parameter :: svp2 = 17.67
92 real(kind_phys), parameter :: svp3 = 29.65
93 real(kind_phys), parameter :: svpt0 = 273.15
94 real(kind_phys), parameter :: vconvc = 1.25
95 real(kind_phys), parameter :: onethird = 1./3.
96 real(kind_phys), parameter :: sqrt3 = 1.7320508075688773
97 real(kind_phys), parameter :: atan1 = 0.785398163397 !in radians
98 real(kind_phys), parameter :: log01 = log(0.01)
99 real(kind_phys), parameter :: log05 = log(0.05)
100 real(kind_phys), parameter :: log07 = log(0.07)
101 real(kind_phys), parameter :: snowz0 = 0.011
102 real(kind_phys), parameter :: coare_opt = 3.0 ! 3.0 or 3.5
103
104 !For debugging purposes:
105 INTEGER, PARAMETER :: debug_code = 0 !0: no extra ouput
106 !1: check input
107 !2: everything - heavy I/O
108
109 REAL(kind_phys), DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab, &
110 psih_stab,psih_unstab
111!$acc declare create(psim_stab, psim_unstab, psih_stab, psih_unstab)
112
113CONTAINS
114
115!-------------------------------------------------------------------
116!-------------------------------------------------------------------
119 SUBROUTINE sfclay_mynn( &
120 U3D,V3D,T3D,QV3D,P3D,dz8w, & !in
121 th3d,pi3d,qc3d, & !in
122 PSFCPA,PBLH,MAVAIL,XLAND,DX, & !in
123 ISFFLX,isftcflx,lsm,lsm_ruc, & !in
124 compute_flux,compute_diag, & !in
125 iz0tlnd,psi_opt, & !in
126 sigmaf,vegtype,shdmax,ivegsrc, & !intent(in)
127 z0pert,ztpert, & !intent(in)
128 redrag,sfc_z0_type, & !intent(in)
129 itimestep,iter,flag_iter, & !in
130 flag_restart, & !in
131 wet, dry, icy, & !intent(in)
132 tskin_wat, tskin_lnd, tskin_ice, & !intent(in)
133 tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in)
134 qsfc_wat, qsfc_lnd, qsfc_ice, & !intent(in)
135 snowh_wat, snowh_lnd, snowh_ice, & !intent(in)
136 znt_wat, znt_lnd, znt_ice, & !intent(inout)
137 ust_wat, ust_lnd, ust_ice, & !intent(inout)
138 cm_wat, cm_lnd, cm_ice, & !intent(inout)
139 ch_wat, ch_lnd, ch_ice, & !intent(inout)
140 rb_wat, rb_lnd, rb_ice, & !intent(inout)
141 stress_wat,stress_lnd,stress_ice, & !intent(inout)
142 fm_wat, fm_lnd, fm_ice, & !intent(inout)
143 fh_wat, fh_lnd, fh_ice, & !intent(inout)
144 fm10_wat, fm10_lnd, fm10_ice, & !intent(inout)
145 fh2_wat, fh2_lnd, fh2_ice, & !intent(inout)
146 hflx_wat, hflx_lnd, hflx_ice, &
147 qflx_wat, qflx_lnd, qflx_ice, &
148 ch,chs,chs2,cqs2,cpm, &
149 znt,ustm,zol,mol,rmol, &
150 psim,psih, &
151 hflx,hfx,qflx,qfx,lh,flhc,flqc, &
152 qgh,qsfc, &
153 u10,v10,th2,t2,q2, &
154 gz1oz0,wspd,wstar, &
155 spp_sfc,pattern_spp_sfc, &
156 ids,ide, jds,jde, kds,kde, &
157 ims,ime, jms,jme, kms,kme, &
158 its,ite, jts,jte, kts,kte, &
159 errmsg, errflg )
160!-------------------------------------------------------------------
161 IMPLICIT NONE
162!-------------------------------------------------------------------
163!-- U3D 3D u-velocity interpolated to theta points (m/s)
164!-- V3D 3D v-velocity interpolated to theta points (m/s)
165!-- T3D 3D temperature (K)
166!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
167!-- P3D 3D pressure (Pa)
168!-- dz8w 3D dz between full levels (m)
169!-- CP heat capacity at constant pressure for dry air (J/kg/K)
170!-- grav acceleration due to gravity (m/s^2)
171!-- ROVCP R/CP
172!-- Rd gas constant for dry air (J/kg/K)
173!-- XLV latent heat of vaporization for water (J/kg)
174!-- PSFCPA surface pressure (Pa)
175!-- ZNT roughness length (m)
176!-- UST u* in similarity theory (m/s)
177!-- USTM u* in similarity theory (m/s) w* added to WSPD. This is
178! used to couple with TKE scheme but not in MYNN.
179! (as of now, USTM = UST in this version)
180!-- PBLH PBL height from previous time (m)
181!-- MAVAIL surface moisture availability (between 0 and 1)
182!-- ZOL z/L height over Monin-Obukhov length
183!-- MOL T* (similarity theory) (K)
184!-- RMOL Reciprocal of M-O length (/m)
185!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
186!-- PSIM similarity stability function for momentum
187!-- PSIH similarity stability function for heat
188!-- XLAND land mask (1 for land, 2 for water)
189!-- HFX upward heat flux at the surface (W/m^2)
190! HFX = HFLX * rho * cp
191!-- HFLX upward temperature flux at the surface (K m s^-1)
192!-- QFX upward moisture flux at the surface (kg/m^2/s)
193! QFX = QFLX * rho
194!-- QFLX upward moisture flux at the surface (kg kg-1 m s-1)
195!-- LH net upward latent heat flux at surface (W/m^2)
196!-- TSK surface temperature (K)
197!-- FLHC exchange coefficient for heat (W/m^2/K)
198!-- FLQC exchange coefficient for moisture (kg/m^2/s)
199!-- CHS heat/moisture exchange coefficient for LSM (m/s)
200!-- QGH lowest-level saturated mixing ratio
201!-- QSFC qv (specific humidity) at the surface
202!-- QSFCMR qv (mixing ratio) at the surface
203!-- U10 diagnostic 10m u wind
204!-- V10 diagnostic 10m v wind
205!-- TH2 diagnostic 2m theta (K)
206!-- T2 diagnostic 2m temperature (K)
207!-- Q2 diagnostic 2m mixing ratio (kg/kg)
208!-- SNOWH Snow height (m)
209!-- GZ1OZ0 log((z1+ZNT)/ZNT) where ZNT is roughness length
210!-- WSPD wind speed at lowest model level (m/s)
211!-- BR bulk Richardson number in surface layer
212!-- ISFFLX isfflx=1 for surface heat and moisture fluxes
213!-- DX horizontal grid size (m)
214!-- SVP1 constant for saturation vapor pressure (=0.6112 kPa)
215!-- SVP2 constant for saturation vapor pressure (=17.67 dimensionless)
216!-- SVP3 constant for saturation vapor pressure (=29.65 K)
217!-- SVPT0 constant for saturation vapor pressure (=273.15 K)
218!-- EP1 constant for virtual temperature (Rv/Rd - 1) (dimensionless)
219!-- EP2 constant for spec. hum. calc (Rd/Rv = 0.622) (dimensionless)
220!-- EP3 constant for spec. hum. calc (1 - Rd/Rv = 0.378 ) (dimensionless)
221!-- KARMAN Von Karman constant
222!-- ck enthalpy exchange coeff at 10 meters
223!-- cd momentum exchange coeff at 10 meters
224!-- cka enthalpy exchange coeff at the lowest model level
225!-- cda momentum exchange coeff at the lowest model level
226!-- isftcflx =0: z0, zt, and zq from COARE3.0/3.5 (Fairall et al 2003/Edson et al 2013)
227! (water =1: z0 from Davis et al (2008), zt & zq from COARE3.0/3.5
228! only) =2: z0 from Davis et al (2008), zt & zq from Garratt (1992)
229! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5
230!-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.095,
231! (land =1: Czil_new (modified according to Chen & Zhang 2008)
232! only) =2: Modified Yang et al (2002, 2008) - generalized for all landuse
233! =3: constant zt = z0/7.4 (Garratt 1992)
234!
235!-- ids start index for i in domain
236!-- ide end index for i in domain
237!-- jds start index for j in domain
238!-- jde end index for j in domain
239!-- kds start index for k in domain
240!-- kde end index for k in domain
241!-- ims start index for i in memory
242!-- ime end index for i in memory
243!-- jms start index for j in memory
244!-- jme end index for j in memory
245!-- kms start index for k in memory
246!-- kme end index for k in memory
247!-- its start index for i in tile
248!-- ite end index for i in tile
249!-- jts start index for j in tile
250!-- jte end index for j in tile
251!-- kts start index for k in tile
252!-- kte end index for k in tile
253!-- errmsg CCPP error message
254!-- errflg CCPP error code
255!=================================================================
256! SCALARS
257!===================================
258 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
259 ims,ime, jms,jme, kms,kme, &
260 its,ite, jts,jte, kts,kte
261 INTEGER, INTENT(IN) :: itimestep,iter
262!NAMELIST/CONFIGURATION OPTIONS:
263 integer, intent(in) :: ISFFLX, LSM, LSM_RUC
264 INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND
265 INTEGER, OPTIONAL, INTENT(IN) :: spp_sfc, psi_opt
266 logical, intent(in) :: compute_flux,compute_diag
267 integer, intent(in) :: ivegsrc
268 integer, intent(in) :: sfc_z0_type ! option for calculating surface roughness length over ocean
269 logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han)
270 logical, intent(in) :: flag_restart
271
272!Input data
273 integer, dimension(ims:ime), intent(in) :: vegtype
274 real(kind_phys), dimension(ims:ime), intent(in) :: &
275 & sigmaf,shdmax,z0pert,ztpert
276!===================================
277! 3D VARIABLES
278!===================================
279 REAL(kind_phys), DIMENSION( ims:ime, kms:kme ) , &
280 INTENT(IN ) :: dz8w, &
281 QV3D, &
282 P3D, &
283 T3D, &
284 QC3D, &
285 U3D,V3D, &
286 th3d,pi3d
287
288 !GJF: This array must be assumed-shape since it is conditionally-allocated
289 REAL(kind_phys), DIMENSION( :,: ), OPTIONAL, &
290 INTENT(IN) :: pattern_spp_sfc
291!===================================
292! 2D VARIABLES
293!===================================
294 REAL(kind_phys), DIMENSION( ims:ime ) , &
295 INTENT(IN ) :: MAVAIL, &
296 PBLH, &
297 XLAND, &
298 PSFCPA, &
299 DX
300
301 REAL(kind_phys), DIMENSION( ims:ime ) , &
302 INTENT(OUT ) :: U10,V10, &
303 TH2,T2,Q2
304
305
306 REAL(kind_phys), DIMENSION( ims:ime ) , &
307 INTENT(INOUT) :: HFLX,HFX, &
308 QFLX,QFX, &
309 RMOL, &
310 QSFC, &
311 QGH, &
312 ZNT, &
313 CPM, &
314 CHS, &
315 CH, &
316 FLHC,FLQC, &
317 GZ1OZ0,WSPD, &
318 PSIM,PSIH
319 REAL(kind_phys), DIMENSION( ims:ime ), OPTIONAL , &
320 INTENT(INOUT) :: USTM, &
321 CHS2, &
322 CQS2, &
323 LH, &
324 ZOL, &
325 MOL, &
326 WSTAR
327
328 LOGICAL, DIMENSION( ims:ime ), INTENT(IN) :: &
329& wet, dry, icy, flag_iter
330
331 REAL(kind_phys), DIMENSION( ims:ime ), INTENT(IN) :: &
332 & tskin_wat, tskin_lnd, tskin_ice, &
333 & tsurf_wat, tsurf_lnd, tsurf_ice, &
334 & snowh_wat, snowh_lnd, snowh_ice
335
336 REAL(kind_phys), DIMENSION( ims:ime), INTENT(INOUT) :: &
337 & ZNT_wat, ZNT_lnd, ZNT_ice, &
338 & UST_wat, UST_lnd, UST_ice, &
339 & cm_wat, cm_lnd, cm_ice, &
340 & ch_wat, ch_lnd, ch_ice, &
341 & rb_wat, rb_lnd, rb_ice, &
342 & stress_wat,stress_lnd,stress_ice, &
343 & fm_wat, fm_lnd, fm_ice, &
344 & fh_wat, fh_lnd, fh_ice, &
345 & fm10_wat, fm10_lnd, fm10_ice, &
346 & fh2_wat, fh2_lnd, fh2_ice, &
347 & HFLX_wat, HFLX_lnd, HFLX_ice, &
348 & QFLX_wat, QFLX_lnd, QFLX_ice, &
349 & qsfc_wat, qsfc_lnd, qsfc_ice
350
351! CCPP error handling
352 character(len=*), intent(inout) :: errmsg
353 integer, intent(inout) :: errflg
354
355!ADDITIONAL OUTPUT
356!JOE-begin
357 REAL(kind_phys), DIMENSION( ims:ime ) :: qstar
358!JOE-end
359!===================================
360! 1D LOCAL ARRAYS
361!===================================
362 REAL(kind_phys), DIMENSION( its:ite ) :: U1D,V1D, & !level1 winds
363 U1D2,V1D2, & !level2 winds
364 QV1D, &
365 P1D, &
366 T1D,QC1D, &
367 dz8w1d, & !level 1 height
368 dz2w1d !level 2 height
369
370 REAL(kind_phys), DIMENSION( its:ite ) :: rstoch1D
371
372 INTEGER :: I,J,K,itf,ktf
373!-----------------------------------------------------------
374
375 ! Initialize error-handling
376 errflg = 0
377 errmsg = ''
378
379!$acc enter data copyin( dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D, &
380!$acc pattern_spp_sfc)
381
382!$acc enter data copyin( UST_WAT(:), UST_LND(:), UST_ICE(:), &
383!$acc MOL(:), QFLX(:), HFLX(:), &
384!$acc QSFC(:), QSFC_WAT(:), QSFC_LND(:), &
385!$acc QSFC_ICE(:))
386
387!$acc enter data create( dz8w1d(:), dz2w1d(:), U1D(:), &
388!$acc V1D(:), U1D2(:), V1D2(:), &
389!$acc QV1D(:), QC1D(:), P1D(:), &
390!$acc T1D(:), rstoch1D(:), qstar(:))
391
392
393 IF (debug_code >= 1) THEN
394 write(*,*)"======= printing of constants:"
395 write(*,*)"cp=", cp," g=", grav
396 write(*,*)"Rd=", rd," ep1=", ep1
397 write(*,*)"xlv=", xlv," xlf=", xlf
398 write(*,*)"ep2=", ep2
399 ENDIF
400
401 itf=ite !MIN0(ite,ide-1)
402 ktf=kte !MIN0(kte,kde-1)
403
404!$acc parallel loop present(dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D, &
405!$acc pattern_spp_sfc,dz8w1d,dz2w1d,U1D, &
406!$acc V1D,U1D2,V1D2,QV1D,QC1D,P1D,T1D, &
407!$acc rstoch1D,qstar)
408 DO i=its,ite
409 dz8w1d(i) = dz8w(i,kts)
410 dz2w1d(i) = dz8w(i,kts+1)
411 u1d(i) =u3d(i,kts)
412 v1d(i) =v3d(i,kts)
413 !2nd model level winds - for diags with high-res grids
414 u1d2(i) =u3d(i,kts+1)
415 v1d2(i) =v3d(i,kts+1)
416 qv1d(i)=qv3d(i,kts)
417 qc1d(i)=qc3d(i,kts)
418 p1d(i) =p3d(i,kts)
419 t1d(i) =t3d(i,kts)
420 if (spp_sfc==1) then
421 rstoch1d(i)=pattern_spp_sfc(i,kts)
422 else
423 rstoch1d(i)=0.0
424 endif
425 qstar(i)=0.0
426 ENDDO
427
428 IF (itimestep==1 .AND. iter==1) THEN
429!$acc parallel loop present(U1D,V1D,UST_WAT,UST_LND,UST_ICE,MOL, &
430!$acc QFLX,HFLX,QV3D,QSFC,QSFC_WAT, &
431!$acc QSFC_LND,QSFC_ICE)
432 DO i=its,ite
433 IF (.not. flag_restart) THEN
434 !Everything here is used before calculated
435 if (ust_wat(i) .lt. 1e-4 .or. ust_wat(i) .gt. 3.0) then
436 ust_wat(i)=max(0.04*sqrt(u1d(i)*u1d(i) + v1d(i)*v1d(i)),0.001_kind_phys)
437 endif
438 if (ust_lnd(i) .lt. 1e-4 .or. ust_lnd(i) .gt. 3.0) then
439 ust_lnd(i)=max(0.04*sqrt(u1d(i)*u1d(i) + v1d(i)*v1d(i)),0.001_kind_phys)
440 endif
441 if (ust_ice(i) .lt. 1e-4 .or. ust_ice(i) .gt. 3.0) then
442 ust_ice(i)=max(0.04*sqrt(u1d(i)*u1d(i) + v1d(i)*v1d(i)),0.001_kind_phys)
443 endif
444 mol(i)=0.0
445 ENDIF ! restart
446 qflx(i)=0.
447 hflx(i)=0.
448 if ( lsm == lsm_ruc ) then
449 !- qsfc_lnd and qsfc_ice are already available
450 qsfc(i)=qv3d(i,kts)/(1.+qv3d(i,kts))
451 qsfc_wat(i)=qsfc(i)
452 else
453 qsfc(i)=qv3d(i,kts)/(1.+qv3d(i,kts))
454 qsfc_wat(i)=qsfc(i)
455 qsfc_lnd(i)=qsfc(i)
456 qsfc_ice(i)=qsfc(i)
457 endif ! lsm==lsm_ruc
458 ENDDO
459 ENDIF
460
461!$acc exit data delete( dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D, &
462!$acc pattern_spp_sfc, QC1D)
463
464 CALL sfclay1d_mynn(flag_iter, &
465 j,u1d,v1d,t1d,qv1d,p1d,dz8w1d, &
466 u1d2,v1d2,dz2w1d, &
467 psfcpa,pblh,mavail,xland,dx, &
468 isfflx,isftcflx,iz0tlnd,psi_opt, &
469 compute_flux,compute_diag, &
470 sigmaf,vegtype,shdmax,ivegsrc, & !intent(in)
471 z0pert,ztpert, & !intent(in)
472 redrag,sfc_z0_type, & !intent(in)
473 itimestep,iter,flag_restart,lsm,lsm_ruc, &
474 wet, dry, icy, & !intent(in)
475 tskin_wat, tskin_lnd, tskin_ice, & !intent(in)
476 tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in)
477 qsfc_wat, qsfc_lnd, qsfc_ice, & !intent(in)
478 snowh_wat, snowh_lnd, snowh_ice, & !intent(in)
479 znt_wat, znt_lnd, znt_ice, & !intent(inout)
480 ust_wat, ust_lnd, ust_ice, & !intent(inout)
481 cm_wat, cm_lnd, cm_ice, & !intent(inout)
482 ch_wat, ch_lnd, ch_ice, & !intent(inout)
483 rb_wat, rb_lnd, rb_ice, & !intent(inout)
484 stress_wat, stress_lnd, stress_ice, & !intent(inout)
485 fm_wat, fm_lnd, fm_ice, & !intent(inout)
486 fh_wat, fh_lnd, fh_ice, & !intent(inout)
487 fm10_wat, fm10_lnd, fm10_ice, & !intent(inout)
488 fh2_wat, fh2_lnd, fh2_ice, &
489 hflx_wat, hflx_lnd, hflx_ice, &
490 qflx_wat, qflx_lnd, qflx_ice, &
491 ch,chs,chs2,cqs2,cpm, &
492 znt,ustm,zol,mol,rmol, &
493 psim,psih, &
494 hflx,hfx,qflx,qfx,lh,flhc,flqc, &
495 qgh,qsfc,u10,v10,th2,t2,q2, &
496 gz1oz0,wspd,wstar,qstar, &
497 spp_sfc,rstoch1d, &
498 ids,ide, jds,jde, kds,kde, &
499 ims,ime, jms,jme, kms,kme, &
500 its,ite, jts,jte, kts,kte, &
501 errmsg, errflg )
502
503!$acc exit data copyout( UST_WAT(:), UST_LND(:), UST_ICE(:), &
504!$acc MOL(:), QFLX(:), HFLX(:), &
505!$acc QSFC(:), QSFC_WAT(:), QSFC_LND(:), &
506!$acc QSFC_ICE(:))
507
508!$acc exit data delete( dz8w1d(:), dz2w1d(:), U1D(:), &
509!$acc V1D(:), U1D2(:), V1D2(:), &
510!$acc QV1D(:), T1D(:), P1D(:), &
511!$acc rstoch1D(:), qstar(:))
512
513 END SUBROUTINE sfclay_mynn
514
515!-------------------------------------------------------------------
521 SUBROUTINE sfclay1d_mynn(flag_iter, &
522 J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,U1D2,V1D2,dz2w1d, &
523 PSFCPA,PBLH,MAVAIL,XLAND,DX, &
524 ISFFLX,isftcflx,iz0tlnd,psi_opt, &
525 compute_flux,compute_diag, &
526 sigmaf,vegtype,shdmax,ivegsrc, & !intent(in)
527 z0pert,ztpert, & !intent(in)
528 redrag,sfc_z0_type, & !intent(in)
529 itimestep,iter,flag_restart,lsm,lsm_ruc, &
530 wet, dry, icy, & !intent(in)
531 tskin_wat, tskin_lnd, tskin_ice, & !intent(in)
532 tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in)
533 qsfc_wat, qsfc_lnd, qsfc_ice, & !intent(in)
534 snowh_wat, snowh_lnd, snowh_ice, & !intent(in)
535 znt_wat, znt_lnd, znt_ice, & !intent(inout)
536 ust_wat, ust_lnd, ust_ice, & !intent(inout)
537 cm_wat, cm_lnd, cm_ice, & !intent(inout)
538 ch_wat, ch_lnd, ch_ice, & !intent(inout)
539 rb_wat, rb_lnd, rb_ice, & !intent(inout)
540 stress_wat, stress_lnd, stress_ice, & !intent(inout)
541 psix_wat, psix_lnd, psix_ice, & !=fm, intent(inout)
542 psit_wat, psit_lnd, psit_ice, & !=fh, intent(inout)
543 psix10_wat, psix10_lnd, psix10_ice, & !=fm10, intent(inout)
544 psit2_wat, psit2_lnd, psit2_ice, & !=fh2, intent(inout)
545 hflx_wat, hflx_lnd, hflx_ice, &
546 qflx_wat, qflx_lnd, qflx_ice, &
547 ch,chs,chs2,cqs2,cpm, &
548 znt,ustm,zol,mol,rmol, &
549 psim,psih, &
550 hflx,hfx,qflx,qfx,lh,flhc,flqc, &
551 qgh,qsfc, &
552 u10,v10,th2,t2,q2, &
553 gz1oz0,wspd,wstar,qstar, &
554 spp_sfc,rstoch1d, &
555 ids,ide, jds,jde, kds,kde, &
556 ims,ime, jms,jme, kms,kme, &
557 its,ite, jts,jte, kts,kte, &
558 errmsg, errflg )
559
560!-------------------------------------------------------------------
561 IMPLICIT NONE
562!-------------------------------------------------------------------
563! SCALARS
564!-----------------------------
565 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
566 ims,ime, jms,jme, kms,kme, &
567 its,ite, jts,jte, kts,kte, &
568 J, itimestep, iter, lsm, lsm_ruc
569 LOGICAL, INTENT(IN) :: flag_restart
570
571 REAL(kind_phys), PARAMETER :: XKA=2.4e-5 !molecular diffusivity
572 REAL(kind_phys), PARAMETER :: PRT=1. !prandlt number
573 REAL(kind_phys), PARAMETER :: snowh_thresh = 50. !mm
574
575!-----------------------------
576! NAMELIST OPTIONS
577!-----------------------------
578 integer, intent(in) :: ISFFLX
579 integer, optional, intent(in) :: ISFTCFLX, IZ0TLND
580 logical, intent(in) :: compute_flux,compute_diag
581 integer, intent(in) :: spp_sfc, psi_opt
582 integer, intent(in) :: ivegsrc
583 integer, intent(in) :: sfc_z0_type ! option for calculating surface roughness length over ocean
584 logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han)
585
586!Input data
587 integer, dimension(ims:ime), intent(in) :: vegtype
588 real(kind_phys), dimension(ims:ime), intent(in) :: &
589 & sigmaf,shdmax,z0pert,ztpert
590
591!-----------------------------
592! 1D ARRAYS
593!-----------------------------
594 REAL(kind_phys), DIMENSION( ims:ime ), &
595 INTENT(IN) :: MAVAIL, &
596 PBLH, &
597 XLAND, &
598 PSFCPA, &
599 DX
600
601 REAL(kind_phys), DIMENSION( its:ite ), &
602 INTENT(IN) :: U1D,V1D, &
603 U1D2,V1D2, &
604 QV1D,P1D, &
605 T1D, &
606 dz8w1d, &
607 dz2w1d
608
609 REAL(kind_phys), DIMENSION( ims:ime ), &
610 INTENT(OUT) :: QFX,HFX, &
611 RMOL
612 REAL(kind_phys), DIMENSION( ims:ime ), &
613 INTENT(INOUT) :: HFLX,QFLX, &
614 QGH,QSFC, &
615 ZNT, &
616 CPM, &
617 CHS,CH, &
618 FLHC,FLQC, &
619 GZ1OZ0, &
620 WSPD, &
621 PSIM, &
622 PSIH
623 REAL(kind_phys), DIMENSION( ims:ime ), OPTIONAL, &
624 INTENT(INOUT) :: MOL, &
625 ZOL, &
626 LH, &
627 CHS2,CQS2, &
628 USTM
629
630 LOGICAL, DIMENSION( ims:ime ), INTENT(IN) :: &
631 & wet, dry, icy, flag_iter
632
633 REAL(kind_phys), DIMENSION( ims:ime ), INTENT(in) :: &
634 & tskin_wat, tskin_lnd, tskin_ice, &
635 & tsurf_wat, tsurf_lnd, tsurf_ice, &
636 & snowh_wat, snowh_lnd, snowh_ice
637
638 REAL(kind_phys), DIMENSION( ims:ime ), INTENT(inout) :: &
639 & ZNT_wat, ZNT_lnd, ZNT_ice, &
640 & UST_wat, UST_lnd, UST_ice, &
641 & cm_wat, cm_lnd, cm_ice, &
642 & ch_wat, ch_lnd, ch_ice, &
643 & rb_wat, rb_lnd, rb_ice, &
644 & stress_wat,stress_lnd,stress_ice, &
645 & psix_wat, psix_lnd, psix_ice, &
646 & psit_wat, psit_lnd, psit_ice, &
647 & psix10_wat,psix10_lnd,psix10_ice, &
648 & psit2_wat, psit2_lnd, psit2_ice, &
649 & HFLX_wat, HFLX_lnd, HFLX_ice, &
650 & QFLX_wat, QFLX_lnd, QFLX_ice, &
651 & qsfc_wat, qsfc_lnd, qsfc_ice
652
653 REAL(kind_phys), DIMENSION( its:ite ), &
654 & INTENT(IN) :: rstoch1D
655
656 ! DIAGNOSTIC OUTPUT
657 REAL(kind_phys), DIMENSION( ims:ime ), &
658 & INTENT(OUT) :: U10, V10, &
659 & TH2, T2, &
660 & Q2
661
662!--------------------------------------------
663!JOE-additinal output
664 REAL(kind_phys), DIMENSION( ims:ime ), &
665 & INTENT(OUT) :: qstar
666 REAL(kind_phys), DIMENSION( ims:ime ), OPTIONAL, &
667 & INTENT(OUT) :: wstar
668
669!JOE-end
670
671! CCPP error handling
672 character(len=*), intent(inout) :: errmsg
673 integer, intent(inout) :: errflg
674
675! Local fixed-size errmsg character array for error messages on accelerator
676! devices distinct from the host (e.g. GPUs). Necessary since OpenACC does
677! not support assumed-size (len=*) arrays like errmsg. Additional
678! device_errflg integer to denote when device_errmsg needs to be synced
679! with errmsg.
680 character(len=512) :: device_errmsg
681 integer :: device_errflg
682
683! Special versions of the fixed-size errmsg character array for error messages
684! on the device and it's errflag counterpart. These are necessary to ensure
685! the return statements at lines 1417 and 2030 are executed only for this
686! special case, and not any and all error messages set on the device.
687 character(len=512) :: device_special_errmsg
688 integer :: device_special_errflg
689
690
691!----------------------------------------------------------------
692! LOCAL VARS
693!----------------------------------------------------------------
694 REAL(kind_phys), DIMENSION(its:ite) :: &
695 ZA, & !Height of lowest 1/2 sigma level(m)
696 ZA2, & !Height of 2nd lowest 1/2 sigma level(m)
697 THV1D, & !Theta-v at lowest 1/2 sigma (K)
698 TH1D, & !Theta at lowest 1/2 sigma (K)
699 TC1D, & !T at lowest 1/2 sigma (Celsius)
700 TV1D, & !Tv at lowest 1/2 sigma (K)
701 RHO1D, & !density at lowest 1/2 sigma level
702 QVSH, & !qv at lowest 1/2 sigma (spec humidity)
703 PSIH2, & !M-O stability functions at z=2 m
704 PSIM10, & !M-O stability functions at z=10 m
705 PSIH10, & !M-O stability functions at z=10 m
706 WSPDI, &
707 GOVRTH, & !grav/theta
708 PSFC, & !press at surface (Pa/1000)
709 QSFCMR, & !qv at surface (mixing ratio, kg/kg)
710 THCON, & !conversion from temp to theta
711 zratio_lnd, zratio_ice, zratio_wat, & !z0/zt
712 TSK_lnd, TSK_ice, TSK_wat, & !absolute temperature
713 THSK_lnd, THSK_ice, THSK_wat, & !theta
714 THVSK_lnd, THVSK_ice, THVSK_wat, & !theta-v
715 GZ1OZ0_lnd, GZ1OZ0_ice, GZ1OZ0_wat, & !LOG((ZA(I)+ZNT(i))/ZNT(i))
716 GZ1OZt_lnd, GZ1OZt_ice, GZ1OZt_wat, & !LOG((ZA(I)+ZT(i))/ZT(i))
717 GZ2OZ0_lnd, GZ2OZ0_ice, GZ2OZ0_wat, & !LOG((2.0+ZNT(I))/ZNT(I))
718 GZ2OZt_lnd, GZ2OZt_ice, GZ2OZt_wat, & !LOG((2.0+ZT(I))/ZT(I))
719 GZ10OZ0_lnd, GZ10OZ0_ice, GZ10OZ0_wat, & !LOG((10.+ZNT(I))/ZNT(I))
720 GZ10OZt_lnd, GZ10OZt_ice, GZ10OZt_wat, & !LOG((10.+ZT(I))/ZT(I))
721 ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat, &
722 ZT_lnd, ZT_ice, ZT_wat, &
723 ZQ_lnd, ZQ_ice, ZQ_wat, &
724 PSIQ_lnd, PSIQ_ice, PSIQ_wat, &
725 PSIQ2_lnd, PSIQ2_ice, PSIQ2_wat, &
726 QSFCMR_lnd, QSFCMR_ice, QSFCMR_wat
727
728 INTEGER :: N,I,K,L,yesno
729
730 REAL(kind_phys) :: PL,E1,TABS
731 REAL(kind_phys) :: WSPD_lnd, WSPD_ice, WSPD_wat
732 REAL(kind_phys) :: DTHVDZ,DTHVM,VCONV,ZOL2,ZOL10,ZOLZA,ZOLZ0,ZOLZT
733 REAL(kind_phys) :: DTG,DTTHX,PSIQ,PSIQ2,PSIQ10,PSIT10
734 REAL(kind_phys) :: FLUXC,VSGD
735 REAL(kind_phys) :: restar,VISC,DQG,OLDUST,OLDTST
736
737 ! Initialize error-handling
738 errflg = 0
739 errmsg = ''
740 device_errflg = errflg
741 device_errmsg = errmsg
742 device_special_errflg = errflg
743 device_special_errmsg = errmsg
744!-------------------------------------------------------------------
745!$acc update device(psim_stab, psim_unstab, psih_stab, psih_unstab)
746
747!$acc enter data create( ZA, ZA2, THV1D, TH1D, TC1D, TV1D, &
748!$acc RHO1D, QVSH, PSIH2, PSIM10, PSIH10, WSPDI, &
749!$acc GOVRTH, PSFC, THCON, &
750!$acc zratio_lnd, zratio_ice, zratio_wat, &
751!$acc TSK_lnd, TSK_ice, TSK_wat, &
752!$acc THSK_lnd, THSK_ice, THSK_wat, &
753!$acc THVSK_lnd, THVSK_ice, THVSK_wat, &
754!$acc GZ1OZ0_lnd, GZ1OZ0_ice, GZ1OZ0_wat, &
755!$acc GZ1OZt_lnd, GZ1OZt_ice, GZ1OZt_wat, &
756!$acc GZ2OZ0_lnd, GZ2OZ0_ice, GZ2OZ0_wat, &
757!$acc GZ2OZt_lnd, GZ2OZt_ice, GZ2OZt_wat, &
758!$acc GZ10OZ0_lnd, GZ10OZ0_ice, GZ10OZ0_wat, &
759!$acc GZ10OZt_lnd, GZ10OZt_ice, GZ10OZt_wat, &
760!$acc ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat, &
761!$acc ZT_lnd, ZT_ice, ZT_wat, &
762!$acc ZQ_lnd, ZQ_ice, ZQ_wat, &
763!$acc PSIQ_lnd, PSIQ_ice, PSIQ_wat, &
764!$acc PSIQ2_lnd, PSIQ2_ice, PSIQ2_wat, &
765!$acc QSFCMR_lnd, QSFCMR_ice, QSFCMR_wat )
766
767!$acc enter data copyin(flag_iter, dry, wet, icy, CPM, MAVAIL, &
768!$acc QFX, FLHC, FLQC, CHS, CH, CHS2, CQS2, USTM, &
769!$acc HFX, LH, wstar, qstar, PBLH, ZOL, MOL, RMOL, &
770!$acc T2, TH2, Q2, QV1D, PSFCPA, &
771!$acc WSPD, U10, V10, U1D, V1D, U1D2, V1D2, &
772!$acc T1D, P1D, rstoch1D, sigmaf, &
773!$acc shdmax, vegtype, z0pert, ztpert, dx, QGH, &
774!$acc dz2w1d, dz8w1d, &
775!$acc stress_wat, stress_lnd, stress_ice, &
776!$acc rb_wat, rb_lnd, rb_ice, &
777!$acc tskin_wat, tskin_lnd, tskin_ice, &
778!$acc tsurf_wat, tsurf_lnd, tsurf_ice, &
779!$acc psim, psih, &
780!$acc UST_wat, UST_lnd, UST_ice, &
781!$acc ZNT_wat, ZNT_lnd, ZNT_ice, &
782!$acc QSFC, QSFC_lnd, QSFC_wat, QSFC_ice, &
783!$acc QFLX, QFLX_lnd, QFLX_wat, QFLX_ice, &
784!$acc HFLX, HFLX_lnd, HFLX_wat, HFLX_ice, &
785!$acc PSIX_wat, PSIX_lnd, PSIX_ice, &
786!$acc PSIX10_wat, PSIX10_lnd, PSIX10_ice, &
787!$acc PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
788!$acc PSIT_lnd, PSIT_wat, PSIT_ice, &
789!$acc ch_lnd, ch_wat, ch_ice, &
790!$acc cm_lnd, cm_wat, cm_ice, &
791!$acc snowh_lnd, snowh_wat, snowh_ice, &
792!$acc device_errmsg, device_errflg, &
793!$acc device_special_errmsg, device_special_errflg)
794
795!$acc parallel loop present(PSFCPA, PSFC, QSFC, T1D, flag_iter, tsurf_lnd, &
796!$acc QSFC_wat, QSFCMR_wat, wet, TSK_wat, tskin_wat, &
797!$acc QSFC_lnd, QSFCMR_lnd, dry, TSK_lnd, tskin_lnd, &
798!$acc QSFC_ice, QSFCMR_ice, icy, TSK_ice, tskin_ice)
799 DO i=its,ite
800
801 ! PSFC ( in cmb) is used later in saturation checks
802 psfc(i)=psfcpa(i)/1000.
803 !tgs - do computations if flag_iter(i) = .true.
804 if ( flag_iter(i) ) then
805
806 IF (itimestep == 1) THEN
807 !initialize surface specific humidity and mixing ratios for land, ice and water
808 IF (wet(i)) THEN
809 tsk_wat(i) = tskin_wat(i)
810 IF (tsk_wat(i) .LT. 273.15) THEN
811 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
812 e1=svp1*exp(4648*(1./273.15 - 1./tsk_wat(i)) - &
813 & 11.64*log(273.15/tsk_wat(i)) + 0.02265*(273.15 - tsk_wat(i)))
814 ELSE
815 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
816 e1=svp1*exp(svp2*(tsk_wat(i)-svpt0)/(tsk_wat(i)-svp3))
817 ENDIF
818 qsfc_wat(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity
819 qsfcmr_wat(i)=ep2*e1/(psfc(i)-e1) !mixing ratio
820 IF(qsfc_wat(i)>1..or.qsfc_wat(i)<0.) print *,' QSFC_wat(I)',itimestep,i,qsfc_wat(i),tsk_wat(i)
821 ENDIF
822 IF (dry(i)) THEN
823 tsk_lnd(i) = tskin_lnd(i)
824 if( lsm == lsm_ruc) then
825 qsfcmr_lnd(i)=qsfc_lnd(i)/(1.-qsfc_lnd(i)) !mixing ratio
826 else
827 tabs = 0.5*(tsk_lnd(i) + t1d(i))
828 IF (tabs .LT. 273.15) THEN
829 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
830 e1=svp1*exp(4648*(1./273.15 - 1./tabs) - &
831 & 11.64*log(273.15/tabs) + 0.02265*(273.15 - tabs))
832 ELSE
833 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
834 e1=svp1*exp(svp2*(tabs-svpt0)/(tabs-svp3))
835 ENDIF
836 qsfc_lnd(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity
837 qsfc_lnd(i)=0.5*(qsfc_lnd(i) + qsfc(i))
838 qsfcmr_lnd(i)=qsfc_lnd(i)/(1.-qsfc_lnd(i)) !mixing ratio
839 endif ! lsm
840 IF(qsfc_lnd(i)>1..or.qsfc_lnd(i)<0.) print *,' QSFC_lnd(I)',itimestep,i,qsfc_lnd(i),tskin_lnd(i),tsurf_lnd(i),qsfc(i)
841 ENDIF
842 IF (icy(i)) THEN
843 tsk_ice(i) = tskin_ice(i)
844 if( lsm == lsm_ruc) then
845 qsfcmr_ice(i)=qsfc_ice(i)/(1.-qsfc_ice(i)) !mixing ratio
846 else
847 IF (tsk_ice(i) .LT. 273.15) THEN
848 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
849 e1=svp1*exp(4648*(1./273.15 - 1./tsk_ice(i)) - &
850 & 11.64*log(273.15/tsk_ice(i)) + 0.02265*(273.15 - tsk_ice(i)))
851 ELSE
852 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
853 e1=svp1*exp(svp2*(tsk_ice(i)-svpt0)/(tsk_ice(i)-svp3))
854 ENDIF
855 qsfc_ice(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity
856 qsfcmr_ice(i)=ep2*e1/(psfc(i)-e1) !mixing ratio
857 endif ! lsm
858 IF(qsfc_ice(i)>1..or.qsfc_ice(i)<0.) print *,' QSFC_ice(I)',itimestep,i,qsfc_ice(i),tsk_ice(i)
859 ENDIF
860
861 ELSE
862
863 ! Use what comes out of the NST, LSM, SICE after check
864 IF (wet(i)) then
865 tsk_wat(i) = tskin_wat(i)
866 IF (tsk_wat(i) .LT. 273.15) THEN
867 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
868 e1=svp1*exp(4648*(1./273.15 - 1./tsk_wat(i)) - &
869 & 11.64*log(273.15/tsk_wat(i)) + 0.02265*(273.15 - tsk_wat(i)))
870 ELSE
871 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
872 e1=svp1*exp(svp2*(tsk_wat(i)-svpt0)/(tsk_wat(i)-svp3))
873 ENDIF
874 qsfc_wat(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity
875 ENDIF
876 IF (dry(i).and.(qsfc_lnd(i)>1..or.qsfc_lnd(i)<0.)) then
877 !print *,'bad QSFC_lnd(I)',itimestep,iter,i,QSFC_lnd(I),TSKin_lnd(I)
878 tabs = 0.5*(tskin_lnd(i) + t1d(i))
879 IF (tabs .LT. 273.15) THEN
880 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
881 e1=svp1*exp(4648*(1./273.15 - 1./tabs) - &
882 & 11.64*log(273.15/tabs) + 0.02265*(273.15 - tabs))
883 ELSE
884 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
885 e1=svp1*exp(svp2*(tabs-svpt0)/(tabs-svp3))
886 ENDIF
887 qsfc_lnd(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity
888 qsfc_lnd(i)=0.5*(qsfc_lnd(i) + qsfc(i))
889 ENDIF
890 IF (icy(i).and.(qsfc_ice(i)>1..or.qsfc_ice(i)<0.)) then
891 !print *,'bad QSFC_ice(I)',itimestep,iter,i,QSFC_ice(I),TSKin_ice(I)
892 IF (tskin_ice(i) .LT. 273.15) THEN
893 !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb)
894 e1=svp1*exp(4648*(1./273.15 - 1./tskin_ice(i)) - &
895 & 11.64*log(273.15/tskin_ice(i)) + 0.02265*(273.15 - tskin_ice(i)))
896 ELSE
897 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
898 e1=svp1*exp(svp2*(tskin_ice(i)-svpt0)/(tskin_ice(i)-svp3))
899 ENDIF
900 qsfc_ice(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity
901 ENDIF
902
903 IF (wet(i)) qsfcmr_wat(i)=qsfc_wat(i)/(1.-qsfc_wat(i))
904 IF (dry(i)) qsfcmr_lnd(i)=qsfc_lnd(i)/(1.-qsfc_lnd(i))
905 IF (icy(i)) qsfcmr_ice(i)=qsfc_ice(i)/(1.-qsfc_ice(i))
906
907 ENDIF
908 endif ! flag_iter
909 ENDDO
910
911!$acc serial present(pblh, PSFCPA, dz8w1d, qflx, hflx, &
912!$acc dry, tskin_lnd, tsurf_lnd, qsfc_lnd, znt_lnd, ust_lnd, snowh_lnd, &
913!$acc icy, tskin_ice, tsurf_ice, qsfc_ice, znt_ice, ust_ice, snowh_ice, &
914!$acc wet, tskin_wat, tsurf_wat, qsfc_wat, znt_wat, ust_wat, snowh_wat)
915 IF (debug_code >= 1) THEN
916 write(0,*)"ITIMESTEP=",itimestep," iter=",iter
917 DO i=its,ite
918 write(0,*)"=== important input to mynnsfclayer, i:", i
919 IF (dry(i)) THEN
920 write(0,*)"dry=",dry(i)," pblh=",pblh(i)," tsk=", tskin_lnd(i),&
921 " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),&
922 " ust=", ust_lnd(i)," snowh=", snowh_lnd(i)," psfcpa=",psfcpa(i), &
923 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
924 ENDIF
925 IF (icy(i)) THEN
926 write(0,*)"icy=",icy(i)," pblh=",pblh(i)," tsk=", tskin_ice(i),&
927 " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),&
928 " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",psfcpa(i), &
929 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
930 ENDIF
931 IF (wet(i)) THEN
932 write(0,*)"wet=",wet(i)," pblh=",pblh(i)," tsk=", tskin_wat(i),&
933 " tsurf=", tsurf_wat(i)," qsfc=", qsfc_wat(i)," znt=", znt_wat(i),&
934 " ust=", ust_wat(i)," snowh=", snowh_wat(i),"psfcpa=",psfcpa(i), &
935 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
936 ENDIF
937 ENDDO
938 ENDIF
939!$acc end serial
940
941!$acc parallel loop present(PSFC, PSFCPA, QVSH, QV1D, THCON, flag_iter, &
942!$acc dry, tskin_lnd, TSK_lnd, tsurf_lnd, THSK_lnd, THVSK_lnd, qsfc_lnd, &
943!$acc icy, tskin_ice, TSK_ice, tsurf_ice, THSK_ice, THVSK_ice, qsfc_ice, &
944!$acc wet, tskin_wat, TSK_wat, tsurf_wat, THSK_wat, THVSK_wat, qsfc_wat)
945 DO i=its,ite
946 ! PSFC ( in cmb) is used later in saturation checks
947 psfc(i)=psfcpa(i)/1000.
948 qvsh(i)=qv1d(i)/(1.+qv1d(i)) !CONVERT TO SPEC HUM (kg/kg)
949 thcon(i)=(100000./psfcpa(i))**rovcp
950 if( flag_iter(i) ) then
951 ! DEFINE SKIN TEMPERATURES FOR LAND/WATER/ICE
952 if(dry(i)) then
953 tsk_lnd(i) = tskin_lnd(i)
954 !TSK_lnd(I) = 0.5 * (tsurf_lnd(i)+tskin_lnd(i))
955 ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE:
956 thsk_lnd(i) = tsk_lnd(i)*thcon(i) !(K)
957 thvsk_lnd(i) = thsk_lnd(i)*(1.+ep1*qsfc_lnd(i))
958 if(thvsk_lnd(i) < 160. .or. thvsk_lnd(i) > 390.) &
959 print *,'THVSK_lnd(I)',itimestep,i,thvsk_lnd(i),thsk_lnd(i),tsurf_lnd(i),tskin_lnd(i),qsfc_lnd(i)
960 endif
961 if(icy(i)) then
962 tsk_ice(i) = tskin_ice(i)
963 !TSK_ice(I) = 0.5 * (tsurf_ice(i)+tskin_ice(i))
964 ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE:
965 thsk_ice(i) = tsk_ice(i)*thcon(i) !(K)
966 thvsk_ice(i) = thsk_ice(i)*(1.+ep1*qsfc_ice(i)) !(K)
967 if(thvsk_ice(i) < 160. .or. thvsk_ice(i) > 390.) &
968 print *,'THVSK_ice(I)',itimestep,i,thvsk_ice(i),thsk_ice(i),tsurf_ice(i),tskin_ice(i),qsfc_ice(i)
969 endif
970 if(wet(i)) then
971 tsk_wat(i) = tskin_wat(i)
972 !TSK_wat(I) = 0.5 * (tsurf_wat(i)+tskin_wat(i))
973 ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE:
974 thsk_wat(i) = tsk_wat(i)*thcon(i) !(K)
975 thvsk_wat(i) = thsk_wat(i)*(1.+ep1*qvsh(i)) !(K)
976 if(thvsk_wat(i) < 160. .or. thvsk_wat(i) > 390.) &
977 print *,'THVSK_wat(I)',i,thvsk_wat(i),thsk_wat(i),tsurf_wat(i),tskin_wat(i),qsfc_wat(i)
978 endif
979 endif ! flag_iter
980 ENDDO
981
982!$acc parallel loop present(TH1D, T1D, P1D, TC1D)
983 DO i=its,ite
984 ! CONVERT LOWEST LAYER TEMPERATURE TO POTENTIAL TEMPERATURE:
985 th1d(i)=t1d(i)*(100000./p1d(i))**rovcp !(Theta, K)
986 tc1d(i)=t1d(i)-273.15 !(T, Celsius)
987 ENDDO
988
989!$acc parallel loop present(THV1D, TH1D, QVSH, TV1D, T1D)
990 DO i=its,ite
991 ! CONVERT TO VIRTUAL TEMPERATURE
992 thv1d(i)=th1d(i)*(1.+ep1*qvsh(i)) !(K)
993 tv1d(i)=t1d(i)*(1.+ep1*qvsh(i)) !(K)
994 ENDDO
995
996!$acc parallel loop present(RHO1D, P1D, TV1D, TH1D, ZA, ZA2, dz2w1d, dz8w1d, GOVRTH)
997 DO i=its,ite
998 rho1d(i)=p1d(i)/(rd*tv1d(i)) !now using value calculated in sfc driver
999 za(i)=0.5*dz8w1d(i) !height of first half-sigma level
1000 za2(i)=dz8w1d(i) + 0.5*dz2w1d(i) !height of 2nd half-sigma level
1001 govrth(i)=grav/th1d(i)
1002 ENDDO
1003
1004 !tgs - should QFX and HFX be separate for land, ice and water?
1005!$acc parallel loop present(QFX, QFLX, RHO1D, HFX, HFLX)
1006 DO i=its,ite
1007 qfx(i)=qflx(i)*rho1d(i)
1008 hfx(i)=hflx(i)*rho1d(i)*cp
1009 ENDDO
1010
1011!$acc serial present(THV1D, TV1D, RHO1D, GOVRTH, &
1012!$acc dry, tsk_lnd, thvsk_lnd, &
1013!$acc icy, tsk_ice, thvsk_ice, &
1014!$acc wet, tsk_wat, thvsk_wat)
1015 IF (debug_code ==2) THEN
1016 !write(*,*)"ITIMESTEP=",ITIMESTEP
1017 DO i=its,ite
1018 write(*,*)"=== derived quantities in mynn sfc layer, i:", i
1019 write(*,*)" land, ice, water"
1020 write(*,*)"dry=",dry(i)," icy=",icy(i)," wet=",wet(i)
1021 write(*,*)"tsk=", tsk_lnd(i),tsk_ice(i),tsk_wat(i)
1022 write(*,*)"thvsk=", thvsk_lnd(i),thvsk_ice(i),thvsk_wat(i)
1023 write(*,*)"THV1D=", thv1d(i)," TV1D=",tv1d(i)
1024 write(*,*)"RHO1D=", rho1d(i)," GOVRTH=",govrth(i)
1025 ENDDO
1026 ENDIF
1027!$acc end serial
1028
1029!$acc parallel loop present(T1D,P1D,QGH,QV1D,CPM)
1030 DO i=its,ite
1031 ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP
1032 ! Q2SAT = QGH IN LSM
1033 IF (t1d(i) .LT. 273.15) THEN
1034 !SATURATION VAPOR PRESSURE WRT ICE
1035 e1=svp1*exp(4648.*(1./273.15 - 1./t1d(i)) - &
1036 & 11.64*log(273.15/t1d(i)) + 0.02265*(273.15 - t1d(i)))
1037 ELSE
1038 !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980)
1039 e1=svp1*exp(svp2*(t1d(i)-svpt0)/(t1d(i)-svp3))
1040 ENDIF
1041 pl=p1d(i)/1000.
1042 !QGH(I)=EP2*E1/(PL-ep3*E1) !specific humidity
1043 qgh(i)=ep2*e1/(pl-e1) !mixing ratio
1044 cpm(i)=cp*(1.+0.84*qv1d(i))
1045 ENDDO
1046
1047!$acc serial present(QGH, &
1048!$acc wet, QSFC_wat, QSFCMR_wat, &
1049!$acc dry, QSFC_lnd, QSFCMR_lnd, &
1050!$acc icy, QSFC_ice, QSFCMR_ice)
1051 IF (debug_code == 2) THEN
1052 write(*,*)"ITIMESTEP=",itimestep
1053 DO i=its,ite
1054 if (wet(i)) then
1055 write(*,*)"==== q-bombs, i:",i," wet"
1056 write(*,*)"QSFC_wat=", qsfc_wat(i)," QSFCMR_wat=", qsfcmr_wat(i)," QGH=",qgh(i)
1057 endif
1058 if(dry(i)) then
1059 write(*,*)"==== q-bombs, i:",i," dry"
1060 write(*,*)"QSFC_lnd=", qsfc_lnd(i)," QSFCMR_lnd=", qsfcmr_lnd(i)," QGH=",qgh(i)
1061 endif
1062 if(icy(i)) then
1063 write(*,*)"==== q-bombs, i:",i," ice"
1064 write(*,*)"QSFC_ice=", qsfc_ice(i)," QSFCMR_ice=", qsfcmr_ice(i)," QGH=",qgh(i)
1065 endif
1066 ENDDO
1067 ENDIF
1068!$acc end serial
1069
1070!$acc parallel loop present(flag_iter,U1D,V1D,WSPD,wet,dry,icy, &
1071!$acc THV1D,THVSK_wat,THVSK_lnd,THVSK_ice, &
1072!$acc hfx,RHO1D,qfx,WSTAR,pblh,dx,GOVRTH,ZA, &
1073!$acc TSK_wat,TSK_lnd,TSK_ice, &
1074!$acc rb_wat,rb_lnd,rb_ice)
1075 DO i=its,ite
1076 if( flag_iter(i) ) then
1077 ! DH* 20200401 - note. A weird bug in Intel 18 on hera prevents using the
1078 ! normal -O2 optimization in Release mode for this file. Not reproducible
1079 ! by every user, the bug manifests itself in the resulting wind speed WSPD(I)
1080 ! being -99.0 despite the assignments in lines 932 and 933. *DH
1081 wspd(i)=sqrt(u1d(i)*u1d(i)+v1d(i)*v1d(i))
1082 wspd_wat = -99.
1083 wspd_ice = -99.
1084 wspd_lnd = -99.
1085
1086 IF (wet(i)) THEN
1087 dthvdz=(thv1d(i)-thvsk_wat(i))
1088 !--------------------------------------------------------
1089 ! Calculate the convective velocity scale (WSTAR) and
1090 ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS)
1091 ! and Mahrt and Sun (1995, MWR), respectively
1092 !-------------------------------------------------------
1093 !tgs - the line below could be used when hflx_wat,qflx_wat are moved from
1094 ! Interstitial to Sfcprop
1095 !fluxc = max(hflx_wat(i) + ep1*THVSK_wat(I)*qflx_wat(i),0.)
1096 fluxc = max(hfx(i)/rho1d(i)/cp &
1097 & + ep1*thvsk_wat(i)*qfx(i)/rho1d(i),0._kind_phys)
1098 !WSTAR(I) = vconvc*(grav/TSK(i)*pblh(i)*fluxc)**onethird
1099 wstar(i) = vconvc*(grav/tsk_wat(i)*pblh(i)*fluxc)**onethird
1100 !--------------------------------------------------------
1101 ! Mahrt and Sun low-res correction - modified for water points (halved)
1102 ! (for 13 km ~ 0.18 m/s; for 3 km == 0 m/s)
1103 !--------------------------------------------------------
1104 vsgd = min( 0.25 * (max(dx(i)/5000.-1.,0._kind_phys))**onethird , 0.5_kind_phys)
1105 wspd_wat=sqrt(wspd(i)*wspd(i)+wstar(i)*wstar(i)+vsgd*vsgd)
1106 wspd_wat=max(wspd_wat,wmin)
1107 !--------------------------------------------------------
1108 ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER,
1109 ! ACCORDING TO AKB(1976), EQ(12).
1110 !--------------------------------------------------------
1111 rb_wat(i)=govrth(i)*za(i)*dthvdz/(wspd_wat*wspd_wat)
1112 rb_wat(i)=max(rb_wat(i),-2.0_kind_phys)
1113 rb_wat(i)=min(rb_wat(i), 2.0_kind_phys)
1114 ENDIF ! end water point
1115
1116 IF (dry(i)) THEN
1117 dthvdz=(thv1d(i)-thvsk_lnd(i))
1118 !--------------------------------------------------------
1119 ! Calculate the convective velocity scale (WSTAR) and
1120 ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS)
1121 ! and Mahrt and Sun (1995, MWR), respectively
1122 !-------------------------------------------------------
1123 !tgs - the line below could be used when hflx_lnd,qflx_wat are moved from
1124 ! Interstitial to Sfcprop
1125 !fluxc = max(hflx_lnd(i) + ep1*THVSK_lnd(I)*qflx_lnd(i),0.)
1126 fluxc = max(hfx(i)/rho1d(i)/cp &
1127 & + ep1*thvsk_lnd(i)*qfx(i)/rho1d(i),0._kind_phys)
1128 ! WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**onethird
1129 ! increase height scale, assuming that the non-local transoport
1130 ! from the mass-flux (plume) mixing exceedsd the PBLH.
1131 wstar(i) = vconvc*(grav/tsk_lnd(i)*min(1.5*pblh(i),4000._kind_phys)*fluxc)**onethird
1132 !--------------------------------------------------------
1133 ! Mahrt and Sun low-res correction
1134 ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s)
1135 !--------------------------------------------------------
1136 vsgd = min( 0.32 * (max(dx(i)/5000.-1.,0._kind_phys))**onethird , 0.5_kind_phys)
1137 wspd_lnd=sqrt(wspd(i)*wspd(i)+wstar(i)*wstar(i)+vsgd*vsgd)
1138 wspd_lnd=max(wspd_lnd,wmin)
1139 !--------------------------------------------------------
1140 ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER,
1141 ! ACCORDING TO AKB(1976), EQ(12).
1142 !--------------------------------------------------------
1143 rb_lnd(i)=govrth(i)*za(i)*dthvdz/(wspd_lnd*wspd_lnd)
1144 !From Tilden Meyers:
1145 !IF (rb_lnd(I) .GE 0.0) THEN
1146 ! ust_lnd(i)=WSPD_lnd*0.1/(1.0 + 10.0*rb_lnd(I))
1147 !ELSE
1148 ! ust_lnd(i)=WSPD_lnd*0.1*(1.0 - 10.0*rb_lnd(I))**onethird
1149 !ENDIF
1150 rb_lnd(i)=max(rb_lnd(i),-2.0_kind_phys)
1151 rb_lnd(i)=min(rb_lnd(i), 2.0_kind_phys)
1152 ENDIF ! end land point
1153
1154 IF (icy(i)) THEN
1155 dthvdz=(thv1d(i)-thvsk_ice(i))
1156 !--------------------------------------------------------
1157 ! Calculate the convective velocity scale (WSTAR) and
1158 ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS)
1159 ! and Mahrt and Sun (1995, MWR), respectively
1160 !-------------------------------------------------------
1161 !tgs - the line below could be used when hflx_ice,qflx_ice are moved from
1162 ! Interstitial to Sfcprop
1163 !fluxc = max(hflx_ice(i) + ep1*THVSK_ice(I)*qflx_ice(i)/RHO1D(i),0.)
1164 fluxc = max(hfx(i)/rho1d(i)/cp &
1165 & + ep1*thvsk_ice(i)*qfx(i)/rho1d(i),0._kind_phys)
1166 ! WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**onethird
1167 ! increase height scale, assuming that the non-local transport
1168 ! from the mass-flux (plume) mixing exceedsd the PBLH.
1169 wstar(i) = vconvc*(grav/tsk_ice(i)*min(1.5*pblh(i),4000._kind_phys)*fluxc)**onethird
1170 !--------------------------------------------------------
1171 ! Mahrt and Sun low-res correction
1172 ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s)
1173 !--------------------------------------------------------
1174 vsgd = min( 0.32 * (max(dx(i)/5000.-1.,0._kind_phys))**onethird , 0.5_kind_phys)
1175 wspd_ice=sqrt(wspd(i)*wspd(i)+wstar(i)*wstar(i)+vsgd*vsgd)
1176 wspd_ice=max(wspd_ice,wmin)
1177 !--------------------------------------------------------
1178 ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER,
1179 ! ACCORDING TO AKB(1976), EQ(12).
1180 !--------------------------------------------------------
1181 rb_ice(i)=govrth(i)*za(i)*dthvdz/(wspd_ice*wspd_ice)
1182 rb_ice(i)=max(rb_ice(i),-2.0_kind_phys)
1183 rb_ice(i)=min(rb_ice(i), 2.0_kind_phys)
1184 ENDIF ! end ice point
1185
1186 !NOW CONDENSE THE POSSIBLE WSPD VALUES BY TAKING THE MAXIMUM
1187 wspd(i) = max(wspd_ice,wspd_wat)
1188 wspd(i) = max(wspd_lnd,wspd(i))
1189
1190 IF (debug_code == 2) THEN
1191 write(*,*)"===== After rb calc in mynn sfc layer:"
1192 write(*,*)"ITIMESTEP=",itimestep
1193 write(*,*)"WSPD=", wspd(i)," WSTAR=", wstar(i)," vsgd=",vsgd
1194 IF (icy(i))write(*,*)"rb_ice=", rb_ice(i)," DTHVDZ=",dthvdz
1195 IF (wet(i))write(*,*)"rb_wat=", rb_wat(i)," DTHVDZ=",dthvdz
1196 IF (dry(i))write(*,*)"rb_lnd=", rb_lnd(i)," DTHVDZ=",dthvdz
1197 ENDIF
1198
1199 ! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2 (STABLE)
1200 !if (itimestep .GT. 1) THEN
1201 ! IF(MOL(I).LT.0.)BR(I)=MIN(BR(I),0.0)
1202 !ENDIF
1203
1204 endif ! flag_iter
1205 ENDDO
1206
1207 1006 format(a,f7.3,a,f9.4,a,f9.5,a,f9.4)
1208 1007 format(a,f2.0,a,f6.2,a,f7.3,a,f7.2)
1209
1210!--------------------------------------------------------------------
1211!--------------------------------------------------------------------
1212!--- BEGIN I-LOOP
1213!--------------------------------------------------------------------
1214!--------------------------------------------------------------------
1215
1216!$acc parallel loop present(flag_iter, PSFCPA, dz8w1d, pblh, &
1217!$acc device_errmsg, device_errflg, &
1218!$acc device_special_errmsg, device_special_errflg, &
1219!$acc wet, dry, icy, &
1220!$acc ZT_wat, ZT_lnd, ZT_ice, &
1221!$acc ZNT_wat, ZNT_lnd, ZNT_ice, &
1222!$acc ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
1223!$acc UST_wat, UST_lnd, UST_ice, &
1224!$acc ZQ_wat, ZQ_lnd, ZQ_ice, &
1225!$acc snowh_wat, snowh_lnd, snowh_ice, &
1226!$acc THVSK_wat, THVSK_lnd, THVSK_ice, &
1227!$acc tskin_wat, tskin_lnd, tskin_ice, &
1228!$acc tsurf_wat, tsurf_lnd, tsurf_ice, &
1229!$acc qsfc_wat, qsfc_lnd, qsfc_ice, &
1230!$acc GZ1OZ0_wat, GZ1OZt_wat, GZ2OZ0_wat, GZ2OZt_wat, GZ10OZ0_wat, GZ10OZt_wat, &
1231!$acc GZ1OZ0_lnd, GZ1OZt_lnd, GZ2OZ0_lnd, GZ2OZt_lnd, GZ10OZ0_lnd, GZ10OZt_lnd, &
1232!$acc GZ1OZ0_ice, GZ1OZt_ice, GZ2OZ0_ice, GZ2OZt_ice, GZ10OZ0_ice, GZ10OZt_ice, &
1233!$acc zratio_wat, zratio_lnd, zratio_ice, &
1234!$acc stress_wat, stress_lnd, stress_ice, &
1235!$acc rb_wat, rb_lnd, rb_ice, &
1236!$acc qflx, qflx_lnd, &
1237!$acc hflx, hflx_lnd, &
1238!$acc psim, psih, psim10, psih10, psih2, &
1239!$acc psix_wat, psix10_wat, psit_wat, psit2_wat, psiq_wat, psiq2_wat, &
1240!$acc psix_lnd, psix10_lnd, psit_lnd, psit2_lnd, psiq_lnd, psiq2_lnd, &
1241!$acc psix_ice, psix10_ice, psit_ice, psit2_ice, psiq_ice, psiq2_ice, &
1242!$acc WSPD, WSPDI, U1D, V1D, TC1D, THV1D, rstoch1D, USTM, ZA, ZOL, QVSH, &
1243!$acc shdmax, vegtype, z0pert, ztpert, mol, rmol, wstar, qstar, sigmaf)
1244
1245 DO i=its,ite
1246 if( flag_iter(i) ) then
1247
1248 !COMPUTE KINEMATIC VISCOSITY (m2/s) Andreas (1989) CRREL Rep. 89-11
1249 !valid between -173 and 277 degrees C.
1250 visc=1.326e-5*(1. + 6.542e-3*tc1d(i) + 8.301e-6*tc1d(i)*tc1d(i) &
1251 - 4.84e-9*tc1d(i)*tc1d(i)*tc1d(i))
1252
1253 IF (wet(i)) THEN
1254 !--------------------------------------
1255 ! WATER
1256 !--------------------------------------
1257 if (sfc_z0_type >= 0) then ! Avoid calculation is using wave model
1258 ! CALCULATE z0 (znt)
1259 !--------------------------------------
1260
1261 IF (debug_code == 2) THEN
1262 write(*,*)"=============Input to ZNT over water:"
1263 write(*,*)"u*:",ust_wat(i)," wspd=",wspd(i)," visc=",visc," za=",za(i)
1264 ENDIF
1265
1266 IF ( PRESENT(isftcflx) ) THEN
1267 IF ( isftcflx .EQ. 0 ) THEN
1268 IF (coare_opt .EQ. 3.0) THEN
1269 !COARE 3.0 (MISLEADING SUBROUTINE NAME)
1270 CALL charnock_1955(znt_wat(i),ust_wat(i),wspd(i),visc,za(i))
1271 ELSE
1272 !COARE 3.5
1273 CALL edson_etal_2013(znt_wat(i),ust_wat(i),wspd(i),visc,za(i))
1274 ENDIF
1275 ELSEIF ( isftcflx .EQ. 1 .OR. isftcflx .EQ. 2 ) THEN
1276 CALL davis_etal_2008(znt_wat(i),ust_wat(i))
1277 ELSEIF ( isftcflx .EQ. 3 ) THEN
1278 CALL taylor_yelland_2001(znt_wat(i),ust_wat(i),wspd(i))
1279 ELSEIF ( isftcflx .EQ. 4 ) THEN
1280 !GFS surface layer scheme
1281 CALL gfs_z0_wat(znt_wat(i),ust_wat(i),wspd(i),za(i),sfc_z0_type,redrag)
1282 ENDIF
1283 ELSE
1284 !DEFAULT TO COARE 3.0/3.5
1285 IF (coare_opt .EQ. 3.0) THEN
1286 !COARE 3.0
1287 CALL charnock_1955(znt_wat(i),ust_wat(i),wspd(i),visc,za(i))
1288 ELSE
1289 !COARE 3.5
1290 CALL edson_etal_2013(znt_wat(i),ust_wat(i),wspd(i),visc,za(i))
1291 ENDIF
1292 ENDIF
1293 endif !-end wave model check
1294
1295 ! add stochastic perturbation of ZNT
1296 if (spp_sfc==1) then
1297 zntstoch_wat(i) = max(znt_wat(i) + znt_wat(i)*1.0*rstoch1d(i), 1e-6_kind_phys)
1298 else
1299 zntstoch_wat(i) = znt_wat(i)
1300 endif
1301
1302 IF (debug_code > 1) THEN
1303 write(*,*)"==========Output ZNT over water:"
1304 write(*,*)"ZNT:",zntstoch_wat(i)
1305 ENDIF
1306
1307 !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING NEW ZNT
1308 ! AHW: Garrattt formula: Calculate roughness Reynolds number
1309 ! Kinematic viscosity of air (linear approx to
1310 ! temp dependence at sea level)
1311 restar=max(ust_wat(i)*zntstoch_wat(i)/visc, 0.1_kind_phys)
1312
1313 !--------------------------------------
1314 !CALCULATE z_t and z_q
1315 !--------------------------------------
1316 IF (debug_code > 1) THEN
1317 write(*,*)"=============Input to ZT over water:"
1318 write(*,*)"u*:",ust_wat(i)," restar=",restar," visc=",visc
1319 ENDIF
1320
1321 IF ( PRESENT(isftcflx) ) THEN
1322 IF ( isftcflx .EQ. 0 ) THEN
1323 IF (coare_opt .EQ. 3.0) THEN
1324 CALL fairall_etal_2003(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1325 rstoch1d(i),spp_sfc)
1326 ELSE
1327 CALL fairall_etal_2014(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1328 rstoch1d(i),spp_sfc)
1329 ENDIF
1330 ELSEIF ( isftcflx .EQ. 1 ) THEN
1331 IF (coare_opt .EQ. 3.0) THEN
1332 CALL fairall_etal_2003(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1333 rstoch1d(i),spp_sfc)
1334 ELSE
1335 CALL fairall_etal_2014(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1336 rstoch1d(i),spp_sfc)
1337 ENDIF
1338 ELSEIF ( isftcflx .EQ. 2 ) THEN
1339 CALL garratt_1992(zt_wat(i),zq_wat(i),zntstoch_wat(i),restar,2.0_kind_phys)
1340 ELSEIF ( isftcflx .EQ. 3 ) THEN
1341 IF (coare_opt .EQ. 3.0) THEN
1342 CALL fairall_etal_2003(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1343 rstoch1d(i),spp_sfc)
1344 ELSE
1345 CALL fairall_etal_2014(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1346 rstoch1d(i),spp_sfc)
1347 ENDIF
1348 ELSEIF ( isftcflx .EQ. 4 ) THEN
1349 !GFS zt formulation
1350 CALL gfs_zt_wat(zt_wat(i),zntstoch_wat(i),restar,wspd(i),za(i),sfc_z0_type,device_errmsg,device_errflg)
1351 zq_wat(i)=zt_wat(i)
1352 ENDIF
1353 ELSE
1354 !DEFAULT TO COARE 3.0/3.5
1355 IF (coare_opt .EQ. 3.0) THEN
1356 CALL fairall_etal_2003(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1357 rstoch1d(i),spp_sfc)
1358 ELSE
1359 CALL fairall_etal_2014(zt_wat(i),zq_wat(i),restar,ust_wat(i),visc,&
1360 rstoch1d(i),spp_sfc)
1361 ENDIF
1362 ENDIF
1363
1364 IF (debug_code > 1) THEN
1365 write(*,*)"=============Output ZT & ZQ over water:"
1366 write(*,*)"ZT:",zt_wat(i)," ZQ:",zq_wat(i)
1367 ENDIF
1368
1369 gz1oz0_wat(i)= log((za(i)+zntstoch_wat(i))/zntstoch_wat(i))
1370 gz1ozt_wat(i)= log((za(i)+zntstoch_wat(i))/zt_wat(i))
1371 gz2oz0_wat(i)= log((2.0+zntstoch_wat(i))/zntstoch_wat(i))
1372 gz2ozt_wat(i)= log((2.0+zntstoch_wat(i))/zt_wat(i))
1373 gz10oz0_wat(i)=log((10.+zntstoch_wat(i))/zntstoch_wat(i))
1374 gz10ozt_wat(i)=log((10.+zntstoch_wat(i))/zt_wat(i))
1375 zratio_wat(i)=zntstoch_wat(i)/zt_wat(i) !need estimate for Li et al.
1376
1377 ENDIF !end water point
1378
1379 IF (dry(i)) THEN
1380
1381 if ( iz0tlnd .EQ. 4 ) then
1382 CALL gfs_z0_lnd(znt_lnd(i),shdmax(i),za(i),vegtype(i),ivegsrc,z0pert(i))
1383 endif
1384
1385 ! add stochastic perturbaction of ZNT
1386 if (spp_sfc==1) then
1387 zntstoch_lnd(i) = max(znt_lnd(i) + znt_lnd(i)*1.0*rstoch1d(i), 1e-6_kind_phys)
1388 else
1389 zntstoch_lnd(i) = znt_lnd(i)
1390 endif
1391 !add limit to prevent ridiculous values of z0 (more than dz/15)
1392 zntstoch_lnd(i) = min(zntstoch_lnd(i), dz8w1d(i)*0.0666_kind_phys)
1393
1394 !--------------------------------------
1395 ! LAND
1396 !--------------------------------------
1397 !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT
1398 restar=max(ust_lnd(i)*zntstoch_lnd(i)/visc, 0.1_kind_phys)
1399
1400 !--------------------------------------
1401 !GET z_t and z_q
1402 !--------------------------------------
1403 IF (snowh_lnd(i) > 50.) THEN ! (mm) Treat as snow cover - use Andreas
1404 CALL andreas_2002(zntstoch_lnd(i),visc,ust_lnd(i),zt_lnd(i),zq_lnd(i))
1405 ELSE
1406 IF ( PRESENT(iz0tlnd) ) THEN
1407 IF ( iz0tlnd .LE. 1 ) THEN
1408 CALL zilitinkevich_1995(zntstoch_lnd(i),zt_lnd(i),zq_lnd(i),restar,&
1409 ust_lnd(i),karman,1.0_kind_phys,iz0tlnd,spp_sfc,rstoch1d(i))
1410 ELSEIF ( iz0tlnd .EQ. 2 ) THEN
1411 ! DH note - at this point, qstar is either not initialized
1412 ! or initialized to zero, but certainly not set correctly
1413 device_special_errmsg = 'Logic error: qstar is not set correctly when calling Yang_2008'
1414 device_special_errflg = 1
1415#ifndef _OPENACC
1416! Necessary since OpenACC does not support branching in parallel code
1417! Must sync errmsg and errflg with device_errmsg and device_errflg, respectively
1418! so that proper error message and error flag codes are returned.
1419 errmsg = device_special_errmsg
1420 errflg = device_special_errflg
1421 return
1422#endif
1423 CALL yang_2008(zntstoch_lnd(i),zt_lnd(i),zq_lnd(i),ust_lnd(i),mol(i),&
1424 qstar(i),restar,visc)
1425 ELSEIF ( iz0tlnd .EQ. 3 ) THEN
1426 !Original MYNN in WRF-ARW used this form:
1427 CALL garratt_1992(zt_lnd(i),zq_lnd(i),zntstoch_lnd(i),restar,1.0_kind_phys)
1428 ELSEIF ( iz0tlnd .EQ. 4 ) THEN
1429 !GFS:
1430 CALL gfs_zt_lnd(zt_lnd(i),zntstoch_lnd(i),sigmaf(i),ztpert(i),ust_lnd(i))
1431 zq_lnd(i)=zt_lnd(i)
1432 ENDIF
1433 ELSE
1434 !DEFAULT TO ZILITINKEVICH
1435 CALL zilitinkevich_1995(zntstoch_lnd(i),zt_lnd(i),zq_lnd(i),restar,&
1436 ust_lnd(i),karman,1.0_kind_phys,0,spp_sfc,rstoch1d(i))
1437 ENDIF
1438 ENDIF
1439
1440 IF (zntstoch_lnd(i) < 1e-8 .OR. zt_lnd(i) < 1e-10) THEN
1441 write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i
1442 write(0,*)" ZNT=", zntstoch_lnd(i)," ZT=",zt_lnd(i)
1443 write(0,*)" tsk=", tskin_lnd(i)," restar=",restar,&
1444 " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),&
1445 " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",psfcpa(i), &
1446 " dz=",dz8w1d(i)," qflx=",qflx_lnd(i)," hflx=",hflx_lnd(i)," hpbl=",pblh(i)
1447 ENDIF
1448
1449 gz1oz0_lnd(i)= log((za(i)+zntstoch_lnd(i))/zntstoch_lnd(i))
1450 gz1ozt_lnd(i)= log((za(i)+zntstoch_lnd(i))/zt_lnd(i))
1451 gz2oz0_lnd(i)= log((2.0+zntstoch_lnd(i))/zntstoch_lnd(i))
1452 gz2ozt_lnd(i)= log((2.0+zntstoch_lnd(i))/zt_lnd(i))
1453 gz10oz0_lnd(i)=log((10.+zntstoch_lnd(i))/zntstoch_lnd(i))
1454 gz10ozt_lnd(i)=log((10.+zntstoch_lnd(i))/zt_lnd(i))
1455 zratio_lnd(i)=zntstoch_lnd(i)/zt_lnd(i) !need estimate for Li et al.
1456
1457 ENDIF !end land point
1458
1459 IF (icy(i)) THEN
1460
1461 ! add stochastic perturbaction of ZNT
1462 if (spp_sfc==1) then
1463 zntstoch_ice(i) = max(znt_ice(i) + znt_ice(i)*1.0*rstoch1d(i), 1e-6_kind_phys)
1464 else
1465 zntstoch_ice(i) = znt_ice(i)
1466 endif
1467
1468 !--------------------------------------
1469 ! ICE
1470 !--------------------------------------
1471 !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT
1472 restar=max(ust_ice(i)*zntstoch_ice(i)/visc, 0.1_kind_phys)
1473 !--------------------------------------
1474 !GET z_t and z_q
1475 !--------------------------------------
1476 CALL andreas_2002(zntstoch_ice(i),visc,ust_ice(i),zt_ice(i),zq_ice(i))
1477
1478 gz1oz0_ice(i)= log((za(i)+zntstoch_ice(i))/zntstoch_ice(i))
1479 gz1ozt_ice(i)= log((za(i)+zntstoch_ice(i))/zt_ice(i))
1480 gz2oz0_ice(i)= log((2.0+zntstoch_ice(i))/zntstoch_ice(i))
1481 gz2ozt_ice(i)= log((2.0+zntstoch_ice(i))/zt_ice(i))
1482 gz10oz0_ice(i)=log((10.+zntstoch_ice(i))/zntstoch_ice(i))
1483 gz10ozt_ice(i)=log((10.+zntstoch_ice(i))/zt_ice(i))
1484 zratio_ice(i)=zntstoch_ice(i)/zt_ice(i) !need estimate for Li et al.
1485
1486 ENDIF !end ice point
1487
1488 !Capture a representative ZNT
1489 !tgs - should this be changed for fractional grid or fractional sea ice?
1490 IF (dry(i)) THEN
1491 znt(i)=zntstoch_lnd(i)
1492 ELSEIF (wet(i)) THEN
1493 znt(i)=zntstoch_wat(i)
1494 ELSEIF (icy(i)) THEN
1495 znt(i)=zntstoch_ice(i)
1496 ENDIF
1497
1498 !--------------------------------------------------------------------
1499 !--- DIAGNOSE STABILITY FUNCTIONS FOR THE APPROPRIATE STABILITY CLASS:
1500 ! THE STABILITY CLASSES ARE DETERMINED BY THE BULK RICHARDSON NUMBER.
1501 !--------------------------------------------------------------------
1502
1503 IF (wet(i)) THEN
1504 IF (rb_wat(i) .GT. 0.0) THEN
1505
1506 IF (.not. flag_restart .or. (flag_restart .and. itimestep > 1) ) THEN
1507 !COMPUTE z/L first guess:
1508 CALL li_etal_2010(zol(i),rb_wat(i),za(i)/zntstoch_wat(i),zratio_wat(i))
1509 !ZOL(I)=ZA(I)*KARMAN*grav*MOL(I)/(TH1D(I)*MAX(UST_wat(I)*UST_wat(I),0.0001))
1510 zol(i)=max(zol(i),0.0_kind_phys)
1511 zol(i)=min(zol(i),20._kind_phys)
1512
1513 IF (debug_code >= 1) THEN
1514 IF (zntstoch_wat(i) < 1e-8 .OR. zt_wat(i) < 1e-10) THEN
1515 write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i
1516 write(0,*)"rb=", rb_wat(i)," ZNT=", zntstoch_wat(i)," ZT=",zt_wat(i)
1517 write(0,*)" tsk=", tskin_wat(i)," prev z/L=",zol(i),&
1518 " tsurf=", tsurf_wat(i)," qsfc=", qsfc_wat(i)," znt=", znt_wat(i),&
1519 " ust=", ust_wat(i)," snowh=", snowh_wat(i),"psfcpa=",psfcpa(i), &
1520 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
1521 ENDIF
1522 ENDIF
1523
1524 !Use Pedros iterative function to find z/L
1525 !zol(I)=zolri(rb_wat(I),ZA(I),ZNTstoch_wat(I),ZT_wat(I),ZOL(I),psi_opt)
1526 !Use brute-force method
1527 zol(i)=zolrib(rb_wat(i),za(i),zntstoch_wat(i),zt_wat(i),gz1oz0_wat(i),gz1ozt_wat(i),zol(i),psi_opt)
1528 ENDIF ! restart
1529 zol(i)=max(zol(i),0.0_kind_phys)
1530 zol(i)=min(zol(i),20._kind_phys)
1531
1532 zolzt = zol(i)*zt_wat(i)/za(i) ! zt/L
1533 zolz0 = zol(i)*zntstoch_wat(i)/za(i) ! z0/L
1534 zolza = zol(i)*(za(i)+zntstoch_wat(i))/za(i) ! (z+z0/L
1535 zol10 = zol(i)*(10.+zntstoch_wat(i))/za(i) ! (10+z0)/L
1536 zol2 = zol(i)*(2.+zntstoch_wat(i))/za(i) ! (2+z0)/L
1537
1538 !COMPUTE PSIM and PSIH
1539 !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I))
1540 !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
1541 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
1542 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_wat(I),ZNTstoch_wat(I),ZA(I))
1543 !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0)
1544 ! or use tables
1545 psim(i)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt)
1546 psih(i)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt)
1547 psim10(i)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt)
1548 psih10(i)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt)
1549 psih2(i)=psih_stable(zol2,psi_opt)-psih_stable(zolzt,psi_opt)
1550
1551 ! 1.0 over Monin-Obukhov length
1552 rmol(i)= zol(i)/za(i)
1553
1554 ELSEIF(rb_wat(i) .EQ. 0.) THEN
1555 !=========================================================
1556 !-----CLASS 3; FORCED CONVECTION/NEUTRAL:
1557 !=========================================================
1558
1559 psim(i)=0.0
1560 psih(i)=psim(i)
1561 psim10(i)=0.
1562 psih10(i)=0.
1563 psih2(i)=0.
1564
1565 zol(i) =0.
1566 rmol(i) =0.
1567
1568 ELSEIF(rb_wat(i) .LT. 0.)THEN
1569 !==========================================================
1570 !-----CLASS 4; FREE CONVECTION:
1571 !==========================================================
1572
1573 !COMPUTE z/L first guess:
1574 IF (.not. flag_restart .or. (flag_restart .and. itimestep > 1) ) THEN
1575 CALL li_etal_2010(zol(i),rb_wat(i),za(i)/zntstoch_wat(i),zratio_wat(i))
1576 !ZOL(I)=ZA(I)*KARMAN*grav*MOL(I)/(TH1D(I)*MAX(UST_wat(I)*UST_wat(I),0.001))
1577 zol(i)=max(zol(i),-20.0_kind_phys)
1578 zol(i)=min(zol(i),0.0_kind_phys)
1579
1580 IF (debug_code >= 1) THEN
1581 IF (zntstoch_wat(i) < 1e-8 .OR. zt_wat(i) < 1e-10) THEN
1582 write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i
1583 write(0,*)"rb=", rb_wat(i)," ZNT=", zntstoch_wat(i)," ZT=",zt_wat(i)
1584 write(0,*)" tsk=", tskin_wat(i)," wstar=",wstar(i)," prev z/L=",zol(i),&
1585 " tsurf=", tsurf_wat(i)," qsfc=", qsfc_wat(i)," znt=", znt_wat(i),&
1586 " ust=", ust_wat(i)," snowh=", snowh_wat(i),"psfcpa=",psfcpa(i), &
1587 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
1588 ENDIF
1589 ENDIF
1590
1591 !Use Pedros iterative function to find z/L
1592 !zol(I)=zolri(rb_wat(I),ZA(I),ZNTstoch_wat(I),ZT_wat(I),ZOL(I),psi_opt)
1593 !Use brute-force method
1594 zol(i)=zolrib(rb_wat(i),za(i),zntstoch_wat(i),zt_wat(i),gz1oz0_wat(i),gz1ozt_wat(i),zol(i),psi_opt)
1595 ENDIF ! restart
1596 zol(i)=max(zol(i),-20.0_kind_phys)
1597 zol(i)=min(zol(i),0.0_kind_phys)
1598
1599 zolzt = zol(i)*zt_wat(i)/za(i) ! zt/L
1600 zolz0 = zol(i)*zntstoch_wat(i)/za(i) ! z0/L
1601 zolza = zol(i)*(za(i)+zntstoch_wat(i))/za(i) ! (z+z0/L
1602 zol10 = zol(i)*(10.+zntstoch_wat(i))/za(i) ! (10+z0)/L
1603 zol2 = zol(i)*(2.+zntstoch_wat(i))/za(i) ! (2+z0)/L
1604
1605 !COMPUTE PSIM and PSIH
1606 !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I))
1607 !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), ZT_wat(I), ZNTstoch_wat(I), ZA(I))
1608 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
1609 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_wat(I),ZNTstoch_wat(I),ZA(I))
1610 ! use tables
1611 psim(i)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt)
1612 psih(i)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt)
1613 psim10(i)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt)
1614 psih10(i)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt)
1615 psih2(i)=psih_unstable(zol2,psi_opt)-psih_unstable(zolzt,psi_opt)
1616
1617 !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND
1618 !---HIGH ROUGHNESS. THIS PREVENTS DENOMINATOR IN FLUXES
1619 !---FROM GETTING TOO SMALL
1620 psih(i)=min(psih(i),0.9*gz1ozt_wat(i))
1621 psim(i)=min(psim(i),0.9*gz1oz0_wat(i))
1622 psih2(i)=min(psih2(i),0.9*gz2ozt_wat(i))
1623 psim10(i)=min(psim10(i),0.9*gz10oz0_wat(i))
1624 psih10(i)=min(psih10(i),0.9*gz10ozt_wat(i))
1625
1626 rmol(i) = zol(i)/za(i)
1627
1628 ENDIF
1629
1630 ! CALCULATE THE RESISTANCE:
1631 psix_wat(i) =max(gz1oz0_wat(i)-psim(i) , 1.0_kind_phys) ! = fm
1632 psix10_wat(i)=max(gz10oz0_wat(i)-psim10(i), 1.0_kind_phys) ! = fm10
1633 psit_wat(i) =max(gz1ozt_wat(i)-psih(i) , 1.0_kind_phys) ! = fh
1634 psit2_wat(i) =max(gz2ozt_wat(i)-psih2(i) , 1.0_kind_phys) ! = fh2
1635 psiq_wat(i) =max(log((za(i)+zq_wat(i))/zq_wat(i))-psih(i) ,1.0_kind_phys)
1636 psiq2_wat(i) =max(log((2.0+zq_wat(i))/zq_wat(i))-psih2(i) ,1.0_kind_phys)
1637
1638 ENDIF ! end water points
1639
1640 IF (dry(i)) THEN
1641 IF (rb_lnd(i) .GT. 0.0) THEN
1642
1643 IF (.not. flag_restart .or. (flag_restart .and. itimestep > 1) ) THEN
1644 !COMPUTE z/L first guess:
1645 CALL li_etal_2010(zol(i),rb_lnd(i),za(i)/zntstoch_lnd(i),zratio_lnd(i))
1646 !ZOL(I)=ZA(I)*KARMAN*grav*MOL(I)/(TH1D(I)*MAX(UST_lnd(I)*UST_lnd(I),0.0001))
1647 zol(i)=max(zol(i),0.0_kind_phys)
1648 zol(i)=min(zol(i),20._kind_phys)
1649
1650 IF (debug_code >= 1) THEN
1651 IF (zntstoch_lnd(i) < 1e-8 .OR. zt_lnd(i) < 1e-10) THEN
1652 write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i
1653 write(0,*)"rb=", rb_lnd(i)," ZNT=", zntstoch_lnd(i)," ZT=",zt_lnd(i)
1654 write(0,*)" tsk=", tskin_lnd(i)," prev z/L=",zol(i),&
1655 " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),&
1656 " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",psfcpa(i), &
1657 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
1658 ENDIF
1659 ENDIF
1660
1661 !Use Pedros iterative function to find z/L
1662 !zol(I)=zolri(rb_lnd(I),ZA(I),ZNTstoch_lnd(I),ZT_lnd(I),ZOL(I),psi_opt)
1663 !Use brute-force method
1664 zol(i)=zolrib(rb_lnd(i),za(i),zntstoch_lnd(i),zt_lnd(i),gz1oz0_lnd(i),gz1ozt_lnd(i),zol(i),psi_opt)
1665 ENDIF ! restart
1666 zol(i)=max(zol(i),0.0_kind_phys)
1667 zol(i)=min(zol(i),20._kind_phys)
1668
1669 zolzt = zol(i)*zt_lnd(i)/za(i) ! zt/L
1670 zolz0 = zol(i)*zntstoch_lnd(i)/za(i) ! z0/L
1671 zolza = zol(i)*(za(i)+zntstoch_lnd(i))/za(i) ! (z+z0/L
1672 zol10 = zol(i)*(10.+zntstoch_lnd(i))/za(i) ! (10+z0)/L
1673 zol2 = zol(i)*(2.+zntstoch_lnd(i))/za(i) ! (2+z0)/L
1674
1675 !COMPUTE PSIM and PSIH
1676 !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
1677 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
1678 !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I))
1679 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_lnd(I),ZNTstoch_lnd(I),ZA(I))
1680 !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0)
1681 psim(i)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt)
1682 psih(i)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt)
1683 psim10(i)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt)
1684 psih10(i)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt)
1685 psih2(i)=psih_stable(zol2,psi_opt)-psih_stable(zolzt,psi_opt)
1686
1687 ! 1.0 over Monin-Obukhov length
1688 rmol(i)= zol(i)/za(i)
1689
1690 ELSEIF(rb_lnd(i) .EQ. 0.) THEN
1691 !=========================================================
1692 !-----CLASS 3; FORCED CONVECTION/NEUTRAL:
1693 !=========================================================
1694
1695 psim(i)=0.0
1696 psih(i)=psim(i)
1697 psim10(i)=0.
1698 psih10(i)=0.
1699 psih2(i)=0.
1700
1701 zol(i) =0.
1702 rmol(i) =0.
1703
1704 ELSEIF(rb_lnd(i) .LT. 0.)THEN
1705 !==========================================================
1706 !-----CLASS 4; FREE CONVECTION:
1707 !==========================================================
1708
1709 IF (.not. flag_restart .or. (flag_restart .and. itimestep > 1) ) THEN
1710 !COMPUTE z/L first guess:
1711 CALL li_etal_2010(zol(i),rb_lnd(i),za(i)/zntstoch_lnd(i),zratio_lnd(i))
1712 !ZOL(I)=ZA(I)*KARMAN*grav*MOL(I)/(TH1D(I)*MAX(UST_lnd(I)*UST_lnd(I),0.001))
1713 zol(i)=max(zol(i),-20.0_kind_phys)
1714 zol(i)=min(zol(i),0.0_kind_phys)
1715
1716 IF (debug_code >= 1) THEN
1717 IF (zntstoch_lnd(i) < 1e-8 .OR. zt_lnd(i) < 1e-10) THEN
1718 write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i
1719 write(0,*)"rb=", rb_lnd(i)," ZNT=", zntstoch_lnd(i)," ZT=",zt_lnd(i)
1720 write(0,*)" tsk=", tskin_lnd(i)," wstar=",wstar(i)," prev z/L=",zol(i),&
1721 " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),&
1722 " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",psfcpa(i), &
1723 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
1724 ENDIF
1725 ENDIF
1726
1727 !Use Pedros iterative function to find z/L
1728 !zol(I)=zolri(rb_lnd(I),ZA(I),ZNTstoch_lnd(I),ZT_lnd(I),ZOL(I),psi_opt)
1729 !Use brute-force method
1730 zol(i)=zolrib(rb_lnd(i),za(i),zntstoch_lnd(i),zt_lnd(i),gz1oz0_lnd(i),gz1ozt_lnd(i),zol(i),psi_opt)
1731 ENDIF ! restart
1732 zol(i)=max(zol(i),-20.0_kind_phys)
1733 zol(i)=min(zol(i),0.0_kind_phys)
1734
1735 zolzt = zol(i)*zt_lnd(i)/za(i) ! zt/L
1736 zolz0 = zol(i)*zntstoch_lnd(i)/za(i) ! z0/L
1737 zolza = zol(i)*(za(i)+zntstoch_lnd(i))/za(i) ! (z+z0/L
1738 zol10 = zol(i)*(10.+zntstoch_lnd(i))/za(i) ! (10+z0)/L
1739 zol2 = zol(i)*(2.+zntstoch_lnd(i))/za(i) ! (2+z0)/L
1740
1741 !COMPUTE PSIM and PSIH
1742 !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), ZT_lnd(I), ZNTstoch_lnd(I), ZA(I))
1743 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
1744 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_lnd(I),ZNTstoch_lnd(I),ZA(I))
1745 ! use tables
1746 psim(i)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt)
1747 psih(i)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt)
1748 psim10(i)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt)
1749 psih10(i)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt)
1750 psih2(i)=psih_unstable(zol2,psi_opt)-psih_unstable(zolzt,psi_opt)
1751
1752 !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND
1753 !---HIGH ROUGHNESS. THIS PREVENTS DENOMINATOR IN FLUXES
1754 !---FROM GETTING TOO SMALL
1755 psih(i)=min(psih(i),0.9*gz1ozt_lnd(i))
1756 psim(i)=min(psim(i),0.9*gz1oz0_lnd(i))
1757 psih2(i)=min(psih2(i),0.9*gz2ozt_lnd(i))
1758 psim10(i)=min(psim10(i),0.9*gz10oz0_lnd(i))
1759 psih10(i)=min(psih10(i),0.9*gz10ozt_lnd(i))
1760
1761 rmol(i) = zol(i)/za(i)
1762
1763 ENDIF
1764
1765 ! CALCULATE THE RESISTANCE:
1766 psix_lnd(i) =max(gz1oz0_lnd(i)-psim(i), 1.0_kind_phys)
1767 psix10_lnd(i)=max(gz10oz0_lnd(i)-psim10(i), 1.0_kind_phys)
1768 psit_lnd(i) =max(gz1ozt_lnd(i)-psih(i) , 1.0_kind_phys)
1769 psit2_lnd(i) =max(gz2ozt_lnd(i)-psih2(i), 1.0_kind_phys)
1770 psiq_lnd(i) =max(log((za(i)+zq_lnd(i))/zq_lnd(i))-psih(i) ,1.0_kind_phys)
1771 psiq2_lnd(i) =max(log((2.0+zq_lnd(i))/zq_lnd(i))-psih2(i) ,1.0_kind_phys)
1772
1773 ENDIF ! end land points
1774
1775 IF (icy(i)) THEN
1776 IF (rb_ice(i) .GT. 0.0) THEN
1777
1778 IF (.not. flag_restart .or. (flag_restart .and. itimestep > 1) ) THEN
1779 !COMPUTE z/L first guess:
1780 CALL li_etal_2010(zol(i),rb_ice(i),za(i)/zntstoch_ice(i),zratio_ice(i))
1781 !ZOL(I)=ZA(I)*KARMAN*grav*MOL(I)/(TH1D(I)*MAX(UST_ice(I)*UST_ice(I),0.0001))
1782 zol(i)=max(zol(i),0.0_kind_phys)
1783 zol(i)=min(zol(i),20._kind_phys)
1784
1785 IF (debug_code >= 1) THEN
1786 IF (zntstoch_ice(i) < 1e-8 .OR. zt_ice(i) < 1e-10) THEN
1787 write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i
1788 write(0,*)"rb=", rb_ice(i)," ZNT=", zntstoch_ice(i)," ZT=",zt_ice(i)
1789 write(0,*)" tsk=", tskin_ice(i)," prev z/L=",zol(i),&
1790 " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),&
1791 " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",psfcpa(i), &
1792 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
1793 ENDIF
1794 ENDIF
1795
1796 !Use Pedros iterative function to find z/L
1797 !zol(I)=zolri(rb_ice(I),ZA(I),ZNTstoch_ice(I),ZT_ice(I),ZOL(I),psi_opt)
1798 !Use brute-force method
1799 zol(i)=zolrib(rb_ice(i),za(i),zntstoch_ice(i),zt_ice(i),gz1oz0_ice(i),gz1ozt_ice(i),zol(i),psi_opt)
1800 ENDIF ! restart
1801 zol(i)=max(zol(i),0.0_kind_phys)
1802 zol(i)=min(zol(i),20._kind_phys)
1803
1804 zolzt = zol(i)*zt_ice(i)/za(i) ! zt/L
1805 zolz0 = zol(i)*zntstoch_ice(i)/za(i) ! z0/L
1806 zolza = zol(i)*(za(i)+zntstoch_ice(i))/za(i) ! (z+z0/L
1807 zol10 = zol(i)*(10.+zntstoch_ice(i))/za(i) ! (10+z0)/L
1808 zol2 = zol(i)*(2.+zntstoch_ice(i))/za(i) ! (2+z0)/L
1809
1810 !COMPUTE PSIM and PSIH
1811 !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
1812 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
1813 !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I))
1814 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_ice(I),ZNTstoch_ice(I),ZA(I))
1815 !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0)
1816 psim(i)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt)
1817 psih(i)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt)
1818 psim10(i)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt)
1819 psih10(i)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt)
1820 psih2(i)=psih_stable(zol2,psi_opt)-psih_stable(zolzt,psi_opt)
1821
1822 ! 1.0 over Monin-Obukhov length
1823 rmol(i)= zol(i)/za(i)
1824
1825 ELSEIF(rb_ice(i) .EQ. 0.) THEN
1826 !=========================================================
1827 !-----CLASS 3; FORCED CONVECTION/NEUTRAL:
1828 !=========================================================
1829
1830 psim(i)=0.0
1831 psih(i)=psim(i)
1832 psim10(i)=0.
1833 psih10(i)=0.
1834 psih2(i)=0.
1835
1836 zol(i) =0.
1837 rmol(i) =0.
1838
1839 ELSEIF(rb_ice(i) .LT. 0.)THEN
1840 !==========================================================
1841 !-----CLASS 4; FREE CONVECTION:
1842 !==========================================================
1843
1844 IF (.not. flag_restart .or. (flag_restart .and. itimestep > 1) ) THEN
1845 !COMPUTE z/L first guess:
1846 CALL li_etal_2010(zol(i),rb_ice(i),za(i)/zntstoch_ice(i),zratio_ice(i))
1847 !ZOL(I)=ZA(I)*KARMAN*grav*MOL(I)/(TH1D(I)*MAX(UST_ice(I)*UST_ice(I),0.001))
1848 zol(i)=max(zol(i),-20.0_kind_phys)
1849 zol(i)=min(zol(i),0.0_kind_phys)
1850
1851 IF (debug_code >= 1) THEN
1852 IF (zntstoch_ice(i) < 1e-8 .OR. zt_ice(i) < 1e-10) THEN
1853 write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i
1854 write(0,*)"rb=", rb_ice(i)," ZNT=", zntstoch_ice(i)," ZT=",zt_ice(i)
1855 write(0,*)" tsk=", tskin_ice(i)," wstar=",wstar(i)," prev z/L=",zol(i),&
1856 " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),&
1857 " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",psfcpa(i), &
1858 " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i)
1859 ENDIF
1860 ENDIF
1861
1862 !Use Pedros iterative function to find z/L
1863 !zol(I)=zolri(rb_ice(I),ZA(I),ZNTstoch_ice(I),ZT_ice(I),ZOL(I),psi_opt)
1864 !Use brute-force method
1865 zol(i)=zolrib(rb_ice(i),za(i),zntstoch_ice(i),zt_ice(i),gz1oz0_ice(i),gz1ozt_ice(i),zol(i),psi_opt)
1866 ENDIF ! restart
1867 zol(i)=max(zol(i),-20.0_kind_phys)
1868 zol(i)=min(zol(i),0.0_kind_phys)
1869
1870 zolzt = zol(i)*zt_ice(i)/za(i) ! zt/L
1871 zolz0 = zol(i)*zntstoch_ice(i)/za(i) ! z0/L
1872 zolza = zol(i)*(za(i)+zntstoch_ice(i))/za(i) ! (z+z0/L
1873 zol10 = zol(i)*(10.+zntstoch_ice(i))/za(i) ! (10+z0)/L
1874 zol2 = zol(i)*(2.+zntstoch_ice(i))/za(i) ! (2+z0)/L
1875
1876 !COMPUTE PSIM and PSIH
1877 !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), ZT_ice(I), ZNTstoch_ice(I), ZA(I))
1878 !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
1879 !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_ice(I),ZNTstoch_ice(I),ZA(I))
1880 ! use tables
1881 psim(i)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt)
1882 psih(i)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt)
1883 psim10(i)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt)
1884 psih10(i)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt)
1885 psih2(i)=psih_unstable(zol2,psi_opt)-psih_unstable(zolzt,psi_opt)
1886
1887 !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND
1888 !---HIGH ROUGHNESS. THIS PREVENTS DENOMINATOR IN FLUXES
1889 !---FROM GETTING TOO SMALL
1890 psih(i)=min(psih(i),0.9*gz1ozt_ice(i))
1891 psim(i)=min(psim(i),0.9*gz1oz0_ice(i))
1892 psih2(i)=min(psih2(i),0.9*gz2ozt_ice(i))
1893 psim10(i)=min(psim10(i),0.9*gz10oz0_ice(i))
1894 psih10(i)=min(psih10(i),0.9*gz10ozt_ice(i))
1895
1896 rmol(i) = zol(i)/za(i)
1897
1898 ENDIF
1899
1900 ! CALCULATE THE RESISTANCE:
1901 psix_ice(i) =max(gz1oz0_ice(i)-psim(i) , 1.0_kind_phys)
1902 psix10_ice(i)=max(gz10oz0_ice(i)-psim10(i), 1.0_kind_phys)
1903 psit_ice(i) =max(gz1ozt_ice(i)-psih(i) , 1.0_kind_phys)
1904 psit2_ice(i) =max(gz2ozt_ice(i)-psih2(i) , 1.0_kind_phys)
1905 psiq_ice(i) =max(log((za(i)+zq_ice(i))/zq_ice(i))-psih(i) ,1.0_kind_phys)
1906 psiq2_ice(i) =max(log((2.0+zq_ice(i))/zq_ice(i))-psih2(i) ,1.0_kind_phys)
1907
1908 ENDIF ! end ice points
1909
1910 !------------------------------------------------------------
1911 !-----COMPUTE THE FRICTIONAL VELOCITY:
1912 !------------------------------------------------------------
1913
1914 IF (wet(i)) THEN
1915 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
1916 oldust = ust_wat(i)
1917 ust_wat(i)=0.5*ust_wat(i)+0.5*karman*wspd(i)/psix_wat(i)
1918 !NON-AVERAGED:
1919 !UST_wat(I)=KARMAN*WSPD(I)/PSIX_wat(I)
1920 stress_wat(i)=ust_wat(i)**2
1921
1922 ! Compute u* without vconv for use in HFX calc when isftcflx > 0
1923 wspdi(i)=max(sqrt(u1d(i)*u1d(i)+v1d(i)*v1d(i)), wmin)
1924 !tgs - should USTM be separater for dry, icy, wet?
1925 ustm(i)=0.5*ustm(i)+0.5*karman*wspdi(i)/psix_wat(i)
1926
1927 ! for possible future changes in sea-ice fraction from 0 to >0:
1928 if (.not. icy(i)) ust_ice(i)=ust_wat(i)
1929 ENDIF ! end water points
1930
1931 IF (dry(i)) THEN
1932 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
1933 oldust = ust_lnd(i)
1934 ust_lnd(i)=0.5*ust_lnd(i)+0.5*karman*wspd(i)/psix_lnd(i)
1935 !NON-AVERAGED:
1936 !UST_lnd(I)=KARMAN*WSPD(I)/PSIX_lnd(I)
1937 !From Tilden Meyers:
1938 !IF (rb_lnd(I) .GE. 0.0) THEN
1939 ! ust_lnd(i)=WSPD_lnd*0.1/(1.0 + 10.0*rb_lnd(I))
1940 !ELSE
1941 ! ust_lnd(i)=WSPD_lnd*0.1*(1.0 - 10.0*rb_lnd(I))**onethird
1942 !ENDIF
1943 ust_lnd(i)=max(ust_lnd(i),0.005_kind_phys)
1944 stress_lnd(i)=ust_lnd(i)**2
1945
1946 !set ustm = ust over land.
1947 !tgs - should USTM be separater for dry, icy, wet?
1948 ustm(i)=ust_lnd(i)
1949 ENDIF ! end water points
1950
1951 IF (icy(i)) THEN
1952 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
1953 oldust = ust_ice(i)
1954 ust_ice(i)=0.5*ust_ice(i)+0.5*karman*wspd(i)/psix_ice(i)
1955 !NON-AVERAGED:
1956 !UST_ice(I)=KARMAN*WSPD(I)/PSIX_ice(I)
1957 ust_ice(i)=max(ust_ice(i),0.005_kind_phys)
1958 stress_ice(i)=ust_ice(i)**2
1959
1960 !Set ustm = ust over ice.
1961 !tgs - should USTM be separate for for dry, icy, wet?
1962 ustm(i)=ust_ice(i)
1963
1964 ! for possible future changes in sea-ice fraction from 1 to <1:
1965 !tgs - sea ice can be <1 now
1966 if (.not. wet(i)) ust_wat(i)=ust_ice(i)
1967 ENDIF ! end ice points
1968
1969 !----------------------------------------------------
1970 !----COMPUTE THE TEMPERATURE SCALE (a.k.a. FRICTION TEMPERATURE, T*, or MOL)
1971 !----AND COMPUTE THE MOISTURE SCALE (or q*)
1972 !----------------------------------------------------
1973
1974 !tgs - should we have MOL and qstar separate for dry, icy and wet?
1975 IF (wet(i)) THEN
1976 dtg=thv1d(i)-thvsk_wat(i)
1977 oldtst=mol(i)
1978 mol(i)=karman*dtg/psit_wat(i)/prt
1979 !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I))
1980 !t_star(I) = MOL(I)
1981 !----------------------------------------------------
1982 dqg=(qvsh(i)-qsfc_wat(i))*1000. !(kg/kg -> g/kg)
1983 qstar(i)=karman*dqg/psiq_wat(i)/prt
1984 ENDIF
1985
1986 IF (dry(i)) THEN
1987 dtg=thv1d(i)-thvsk_lnd(i)
1988 oldtst=mol(i)
1989 mol(i)=karman*dtg/psit_lnd(i)/prt
1990 !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I))
1991 !t_star(I) = MOL(I)
1992 !----------------------------------------------------
1993 dqg=(qvsh(i)-qsfc_lnd(i))*1000. !(kg/kg -> g/kg)
1994 qstar(i)=karman*dqg/psiq_lnd(i)/prt
1995 ENDIF
1996
1997 IF (icy(i)) THEN
1998 dtg=thv1d(i)-thvsk_ice(i)
1999 oldtst=mol(i)
2000 mol(i)=karman*dtg/psit_ice(i)/prt
2001 !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I))
2002 !t_star(I) = MOL(I)
2003 !----------------------------------------------------
2004 dqg=(qvsh(i)-qsfc_ice(i))*1000. !(kg/kg -> g/kg)
2005 qstar(i)=karman*dqg/psiq_ice(i)/prt
2006 ENDIF
2007
2008 endif ! flag_iter
2009 ENDDO ! end i-loop
2010
2011#ifdef _OPENACC
2012! Necessary since OpenACC does not support branching in parallel code.
2013! Must sync host errflg, errmsg to determine if return must be triggered
2014! and correct error message and error flag code returned.
2015! This code is being executed on the HOST side only, pulling data from DEVICE.
2016!$acc exit data copyout(device_special_errflg, device_special_errmsg)
2017 IF (device_special_errflg /= 0) THEN
2018 errflg = device_special_errflg
2019 errmsg = device_special_errmsg
2020 return
2021 ENDIF
2022#endif
2023
2024!$acc serial present(wet, dry, icy, &
2025!$acc PSIM, PSIH, CPM, RHO1D, ZOL, wspd, MOL, &
2026!$acc wstar, qstar, THV1D, HFX, MAVAIL, QVSH, &
2027!$acc THVSK_wat, THVSK_lnd, THVSK_ice, &
2028!$acc UST_wat, UST_lnd, UST_ice, &
2029!$acc ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
2030!$acc zt_wat, zt_lnd, zt_ice)
2031 IF (debug_code == 2) THEN
2032 DO i=its,ite
2033 IF(wet(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(wet)"
2034 IF(dry(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(land)"
2035 IF(icy(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(ice)"
2036 write(*,*)"z/L:",zol(i)," wspd:",wspd(i)," Tstar:",mol(i)
2037 IF(wet(i))write(*,*)"PSIM:",psim(i)," PSIH:",psih(i)," W*:",wstar(i),&
2038 " DTHV:",thv1d(i)-thvsk_wat(i)
2039 IF(dry(i))write(*,*)"PSIM:",psim(i)," PSIH:",psih(i)," W*:",wstar(i),&
2040 " DTHV:",thv1d(i)-thvsk_lnd(i)
2041 IF(icy(i))write(*,*)"PSIM:",psim(i)," PSIH:",psih(i)," W*:",wstar(i),&
2042 " DTHV:",thv1d(i)-thvsk_ice(i)
2043 write(*,*)"CPM:",cpm(i)," RHO1D:",rho1d(i)," q*:",qstar(i)," T*:",mol(i)
2044 IF(wet(i))write(*,*)"U*:",ust_wat(i)," Z0:",zntstoch_wat(i)," Zt:",zt_wat(i)
2045 IF(dry(i))write(*,*)"U*:",ust_lnd(i)," Z0:",zntstoch_lnd(i)," Zt:",zt_lnd(i)
2046 IF(icy(i))write(*,*)"U*:",ust_ice(i)," Z0:",zntstoch_ice(i)," Zt:",zt_ice(i)
2047 write(*,*)"hfx:",hfx(i)," MAVAIL:",mavail(i)," QVSH(I):",qvsh(i)
2048 write(*,*)"============================================="
2049 ENDDO ! end i-loop
2050 ENDIF
2051!$acc end serial
2052
2053 !----------------------------------------------------------
2054 ! COMPUTE SURFACE HEAT AND MOISTURE FLUXES
2055 !----------------------------------------------------------
2056!$acc parallel loop present(flag_iter, dry, wet, icy, &
2057!$acc QFX, HFX, FLHC, FLQC, LH, CHS, CH, CHS2, CQS2, &
2058!$acc RHO1D, MAVAIL, USTM, &
2059!$acc UST_lnd, UST_wat, UST_ice, &
2060!$acc PSIQ_lnd, PSIT_lnd, PSIX_lnd, &
2061!$acc PSIQ_wat, PSIT_wat, PSIX_wat, &
2062!$acc PSIQ_ice, PSIT_ice, PSIX_ice, &
2063!$acc PSIQ2_lnd, PSIT2_lnd, &
2064!$acc PSIQ2_wat, PSIT2_wat, &
2065!$acc PSIQ2_ice, PSIT2_ice, &
2066!$acc QSFC, QSFC_lnd, QSFC_wat, QSFC_ice, &
2067!$acc QFLX, QFLX_lnd, QFLX_wat, QFLX_ice, &
2068!$acc HFLX, HFLX_lnd, HFLX_wat, HFLX_ice, &
2069!$acc QSFCMR_lnd, QSFCMR_wat, QSFCMR_ice, &
2070!$acc QV1D, WSPD, WSPDI, CPM, TH1D, &
2071!$acc THSK_lnd, THSK_wat, THSK_ice, &
2072!$acc ch_lnd, ch_wat, ch_ice, &
2073!$acc cm_lnd, cm_wat, cm_ice)
2074 DO i=its,ite
2075 if( flag_iter(i) ) then
2076
2077 IF (isfflx .LT. 1) THEN
2078
2079 qfx(i) = 0.
2080 hfx(i) = 0.
2081 hflx(i) = 0.
2082 flhc(i) = 0.
2083 flqc(i) = 0.
2084 lh(i) = 0.
2085 chs(i) = 0.
2086 ch(i) = 0.
2087 chs2(i) = 0.
2088 cqs2(i) = 0.
2089 ch_wat(i)= 0.
2090 cm_wat(i)= 0.
2091 ch_lnd(i)= 0.
2092 cm_lnd(i)= 0.
2093 ch_ice(i)= 0.
2094 cm_ice(i)= 0.
2095
2096 ELSE
2097
2098 IF (dry(i)) THEN
2099
2100 !------------------------------------------
2101 ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC)
2102 ! AND MOISTURE (FLQC)
2103 !------------------------------------------
2104 !tgs - should FLQC, FLHC be separate for dry, icy and wet?
2105 flqc(i)=rho1d(i)*mavail(i)*ust_lnd(i)*karman/psiq_lnd(i)
2106 flhc(i)=rho1d(i)*cpm(i)*ust_lnd(i)*karman/psit_lnd(i)
2107
2108 IF (compute_flux) THEN
2109 !----------------------------------
2110 ! COMPUTE SURFACE MOISTURE FLUX:
2111 !----------------------------------
2112 qfx(i)=flqc(i)*(qsfcmr_lnd(i)-qv1d(i))
2113 !QFX(I)=FLQC(I)*(QSFC_lnd(I)-QV1D(I))
2114 qfx(i)=max(qfx(i),-0.02_kind_phys) !allows small neg QFX
2115 lh(i)=xlv*qfx(i)
2116 ! BWG, 2020-06-17: Mod next 2 lines for fractional
2117 qflx_lnd(i)=qfx(i)/rho1d(i)
2118 qflx(i)=qflx_lnd(i)
2119
2120 !----------------------------------
2121 ! COMPUTE SURFACE HEAT FLUX:
2122 !----------------------------------
2123 !HFX(I)=FLHC(I)*(THSK_lnd(I)-TH1D(I))
2124 hfx(i)=rho1d(i)*cpm(i)*karman*wspd(i)/psix_lnd(i)*karman/psit_lnd(i)*(thsk_lnd(i)-th1d(i))
2125 hfx(i)=max(hfx(i),-250._kind_phys)
2126 ! BWG, 2020-06-17: Mod next 2 lines for fractional
2127 hflx_lnd(i)=hfx(i)/(rho1d(i)*cpm(i))
2128 hflx(i)=hflx_lnd(i)
2129 ENDIF
2130
2131 !TRANSFER COEFF FOR SOME LSMs:
2132 !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
2133 ! /XKA+ZA(I)/ZL)-PSIH(I))
2134
2135 !tgs - should QSFC, CHS, CHS2 and CQS2 be separate for dry, icy and wet?
2136 chs(i)=ust_lnd(i)*karman/psit_lnd(i)
2137
2138 !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY
2139 cqs2(i)=ust_lnd(i)*karman/psiq2_lnd(i)
2140 chs2(i)=ust_lnd(i)*karman/psit2_lnd(i)
2141
2142 qsfc(i)=qsfc_lnd(i)
2143
2144 ELSEIF (wet(i)) THEN
2145
2146 !------------------------------------------
2147 ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC)
2148 ! AND MOISTURE (FLQC)
2149 !------------------------------------------
2150 flqc(i)=rho1d(i)*mavail(i)*ust_wat(i)*karman/psiq_wat(i)
2151 flhc(i)=rho1d(i)*cpm(i)*ust_wat(i)*karman/psit_wat(i)
2152
2153 IF (compute_flux) THEN
2154 !----------------------------------
2155 ! COMPUTE SURFACE MOISTURE FLUX:
2156 !----------------------------------
2157 qfx(i)=flqc(i)*(qsfcmr_wat(i)-qv1d(i))
2158 !QFX(I)=FLQC(I)*(QSFC_wat(I)-QV1D(I))
2159 qfx(i)=max(qfx(i),-0.02_kind_phys) !allows small neg QFX
2160 lh(i)=xlv*qfx(i)
2161 ! BWG, 2020-06-17: Mod next 2 lines for fractional
2162 qflx_wat(i)=qfx(i)/rho1d(i)
2163 qflx(i)=qflx_wat(i)
2164
2165 !----------------------------------
2166 ! COMPUTE SURFACE HEAT FLUX:
2167 !----------------------------------
2168 !HFX(I)=FLHC(I)*(THSK_wat(I)-TH1D(I))
2169 hfx(i)=rho1d(i)*cpm(i)*karman*wspd(i)/psix_wat(i)*karman/psit_wat(i)*(thsk_wat(i)-th1d(i))
2170 IF ( PRESENT(isftcflx) ) THEN
2171 IF ( isftcflx.NE.0 ) THEN
2172 ! AHW: add dissipative heating term
2173 hfx(i)=hfx(i)+rho1d(i)*ustm(i)*ustm(i)*wspdi(i)
2174 ENDIF
2175 ENDIF
2176 ! BWG, 2020-06-17: Mod next 2 lines for fractional
2177 hflx_wat(i)=hfx(i)/(rho1d(i)*cpm(i))
2178 hflx(i)=hflx_wat(i)
2179 ENDIF
2180
2181 !TRANSFER COEFF FOR SOME LSMs:
2182 !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
2183 ! /XKA+ZA(I)/ZL)-PSIH(I))
2184 chs(i)=ust_wat(i)*karman/psit_wat(i)
2185
2186 !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY
2187 cqs2(i)=ust_wat(i)*karman/psiq2_wat(i)
2188 chs2(i)=ust_wat(i)*karman/psit2_wat(i)
2189
2190 qsfc(i)=qsfc_wat(i)
2191
2192 ELSEIF (icy(i)) THEN
2193
2194 !------------------------------------------
2195 ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC)
2196 ! AND MOISTURE (FLQC)
2197 !------------------------------------------
2198 flqc(i)=rho1d(i)*mavail(i)*ust_ice(i)*karman/psiq_ice(i)
2199 flhc(i)=rho1d(i)*cpm(i)*ust_ice(i)*karman/psit_ice(i)
2200
2201 IF (compute_flux) THEN
2202 !----------------------------------
2203 ! COMPUTE SURFACE MOISTURE FLUX:
2204 !----------------------------------
2205 qfx(i)=flqc(i)*(qsfcmr_ice(i)-qv1d(i))
2206 !QFX(I)=FLQC(I)*(QSFC_ice(I)-QV1D(I))
2207 qfx(i)=max(qfx(i),-0.02_kind_phys) !allows small neg QFX
2208 lh(i)=xlf*qfx(i)
2209 ! BWG, 2020-06-17: Mod next 2 lines for fractional
2210 qflx_ice(i)=qfx(i)/rho1d(i)
2211 qflx(i)=qflx_ice(i)
2212
2213 !----------------------------------
2214 ! COMPUTE SURFACE HEAT FLUX:
2215 !----------------------------------
2216 !HFX(I)=FLHC(I)*(THSK_ice(I)-TH1D(I))
2217 hfx(i)=rho1d(i)*cpm(i)*karman*wspd(i)/psix_ice(i)*karman/psit_ice(i)*(thsk_ice(i)-th1d(i))
2218 hfx(i)=max(hfx(i),-250._kind_phys)
2219 ! BWG, 2020-06-17: Mod next 2 lines for fractional
2220 hflx_ice(i)=hfx(i)/(rho1d(i)*cpm(i))
2221 hflx(i)=hflx_ice(i)
2222 ENDIF
2223
2224 !TRANSFER COEFF FOR SOME LSMs:
2225 !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
2226 ! /XKA+ZA(I)/ZL)-PSIH(I))
2227 chs(i)=ust_ice(i)*karman/psit_ice(i)
2228
2229 !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY
2230 cqs2(i)=ust_ice(i)*karman/psiq2_ice(i)
2231 chs2(i)=ust_ice(i)*karman/psit2_ice(i)
2232
2233 qsfc(i)=qsfc_ice(i)
2234
2235 ENDIF
2236
2237 IF (debug_code > 1) THEN
2238 write(*,*)"QFX=",qfx(i),"FLQC=",flqc(i)
2239 if(icy(i))write(*,*)"ice, MAVAIL:",mavail(i)," u*=",ust_ice(i)," psiq=",psiq_ice(i)
2240 if(dry(i))write(*,*)"lnd, MAVAIL:",mavail(i)," u*=",ust_lnd(i)," psiq=",psiq_lnd(i)
2241 if(wet(i))write(*,*)"ocn, MAVAIL:",mavail(i)," u*=",ust_wat(i)," psiq=",psiq_wat(i)
2242 ENDIF
2243
2244 ! The exchange coefficient for cloud water is assumed to be the
2245 ! same as that for heat. CH is multiplied by WSPD.
2246 ch(i)=flhc(i)/( cpm(i)*rho1d(i) )
2247
2248 !-----------------------------------------
2249 !--- COMPUTE EXCHANGE COEFFICIENTS FOR FV3
2250 !-----------------------------------------
2251 IF (wet(i)) THEN
2252 ch_wat(i)=(karman/psix_wat(i))*(karman/psit_wat(i))
2253 cm_wat(i)=(karman/psix_wat(i))*(karman/psix_wat(i))
2254 ENDIF
2255 IF (dry(i)) THEN
2256 ch_lnd(i)=(karman/psix_lnd(i))*(karman/psit_lnd(i))
2257 cm_lnd(i)=(karman/psix_lnd(i))*(karman/psix_lnd(i))
2258 ENDIF
2259 IF (icy(i)) THEN
2260 ch_ice(i)=(karman/psix_ice(i))*(karman/psit_ice(i))
2261 cm_ice(i)=(karman/psix_ice(i))*(karman/psix_ice(i))
2262 ENDIF
2263
2264 ENDIF !end ISFFLX option
2265 endif ! flag_iter
2266ENDDO ! end i-loop
2267
2268IF (compute_diag) then
2269 !$acc parallel loop present(flag_iter, dry, wet, icy, &
2270 !$acc ZA, ZA2, T2, TH2, TH1D, Q2, QV1D, PSFCPA, &
2271 !$acc THSK_lnd, THSK_wat, THSK_ice, &
2272 !$acc QSFC_lnd, QSFC_wat, QSFC_ice, &
2273 !$acc U10, V10, U1D, V1D, U1D2, V1D2, &
2274 !$acc ZNTstoch_lnd, ZNTstoch_lnd, ZNTstoch_ice, &
2275 !$acc PSIX_lnd, PSIX_wat, PSIX_ice, &
2276 !$acc PSIX10_lnd, PSIX10_wat, PSIX10_ice, &
2277 !$acc PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
2278 !$acc PSIT_lnd, PSIT_wat, PSIT_ice, &
2279 !$acc PSIQ2_lnd, PSIQ2_wat, PSIQ2_ice, &
2280 !$acc PSIQ_lnd, PSIQ_wat, PSIQ_ice)
2281 DO i=its,ite
2282 if( flag_iter(i) ) then
2283 !-----------------------------------------------------
2284 !COMPUTE DIAGNOSTICS
2285 !-----------------------------------------------------
2286 !COMPUTE 10 M WNDS
2287 !-----------------------------------------------------
2288 ! If the lowest model level is close to 10-m, use it
2289 ! instead of the flux-based diagnostic formula.
2290 if (za(i) .le. 7.0) then
2291 ! high vertical resolution
2292 if(za2(i) .gt. 7.0 .and. za2(i) .lt. 13.0) then
2293 !use 2nd model level
2294 u10(i)=u1d2(i)
2295 v10(i)=v1d2(i)
2296 else
2297 IF (dry(i)) THEN
2298 !U10(I)=U1D(I)*PSIX10_lnd(I)/PSIX_lnd(I)
2299 !V10(I)=V1D(I)*PSIX10_lnd(I)/PSIX_lnd(I)
2300 !use neutral-log:
2301 u10(i)=u1d(i)*log(10./zntstoch_lnd(i))/log(za(i)/zntstoch_lnd(i))
2302 v10(i)=v1d(i)*log(10./zntstoch_lnd(i))/log(za(i)/zntstoch_lnd(i))
2303 ELSEIF (wet(i)) THEN
2304 u10(i)=u1d(i)*log(10./zntstoch_wat(i))/log(za(i)/zntstoch_wat(i))
2305 v10(i)=v1d(i)*log(10./zntstoch_wat(i))/log(za(i)/zntstoch_wat(i))
2306 ELSEIF (icy(i)) THEN
2307 u10(i)=u1d(i)*log(10./zntstoch_ice(i))/log(za(i)/zntstoch_ice(i))
2308 v10(i)=v1d(i)*log(10./zntstoch_ice(i))/log(za(i)/zntstoch_ice(i))
2309 ENDIF
2310 endif
2311 elseif (za(i) .gt. 7.0 .and. za(i) .lt. 13.0) then
2312 !moderate vertical resolution
2313 IF (dry(i)) THEN
2314 !U10(I)=U1D(I)*PSIX10_lnd(I)/PSIX_lnd(I)
2315 !V10(I)=V1D(I)*PSIX10_lnd(I)/PSIX_lnd(I)
2316 !use neutral-log:
2317 u10(i)=u1d(i)*log(10./zntstoch_lnd(i))/log(za(i)/zntstoch_lnd(i))
2318 v10(i)=v1d(i)*log(10./zntstoch_lnd(i))/log(za(i)/zntstoch_lnd(i))
2319 ELSEIF (wet(i)) THEN
2320 u10(i)=u1d(i)*log(10./zntstoch_wat(i))/log(za(i)/zntstoch_wat(i))
2321 v10(i)=v1d(i)*log(10./zntstoch_wat(i))/log(za(i)/zntstoch_wat(i))
2322 ELSEIF (icy(i)) THEN
2323 u10(i)=u1d(i)*log(10./zntstoch_ice(i))/log(za(i)/zntstoch_ice(i))
2324 v10(i)=v1d(i)*log(10./zntstoch_ice(i))/log(za(i)/zntstoch_ice(i))
2325 ENDIF
2326 else
2327 ! very coarse vertical resolution
2328 IF (dry(i)) THEN
2329 u10(i)=u1d(i)*psix10_lnd(i)/psix_lnd(i)
2330 v10(i)=v1d(i)*psix10_lnd(i)/psix_lnd(i)
2331 ELSEIF (wet(i)) THEN
2332 u10(i)=u1d(i)*psix10_wat(i)/psix_wat(i)
2333 v10(i)=v1d(i)*psix10_wat(i)/psix_wat(i)
2334 ELSEIF (icy(i)) THEN
2335 u10(i)=u1d(i)*psix10_ice(i)/psix_ice(i)
2336 v10(i)=v1d(i)*psix10_ice(i)/psix_ice(i)
2337 ENDIF
2338 endif
2339
2340 !-----------------------------------------------------
2341 !COMPUTE 2m T, TH, AND Q
2342 !THESE WILL BE OVERWRITTEN FOR LAND POINTS IN THE LSM
2343 !-----------------------------------------------------
2344 IF (dry(i)) THEN
2345 dtg=th1d(i)-thsk_lnd(i)
2346 th2(i)=thsk_lnd(i)+dtg*psit2_lnd(i)/psit_lnd(i)
2347 !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY
2348 !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL.
2349 IF ((th1d(i)>thsk_lnd(i) .AND. (th2(i)<thsk_lnd(i) .OR. th2(i)>th1d(i))) .OR. &
2350 (th1d(i)<thsk_lnd(i) .AND. (th2(i)>thsk_lnd(i) .OR. th2(i)<th1d(i)))) THEN
2351 th2(i)=thsk_lnd(i) + 2.*(th1d(i)-thsk_lnd(i))/za(i)
2352 ENDIF
2353 t2(i)=th2(i)*(psfcpa(i)/100000.)**rovcp
2354
2355 q2(i)=qsfc_lnd(i)+(qv1d(i)-qsfc_lnd(i))*psiq2_lnd(i)/psiq_lnd(i)
2356 q2(i)= max(q2(i), min(qsfc_lnd(i), qv1d(i)))
2357 q2(i)= min(q2(i), 1.05*qv1d(i))
2358 ELSEIF (wet(i)) THEN
2359 dtg=th1d(i)-thsk_wat(i)
2360 th2(i)=thsk_wat(i)+dtg*psit2_wat(i)/psit_wat(i)
2361 !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY
2362 !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL.
2363 IF ((th1d(i)>thsk_wat(i) .AND. (th2(i)<thsk_wat(i) .OR. th2(i)>th1d(i))) .OR. &
2364 (th1d(i)<thsk_wat(i) .AND. (th2(i)>thsk_wat(i) .OR. th2(i)<th1d(i)))) THEN
2365 th2(i)=thsk_wat(i) + 2.*(th1d(i)-thsk_wat(i))/za(i)
2366 ENDIF
2367 t2(i)=th2(i)*(psfcpa(i)/100000.)**rovcp
2368
2369 q2(i)=qsfc_wat(i)+(qv1d(i)-qsfc_wat(i))*psiq2_wat(i)/psiq_wat(i)
2370 q2(i)= max(q2(i), min(qsfc_wat(i), qv1d(i)))
2371 q2(i)= min(q2(i), 1.05*qv1d(i))
2372 ELSEIF (icy(i)) THEN
2373 dtg=th1d(i)-thsk_ice(i)
2374 th2(i)=thsk_ice(i)+dtg*psit2_ice(i)/psit_ice(i)
2375 !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY
2376 !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL.
2377 IF ((th1d(i)>thsk_ice(i) .AND. (th2(i)<thsk_ice(i) .OR. th2(i)>th1d(i))) .OR. &
2378 (th1d(i)<thsk_ice(i) .AND. (th2(i)>thsk_ice(i) .OR. th2(i)<th1d(i)))) THEN
2379 th2(i)=thsk_ice(i) + 2.*(th1d(i)-thsk_ice(i))/za(i)
2380 ENDIF
2381 t2(i)=th2(i)*(psfcpa(i)/100000.)**rovcp
2382
2383 q2(i)=qsfc_ice(i)+(qv1d(i)-qsfc_ice(i))*psiq2_ice(i)/psiq_ice(i)
2384 q2(i)= max(q2(i), min(qsfc_ice(i), qv1d(i)))
2385 q2(i)= min(q2(i), 1.05*qv1d(i))
2386 ENDIF
2387 endif ! flag_iter
2388 ENDDO
2389ENDIF ! end compute_diag
2390
2391!-----------------------------------------------------
2392! DEBUG - SUSPICIOUS VALUES
2393!-----------------------------------------------------
2394!$acc serial present(dry, wet, icy, CPM, MAVAIL, &
2395!$acc HFX, LH, wstar, RHO1D, PBLH, ZOL, ZA, MOL, &
2396!$acc PSIM, PSIH, WSTAR, T1D, TH1D, THV1D, QVSH, &
2397!$acc UST_wat, UST_lnd, UST_ice, &
2398!$acc THSK_wat, THSK_lnd, THSK_ice, &
2399!$acc THVSK_wat, THVSK_lnd, THVSK_ice, &
2400!$acc ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
2401!$acc ZT_wat, ZT_lnd, ZT_ice, &
2402!$acc QSFC_wat, QSFC_lnd, QSFC_ice, &
2403!$acc PSIX_wat, PSIX_lnd, PSIX_ice)
2404IF ( debug_code == 2) THEN
2405 DO i=its,ite
2406 yesno = 0
2407 IF (compute_flux) THEN
2408 IF (hfx(i) > 1200. .OR. hfx(i) < -700.)THEN
2409 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2410 i,j, "HFX: ",hfx(i)
2411 yesno = 1
2412 ENDIF
2413 IF (lh(i) > 1200. .OR. lh(i) < -700.)THEN
2414 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2415 i,j, "LH: ",lh(i)
2416 yesno = 1
2417 ENDIF
2418 ENDIF
2419 IF (wet(i)) THEN
2420 IF (ust_wat(i) < 0.0 .OR. ust_wat(i) > 4.0 )THEN
2421 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2422 i,j, "UST_wat: ",ust_wat(i)
2423 yesno = 1
2424 ENDIF
2425 ENDIF
2426 IF (dry(i)) THEN
2427 IF (ust_lnd(i) < 0.0 .OR. ust_lnd(i) > 4.0 )THEN
2428 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2429 i,j, "UST_lnd: ",ust_lnd(i)
2430 yesno = 1
2431 ENDIF
2432 ENDIF
2433 IF (icy(i)) THEN
2434 IF (ust_ice(i) < 0.0 .OR. ust_ice(i) > 4.0 )THEN
2435 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2436 i,j, "UST_ice: ",ust_ice(i)
2437 yesno = 1
2438 ENDIF
2439 ENDIF
2440 IF (wstar(i)<0.0 .OR. wstar(i) > 6.0)THEN
2441 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2442 i,j, "WSTAR: ",wstar(i)
2443 yesno = 1
2444 ENDIF
2445 IF (rho1d(i)<0.0 .OR. rho1d(i) > 1.6 )THEN
2446 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2447 i,j, "rho: ",rho1d(i)
2448 yesno = 1
2449 ENDIF
2450 IF (dry(i)) THEN
2451 IF (qsfc_lnd(i)*1000. <0.0 .OR. qsfc_lnd(i)*1000. >40.)THEN
2452 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2453 i,j, "QSFC_lnd: ",qsfc_lnd(i)
2454 yesno = 1
2455 ENDIF
2456 ENDIF
2457 IF (pblh(i)<0. .OR. pblh(i)>6000.)THEN
2458 print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",&
2459 i,j, "PBLH: ",pblh(i)
2460 yesno = 1
2461 ENDIF
2462
2463 IF (yesno == 1) THEN
2464 IF (wet(i)) THEN
2465 print*," OTHER INFO over water:"
2466 print*,"z/L:",zol(i)," U*:",ust_wat(i)," Tstar:",mol(i)
2467 print*,"PSIM:",psim(i)," PSIH:",psih(i)," W*:",wstar(i),&
2468 " DTHV:",thv1d(i)-thvsk_wat(i)
2469 print*,"CPM:",cpm(i)," RHO1D:",rho1d(i)," L:",&
2470 zol(i)/za(i)," DTH:",th1d(i)-thsk_wat(i)
2471 print*," Z0:",zntstoch_wat(i)," Zt:",zt_wat(i)," za:",za(i)
2472 print*,"MAVAIL:",mavail(i)," QSFC_wat(I):",&
2473 qsfc_wat(i)," QVSH(I):",qvsh(i)
2474 print*,"PSIX=",psix_wat(i)," T1D(i):",t1d(i)
2475 write(*,*)"============================================="
2476 ENDIF
2477 IF (dry(i)) THEN
2478 print*," OTHER INFO over land:"
2479 print*,"z/L:",zol(i)," U*:",ust_lnd(i),&
2480 " Tstar:",mol(i)
2481 print*,"PSIM:",psim(i)," PSIH:",psih(i)," W*:",wstar(i),&
2482 " DTHV:",thv1d(i)-thvsk_lnd(i)
2483 print*,"CPM:",cpm(i)," RHO1D:",rho1d(i)," L:",&
2484 zol(i)/za(i)," DTH:",th1d(i)-thsk_lnd(i)
2485 print*," Z0:",zntstoch_lnd(i)," Zt:",zt_lnd(i)," za:",za(i)
2486 print*," MAVAIL:",mavail(i)," QSFC_lnd(I):",&
2487 qsfc_lnd(i)," QVSH(I):",qvsh(i)
2488 print*,"PSIX=",psix_lnd(i)," T1D(i):",t1d(i)
2489 write(*,*)"============================================="
2490 ENDIF
2491 IF (icy(i)) THEN
2492 print*," OTHER INFO:"
2493 print*,"z/L:",zol(i)," U*:",ust_ice(i),&
2494 " Tstar:",mol(i)
2495 print*,"PSIM:",psim(i)," PSIH:",psih(i)," W*:",wstar(i),&
2496 " DTHV:",thv1d(i)-thvsk_ice(i)
2497 print*,"CPM:",cpm(i)," RHO1D:",rho1d(i)," L:",&
2498 zol(i)/za(i)," DTH:",th1d(i)-thsk_ice(i)
2499 print*," Z0:",zntstoch_ice(i)," Zt:",zt_ice(i)," za:",za(i)
2500 print*," MAVAIL:",mavail(i)," QSFC_ice(I):",&
2501 qsfc_ice(i)," QVSH(I):",qvsh(i)
2502 print*,"PSIX=",psix_ice(i)," T1D(i):",t1d(i)
2503 write(*,*)"============================================="
2504 ENDIF
2505 ENDIF
2506 ENDDO ! end i-loop
2507 ENDIF ! end debug option
2508!$acc end serial
2509
2510!$acc exit data copyout(CPM, FLHC, FLQC, CHS, CH, CHS2, CQS2,&
2511!$acc USTM, wstar, qstar, ZOL, MOL, RMOL, &
2512!$acc HFX, QFX, LH, QSFC, QFLX, HFLX, &
2513!$acc T2, TH2, Q2, WSPD, U10, V10, &
2514!$acc QGH, psim, psih, &
2515!$acc stress_wat, stress_lnd, stress_ice, &
2516!$acc rb_wat, rb_lnd, rb_ice, &
2517!$acc UST_wat, UST_lnd, UST_ice, &
2518!$acc ZNT_wat, ZNT_lnd, ZNT_ice, &
2519!$acc QSFC_lnd, QSFC_wat, QSFC_ice, &
2520!$acc QFLX_lnd, QFLX_wat, QFLX_ice, &
2521!$acc HFLX_lnd, HFLX_wat, HFLX_ice, &
2522!$acc PSIX_wat, PSIX_lnd, PSIX_ice, &
2523!$acc PSIX10_wat, PSIX10_lnd, PSIX10_ice, &
2524!$acc PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
2525!$acc PSIT_lnd, PSIT_wat, PSIT_ice, &
2526!$acc ch_lnd, ch_wat, ch_ice, &
2527!$acc cm_lnd, cm_wat, cm_ice, &
2528!$acc device_errmsg, device_errflg)
2529
2530! Final sync of device and host error flags and messages
2531IF (device_errflg /= 0) THEN
2532 errflg = device_errflg
2533 errmsg = device_errmsg
2534ENDIF
2535
2536!$acc exit data delete( flag_iter, dry, wet, icy, dx, &
2537!$acc MAVAIL, PBLH, PSFCPA, z0pert, ztpert, &
2538!$acc QV1D, U1D, V1D, U1D2, V1D2, T1D, P1D, &
2539!$acc rstoch1D, sigmaf, shdmax, vegtype, &
2540!$acc dz2w1d, dz8w1d, &
2541!$acc snowh_wat, snowh_lnd, snowh_ice, &
2542!$acc tskin_wat, tskin_lnd, tskin_ice, &
2543!$acc tsurf_wat, tsurf_lnd, tsurf_ice)
2544
2545!$acc exit data delete( ZA, ZA2, THV1D, TH1D, TC1D, TV1D, &
2546!$acc RHO1D, QVSH, PSIH2, PSIM10, PSIH10, WSPDI, &
2547!$acc GOVRTH, PSFC, THCON, &
2548!$acc zratio_lnd, zratio_ice, zratio_wat, &
2549!$acc TSK_lnd, TSK_ice, TSK_wat, &
2550!$acc THSK_lnd, THSK_ice, THSK_wat, &
2551!$acc THVSK_lnd, THVSK_ice, THVSK_wat, &
2552!$acc GZ1OZ0_lnd, GZ1OZ0_ice, GZ1OZ0_wat, &
2553!$acc GZ1OZt_lnd, GZ1OZt_ice, GZ1OZt_wat, &
2554!$acc GZ2OZ0_lnd, GZ2OZ0_ice, GZ2OZ0_wat, &
2555!$acc GZ2OZt_lnd, GZ2OZt_ice, GZ2OZt_wat, &
2556!$acc GZ10OZ0_lnd, GZ10OZ0_ice, GZ10OZ0_wat, &
2557!$acc GZ10OZt_lnd, GZ10OZt_ice, GZ10OZt_wat, &
2558!$acc ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat, &
2559!$acc ZT_lnd, ZT_ice, ZT_wat, &
2560!$acc ZQ_lnd, ZQ_ice, ZQ_wat, &
2561!$acc PSIQ_lnd, PSIQ_ice, PSIQ_wat, &
2562!$acc PSIQ2_lnd, PSIQ2_ice, PSIQ2_wat, &
2563!$acc QSFCMR_lnd, QSFCMR_ice, QSFCMR_wat )
2564
2565END SUBROUTINE sfclay1d_mynn
2566!-------------------------------------------------------------------
2576 SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,&
2577 & landsea,IZ0TLND2,spp_sfc,rstoch)
2578
2579 !$acc routine seq
2580 IMPLICIT NONE
2581 REAL(kind_phys), INTENT(IN) :: Z_0,restar,ustar,KARMAN,landsea
2582 INTEGER, OPTIONAL, INTENT(IN) :: IZ0TLND2
2583 REAL(kind_phys), INTENT(OUT) :: Zt,Zq
2584 REAL(kind_phys) :: CZIL !=0.100 in Chen et al. (1997)
2585 !=0.075 in Zilitinkevich (1995)
2586 !=0.500 in Lemone et al. (2008)
2587 INTEGER, INTENT(IN) :: spp_sfc
2588 REAL(kind_phys), INTENT(IN) :: rstoch
2589
2590
2591 IF (landsea-1.5 .GT. 0) THEN !WATER
2592
2593 !THIS IS BASED ON Zilitinkevich, Grachev, and Fairall (2001;
2594 !Their equations 15 and 16).
2595 IF (restar .LT. 0.1) THEN
2596 zt = z_0*exp(karman*2.0)
2597 zt = min( zt, 6.0e-5_kind_phys)
2598 zt = max( zt, 2.0e-9_kind_phys)
2599 zq = z_0*exp(karman*3.0)
2600 zq = min( zq, 6.0e-5_kind_phys)
2601 zq = max( zq, 2.0e-9_kind_phys)
2602 ELSE
2603 zt = z_0*exp(-karman*(4.0*sqrt(restar)-3.2))
2604 zt = min( zt, 6.0e-5_kind_phys)
2605 zt = max( zt, 2.0e-9_kind_phys)
2606 zq = z_0*exp(-karman*(4.0*sqrt(restar)-4.2))
2607 zq = min( zt, 6.0e-5_kind_phys)
2608 zq = max( zt, 2.0e-9_kind_phys)
2609 ENDIF
2610
2611 ELSE !LAND
2612
2613 !Option to modify CZIL according to Chen & Zhang, 2009
2614 IF ( iz0tlnd2 .EQ. 1 ) THEN
2615 czil = 10.0 ** ( -0.40 * ( z_0 / 0.07 ) )
2616 ELSE
2617 czil = 0.095 !0.075 !0.10
2618 END IF
2619
2620 zt = z_0*exp(-karman*czil*sqrt(restar))
2621 zt = min( zt, 0.75*z_0)
2622
2623 zq = z_0*exp(-karman*czil*sqrt(restar))
2624 zq = min( zq, 0.75*z_0)
2625
2626! stochastically perturb thermal and moisture roughness length.
2627! currently set to half the amplitude:
2628 if (spp_sfc==1) then
2629 zt = zt + zt * 0.5 * rstoch
2630 zt = max(zt, 0.0001_kind_phys)
2631 zq = zt
2632 endif
2633
2634 ENDIF
2635
2636 return
2637
2638 END SUBROUTINE zilitinkevich_1995
2639!--------------------------------------------------------------------
2641 SUBROUTINE davis_etal_2008(Z_0,ustar)
2642
2643 !a.k.a. : Donelan et al. (2004)
2644 !This formulation for roughness length was designed to match
2645 !the labratory experiments of Donelan et al. (2004).
2646 !This is an update version from Davis et al. 2008, which
2647 !corrects a small-bias in Z_0 (AHW real-time 2012).
2648
2649 !$acc routine seq
2650 IMPLICIT NONE
2651 REAL(kind_phys), INTENT(IN) :: ustar
2652 REAL(kind_phys), INTENT(OUT) :: Z_0
2653 REAL(kind_phys) :: ZW, ZN1, ZN2
2654 REAL(kind_phys), PARAMETER :: OZO=1.59e-5
2655
2656 !OLD FORM: Z_0 = 10.*EXP(-10./(ustar**onethird))
2657 !NEW FORM:
2658
2659 zw = min((ustar/1.06)**(0.3),1.0_kind_phys)
2660 zn1 = 0.011*ustar*ustar*g_inv + ozo
2661 zn2 = 10.*exp(-9.5*ustar**(-onethird)) + &
2662 0.11*1.5e-5/max(ustar,0.01_kind_phys)
2663 !0.11*1.5E-5/AMAX1(ustar,0.01)
2664 z_0 = (1.0-zw) * zn1 + zw * zn2
2665
2666 z_0 = max( z_0, 1.27e-7_kind_phys) !These max/mins were suggested by
2667 z_0 = min( z_0, 2.85e-3_kind_phys) !Davis et al. (2008)
2668
2669 return
2670
2671 END SUBROUTINE davis_etal_2008
2672!--------------------------------------------------------------------
2676 SUBROUTINE taylor_yelland_2001(Z_0,ustar,wsp10)
2677 !$acc routine seq
2678 IMPLICIT NONE
2679 REAL(kind_phys), INTENT(IN) :: ustar,wsp10
2680 REAL(kind_phys), INTENT(OUT) :: Z_0
2681 REAL(kind_phys), parameter :: pi=3.14159265
2682 REAL(kind_phys) :: hs, Tp, Lp
2683
2684 !hs is the significant wave height
2685 hs = 0.0248*(wsp10**2.)
2686 !Tp dominant wave period
2687 tp = 0.729*max(wsp10,0.1_kind_phys)
2688 !Lp is the wavelength of the dominant wave
2689 lp = grav*tp**2/(2*pi)
2690
2691 z_0 = 1200.*hs*(hs/lp)**4.5
2692 z_0 = max( z_0, 1.27e-7_kind_phys) !These max/mins were suggested by
2693 z_0 = min( z_0, 2.85e-3_kind_phys) !Davis et al. (2008)
2694
2695 return
2696
2697 END SUBROUTINE taylor_yelland_2001
2698!--------------------------------------------------------------------
2704 SUBROUTINE charnock_1955(Z_0,ustar,wsp10,visc,zu)
2705 !$acc routine seq
2706 IMPLICIT NONE
2707 REAL(kind_phys), INTENT(IN) :: ustar, visc, wsp10, zu
2708 REAL(kind_phys), INTENT(OUT) :: Z_0
2709 REAL(kind_phys), PARAMETER :: CZO2=0.011
2710 REAL(kind_phys) :: CZC ! variable charnock "constant"
2711 REAL(kind_phys) :: wsp10m ! logarithmically calculated 10 m
2712
2713 wsp10m = wsp10*log(10./1e-4)/log(zu/1e-4)
2714 czc = czo2 + 0.007*min(max((wsp10m-10.)/8., 0._kind_phys), 1.0_kind_phys)
2715
2716 z_0 = czc*ustar*ustar*g_inv + (0.11*visc/max(ustar,0.05_kind_phys))
2717 z_0 = max( z_0, 1.27e-7_kind_phys) !These max/mins were suggested by
2718 z_0 = min( z_0, 2.85e-3_kind_phys) !Davis et al. (2008)
2719
2720 return
2721
2722 END SUBROUTINE charnock_1955
2723!--------------------------------------------------------------------
2729 SUBROUTINE edson_etal_2013(Z_0,ustar,wsp10,visc,zu)
2730 !$acc routine seq
2731 IMPLICIT NONE
2732 REAL(kind_phys), INTENT(IN) :: ustar, visc, wsp10, zu
2733 REAL(kind_phys), INTENT(OUT) :: Z_0
2734 REAL(kind_phys), PARAMETER :: m=0.0017, b=-0.005
2735 REAL(kind_phys) :: CZC ! variable charnock "constant"
2736 REAL(kind_phys) :: wsp10m ! logarithmically calculated 10 m
2737
2738 wsp10m = wsp10*log(10/1e-4)/log(zu/1e-4)
2739 wsp10m = min(19._kind_phys, wsp10m)
2740 czc = m*wsp10m + b
2741 czc = max(czc, 0.0_kind_phys)
2742
2743 z_0 = czc*ustar*ustar*g_inv + (0.11*visc/max(ustar,0.07_kind_phys))
2744 z_0 = max( z_0, 1.27e-7_kind_phys) !These max/mins were suggested by
2745 z_0 = min( z_0, 2.85e-3_kind_phys) !Davis et al. (2008)
2746
2747 return
2748
2749 END SUBROUTINE edson_etal_2013
2750!--------------------------------------------------------------------
2758 SUBROUTINE garratt_1992(Zt,Zq,Z_0,Ren,landsea)
2759 !$acc routine seq
2760 IMPLICIT NONE
2761 REAL(kind_phys), INTENT(IN) :: Ren, Z_0,landsea
2762 REAL(kind_phys), INTENT(OUT) :: Zt,Zq
2763 REAL(kind_phys) :: Rq
2764 REAL(kind_phys), PARAMETER :: e=2.71828183
2765
2766 IF (landsea-1.5 .GT. 0) THEN !WATER
2767
2768 zt = z_0*exp(2.0 - (2.48*(ren**0.25)))
2769 zq = z_0*exp(2.0 - (2.28*(ren**0.25)))
2770
2771 zq = min( zq, 5.5e-5_kind_phys)
2772 zq = max( zq, 2.0e-9_kind_phys)
2773 zt = min( zt, 5.5e-5_kind_phys)
2774 zt = max( zt, 2.0e-9_kind_phys) !same lower limit as ECMWF
2775 ELSE !LAND
2776 zq = z_0/(e**2.) !taken from Garratt (1980,1992)
2777 zt = zq
2778 ENDIF
2779
2780 return
2781
2782 END SUBROUTINE garratt_1992
2783!--------------------------------------------------------------------
2794 SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc)
2795 !$acc routine seq
2796 IMPLICIT NONE
2797 REAL(kind_phys), INTENT(IN) :: Ren,ustar,visc,rstoch
2798 INTEGER, INTENT(IN) :: spp_sfc
2799 REAL(kind_phys), INTENT(OUT) :: Zt,Zq
2800
2801 IF (ren .le. 2.) then
2802
2803 zt = (5.5e-5)*(ren**(-0.60))
2804 zq = zt
2805 !FOR SMOOTH SEAS, CAN USE GARRATT
2806 !Zq = 0.2*visc/MAX(ustar,0.1)
2807 !Zq = 0.3*visc/MAX(ustar,0.1)
2808
2809 ELSE
2810
2811 !FOR ROUGH SEAS, USE COARE
2812 zt = (5.5e-5)*(ren**(-0.60))
2813 zq = zt
2814
2815 ENDIF
2816
2817 if (spp_sfc==1) then
2818 zt = zt + zt * 0.5 * rstoch
2819 zq = zt
2820 endif
2821
2822 zt = min(zt,1.0e-4_kind_phys)
2823 zt = max(zt,2.0e-9_kind_phys)
2824
2825 zq = min(zt,1.0e-4_kind_phys)
2826 zq = max(zt,2.0e-9_kind_phys)
2827
2828 return
2829
2830 END SUBROUTINE fairall_etal_2003
2831!--------------------------------------------------------------------
2838 SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc)
2839 !$acc routine seq
2840 IMPLICIT NONE
2841 REAL(kind_phys), INTENT(IN) :: Ren,ustar,visc,rstoch
2842 INTEGER, INTENT(IN) :: spp_sfc
2843 REAL(kind_phys), INTENT(OUT) :: Zt,Zq
2844
2845 !Zt = (5.5e-5)*(Ren**(-0.60))
2846 zt = min(1.6e-4_kind_phys, 5.8e-5/(ren**0.72))
2847 zq = zt
2848
2849 IF (spp_sfc ==1) THEN
2850 zt = max(zt + zt*0.5*rstoch,2.0e-9_kind_phys)
2851 zq = max(zt + zt*0.5*rstoch,2.0e-9_kind_phys)
2852 ELSE
2853 zt = max(zt,2.0e-9_kind_phys)
2854 zq = max(zt,2.0e-9_kind_phys)
2855 ENDIF
2856
2857 return
2858
2859 END SUBROUTINE fairall_etal_2014
2860!--------------------------------------------------------------------
2885 SUBROUTINE yang_2008(Z_0,Zt,Zq,ustar,tstar,qst,Ren,visc)
2886
2887 !$acc routine seq
2888 IMPLICIT NONE
2889 REAL(kind_phys), INTENT(IN) :: Z_0, Ren, ustar, tstar, qst, visc
2890 REAL(kind_phys) :: ht, &! roughness height at critical Reynolds number
2891 tstar2, &! bounded T*, forced to be non-positive
2892 qstar2, &! bounded q*, forced to be non-positive
2893 Z_02, &! bounded Z_0 for variable Renc2 calc
2894 Renc2 ! variable Renc, function of Z_0
2895 REAL(kind_phys), INTENT(OUT) :: Zt,Zq
2896 REAL(kind_phys), PARAMETER :: Renc=300., & !old constant Renc
2897 beta=1.5, & !important for diurnal variation
2898 m=170., & !slope for Renc2 function
2899 b=691. !y-intercept for Renc2 function
2900
2901 z_02 = min(z_0,0.5_kind_phys)
2902 z_02 = max(z_02,0.04_kind_phys)
2903 renc2= b + m*log(z_02)
2904 ht = renc2*visc/max(ustar,0.01_kind_phys)
2905 tstar2 = min(tstar, 0.0_kind_phys)
2906 qstar2 = min(qst,0.0_kind_phys)
2907
2908 zt = ht * exp(-beta*(ustar**0.5)*(abs(tstar2)**1.0))
2909 zq = ht * exp(-beta*(ustar**0.5)*(abs(qstar2)**1.0))
2910 !Zq = Zt
2911
2912 zt = min(zt, z_0/2.0)
2913 zq = min(zq, z_0/2.0)
2914
2915 return
2916
2917 END SUBROUTINE yang_2008
2918!--------------------------------------------------------------------
2919! Taken from the GFS (sfc_diff.f) for comparison
2921 SUBROUTINE gfs_z0_lnd(z0max,shdmax,z1,vegtype,ivegsrc,z0pert)
2922
2923 !$acc routine seq
2924 REAL(kind_phys), INTENT(OUT) :: z0max
2925 REAL(kind_phys), INTENT(IN) :: shdmax,z1,z0pert
2926 INTEGER, INTENT(IN) :: vegtype,ivegsrc
2927 REAL(kind_phys) :: tem1, tem2
2928
2929! z0max = max(1.0e-6, min(0.01 * z0max, z1))
2930!already converted into meters in the wrapper
2931 z0max = max(1.0e-6_kind_phys, min(z0max, z1))
2932!** xubin's new z0 over land
2933 tem1 = 1.0 - shdmax
2934 tem2 = tem1 * tem1
2935 tem1 = 1.0 - tem2
2936
2937 if( ivegsrc == 1 ) then
2938
2939 if (vegtype == 10) then
2940 z0max = exp( tem2*log01 + tem1*log07 )
2941 elseif (vegtype == 6) then
2942 z0max = exp( tem2*log01 + tem1*log05 )
2943 elseif (vegtype == 7) then
2944! z0max = exp( tem2*log01 + tem1*log01 )
2945 z0max = 0.01
2946 elseif (vegtype == 16) then
2947! z0max = exp( tem2*log01 + tem1*log01 )
2948 z0max = 0.01
2949 else
2950 z0max = exp( tem2*log01 + tem1*log(z0max) )
2951 endif
2952
2953 elseif (ivegsrc == 2 ) then
2954
2955 if (vegtype == 7) then
2956 z0max = exp( tem2*log01 + tem1*log07 )
2957 elseif (vegtype == 8) then
2958 z0max = exp( tem2*log01 + tem1*log05 )
2959 elseif (vegtype == 9) then
2960! z0max = exp( tem2*log01 + tem1*log01 )
2961 z0max = 0.01
2962 elseif (vegtype == 11) then
2963! z0max = exp( tem2*log01 + tem1*log01 )
2964 z0max = 0.01
2965 else
2966 z0max = exp( tem2*log01 + tem1*log(z0max) )
2967 endif
2968
2969 endif
2970
2971! mg, sfc-perts: add surface perturbations to z0max over land
2972 if (z0pert /= 0.0 ) then
2973 z0max = z0max * (10.**z0pert)
2974 endif
2975
2976 z0max = max(z0max, 1.0e-6_kind_phys)
2977
2978 END SUBROUTINE gfs_z0_lnd
2979!--------------------------------------------------------------------
2980! Taken from the GFS (sfc_diff.f) for comparison
2982 SUBROUTINE gfs_zt_lnd(ztmax,z0max,sigmaf,ztpert,ustar_lnd)
2983
2984 !$acc routine seq
2985 REAL(kind_phys), INTENT(OUT) :: ztmax
2986 REAL(kind_phys), INTENT(IN) :: z0max,sigmaf,ztpert,ustar_lnd
2987 REAL(kind_phys) :: czilc, tem1, tem2
2988 REAL(kind_phys), PARAMETER :: ca = 0.4
2989
2990! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil
2991 czilc = 0.8
2992
2993 tem1 = 1.0 - sigmaf
2994 ztmax = z0max*exp( - tem1*tem1 &
2995 & * czilc*ca*sqrt(ustar_lnd*(0.01/1.5e-05)))
2996!
2997! czilc = 10.0 ** (- 4. * z0max) ! Trier et al. (2011, WAF)
2998! ztmax = z0max * exp( - czilc * ca &
2999! & * 258.2 * sqrt(ustar_lnd*z0max) )
3000
3001
3002! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land
3003 if (ztpert /= 0.0) then
3004 ztmax = ztmax * (10.**ztpert)
3005 endif
3006 ztmax = max(ztmax, 1.0e-6_kind_phys)
3007
3008 END SUBROUTINE gfs_zt_lnd
3009!--------------------------------------------------------------------
3011 SUBROUTINE gfs_z0_wat(z0rl_wat,ustar_wat,WSPD,z1,sfc_z0_type,redrag)
3012
3013 !$acc routine seq
3014 REAL(kind_phys), INTENT(OUT) :: z0rl_wat
3015 REAL(kind_phys), INTENT(INOUT):: ustar_wat
3016 REAL(kind_phys), INTENT(IN) :: wspd,z1
3017 LOGICAL, INTENT(IN) :: redrag
3018 INTEGER, INTENT(IN) :: sfc_z0_type
3019 REAL(kind_phys) :: z0,z0max,wind10m
3020 REAL(kind_phys), PARAMETER :: charnock = 0.014, z0s_max=.317e-2
3021
3022! z0 = 0.01 * z0rl_wat
3023!Already converted to meters in the wrapper
3024 z0 = z0rl_wat
3025 z0max = max(1.0e-6_kind_phys, min(z0,z1))
3026 ustar_wat = sqrt(grav * z0 / charnock)
3027 wind10m = wspd*log(10./1e-4)/log(z1/1e-4)
3028 !wind10m = sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i))
3029!
3030 if (sfc_z0_type >= 0) then
3031 if (sfc_z0_type == 0) then
3032 z0 = (charnock / grav) * ustar_wat * ustar_wat
3033
3034! mbek -- toga-coare flux algorithm
3035! z0 = (charnock / g) * ustar(i)*ustar(i) + arnu/ustar(i)
3036! new implementation of z0
3037! cc = ustar(i) * z0 / rnu
3038! pp = cc / (1. + cc)
3039! ff = g * arnu / (charnock * ustar(i) ** 3)
3040! z0 = arnu / (ustar(i) * ff ** pp)
3041
3042 if (redrag) then
3043 !z0rl_wat = 100.0 * max(min(z0, z0s_max), 1.e-7)
3044 z0rl_wat = max(min(z0, z0s_max), 1.e-7_kind_phys)
3045 else
3046 !z0rl_wat = 100.0 * max(min(z0,.1), 1.e-7)
3047 z0rl_wat = max(min(z0,.1_kind_phys), 1.e-7_kind_phys)
3048 endif
3049
3050 elseif (sfc_z0_type == 6) then ! wang
3051 call znot_m_v6(wind10m, z0) ! wind, m/s, z0, m
3052 !z0rl_wat = 100.0 * z0 ! cm
3053 elseif (sfc_z0_type == 7) then ! wang
3054 call znot_m_v7(wind10m, z0) ! wind, m/s, z0, m
3055 !z0rl_wat = 100.0 * z0 ! cm
3056 else
3057 z0rl_wat = 1.0e-6
3058 endif
3059
3060 endif
3061
3062 END SUBROUTINE gfs_z0_wat
3063!--------------------------------------------------------------------
3065 SUBROUTINE gfs_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,device_errmsg,device_errflg)
3066 !$acc routine seq
3067 real(kind_phys), INTENT(OUT) :: ztmax
3068 real(kind_phys), INTENT(IN) :: wspd,z1,z0rl_wat,restar
3069 INTEGER, INTENT(IN) :: sfc_z0_type
3070
3071! Using device_errmsg and device_errflg rather than the CCPP errmsg and errflg
3072! so that this subroutine can be run on an accelerator device with OpenACC.
3073! character(len=*), intent(out) :: errmsg
3074! integer, intent(out) :: errflg
3075 character(len=512), intent(out) :: device_errmsg
3076 integer, intent(out) :: device_errflg
3077
3078 real(kind_phys) :: z0,z0max,wind10m,rat,ustar_wat
3079 real(kind_phys), PARAMETER :: charnock = 0.014, z0s_max=.317e-2
3080
3081 ! Initialize error-handling
3082! errflg = 0
3083! errmsg = ''
3084 device_errflg = 0
3085 device_errmsg = ''
3086
3087! z0 = 0.01 * z0rl_wat
3088!Already converted to meters in the wrapper
3089 z0 = z0rl_wat
3090 z0max = max(1.0e-6_kind_phys, min(z0,z1))
3091 ustar_wat = sqrt(grav * z0 / charnock)
3092 wind10m = wspd*log(10./1e-4)/log(z1/1e-4)
3093
3094!** test xubin's new z0
3095
3096! ztmax = z0max
3097
3098!input restar = max(ustar_wat(i)*z0max*visi, 0.000001)
3099
3100! restar = log(restar)
3101! restar = min(restar,5.)
3102! restar = max(restar,-5.)
3103! rat = aa1 + (bb1 + cc1*restar) * restar
3104! rat = rat / (1. + (bb2 + cc2*restar) * restar))
3105! rat taken from zeng, zhao and dickinson 1997
3106
3107 rat = min(7.0_kind_phys, 2.67 * sqrt(sqrt(restar)) - 2.57)
3108 ztmax = max(z0max * exp(-rat), 1.0e-6_kind_phys)
3109!
3110 if (sfc_z0_type == 6) then
3111 call znot_t_v6(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m)
3112 else if (sfc_z0_type == 7) then
3113 call znot_t_v7(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m)
3114 else if (sfc_z0_type > 0) then
3115 write(0,*)'no option for sfc_z0_type=',sfc_z0_type
3116! errflg = 1
3117! errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
3118 device_errflg = 1
3119 device_errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
3120 return
3121
3122 endif
3123
3124 END SUBROUTINE gfs_zt_wat
3125!--------------------------------------------------------------------
3129
3130 SUBROUTINE znot_m_v6(uref, znotm)
3131 !$acc routine seq
3132 use machine , only : kind_phys
3133 IMPLICIT NONE
3134! Calculate areodynamical roughness over water with input 10-m wind
3135! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
3136! For high winds, try to fit available observational data
3137!
3138! Bin Liu, NOAA/NCEP/EMC 2017
3139!
3140! uref(m/s) : wind speed at 10-m height
3141! znotm(meter): areodynamical roughness scale over water
3142!
3143
3144 REAL(kind_phys), INTENT(IN) :: uref
3145 REAL(kind_phys), INTENT(OUT):: znotm
3146 REAL(kind_phys), PARAMETER :: p13 = -1.296521881682694e-02, &
3147 & p12 = 2.855780863283819e-01, p11 = -1.597898515251717e+00,&
3148 & p10 = -8.396975715683501e+00, &
3149
3150 & p25 = 3.790846746036765e-10, p24 = 3.281964357650687e-09,&
3151 & p23 = 1.962282433562894e-07, p22 = -1.240239171056262e-06,&
3152 & p21 = 1.739759082358234e-07, p20 = 2.147264020369413e-05,&
3153
3154 & p35 = 1.840430200185075e-07, p34 = -2.793849676757154e-05,&
3155 & p33 = 1.735308193700643e-03, p32 = -6.139315534216305e-02,&
3156 & p31 = 1.255457892775006e+00, p30 = -1.663993561652530e+01,&
3157
3158 & p40 = 4.579369142033410e-04
3159
3160
3161 if (uref >= 0.0 .and. uref <= 6.5 ) then
3162 znotm = exp(p10 + uref * (p11 + uref * (p12 + uref*p13)))
3163 elseif (uref > 6.5 .and. uref <= 15.7) then
3164 znotm = p20 + uref * (p21 + uref * (p22 + uref * (p23 &
3165 & + uref * (p24 + uref * p25))))
3166 elseif (uref > 15.7 .and. uref <= 53.0) then
3167 znotm = exp( p30 + uref * (p31 + uref * (p32 + uref * (p33 &
3168 & + uref * (p34 + uref * p35)))))
3169 elseif ( uref > 53.0) then
3170 znotm = p40
3171 else
3172 print*, 'Wrong input uref value:',uref
3173 endif
3174
3175 END SUBROUTINE znot_m_v6
3176!--------------------------------------------------------------------
3183!
3184! uref(m/s) : wind speed at 10-m height
3185! znott(meter): scalar roughness scale over water
3186 SUBROUTINE znot_t_v6(uref, znott)
3187
3188 !$acc routine seq
3189 IMPLICIT NONE
3190!
3191 REAL(kind_phys), INTENT(IN) :: uref
3192 REAL(kind_phys), INTENT(OUT):: znott
3193 REAL(kind_phys), PARAMETER :: p00 = 1.100000000000000e-04,&
3194 & p15 = -9.144581627678278e-10, p14 = 7.020346616456421e-08,&
3195 & p13 = -2.155602086883837e-06, p12 = 3.333848806567684e-05,&
3196 & p11 = -2.628501274963990e-04, p10 = 8.634221567969181e-04,&
3197
3198 & p25 = -8.654513012535990e-12, p24 = 1.232380050058077e-09,&
3199 & p23 = -6.837922749505057e-08, p22 = 1.871407733439947e-06,&
3200 & p21 = -2.552246987137160e-05, p20 = 1.428968311457630e-04,&
3201
3202 & p35 = 3.207515102100162e-12, p34 = -2.945761895342535e-10,&
3203 & p33 = 8.788972147364181e-09, p32 = -3.814457439412957e-08,&
3204 & p31 = -2.448983648874671e-06, p30 = 3.436721779020359e-05,&
3205
3206 & p45 = -3.530687797132211e-11, p44 = 3.939867958963747e-09,&
3207 & p43 = -1.227668406985956e-08, p42 = -1.367469811838390e-05,&
3208 & p41 = 5.988240863928883e-04, p40 = -7.746288511324971e-03,&
3209
3210 & p56 = -1.187982453329086e-13, p55 = 4.801984186231693e-11,&
3211 & p54 = -8.049200462388188e-09, p53 = 7.169872601310186e-07,&
3212 & p52 = -3.581694433758150e-05, p51 = 9.503919224192534e-04,&
3213 & p50 = -1.036679430885215e-02, &
3214
3215 & p60 = 4.751256171799112e-05
3216
3217 if (uref >= 0.0 .and. uref < 5.9 ) then
3218 znott = p00
3219 elseif (uref >= 5.9 .and. uref <= 15.4) then
3220 znott = p10 + uref * (p11 + uref * (p12 + uref * (p13 &
3221 & + uref * (p14 + uref * p15))))
3222 elseif (uref > 15.4 .and. uref <= 21.6) then
3223 znott = p20 + uref * (p21 + uref * (p22 + uref * (p23 &
3224 & + uref * (p24 + uref * p25))))
3225 elseif (uref > 21.6 .and. uref <= 42.2) then
3226 znott = p30 + uref * (p31 + uref * (p32 + uref * (p33 &
3227 & + uref * (p34 + uref * p35))))
3228 elseif ( uref > 42.2 .and. uref <= 53.3) then
3229 znott = p40 + uref * (p41 + uref * (p42 + uref * (p43 &
3230 & + uref * (p44 + uref * p45))))
3231 elseif ( uref > 53.3 .and. uref <= 80.0) then
3232 znott = p50 + uref * (p51 + uref * (p52 + uref * (p53 &
3233 & + uref * (p54 + uref * (p55 + uref * p56)))))
3234 elseif ( uref > 80.0) then
3235 znott = p60
3236 else
3237 print*, 'Wrong input uref value:',uref
3238 endif
3239
3240 END SUBROUTINE znot_t_v6
3241
3242!-------------------------------------------------------------------
3250 SUBROUTINE znot_m_v7(uref, znotm)
3251
3252 !$acc routine seq
3253 IMPLICIT NONE
3254!
3255! uref(m/s) : wind speed at 10-m height
3256! znotm(meter): areodynamical roughness scale over water
3257!
3258
3259 REAL(kind_phys), INTENT(IN) :: uref
3260 REAL(kind_phys), INTENT(OUT):: znotm
3261
3262 REAL(kind_phys), PARAMETER :: p13 = -1.296521881682694e-02,&
3263 & p12 = 2.855780863283819e-01, p11 = -1.597898515251717e+00,&
3264 & p10 = -8.396975715683501e+00, &
3265
3266 & p25 = 3.790846746036765e-10, p24 = 3.281964357650687e-09,&
3267 & p23 = 1.962282433562894e-07, p22 = -1.240239171056262e-06,&
3268 & p21 = 1.739759082358234e-07, p20 = 2.147264020369413e-05,&
3269
3270 & p35 = 1.897534489606422e-07, p34 = -3.019495980684978e-05,&
3271 & p33 = 1.931392924987349e-03, p32 = -6.797293095862357e-02,&
3272 & p31 = 1.346757797103756e+00, p30 = -1.707846930193362e+01,&
3273
3274 & p40 = 3.371427455376717e-04
3275
3276 if (uref >= 0.0 .and. uref <= 6.5 ) then
3277 znotm = exp( p10 + uref * (p11 + uref * (p12 + uref * p13)))
3278 elseif (uref > 6.5 .and. uref <= 15.7) then
3279 znotm = p20 + uref * (p21 + uref * (p22 + uref * (p23 &
3280 & + uref * (p24 + uref * p25))))
3281 elseif (uref > 15.7 .and. uref <= 53.0) then
3282 znotm = exp( p30 + uref * (p31 + uref * (p32 + uref * (p33 &
3283 & + uref * (p34 + uref * p35)))))
3284 elseif ( uref > 53.0) then
3285 znotm = p40
3286 else
3287 print*, 'Wrong input uref value:',uref
3288 endif
3289
3290 END SUBROUTINE znot_m_v7
3291!--------------------------------------------------------------------
3299 SUBROUTINE znot_t_v7(uref, znott)
3300
3301 !$acc routine seq
3302 IMPLICIT NONE
3303!
3304! uref(m/s) : wind speed at 10-m height
3305! znott(meter): scalar roughness scale over water
3306!
3307
3308 REAL(kind_phys), INTENT(IN) :: uref
3309 REAL(kind_phys), INTENT(OUT):: znott
3310 REAL(kind_phys), PARAMETER :: p00 = 1.100000000000000e-04,&
3311 & p15 = -9.193764479895316e-10, p14 = 7.052217518653943e-08,&
3312 & p13 = -2.163419217747114e-06, p12 = 3.342963077911962e-05,&
3313 & p11 = -2.633566691328004e-04, p10 = 8.644979973037803e-04,&
3314
3315 & p25 = -9.402722450219142e-12, p24 = 1.325396583616614e-09,&
3316 & p23 = -7.299148051141852e-08, p22 = 1.982901461144764e-06,&
3317 & p21 = -2.680293455916390e-05, p20 = 1.484341646128200e-04,&
3318
3319 & p35 = 7.921446674311864e-12, p34 = -1.019028029546602e-09,&
3320 & p33 = 5.251986927351103e-08, p32 = -1.337841892062716e-06,&
3321 & p31 = 1.659454106237737e-05, p30 = -7.558911792344770e-05,&
3322
3323 & p45 = -2.694370426850801e-10, p44 = 5.817362913967911e-08,&
3324 & p43 = -5.000813324746342e-06, p42 = 2.143803523428029e-04,&
3325 & p41 = -4.588070983722060e-03, p40 = 3.924356617245624e-02,&
3326
3327 & p56 = -1.663918773476178e-13, p55 = 6.724854483077447e-11,&
3328 & p54 = -1.127030176632823e-08, p53 = 1.003683177025925e-06,&
3329 & p52 = -5.012618091180904e-05, p51 = 1.329762020689302e-03,&
3330 & p50 = -1.450062148367566e-02, p60 = 6.840803042788488e-05
3331
3332 if (uref >= 0.0 .and. uref < 5.9 ) then
3333 znott = p00
3334 elseif (uref >= 5.9 .and. uref <= 15.4) then
3335 znott = p10 + uref * (p11 + uref * (p12 + uref * (p13 &
3336 & + uref * (p14 + uref * p15))))
3337 elseif (uref > 15.4 .and. uref <= 21.6) then
3338 znott = p20 + uref * (p21 + uref * (p22 + uref * (p23 &
3339 & + uref * (p24 + uref * p25))))
3340 elseif (uref > 21.6 .and. uref <= 42.6) then
3341 znott = p30 + uref * (p31 + uref * (p32 + uref * (p33 &
3342 & + uref * (p34 + uref * p35))))
3343 elseif ( uref > 42.6 .and. uref <= 53.0) then
3344 znott = p40 + uref * (p41 + uref * (p42 + uref * (p43 &
3345 & + uref * (p44 + uref * p45))))
3346 elseif ( uref > 53.0 .and. uref <= 80.0) then
3347 znott = p50 + uref * (p51 + uref * (p52 + uref * (p53 &
3348 & + uref * (p54 + uref * (p55 + uref * p56)))))
3349 elseif ( uref > 80.0) then
3350 znott = p60
3351 else
3352 print*, 'Wrong input uref value:',uref
3353 endif
3354
3355 END SUBROUTINE znot_t_v7
3356
3357!--------------------------------------------------------------------
3363 SUBROUTINE andreas_2002(Z_0,bvisc,ustar,Zt,Zq)
3364
3365 !$acc routine seq
3366 IMPLICIT NONE
3367 REAL(kind_phys), INTENT(IN) :: Z_0, bvisc, ustar
3368 REAL(kind_phys), INTENT(OUT) :: Zt, Zq
3369 REAL(kind_phys) :: Ren2, zntsno
3370
3371 REAL(kind_phys), PARAMETER :: &
3372 bt0_s=1.25, bt0_t=0.149, bt0_r=0.317, &
3373 bt1_s=0.0, bt1_t=-0.55, bt1_r=-0.565, &
3374 bt2_s=0.0, bt2_t=0.0, bt2_r=-0.183
3375
3376 REAL(kind_phys), PARAMETER :: &
3377 bq0_s=1.61, bq0_t=0.351, bq0_r=0.396, &
3378 bq1_s=0.0, bq1_t=-0.628, bq1_r=-0.512, &
3379 bq2_s=0.0, bq2_t=0.0, bq2_r=-0.180
3380
3381 !Calculate zo for snow (Andreas et al. 2005, BLM)
3382 zntsno = 0.135*bvisc/ustar + &
3383 (0.035*(ustar*ustar)*g_inv) * &
3384 (5.*exp(-1.*(((ustar - 0.18)/0.1)*((ustar - 0.18)/0.1))) + 1.)
3385 ren2 = ustar*zntsno/bvisc
3386
3387 ! Make sure that Re is not outside of the range of validity
3388 ! for using their equations
3389 IF (ren2 .gt. 1000.) ren2 = 1000.
3390
3391 IF (ren2 .le. 0.135) then
3392
3393 zt = zntsno*exp(bt0_s + bt1_s*log(ren2) + bt2_s*log(ren2)**2)
3394 zq = zntsno*exp(bq0_s + bq1_s*log(ren2) + bq2_s*log(ren2)**2)
3395
3396 ELSE IF (ren2 .gt. 0.135 .AND. ren2 .lt. 2.5) then
3397
3398 zt = zntsno*exp(bt0_t + bt1_t*log(ren2) + bt2_t*log(ren2)**2)
3399 zq = zntsno*exp(bq0_t + bq1_t*log(ren2) + bq2_t*log(ren2)**2)
3400
3401 ELSE
3402
3403 zt = zntsno*exp(bt0_r + bt1_r*log(ren2) + bt2_r*log(ren2)**2)
3404 zq = zntsno*exp(bq0_r + bq1_r*log(ren2) + bq2_r*log(ren2)**2)
3405
3406 ENDIF
3407
3408 return
3409
3410 END SUBROUTINE andreas_2002
3411!--------------------------------------------------------------------
3415 SUBROUTINE psi_hogstrom_1996(psi_m, psi_h, zL, Zt, Z_0, Za)
3416
3417 IMPLICIT NONE
3418 REAL(kind_phys), INTENT(IN) :: zL, Zt, Z_0, Za
3419 REAL(kind_phys), INTENT(OUT) :: psi_m, psi_h
3420 REAL(kind_phys) :: x, x0, y, y0, zmL, zhL
3421
3422 zml = z_0*zl/za
3423 zhl = zt*zl/za
3424
3425 IF (zl .gt. 0.) THEN !STABLE (not well tested - seem large)
3426
3427 psi_m = -5.3*(zl - zml)
3428 psi_h = -8.0*(zl - zhl)
3429
3430 ELSE !UNSTABLE
3431
3432 x = (1.-19.0*zl)**0.25
3433 x0= (1.-19.0*zml)**0.25
3434 y = (1.-11.6*zl)**0.5
3435 y0= (1.-11.6*zhl)**0.5
3436
3437 psi_m = 2.*log((1.+x)/(1.+x0)) + &
3438 &log((1.+x**2.)/(1.+x0**2.)) - &
3439 &2.0*atan(x) + 2.0*atan(x0)
3440 psi_h = 2.*log((1.+y)/(1.+y0))
3441
3442 ENDIF
3443
3444 return
3445
3446 END SUBROUTINE psi_hogstrom_1996
3447!--------------------------------------------------------------------
3453 SUBROUTINE psi_dyerhicks(psi_m, psi_h, zL, Zt, Z_0, Za)
3454
3455 IMPLICIT NONE
3456 REAL(kind_phys), INTENT(IN) :: zL, Zt, Z_0, Za
3457 REAL(kind_phys), INTENT(OUT) :: psi_m, psi_h
3458 REAL(kind_phys) :: x, x0, y, y0, zmL, zhL
3459
3460 zml = z_0*zl/za !Zo/L
3461 zhl = zt*zl/za !Zt/L
3462
3463 IF (zl .gt. 0.) THEN !STABLE
3464
3465 psi_m = -5.0*(zl - zml)
3466 psi_h = -5.0*(zl - zhl)
3467
3468 ELSE !UNSTABLE
3469
3470 x = (1.-16.*zl)**0.25
3471 x0= (1.-16.*zml)**0.25
3472
3473 y = (1.-16.*zl)**0.5
3474 y0= (1.-16.*zhl)**0.5
3475
3476 psi_m = 2.*log((1.+x)/(1.+x0)) + &
3477 &log((1.+x**2.)/(1.+x0**2.)) - &
3478 &2.0*atan(x) + 2.0*atan(x0)
3479 psi_h = 2.*log((1.+y)/(1.+y0))
3480
3481 ENDIF
3482
3483 return
3484
3485 END SUBROUTINE psi_dyerhicks
3486!--------------------------------------------------------------------
3491 SUBROUTINE psi_beljaars_holtslag_1991(psi_m, psi_h, zL)
3492
3493 IMPLICIT NONE
3494 REAL(kind_phys), INTENT(IN) :: zL
3495 REAL(kind_phys), INTENT(OUT) :: psi_m, psi_h
3496 REAL(kind_phys), PARAMETER :: a=1., b=0.666, c=5., d=0.35
3497
3498 IF (zl .lt. 0.) THEN !UNSTABLE
3499
3500 WRITE(*,*)"WARNING: Universal stability functions from"
3501 WRITE(*,*)" Beljaars and Holtslag (1991) should only"
3502 WRITE(*,*)" be used in the stable regime!"
3503 psi_m = 0.
3504 psi_h = 0.
3505
3506 ELSE !STABLE
3507
3508 psi_m = -(a*zl + b*(zl -(c/d))*exp(-d*zl) + (b*c/d))
3509 psi_h = -((1.+.666*a*zl)**1.5 + &
3510 b*(zl - (c/d))*exp(-d*zl) + (b*c/d) -1.)
3511
3512 ENDIF
3513
3514 return
3515
3516 END SUBROUTINE psi_beljaars_holtslag_1991
3517!--------------------------------------------------------------------
3523 SUBROUTINE psi_zilitinkevich_esau_2007(psi_m, psi_h, zL)
3524
3525 IMPLICIT NONE
3526 REAL(kind_phys), INTENT(IN) :: zL
3527 REAL(kind_phys), INTENT(OUT) :: psi_m, psi_h
3528 REAL(kind_phys), PARAMETER :: Cm=3.0, ct=2.5
3529
3530 IF (zl .lt. 0.) THEN !UNSTABLE
3531
3532 WRITE(*,*)"WARNING: Universal stability function from"
3533 WRITE(*,*)" Zilitinkevich and Esau (2007) should only"
3534 WRITE(*,*)" be used in the stable regime!"
3535 psi_m = 0.
3536 psi_h = 0.
3537
3538 ELSE !STABLE
3539
3540 psi_m = -cm*(zl**(5./6.))
3541 psi_h = -ct*(zl**(4./5.))
3542
3543 ENDIF
3544
3545 return
3546
3547 END SUBROUTINE psi_zilitinkevich_esau_2007
3548!--------------------------------------------------------------------
3552 SUBROUTINE psi_businger_1971(psi_m, psi_h, zL)
3553
3554 IMPLICIT NONE
3555 REAL(kind_phys), INTENT(IN) :: zL
3556 REAL(kind_phys), INTENT(OUT) :: psi_m, psi_h
3557 REAL(kind_phys) :: x, y
3558 REAL(kind_phys), PARAMETER :: Pi180 = 3.14159265/180.
3559
3560 IF (zl .lt. 0.) THEN !UNSTABLE
3561
3562 x = (1. - 15.0*zl)**0.25
3563 y = (1. - 9.0*zl)**0.5
3564
3565 psi_m = log(((1.+x)/2.)**2.) + &
3566 &log((1.+x**2.)/2.) - &
3567 &2.0*atan(x) + pi180*90.
3568 psi_h = 2.*log((1.+y)/2.)
3569
3570 ELSE !STABLE
3571
3572 psi_m = -4.7*zl
3573 psi_h = -(4.7/0.74)*zl
3574
3575 ENDIF
3576
3577 return
3578
3579 END SUBROUTINE psi_businger_1971
3580!--------------------------------------------------------------------
3588 SUBROUTINE psi_suselj_sood_2010(psi_m, psi_h, zL)
3589
3590 IMPLICIT NONE
3591 REAL(kind_phys), INTENT(IN) :: zL
3592 REAL(kind_phys), INTENT(OUT) :: psi_m, psi_h
3593 REAL(kind_phys), PARAMETER :: Rfc=0.19, ric=0.183, phit=0.8
3594
3595 IF (zl .gt. 0.) THEN !STABLE
3596
3597 psi_m = -(zl/rfc + 1.1223*exp(1.-1.6666/zl))
3598 !psi_h = -zL*Ric/((Rfc**2.)*PHIT) + 8.209*(zL**1.1091)
3599 !THEIR EQ FOR PSI_H CRASHES THE MODEL AND DOES NOT MATCH
3600 !THEIR FIG 1. THIS EQ (BELOW) MATCHES THEIR FIG 1 BETTER:
3601 psi_h = -(zl*ric/((rfc**2.)*5.) + 7.09*(zl**1.1091))
3602
3603 ELSE !UNSTABLE
3604
3605 psi_m = 0.9904*log(1. - 14.264*zl)
3606 psi_h = 1.0103*log(1. - 16.3066*zl)
3607
3608 ENDIF
3609
3610 return
3611
3612 END SUBROUTINE psi_suselj_sood_2010
3613!--------------------------------------------------------------------
3618 SUBROUTINE psi_cb2005(psim1,psih1,zL,z0L)
3619
3620 IMPLICIT NONE
3621 REAL(kind_phys), INTENT(IN) :: zL,z0L
3622 REAL(kind_phys), INTENT(OUT) :: psim1,psih1
3623
3624 psim1 = -6.1*log(zl + (1.+ zl**2.5)**0.4) &
3625 -6.1*log(z0l + (1.+ z0l**2.5)**0.4)
3626 psih1 = -5.5*log(zl + (1.+ zl**1.1)**0.90909090909) &
3627 -5.5*log(z0l + (1.+ z0l**1.1)**0.90909090909)
3628
3629 return
3630
3631 END SUBROUTINE psi_cb2005
3632!--------------------------------------------------------------------
3637 SUBROUTINE li_etal_2010(zL, Rib, zaz0, z0zt)
3638
3639 !$acc routine seq
3640 IMPLICIT NONE
3641 REAL(kind_phys), INTENT(OUT) :: zL
3642 REAL(kind_phys), INTENT(IN) :: Rib, zaz0, z0zt
3643 REAL(kind_phys) :: alfa, beta, zaz02, z0zt2
3644 REAL(kind_phys), PARAMETER :: &
3645 & au11=0.045, bu11=0.003, bu12=0.0059, &
3646 & bu21=-0.0828, bu22=0.8845, bu31=0.1739, &
3647 & bu32=-0.9213, bu33=-0.1057
3648 REAL(kind_phys), PARAMETER :: &
3649 & aw11=0.5738, aw12=-0.4399, aw21=-4.901, &
3650 & aw22=52.50, bw11=-0.0539, bw12=1.540, &
3651 & bw21=-0.669, bw22=-3.282
3652 REAL(kind_phys), PARAMETER :: &
3653 & as11=0.7529, as21=14.94, bs11=0.1569, &
3654 & bs21=-0.3091, bs22=-1.303
3655
3656 !set limits according to Li et al (2010), p 157.
3657 zaz02=zaz0
3658 IF (zaz0 .lt. 100.0) zaz02=100.
3659 IF (zaz0 .gt. 100000.0) zaz02=100000.
3660
3661 !set more limits according to Li et al (2010)
3662 z0zt2=z0zt
3663 IF (z0zt .lt. 0.5) z0zt2=0.5
3664 IF (z0zt .gt. 100.0) z0zt2=100.
3665
3666 alfa = log(zaz02)
3667 beta = log(z0zt2)
3668
3669 IF (rib .le. 0.0) THEN
3670 zl = au11*alfa*rib**2 + ( &
3671 & (bu11*beta + bu12)*alfa**2 + &
3672 & (bu21*beta + bu22)*alfa + &
3673 & (bu31*beta**2 + bu32*beta + bu33))*rib
3674 !if(zL .LT. -15 .OR. zl .GT. 0.)print*,"VIOLATION Rib<0:",zL
3675 zl = max(zl,-15._kind_phys) !LIMITS SET ACCORDING TO Li et al (2010)
3676 zl = min(zl,0._kind_phys) !Figure 1.
3677 ELSEIF (rib .gt. 0.0 .AND. rib .le. 0.2) THEN
3678 zl = ((aw11*beta + aw12)*alfa + &
3679 & (aw21*beta + aw22))*rib**2 + &
3680 & ((bw11*beta + bw12)*alfa + &
3681 & (bw21*beta + bw22))*rib
3682 !if(zL .LT. 0 .OR. zl .GT. 4)print*,"VIOLATION 0<Rib<0.2:",zL
3683 zl = min(zl,4._kind_phys) !LIMITS APPROX SET ACCORDING TO Li et al (2010)
3684 zl = max(zl,0._kind_phys) !THEIR FIGURE 1B.
3685 ELSE
3686 zl = (as11*alfa + as21)*rib + bs11*alfa + &
3687 & bs21*beta + bs22
3688 !if(zL .LE. 1 .OR. zl .GT. 23)print*,"VIOLATION Rib>0.2:",zL
3689 zl = min(zl,20._kind_phys) !LIMITS ACCORDING TO Li et al (2010), THIER
3690 !FIGUE 1C.
3691 zl = max(zl,1._kind_phys)
3692 ENDIF
3693
3694 return
3695
3696 END SUBROUTINE li_etal_2010
3697!-------------------------------------------------------------------
3699 REAL(kind_phys) function zolri(ri,za,z0,zt,zol1,psi_opt)
3700
3706
3707 IMPLICIT NONE
3708 REAL(kind_phys), INTENT(IN) :: ri,za,z0,zt,zol1
3709 INTEGER, INTENT(IN) :: psi_opt
3710 REAL(kind_phys) :: x1,x2,fx1,fx2
3711 INTEGER :: n
3712 INTEGER, PARAMETER :: nmax = 20
3713 !REAL(kind_phys), DIMENSION(nmax):: zLhux
3714
3715 if (ri.lt.0.)then
3716 x1=zol1 - 0.02 !-5.
3717 x2=0.
3718 else
3719 x1=0.
3720 x2=zol1 + 0.02 !5.
3721 endif
3722
3723 n=1
3724 fx1=zolri2(x1,ri,za,z0,zt,psi_opt)
3725 fx2=zolri2(x2,ri,za,z0,zt,psi_opt)
3726
3727 Do While (abs(x1 - x2) > 0.01 .and. n < nmax)
3728 if(abs(fx2).lt.abs(fx1))then
3729 x1=x1-fx1/(fx2-fx1)*(x2-x1)
3730 fx1=zolri2(x1,ri,za,z0,zt,psi_opt)
3731 zolri=x1
3732 else
3733 x2=x2-fx2/(fx2-fx1)*(x2-x1)
3734 fx2=zolri2(x2,ri,za,z0,zt,psi_opt)
3735 zolri=x2
3736 endif
3737 n=n+1
3738 !print*," n=",n," x1=",x1," x2=",x2
3739 !zLhux(n)=zolri
3740 enddo
3741
3742 if (n==nmax .and. abs(x1 - x2) >= 0.01) then
3743 !if convergence fails, use approximate values:
3744 CALL li_etal_2010(zolri, ri, za/z0, z0/zt)
3745 !zLhux(n)=zolri
3746 !print*,"iter FAIL, n=",n," Ri=",ri," z0=",z0
3747 else
3748 !print*,"SUCCESS,n=",n," Ri=",ri," z0=",z0
3749 endif
3750
3751 return
3752 end function
3753!-------------------------------------------------------------------
3754 REAL(kind_phys) function zolri2(zol2,ri2,za,z0,zt,psi_opt)
3755
3756 ! INPUT: =================================
3757 ! zol2 - estimated z/L
3758 ! ri2 - calculated bulk Richardson number
3759 ! za - 1/2 depth of first model layer
3760 ! z0 - aerodynamic roughness length
3761 ! zt - thermal roughness length
3762 ! OUTPUT: ================================
3763 ! zolri2 - delta Ri
3764
3765 IMPLICIT NONE
3766 INTEGER, INTENT(IN) :: psi_opt
3767 REAL(kind_phys), INTENT(IN) :: ri2,za,z0,zt
3768 REAL(kind_phys), INTENT(INOUT) :: zol2
3769 REAL(kind_phys) :: zol20,zol3,psim1,psih1,psix2,psit2,zolt
3770
3771 if(zol2*ri2 .lt. 0.)zol2=0. ! limit zol2 - must be same sign as ri2
3772
3773 zol20=zol2*z0/za ! z0/L
3774 zol3=zol2+zol20 ! (z+z0)/L
3775 zolt=zol2*zt/za ! zt/L
3776
3777 if (ri2.lt.0) then
3778 !psix2=log((za+z0)/z0)-(psim_unstable(zol3)-psim_unstable(zol20))
3779 !psit2=log((za+zt)/zt)-(psih_unstable(zol3)-psih_unstable(zol20))
3780 psit2=max(log((za+z0)/zt)-(psih_unstable(zol3,psi_opt)-psih_unstable(zolt,psi_opt)), 1.0_kind_phys)
3781 psix2=max(log((za+z0)/z0)-(psim_unstable(zol3,psi_opt)-psim_unstable(zol20,psi_opt)), 1.0_kind_phys)
3782 else
3783 !psix2=log((za+z0)/z0)-(psim_stable(zol3)-psim_stable(zol20))
3784 !psit2=log((za+zt)/zt)-(psih_stable(zol3)-psih_stable(zol20))
3785 psit2=max(log((za+z0)/zt)-(psih_stable(zol3,psi_opt)-psih_stable(zolt,psi_opt)), 1.0_kind_phys)
3786 psix2=max(log((za+z0)/z0)-(psim_stable(zol3,psi_opt)-psim_stable(zol20,psi_opt)),1.0_kind_phys)
3787 endif
3788
3789 zolri2=zol2*psit2/psix2**2 - ri2
3790 !print*," target ri=",ri2," est ri=",zol2*psit2/psix2**2
3791
3792 return
3793 end function
3794!====================================================================
3795
3796 REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
3797
3798 !$acc routine seq
3799 ! This iterative algorithm to compute z/L from bulk-Ri
3800
3801 IMPLICIT NONE
3802 REAL(kind_phys), INTENT(IN) :: ri,za,z0,zt,logz0,logzt
3803 INTEGER, INTENT(IN) :: psi_opt
3804 REAL(kind_phys), INTENT(INOUT) :: zol1
3805 REAL(kind_phys) :: zol20,zol3,zolt,zolold
3806 INTEGER :: n
3807 INTEGER, PARAMETER :: nmax = 20
3808 !REAL(kind_phys), DIMENSION(nmax):: zLhux
3809 REAL(kind_phys) :: psit2,psix2
3810
3811 !print*,"+++++++INCOMING: z/L=",zol1," ri=",ri
3812 if (zol1*ri .lt. 0.) THEN
3813 !print*,"begin: WRONG QUADRANTS: z/L=",zol1," ri=",ri
3814 zol1=0.
3815 endif
3816
3817 if (ri .lt. 0.) then
3818 zolold=-99999.
3819 zolrib=-66666.
3820 else
3821 zolold=99999.
3822 zolrib=66666.
3823 endif
3824 n=1
3825
3826 DO While (abs(zolold - zolrib) > 0.01 .and. n < nmax)
3827
3828 if(n==1)then
3829 zolold=zol1
3830 else
3831 zolold=zolrib
3832 endif
3833 zol20=zolold*z0/za ! z0/L
3834 zol3=zolold+zol20 ! (z+z0)/L
3835 zolt=zolold*zt/za ! zt/L
3836 !print*,"z0/L=",zol20," (z+z0)/L=",zol3," zt/L=",zolt
3837 if (ri.lt.0) then
3838 !psit2=log((za+zt)/zt)-(psih_unstable(zol3)-psih_unstable(zol20))
3839 !psit2=log((za+z0)/zt)-(psih_unstable(zol3)-psih_unstable(zol20))
3840 psit2=max(logzt-(psih_unstable(zol3,psi_opt)-psih_unstable(zolt,psi_opt)), 1.0_kind_phys)
3841 psix2=max(logz0-(psim_unstable(zol3,psi_opt)-psim_unstable(zol20,psi_opt)), 1.0_kind_phys)
3842 else
3843 !psit2=log((za+zt)/zt)-(psih_stable(zol3)-psih_stable(zol20))
3844 !psit2=log((za+z0)/zt)-(psih_stable(zol3)-psih_stable(zol20))
3845 psit2=max(logzt-(psih_stable(zol3,psi_opt)-psih_stable(zolt,psi_opt)), 1.0_kind_phys)
3846 psix2=max(logz0-(psim_stable(zol3,psi_opt)-psim_stable(zol20,psi_opt)), 1.0_kind_phys)
3847 endif
3848 !print*,"n=",n," psit2=",psit2," psix2=",psix2
3849 zolrib=ri*psix2**2/psit2
3850 !zLhux(n)=zolrib
3851 n=n+1
3852 enddo
3853
3854 if (n==nmax .and. abs(zolold - zolrib) > 0.01 ) then
3855 !print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri
3856 !if convergence fails, use approximate values:
3857 CALL li_etal_2010(zolrib, ri, za/z0, z0/zt)
3858 !zLhux(n)=zolrib
3859 !print*,"FAILED, n=",n," Ri=",ri," z0=",z0
3860 !print*,"z/L=",zLhux(1:nmax)
3861 else
3862 !if(zolrib*ri .lt. 0.) THEN
3863 ! !print*,"end: WRONG QUADRANTS: z/L=",zolrib," ri=",ri
3864 ! !CALL Li_etal_2010(zolrib, ri, za/z0, z0/zt)
3865 !endif
3866 !print*,"SUCCESS,n=",n," Ri=",ri," z0=",z0
3867 endif
3868
3869 return
3870 end function
3871!====================================================================
3874 SUBROUTINE psi_init(psi_opt,errmsg,errflg)
3875
3876 integer :: N,psi_opt
3877 real(kind_phys) :: zolf
3878 character(len=*), intent(out) :: errmsg
3879 integer, intent(out) :: errflg
3880
3881 if (psi_opt == 0) then
3882 DO n=0,1000
3883 ! stable function tables
3884 zolf = float(n)*0.01
3885 psim_stab(n)=psim_stable_full(zolf)
3886 psih_stab(n)=psih_stable_full(zolf)
3887
3888 ! unstable function tables
3889 zolf = -float(n)*0.01
3890 psim_unstab(n)=psim_unstable_full(zolf)
3891 psih_unstab(n)=psih_unstable_full(zolf)
3892 ENDDO
3893 else
3894 DO n=0,1000
3895 ! stable function tables
3896 zolf = float(n)*0.01
3897 psim_stab(n)=psim_stable_full_gfs(zolf)
3898 psih_stab(n)=psih_stable_full_gfs(zolf)
3899
3900 ! unstable function tables
3901 zolf = -float(n)*0.01
3902 psim_unstab(n)=psim_unstable_full_gfs(zolf)
3903 psih_unstab(n)=psih_unstable_full_gfs(zolf)
3904 ENDDO
3905 endif
3906
3907 !Simple test to see if initialization worked:
3908 if (psim_stab(1) < 0. .AND. psih_stab(1) < 0. .AND. &
3909 psim_unstab(1) > 0. .AND. psih_unstab(1) > 0.) then
3910 errmsg = 'In MYNN SFC, Psi tables have been initialized'
3911 errflg = 0
3912 else
3913 errmsg = 'Error in MYNN SFC: Problem initializing psi tables'
3914 errflg = 1
3915 endif
3916
3917 END SUBROUTINE psi_init
3918! ==================================================================
3919! ... integrated similarity functions from MYNN...
3920!
3922 real(kind_phys) function psim_stable_full(zolf)
3923 !$acc routine seq
3924 real(kind_phys) :: zolf
3925
3926 !psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
3927 psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**0.4)
3928
3929 return
3930 end function
3931
3933 real(kind_phys) function psih_stable_full(zolf)
3934 !$acc routine seq
3935 real(kind_phys) :: zolf
3936
3937 !psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
3938 psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**0.9090909090909090909)
3939
3940 return
3941 end function
3942
3944 real(kind_phys) function psim_unstable_full(zolf)
3945 !$acc routine seq
3946 real(kind_phys) :: zolf,x,ym,psimc,psimk
3947
3948 x=(1.-16.*zolf)**.25
3949 !psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*ATAN(1.)
3950 !psimk=2.*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*atan1
3951 psimk=2.*log(0.5*(1+x))+log(0.5*(1+x*x))-2.*atan(x)+2.*atan1
3952
3953 ym=(1.-10.*zolf)**onethird
3954 !psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
3955 psimc=1.5*log((ym**2 + ym+1.)*onethird)-sqrt3*atan((2.*ym+1)/sqrt3)+4.*atan1/sqrt3
3956
3957 psim_unstable_full=(psimk+zolf**2*(psimc))/(1+zolf**2.)
3958
3959 return
3960 end function
3961
3963 real(kind_phys) function psih_unstable_full(zolf)
3964 !$acc routine seq
3965 real(kind_phys) :: zolf,y,yh,psihc,psihk
3966
3967 y=(1.-16.*zolf)**.5
3968 !psihk=2.*log((1+y)/2.)
3969 psihk=2.*log((1+y)*0.5)
3970
3971 yh=(1.-34.*zolf)**onethird
3972 !psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
3973 psihc=1.5*log((yh**2.+yh+1.)*onethird)-sqrt3*atan((2.*yh+1)/sqrt3)+4.*atan1/sqrt3
3974
3975 psih_unstable_full=(psihk+zolf**2*(psihc))/(1+zolf**2)
3976
3977 return
3978 end function
3979
3980! ==================================================================
3981! ... integrated similarity functions from GFS...
3982!
3985 REAL(kind_phys) function psim_stable_full_gfs(zolf)
3986 !$acc routine seq
3987 REAL(kind_phys) :: zolf
3988 REAL(kind_phys), PARAMETER :: alpha4 = 20.
3989 REAL(kind_phys) :: aa
3990
3991 aa = sqrt(1. + alpha4 * zolf)
3992 psim_stable_full_gfs = -1.*aa + log(aa + 1.)
3993
3994 return
3995 end function
3996
3999 real(kind_phys) function psih_stable_full_gfs(zolf)
4000 !$acc routine seq
4001 real(kind_phys) :: zolf
4002 real(kind_phys), PARAMETER :: alpha4 = 20.
4003 real(kind_phys) :: bb
4004
4005 bb = sqrt(1. + alpha4 * zolf)
4006 psih_stable_full_gfs = -1.*bb + log(bb + 1.)
4007
4008 return
4009 end function
4010
4013 real(kind_phys) function psim_unstable_full_gfs(zolf)
4014 !$acc routine seq
4015 real(kind_phys) :: zolf
4016 real(kind_phys) :: hl1,tem1
4017 real(kind_phys), PARAMETER :: a0=-3.975, a1=12.32, &
4018 b1=-7.755, b2=6.041
4019
4020 if (zolf .ge. -0.5) then
4021 hl1 = zolf
4022 psim_unstable_full_gfs = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1)
4023 else
4024 hl1 = -zolf
4025 tem1 = 1.0 / sqrt(hl1)
4026 psim_unstable_full_gfs = log(hl1) + 2. * sqrt(tem1) - .8776
4027 end if
4028
4029 return
4030 end function
4031
4034 real(kind_phys) function psih_unstable_full_gfs(zolf)
4035 !$acc routine seq
4036 real(kind_phys) :: zolf
4037 real(kind_phys) :: hl1,tem1
4038 real(kind_phys), PARAMETER :: a0p=-7.941, a1p=24.75, &
4039 b1p=-8.705, b2p=7.899
4040
4041 if (zolf .ge. -0.5) then
4042 hl1 = zolf
4043 psih_unstable_full_gfs = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1)
4044 else
4045 hl1 = -zolf
4046 tem1 = 1.0 / sqrt(hl1)
4047 psih_unstable_full_gfs = log(hl1) + .5 * tem1 + 1.386
4048 end if
4049
4050 return
4051 end function
4052
4055 real(kind_phys) function psim_stable(zolf,psi_opt)
4056 !$acc routine seq
4057 integer :: nzol,psi_opt
4058 real(kind_phys) :: rzol,zolf
4059
4060 nzol = int(zolf*100.)
4061 rzol = zolf*100. - nzol
4062 if(nzol+1 .lt. 1000)then
4063 psim_stable = psim_stab(nzol) + rzol*(psim_stab(nzol+1)-psim_stab(nzol))
4064 else
4065 if (psi_opt == 0) then
4066 psim_stable = psim_stable_full(zolf)
4067 else
4068 psim_stable = psim_stable_full_gfs(zolf)
4069 endif
4070 endif
4071
4072 return
4073 end function
4074
4076 real(kind_phys) function psih_stable(zolf,psi_opt)
4077 !$acc routine seq
4078 integer :: nzol,psi_opt
4079 real(kind_phys) :: rzol,zolf
4080
4081 nzol = int(zolf*100.)
4082 rzol = zolf*100. - nzol
4083 if(nzol+1 .lt. 1000)then
4084 psih_stable = psih_stab(nzol) + rzol*(psih_stab(nzol+1)-psih_stab(nzol))
4085 else
4086 if (psi_opt == 0) then
4087 psih_stable = psih_stable_full(zolf)
4088 else
4089 psih_stable = psih_stable_full_gfs(zolf)
4090 endif
4091 endif
4092
4093 return
4094 end function
4095
4097 real(kind_phys) function psim_unstable(zolf,psi_opt)
4098 !$acc routine seq
4099 integer :: nzol,psi_opt
4100 real(kind_phys) :: rzol,zolf
4101
4102 nzol = int(-zolf*100.)
4103 rzol = -zolf*100. - nzol
4104 if(nzol+1 .lt. 1000)then
4105 psim_unstable = psim_unstab(nzol) + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol))
4106 else
4107 if (psi_opt == 0) then
4108 psim_unstable = psim_unstable_full(zolf)
4109 else
4110 psim_unstable = psim_unstable_full_gfs(zolf)
4111 endif
4112 endif
4113
4114 return
4115 end function
4116
4118 real(kind_phys) function psih_unstable(zolf,psi_opt)
4119 !$acc routine seq
4120 integer :: nzol,psi_opt
4121 real(kind_phys) :: rzol,zolf
4122
4123 nzol = int(-zolf*100.)
4124 rzol = -zolf*100. - nzol
4125 if(nzol+1 .lt. 1000)then
4126 psih_unstable = psih_unstab(nzol) + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol))
4127 else
4128 if (psi_opt == 0) then
4129 psih_unstable = psih_unstable_full(zolf)
4130 else
4131 psih_unstable = psih_unstable_full_gfs(zolf)
4132 endif
4133 endif
4134
4135 return
4136 end function
4137!========================================================================
4138
4139END MODULE module_sf_mynn
snowz0
subroutine snowz0
This subroutine calculates total roughness length over snow.
Definition sflx.f:2957
module_sf_mynn::psim_stable_full_gfs
real(kind_phys) function psim_stable_full_gfs(zolf)
Definition module_sf_mynn.F90:3986
module_sf_mynn::znot_t_v7
subroutine znot_t_v7(uref, znott)
Calculate scalar roughness over water with input 10-m wind For low-to-moderate winds,...
Definition module_sf_mynn.F90:3300
module_sf_mynn::gfs_zt_wat
subroutine gfs_zt_wat(ztmax, z0rl_wat, restar, wspd, z1, sfc_z0_type, device_errmsg, device_errflg)
Definition module_sf_mynn.F90:3066
module_sf_mynn::psih_stable_full
real(kind_phys) function psih_stable_full(zolf)
Definition module_sf_mynn.F90:3934
module_sf_mynn::psih_stable_full_gfs
real(kind_phys) function psih_stable_full_gfs(zolf)
Definition module_sf_mynn.F90:4000
module_sf_mynn::charnock_1955
subroutine charnock_1955(z_0, ustar, wsp10, visc, zu)
This version of Charnock's relation employs a varying Charnock parameter, similar to COARE3....
Definition module_sf_mynn.F90:2705
module_sf_mynn::yang_2008
subroutine yang_2008(z_0, zt, zq, ustar, tstar, qst, ren, visc)
This is a modified version of Yang et al (2002 QJRMS, 2008 JAMC) and Chen et al (2010,...
Definition module_sf_mynn.F90:2886
module_sf_mynn::psi_suselj_sood_2010
subroutine psi_suselj_sood_2010(psi_m, psi_h, zl)
This subroutine returns flux-profile relatioships based off of Lobocki (1993), which is derived from ...
Definition module_sf_mynn.F90:3589
module_sf_mynn::psih_stable
real(kind_phys) function psih_stable(zolf, psi_opt)
Definition module_sf_mynn.F90:4077
module_sf_mynn::psi_cb2005
subroutine psi_cb2005(psim1, psih1, zl, z0l)
This subroutine returns the stability functions based off of Cheng and Brutseart (2005,...
Definition module_sf_mynn.F90:3619
module_sf_mynn::li_etal_2010
subroutine li_etal_2010(zl, rib, zaz0, z0zt)
This subroutine returns a more robust z/L that best matches the z/L from Hogstrom (1996) for unstable...
Definition module_sf_mynn.F90:3638
module_sf_mynn::psim_stable
real(kind_phys) function psim_stable(zolf, psi_opt)
look-up table functions - or, if beyond -10 < z/L < 10, recalculate
Definition module_sf_mynn.F90:4056
module_sf_mynn::andreas_2002
subroutine andreas_2002(z_0, bvisc, ustar, zt, zq)
This is taken from Andreas (2002; J. of Hydromet) and Andreas et al. (2005; BLM).
Definition module_sf_mynn.F90:3364
module_sf_mynn::fairall_etal_2014
subroutine fairall_etal_2014(zt, zq, ren, ustar, visc, rstoch, spp_sfc)
This formulation for thermal and moisture roughness length (Zt and Zq) as a function of the roughness...
Definition module_sf_mynn.F90:2839
module_sf_mynn::sfclay1d_mynn
subroutine sfclay1d_mynn(flag_iter, j, u1d, v1d, t1d, qv1d, p1d, dz8w1d, u1d2, v1d2, dz2w1d, psfcpa, pblh, mavail, xland, dx, isfflx, isftcflx, iz0tlnd, psi_opt, compute_flux, compute_diag, sigmaf, vegtype, shdmax, ivegsrc, z0pert, ztpert, redrag, sfc_z0_type, itimestep, iter, flag_restart, lsm, lsm_ruc, wet, dry, icy, tskin_wat, tskin_lnd, tskin_ice, tsurf_wat, tsurf_lnd, tsurf_ice, qsfc_wat, qsfc_lnd, qsfc_ice, snowh_wat, snowh_lnd, snowh_ice, znt_wat, znt_lnd, znt_ice, ust_wat, ust_lnd, ust_ice, cm_wat, cm_lnd, cm_ice, ch_wat, ch_lnd, ch_ice, rb_wat, rb_lnd, rb_ice, stress_wat, stress_lnd, stress_ice, psix_wat, psix_lnd, psix_ice, psit_wat, psit_lnd, psit_ice, psix10_wat, psix10_lnd, psix10_ice, psit2_wat, psit2_lnd, psit2_ice, hflx_wat, hflx_lnd, hflx_ice, qflx_wat, qflx_lnd, qflx_ice, ch, chs, chs2, cqs2, cpm, znt, ustm, zol, mol, rmol, psim, psih, hflx, hfx, qflx, qfx, lh, flhc, flqc, qgh, qsfc, u10, v10, th2, t2, q2, gz1oz0, wspd, wstar, qstar, spp_sfc, rstoch1d, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte, errmsg, errflg)
This subroutine calculates u*, z/L, and the exchange coefficients which are passed to subsequent sche...
Definition module_sf_mynn.F90:559
module_sf_mynn::psi_hogstrom_1996
subroutine psi_hogstrom_1996(psi_m, psi_h, zl, zt, z_0, za)
This subroutine returns the stability functions based off of Hogstrom (1996).
Definition module_sf_mynn.F90:3416
module_sf_mynn::psih_unstable
real(kind_phys) function psih_unstable(zolf, psi_opt)
Definition module_sf_mynn.F90:4119
module_sf_mynn::davis_etal_2008
subroutine davis_etal_2008(z_0, ustar)
Definition module_sf_mynn.F90:2642
module_sf_mynn::znot_m_v7
subroutine znot_m_v7(uref, znotm)
Calculate areodynamical roughness over water with input 10-m wind For low-to-moderate winds,...
Definition module_sf_mynn.F90:3251
module_sf_mynn::psi_init
subroutine psi_init(psi_opt, errmsg, errflg)
Definition module_sf_mynn.F90:3875
module_sf_mynn::zolri
real(kind_phys) function zolri(ri, za, z0, zt, zol1, psi_opt)
Definition module_sf_mynn.F90:3700
module_sf_mynn::gfs_z0_lnd
subroutine gfs_z0_lnd(z0max, shdmax, z1, vegtype, ivegsrc, z0pert)
Definition module_sf_mynn.F90:2922
module_sf_mynn::psim_unstable
real(kind_phys) function psim_unstable(zolf, psi_opt)
Definition module_sf_mynn.F90:4098
module_sf_mynn::znot_t_v6
subroutine znot_t_v6(uref, znott)
Calculate scalar roughness over water with input 10-m wind For low-to-moderate winds,...
Definition module_sf_mynn.F90:3187
module_sf_mynn::psim_unstable_full_gfs
real(kind_phys) function psim_unstable_full_gfs(zolf)
Definition module_sf_mynn.F90:4014
module_sf_mynn::taylor_yelland_2001
subroutine taylor_yelland_2001(z_0, ustar, wsp10)
This formulation for roughness length was designed account for. wave steepness.
Definition module_sf_mynn.F90:2677
module_sf_mynn::psih_unstable_full_gfs
real(kind_phys) function psih_unstable_full_gfs(zolf)
Definition module_sf_mynn.F90:4035
module_sf_mynn::gfs_z0_wat
subroutine gfs_z0_wat(z0rl_wat, ustar_wat, wspd, z1, sfc_z0_type, redrag)
Definition module_sf_mynn.F90:3012
module_sf_mynn::psim_stable_full
real(kind_phys) function psim_stable_full(zolf)
Definition module_sf_mynn.F90:3923
module_sf_mynn::psi_dyerhicks
subroutine psi_dyerhicks(psi_m, psi_h, zl, zt, z_0, za)
This subroutine returns the stability functions based off of Hogstrom (1996), but with different cons...
Definition module_sf_mynn.F90:3454
module_sf_mynn::psi_zilitinkevich_esau_2007
subroutine psi_zilitinkevich_esau_2007(psi_m, psi_h, zl)
This subroutine returns the stability functions come from Zilitinkevich and Esau (2007,...
Definition module_sf_mynn.F90:3524
module_sf_mynn::garratt_1992
subroutine garratt_1992(zt, zq, z_0, ren, landsea)
This formulation for the thermal and moisture roughness lengths (Zt and Zq) relates them to Z0 via th...
Definition module_sf_mynn.F90:2759
module_sf_mynn::sfclay_mynn
subroutine sfclay_mynn(u3d, v3d, t3d, qv3d, p3d, dz8w, th3d, pi3d, qc3d, psfcpa, pblh, mavail, xland, dx, isfflx, isftcflx, lsm, lsm_ruc, compute_flux, compute_diag, iz0tlnd, psi_opt, sigmaf, vegtype, shdmax, ivegsrc, z0pert, ztpert, redrag, sfc_z0_type, itimestep, iter, flag_iter, flag_restart, wet, dry, icy, tskin_wat, tskin_lnd, tskin_ice, tsurf_wat, tsurf_lnd, tsurf_ice, qsfc_wat, qsfc_lnd, qsfc_ice, snowh_wat, snowh_lnd, snowh_ice, znt_wat, znt_lnd, znt_ice, ust_wat, ust_lnd, ust_ice, cm_wat, cm_lnd, cm_ice, ch_wat, ch_lnd, ch_ice, rb_wat, rb_lnd, rb_ice, stress_wat, stress_lnd, stress_ice, fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fm10_wat, fm10_lnd, fm10_ice, fh2_wat, fh2_lnd, fh2_ice, hflx_wat, hflx_lnd, hflx_ice, qflx_wat, qflx_lnd, qflx_ice, ch, chs, chs2, cqs2, cpm, znt, ustm, zol, mol, rmol, psim, psih, hflx, hfx, qflx, qfx, lh, flhc, flqc, qgh, qsfc, u10, v10, th2, t2, q2, gz1oz0, wspd, wstar, spp_sfc, pattern_spp_sfc, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte, errmsg, errflg)
This subroutine.
Definition module_sf_mynn.F90:160
module_sf_mynn::psi_beljaars_holtslag_1991
subroutine psi_beljaars_holtslag_1991(psi_m, psi_h, zl)
This subroutine returns the stability functions based off of Beljaar and Holtslag 1991,...
Definition module_sf_mynn.F90:3492
module_sf_mynn::zilitinkevich_1995
subroutine zilitinkevich_1995(z_0, zt, zq, restar, ustar, karman, landsea, iz0tlnd2, spp_sfc, rstoch)
This subroutine returns the thermal and moisture roughness lengths from Zilitinkevich (1995) and Zili...
Definition module_sf_mynn.F90:2578
module_sf_mynn::psi_businger_1971
subroutine psi_businger_1971(psi_m, psi_h, zl)
This subroutine returns the flux-profile relationships of Businger el al. 1971.
Definition module_sf_mynn.F90:3553
module_sf_mynn::psim_unstable_full
real(kind_phys) function psim_unstable_full(zolf)
Definition module_sf_mynn.F90:3945
module_sf_mynn::psih_unstable_full
real(kind_phys) function psih_unstable_full(zolf)
Definition module_sf_mynn.F90:3964
module_sf_mynn::edson_etal_2013
subroutine edson_etal_2013(z_0, ustar, wsp10, visc, zu)
This version of Charnock's relation employs a varying Charnock parameter, taken from COARE 3....
Definition module_sf_mynn.F90:2730
module_sf_mynn::znot_m_v6
subroutine znot_m_v6(uref, znotm)
add fitted z0,zt curves for hurricane application (used in HWRF/HMON) Weiguo Wang,...
Definition module_sf_mynn.F90:3131
module_sf_mynn::gfs_zt_lnd
subroutine gfs_zt_lnd(ztmax, z0max, sigmaf, ztpert, ustar_lnd)
Definition module_sf_mynn.F90:2983
module_sf_mynn::fairall_etal_2003
subroutine fairall_etal_2003(zt, zq, ren, ustar, visc, rstoch, spp_sfc)
This formulation for thermal and moisture roughness length (Zt and Zq) as a function of the roughness...
Definition module_sf_mynn.F90:2795
lsm_ruc
This module contain the RUC land surface model driver.
Definition lsm_ruc.F90:5
module_sf_mynn
This module contain routines to calculate stability parameters, kinematic siscosity in MYNN surface l...
Definition module_sf_mynn.F90:8