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