CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
gfdl_sfc_layer.F90
1
3
6
7 use machine , only : kind_phys
8
9 implicit none
10
11 public :: gfdl_sfc_layer_init, gfdl_sfc_layer_run
12
13 private
14
15 contains
16
20 subroutine gfdl_sfc_layer_init (icoef_sf, cplwav, cplwav2atm, lcurr_sf, &
21 pert_cd, ntsflg, errmsg, errflg)
22
23 implicit none
24
25 integer, intent(in) :: icoef_sf, ntsflg
26 logical, intent(in) :: cplwav, cplwav2atm, lcurr_sf, pert_cd
27
28 character(len=*), intent(out) :: errmsg
29 integer, intent(out) :: errflg
30
31 ! Initialize CCPP error handling variables
32 errmsg = ''
33 errflg = 0
34
35#if HWRF==1
36 write(errmsg,'(*(a))') 'The GFDL surface layer scheme does not support '&
37 //'use of the HWRF preprocessor flag in gfdl_sfc_layer.F90'
38 errflg = 1
39 return
40#endif
41
42 if (icoef_sf < 0 .or. icoef_sf > 8) then
43 write(errmsg,'(*(a))') 'The value of icoef_sf is outside of the ' &
44 //'supported range (0-8) in gfdl_sfc_layer.F90'
45 errflg = 1
46 return
47 end if
48
49 if (cplwav .or. cplwav2atm) then
50 write(errmsg,'(*(a))') 'The GFDL surface layer scheme is not set up ' &
51 //'to be coupled to waves in gfdl_sfc_layer.F90'
52 errflg = 1
53 return
54 end if
55
56 if (lcurr_sf) then
57 write(errmsg,'(*(a))') 'The GFDL surface layer scheme is not set up ' &
58 //'to be used with the lcurr_sf option in gfdl_sfc_layer.F90'
59 errflg = 1
60 return
61 end if
62
63 if (pert_cd) then
64 write(errmsg,'(*(a))') 'The GFDL surface layer scheme is not set up ' &
65 //'to be used with the pert_cd option in gfdl_sfc_layer.F90'
66 errflg = 1
67 return
68 end if
69
70 if (ntsflg > 0) then
71 !GJF: In order to enable ntsflg > 0, the variable 'tstrc' passed into MFLUX2 should be set
72 ! to the surface_skin_temperature_over_X_interstitial rather than the average of it and
73 ! surface_skin_temperature_after_iteration_over_X
74 write(errmsg,'(*(a))') 'Setting ntsflg > 0 is currently not supported'&
75 //' in gfdl_sfc_layer.F90'
76 errflg = 1
77 return
78 end if
79
80 !GJF: Initialization notes: In WRF, the subroutine module_sf_myjsfc/myjsfcinit
81 ! is called for initialization of the GFDL surface layer scheme from
82 ! the module_physics_init subroutine. It contains the following
83 ! initializations which should already have been done by other
84 ! code in UFS-related host models:
85 ! IF(.NOT.RESTART)THEN
86 ! DO J=JTS,JTE
87 ! DO I=ITS,ITF
88 ! USTAR(I,J)=0.1
89 ! ENDDO
90 ! ENDDO
91 ! ENDIF
92 !also initialize surface roughness length
93
94 end subroutine gfdl_sfc_layer_init
95
99 subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, &
100 lsm_noah, lsm_noahmp, lsm_ruc, icoef_sf, cplwav, karman, &
101 cplwav2atm, lcurr_sf, pert_Cd, ntsflg, sfenth, z1, shdmax, ivegsrc, &
102 vegtype, sigmaf, dt, wet, dry, icy, isltyp, rd, grav, ep1, ep2, smois, &
103 psfc, prsl1, q1, t1, u1, v1, wspd, u10, v10, gsw, glw, tsurf_wat, &
104 tsurf_lnd, tsurf_ice, tskin_wat, tskin_lnd, tskin_ice, ustar_wat, &
105 ustar_lnd, ustar_ice, znt_wat, znt_lnd, znt_ice, cdm_wat, cdm_lnd, &
106 cdm_ice, stress_wat, stress_lnd, stress_ice, rib_wat, rib_lnd, rib_ice, &
107 fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fh2_wat, fh2_lnd, &
108 fh2_ice, ch_wat, ch_lnd, ch_ice, fm10_wat, fm10_lnd, fm10_ice, qss_wat, &
109 qss_lnd, qss_ice, errmsg, errflg)
110
111 use funcphys, only: fpvs
112
113 !#### GJF: temporarily grab parameters from LSM-specific modules -- should go through CCPP ####
114 ! (fixing this involves replacing the functionality of set_soilveg and namelist_soilveg)
115 use namelist_soilveg, only: maxsmc_noah => maxsmc, drysmc_noah => drysmc
116 use namelist_soilveg_ruc, only: maxsmc_ruc => maxsmc, drysmc_ruc => drysmc
117 use noahmp_tables, only: maxsmc_noahmp => smcmax_table, drysmc_noahmp => smcdry_table
118 !################################################################################################
119
120 implicit none
121
122 integer, intent(in) :: im, nsoil, km, ivegsrc
123 integer, intent(in) :: lsm, lsm_noah, lsm_noahmp, &
124 lsm_ruc, icoef_sf, ntsflg
125 logical, intent(in) :: cplwav, cplwav2atm !GJF: this scheme has not been tested with these on
126 logical, intent(in) :: lcurr_sf !GJF: this scheme has not been tested with this option turned on; the variables scurx and scury need to be input in order to use this
127 logical, intent(in) :: pert_cd !GJF: this scheme has not been tested with this option turned on; the variables ens_random_seed and ens_Cdamp need to be input in order to use this
128 logical, dimension(:), intent(in) :: flag_iter, wet, dry, icy
129 integer, dimension(:), intent(in) :: isltyp, vegtype
130 real(kind=kind_phys), intent(in) :: dt, sfenth
131 real(kind=kind_phys), intent(in) :: rd,grav,ep1,ep2
132 real(kind=kind_phys), dimension(:,:), intent(in) :: smois
133 real(kind=kind_phys), dimension(:), intent(in) :: psfc, prsl1, &
134 q1, t1, u1, v1, wspd, u10, v10, gsw, glw, z1, shdmax, sigmaf, xlat, &
135 xlon, tsurf_wat, tsurf_lnd, tsurf_ice
136
137 real(kind=kind_phys), intent(inout), dimension(:) :: tskin_wat, &
138 tskin_lnd, tskin_ice, ustar_wat, ustar_lnd, ustar_ice, &
139 znt_wat, znt_lnd, znt_ice, cdm_wat, cdm_lnd, cdm_ice, &
140 stress_wat, stress_lnd, stress_ice, rib_wat, rib_lnd, rib_ice, &
141 fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fh2_wat, fh2_lnd, &
142 fh2_ice, ch_wat, ch_lnd, ch_ice, fm10_wat, fm10_lnd, fm10_ice, &
143 qss_wat, qss_lnd, qss_ice
144
145 character(len=*), intent(out) :: errmsg
146 integer, intent(out) :: errflg
147
148 !local variables
149
150 integer :: i, its, ite, ims, ime
151
152 logical :: ch_bound_excursion
153
154 !GJF: the vonKarman constant should come in through the CCPP and be defined by the host model
155 real (kind=kind_phys), intent(in) :: karman
156 real (kind=kind_phys), parameter :: log01=log(0.01), log05=log(0.05), &
157 log07=log(0.07)
158
159 !GJF: if the following variables will be used, they should be turned into intent(in) namelist options
160 integer :: iwavecpl, ens_random_seed, issflx
161 logical :: diag_wind10m, diag_qss
162 real(kind=kind_phys) :: ens_cdamp
163
164 real(kind=kind_phys), dimension(im) :: wetc, pspc, pkmax, tstrc, upc, &
165 vpc, mznt, slwdc, wind10, qfx, qgh, zkmax, z1_cm, z0max, ztmax
166 real(kind=kind_phys), dimension(im) :: u10_lnd, u10_ocn, u10_ice, &
167 v10_lnd, v10_ocn, v10_ice
168
169 !GJF: the following variables are identified as:
170 !"SCURX" "Surface Currents(X)" "m s-1"
171 !"SCURY" "Surface Currents(Y)" "m s-1
172 !"CHARN" "Charnock Coeff" " "
173 !"MSANG" "Wind/Stress Angle" "Radian"
174 real(kind=kind_phys), dimension(im) :: charn, msang, scurx, scury
175
176 real(kind=kind_phys), dimension(im) :: fxh, fxe, fxmx, fxmy, xxfh, &
177 xxfh2, tzot
178 real(kind=kind_phys), dimension(1:30) :: maxsmc, drysmc
179 real(kind=kind_phys) :: smcmax, smcdry, zhalf, cd10, &
180 esat, fm_lnd_old, fh_lnd_old, tem1, tem2, czilc, cd_low_limit, &
181 cd_high_limit, ch_low_limit, ch_high_limit, fh2_fh_ratio
182
183 !#### This block will become unnecessary when maxsmc and drysmc come through the CCPP ####
184 if (lsm == lsm_noah) then
185 maxsmc = maxsmc_noah
186 drysmc = drysmc_noah
187 else if (lsm == lsm_noahmp) then
188 maxsmc = maxsmc_noahmp
189 drysmc = drysmc_noahmp
190 else if (lsm == lsm_ruc) then
191 maxsmc = maxsmc_ruc
192 drysmc = drysmc_ruc
193 else
194 !GJF: These data were from the original GFDL surface layer scheme, but
195 ! rather than being hard-coded here, they should be shared with the
196 ! LSM. These data are kept for legacy purposes. Note that these only
197 ! have nonzero values for 16 soil types vs 19 for other STAS datasets
198 data maxsmc/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
199 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
200 0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
201 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
202 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
203 data drysmc/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
204 0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
205 0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
206 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
207 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
208 end if
209 !########################################################################
210
211 !GJF: This code has not been tested with iwavecpl = 1; the variables 'charn' and 'msang' (and others?) need to be input in order to use this
212 ! if (cplwav .or. cplwav2atm) then
213 ! iwavecpl = 1
214 ! else
215 ! iwavecpl = 0
216 ! end if
217 iwavecpl = 0
218
219 !GJF: temporary setting of variables that should be moved to namelist is they are used
220 ens_random_seed = 0 !used for HWRF ensemble?
221 ens_cdamp = 0.0 !used for HWRF ensemble?
222
223 issflx = 0 !GJF: 1 = calculate surface fluxes, 0 = don't
224 diag_wind10m = .false. !GJF: if one wants 10m wind speeds to come from this scheme, set this to True,
225 ! put [u,v]10_[lnd/ocn/ice] in the scheme argument list (and metadata), and modify
226 ! GFS_surface_compsites to receive the individual components and calculate an all-grid value
227 diag_qss = .false. !GJF: saturation specific humidities are calculated by LSM, sea surface, and sea ice schemes in
228 ! GFS-based suites
229
230 ! Initialize CCPP error handling variables
231 errmsg = ''
232 errflg = 0
233
234 its = 1
235 ims = 1
236 ite = im
237 ime = im
238
239 do i=its, ite
240 if (flag_iter(i)) then
241 !GJF: Perform data preparation that is the same for all surface types
242
243 pspc(i) = psfc(i)*10. ! convert from Pa to cgs
244 pkmax(i) = prsl1(i)*10. ! convert from Pa to cgs
245
246 upc(i) = u1(i)*100. ! convert from m s-1 to cm s-1
247 vpc(i) = v1(i)*100. ! convert from m s-1 to cm s-1
248
249 !Wang: use previous u10 v10 to compute wind10, input to MFLUX2 to compute z0 (for first time step, u10 and v10 may be zero)
250 wind10(i)=sqrt(u10(i)*u10(i)+v10(i)*v10(i)) !m s-1
251
252 !Wang: calulate height of the first half level
253 ! if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then
254 ! zhalf = -rd*t1(i)*alog(pkmax(i)/pspc(i))/grav !m
255 ! endif
256
257 !GJF: rather than calculate the height of the first half level, if it is precalculated
258 ! in a different scheme, pass it in and use it; note that in FV3, calculating via the hypsometric equation
259 ! occasionally produced values much shallower than those passed in
260 !zkmax(i) = -rd*t1(i)*alog(pkmax(i)/pspc(i))/grav !m
261 zkmax(i) = z1(i)
262 z1_cm(i) = 100.0*z1(i)
263
264 !GJF: these drag coefficient limits were suggested by Chunxi Zhang via his module_sf_sfclayrev.f90
265 cd_low_limit = 1.0e-5/zkmax(i)
266 cd_high_limit = 0.1
267 !GJF: use the lower of 0.1 from Chunxi Zhang or 0.05/wspd from WRF's module_sf_gfdl.F
268 ! (this will always be the latter if wspd has a minimum of 1.0 m s-1 from above)
269 ch_low_limit = cd_low_limit
270 ch_high_limit = min(0.1,0.05/wspd(i))
271
272 !slwdc... GFDL downward net flux in units of cal/(cm**2/min)
273 !also divide by 10**4 to convert from /m**2 to /cm**2
274 slwdc(i)=gsw(i)+glw(i)
275 slwdc(i)=0.239*60.*slwdc(i)*1.e-4
276
277 !GJF: these variables should be passed in if these options are used
278 charn(i) = 0.0 !used with wave coupling (iwavecpl == 1)
279 msang(i) = 0.0 !used with wave coupling (iwavecpl == 1)
280 scurx(i) = 0.0 !used with ocean currents? (lcurr_sf == T)
281 scury(i) = 0.0 !used with ocean currents? (lcurr_sf == T)
282
283 if (diag_qss) then
284 esat = fpvs(t1(i))
285 qgh(i) = ep2*esat/(psfc(i)-esat)
286 end if
287
288 !GJF: these vars are not needed in a GFS-based suite
289 !rho1(i)=prsl1(i)/(rd*t1(i)*(1.+ep1*q1(i)))
290 !cpm(i)=cp*(1.+0.8*q1(i))
291
292 !GJF: perform data preparation that depends on surface types and call the mflux2 subroutine for each surface type
293 ! Note that this is different than the original WRF module_sf_gfdl.F where mflux2 is called once for all surface
294 ! types, with negative roughness lengths denoting open ocean.
295 if (dry(i)) then
296 !GJF: from WRF's module_sf_gfdl.F
297 smcdry=drysmc(isltyp(i))
298 smcmax=maxsmc(isltyp(i))
299 wetc(i)=(smois(i,1)-smcdry)/(smcmax-smcdry)
300 wetc(i)=amin1(1.,amax1(wetc(i),0.))
301
302 !GJF: the lower boundary temperature passed in to MFLUX2 either follows GFS:
303 tstrc(i) = 0.5*(tskin_lnd(i) + tsurf_lnd(i)) !averaging tskin_lnd and tsurf_lnd as in GFS surface layer breaks ntsflg functionality
304 !GJF: or WRF module_sf_gfdl.F:
305 !tstrc(i) = tskin_lnd(i)
306
307 !GJF: Roughness Length Limitation section
308 ! The WRF version of module_sf_gfdl.F has no checks on the roughness lengths prior to entering MFLUX2.
309 ! The following limits were placed on roughness lengths from the GFS surface layer scheme at the suggestion
310 ! of Chunxi Zhang. Using the GFDL surface layer without such checks can lead to instability in the UFS.
311
312 !znt_lnd is in cm, z0max/ztmax are in m at this point
313 z0max(i) = max(1.0e-6, min(0.01 * znt_lnd(i), zkmax(i)))
314
315 tem1 = 1.0 - shdmax(i)
316 tem2 = tem1 * tem1
317 tem1 = 1.0 - tem2
318
319 if( ivegsrc == 1 ) then
320 if (vegtype(i) == 10) then
321 z0max(i) = exp( tem2*log01 + tem1*log07 )
322 elseif (vegtype(i) == 6) then
323 z0max(i) = exp( tem2*log01 + tem1*log05 )
324 elseif (vegtype(i) == 7) then
325 ! z0max(i) = exp( tem2*log01 + tem1*log01 )
326 z0max(i) = 0.01
327 elseif (vegtype(i) == 16) then
328 ! z0max(i) = exp( tem2*log01 + tem1*log01 )
329 z0max(i) = 0.01
330 else
331 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
332 endif
333 elseif (ivegsrc == 2 ) then
334 if (vegtype(i) == 7) then
335 z0max(i) = exp( tem2*log01 + tem1*log07 )
336 elseif (vegtype(i) == 8) then
337 z0max(i) = exp( tem2*log01 + tem1*log05 )
338 elseif (vegtype(i) == 9) then
339 ! z0max(i) = exp( tem2*log01 + tem1*log01 )
340 z0max(i) = 0.01
341 elseif (vegtype(i) == 11) then
342 ! z0max(i) = exp( tem2*log01 + tem1*log01 )
343 z0max(i) = 0.01
344 else
345 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
346 endif
347 endif
348
349 z0max(i) = max(z0max(i), 1.0e-6)
350
351 ! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil
352 czilc = 0.8
353
354 tem1 = 1.0 - sigmaf(i)
355 ztmax(i) = z0max(i)*exp( - tem1*tem1 &
356 & * czilc*karman*sqrt(ustar_lnd(i)*(0.01/1.5e-05)))
357 ztmax(i) = max(ztmax(i), 1.0e-6)
358
359 !GJF: from WRF's module_sf_gfdl.F
360 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then
361 wind10(i)=wspd(i)*alog(10.0/z0max(i))/alog(z1(i)/z0max(i)) !m s-1
362 end if
363 wind10(i)=wind10(i)*100.0 !convert from m/s to cm/s
364
365 ztmax(i) = ztmax(i)*100.0 !convert from m to cm
366 z0max(i) = z0max(i)*100.0 !convert from m to cm
367
368 call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_lnd(i:i), rib_lnd(i:i), &
369 xxfh(i), ztmax(i), z0max(i), tstrc(i), &
370 pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), &
371 scurx(i), scury(i), pert_cd, ens_random_seed, ens_cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), &
372 dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, &
373 errflg)
374 if (errflg /= 0) return
375
376 !GJF: this is broken when tstrc is set to an average of two variables
377 if (ntsflg==1) then
378 tskin_lnd(i) = tstrc(i) ! gopal's doing
379 end if
380
381 if (diag_wind10m) then
382 u10_lnd(i) = u1(i)*(0.01*wind10(i)/wspd(i))
383 v10_lnd(i) = v1(i)*(0.01*wind10(i)/wspd(i))
384 end if
385
386 !GJF: these variables are not needed in a GFS-based suite, but are found in WRF's module_sf_gfdl.F and kept in comments for legacy
387 !gz1oz0(i) = alog(zkmax(i)/(0.01*znt_lnd(i)))
388 !taux(i) = fxmx(i)/10. ! gopal's doing for Ocean coupling
389 !tauy(i) = fxmy(i)/10. ! gopal's doing for Ocean coupling
390
391 cdm_lnd(i) = max(cdm_lnd(i), cd_low_limit)
392 cdm_lnd(i) = min(cdm_lnd(i), cd_high_limit)
393 fm_lnd(i) = karman/sqrt(cdm_lnd(i))
394
395 !1) try fh_lnd from MFLUX2
396 fh_lnd(i) = karman*xxfh(i)
397
398 !2) calc ch_lnd from fm_lnd and fh_lnd
399 ch_lnd(i) = karman*karman/(fm_lnd(i) * fh_lnd(i))
400
401 !3) check if ch_lnd is out of bounds (if so, recalculate fh_lnd from bounded value)
402 ch_bound_excursion = .false.
403 if (ch_lnd(i) < ch_low_limit) then
404 ch_bound_excursion = .true.
405 ch_lnd(i) = ch_low_limit
406 else if (ch_lnd(i) > ch_high_limit) then
407 ch_bound_excursion = .true.
408 ch_lnd(i) = ch_high_limit
409 end if
410
411 fh2_lnd(i) = karman*xxfh2(i)
412
413 if (ch_bound_excursion) then
414 fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0)
415 fh_lnd(i) = karman*karman/(fm_lnd(i)*ch_lnd(i))
416 fh2_lnd(i) = fh2_fh_ratio*fh_lnd(i)
417 end if
418
419 !GJF: Other CCPP schemes (PBL) ask for fm/fh instead of psim/psih
420 !psim_lnd(i)=gz1oz0(i)-fm_lnd(i)
421 !psih_lnd(i)=gz1oz0(i)-fh_lnd(i)
422
423 !GJF: from WRF's module_sf_gfdl.F
424 ustar_lnd(i) = 0.01*sqrt(cdm_lnd(i)* &
425 (upc(i)*upc(i) + vpc(i)*vpc(i)))
426 !GJF: from Chunxi Zhang's module_sf_sfclayrev.f90 (I'm not sure it's necessary.)
427 ustar_lnd(i) = amax1(ustar_lnd(i),0.001)
428
429 stress_lnd(i) = cdm_lnd(i)*wspd(i)*wspd(i)
430
431 !GJF: from WRF's module_sf_gfdl.F
432 ! convert cd, ch to values at 10m, for output
433 cd10 = cdm_lnd(i)
434 if ( wind10(i) .ge. 0.1 ) then
435 cd10=cdm_lnd(i)* (wspd(i)/(0.01*wind10(i)) )**2
436 !tmp9=0.01*abs(tzot(i))
437 !ch_out(i)=ch_lnd(i)*(wspd(i)/(0.01*wind10(i)) ) * &
438 ! (alog(zkmax(i)/tmp9)/alog(10.0/tmp9))
439 end if
440 fm10_lnd(i) = karman/sqrt(cd10)
441
442 !GJF: conductances aren't used in other CCPP schemes, but this limit
443 ! might be able to replace the limits on drag coefficients above
444
445 !chs_lnd(i)=ch_lnd(i)*wspd (i) !conductance
446 !chs2_lnd(i)=ustar_lnd(i)*karman/fh2_lnd(i) !2m conductance
447
448 !!!2014-0922 cap CHS over land points
449 ! chs_lnd(i)=amin1(chs_lnd(i), 0.05)
450 ! chs2_lnd(i)=amin1(chs2_lnd(i), 0.05)
451 ! if (chs2_lnd(i) < 0) chs2_lnd(i)=1.0e-6
452
453 if (diag_qss) then
454 esat = fpvs(tskin_lnd(i))
455 qss_lnd(i) = ep2*esat/(psfc(i)-esat)
456 end if
457
458 !GJF: not used in CCPP
459 !flhc_lnd(i)=cpm(i)*rho1(i)*chs_lnd(i)
460 !flqc_lnd(i)=rho1(i)*chs_lnd(i)
461 !cqs2_lnd(i)=chs2_lnd(i)
462 end if !dry
463
464 if (icy(i)) then
465 !GJF: from WRF's module_sf_gfdl.F
466 smcdry=drysmc(isltyp(i))
467 smcmax=maxsmc(isltyp(i))
468 wetc(i)=(smois(i,1)-smcdry)/(smcmax-smcdry)
469 wetc(i)=amin1(1.,amax1(wetc(i),0.))
470
471
472 !GJF: the lower boundary temperature passed in to MFLUX2 either follows GFS:
473 tstrc(i) = 0.5*(tskin_ice(i) + tsurf_ice(i)) !averaging tskin_ice and tsurf_ice as in GFS surface layer breaks ntsflg functionality
474 !GJF: or WRF module_sf_gfdl.F:
475 !tstrc(i) = tskin_ice(i)
476 !averaging tskin_ice and tsurf_ice as in GFS surface layer breaks ntsflg functionality
477
478 !GJF: Roughness Length Limitation section
479 ! The WRF version of module_sf_gfdl.F has no checks on the roughness lengths prior to entering MFLUX2.
480 ! The following limits were placed on roughness lengths from the GFS surface layer scheme at the suggestion
481 ! of Chunxi Zhang. Using the GFDL surface layer without such checks can lead to instability in the UFS.
482
483 !znt_ice is in cm, z0max/ztmax are in m at this point
484 z0max(i) = max(1.0e-6, min(0.01 * znt_ice(i), zkmax(i)))
485 !** xubin's new z0 over land and sea ice
486 tem1 = 1.0 - shdmax(i)
487 tem2 = tem1 * tem1
488 tem1 = 1.0 - tem2
489
490 if( ivegsrc == 1 ) then
491 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
492 elseif (ivegsrc == 2 ) then
493 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
494 endif
495
496 z0max(i) = max(z0max(i), 1.0e-6)
497
498 ! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height
499 ! dependance of czil
500 czilc = 0.8
501
502 tem1 = 1.0 - sigmaf(i)
503 ztmax(i) = z0max(i)*exp( - tem1*tem1 &
504 & * czilc*karman*sqrt(ustar_ice(i)*(0.01/1.5e-05)))
505 ztmax(i) = max(ztmax(i), 1.0e-6)
506
507
508 !GJF: from WRF's module_sf_gfdl.F
509 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then
510 wind10(i)=wspd(i)*alog(10.0/z0max(i))/alog(z1(i)/z0max(i))
511 end if
512 wind10(i)=wind10(i)*100.0 !! m/s to cm/s
513
514 ztmax(i) = ztmax(i)*100.0 !m to cm
515 z0max(i) = z0max(i)*100.0 !m to cm
516
517 call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_ice(i:i), rib_ice(i:i), &
518 xxfh(i), ztmax(i), z0max(i), tstrc(i), &
519 pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), &
520 scurx(i), scury(i), pert_cd, ens_random_seed, ens_cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), &
521 dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, &
522 errflg)
523 if (errflg /= 0) return
524
525 !GJF: this is broken when tstrc is set to an average of two variables
526 if (ntsflg==1) then
527 tskin_ice(i) = tstrc(i) ! gopal's doing
528 end if
529
530 if (diag_wind10m) then
531 u10_ice(i) = u1(i)*(0.01*wind10(i)/wspd(i))
532 v10_ice(i) = v1(i)*(0.01*wind10(i)/wspd(i))
533 end if
534
535 !GJF: these variables are not needed in a GFS-based suite, but are found in WRF's module_sf_gfdl.F and kept in comments for legacy
536 !gz1oz0(i) = alog(zkmax(i)/znt_ice(i))
537 !taux(i) = fxmx(i)/10. ! gopal's doing for Ocean coupling
538 !tauy(i) = fxmy(i)/10. ! gopal's doing for Ocean coupling
539
540 cdm_ice(i) = max(cdm_ice(i), cd_low_limit)
541 cdm_ice(i) = min(cdm_ice(i), cd_high_limit)
542 fm_ice(i) = karman/sqrt(cdm_ice(i))
543
544 !1) try fh_ice from MFLUX2
545 fh_ice(i) = karman*xxfh(i)
546
547 !2) calc ch_ice from fm_ice and fh_ice
548 ch_ice(i) = karman*karman/(fm_ice(i) * fh_ice(i))
549
550 !3) check if ch_ice is out of bounds (if so, recalculate fh_ice from bounded value)
551 ch_bound_excursion = .false.
552 if (ch_ice(i) < ch_low_limit) then
553 ch_bound_excursion = .true.
554 ch_ice(i) = ch_low_limit
555 else if (ch_ice(i) > ch_high_limit) then
556 ch_bound_excursion = .true.
557 ch_ice(i) = ch_high_limit
558 end if
559
560 fh2_ice(i) = karman*xxfh2(i)
561
562 if (ch_bound_excursion) then
563 fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0)
564 fh_ice(i) = karman*karman/(fm_ice(i)*ch_ice(i))
565 fh2_ice(i) = fh2_fh_ratio*fh_ice(i)
566 end if
567
568 !Other CCPP schemes (PBL) ask for fm/fh instead of psim/psih
569 !psim_ice(i)=gz1oz0(i)-fm_ice(i)
570 !psih_ice(i)=gz1oz0(i)-fh_ice(i)
571
572 ustar_ice(i) = 0.01*sqrt(cdm_ice(i)* &
573 (upc(i)*upc(i) + vpc(i)*vpc(i)))
574 !GJF: from Chunxi Zhang's module_sf_sfclayrev.f90 (I'm not sure it's necessary.)
575 ustar_ice(i) = amax1(ustar_ice(i),0.001)
576
577 stress_ice(i) = cdm_ice(i)*wspd(i)*wspd(i)
578
579 !GJF: from WRF's module_sf_gfdl.F
580 !!! convert cd, ch to values at 10m, for output
581 cd10 = cdm_ice(i)
582 if ( wind10(i) .ge. 0.1 ) then
583 cd10=cdm_ice(i)* (wspd(i)/(0.01*wind10(i)) )**2
584 !tmp9=0.01*abs(tzot(i))
585 !ch_out(i)=ch_ice(i)*(wspd(i)/(0.01*wind10(i)) ) * &
586 ! (alog(zkmax(i)/tmp9)/alog(10.0/tmp9))
587 end if
588 fm10_ice(i) = karman/sqrt(cd10)
589
590 !GJF: conductances aren't used in other CCPP schemes
591 !chs_ice(i)=ch_ice(i)*wspd (i) !conductance
592 !chs2_ice(i)=ustar_ice(i)*karman/fh2_ice(i) !2m conductance
593
594 if (diag_qss) then
595 esat = fpvs(tskin_ice(i))
596 qss_ice(i) = ep2*esat/(psfc(i)-esat)
597 end if
598
599 !flhc_ice(i)=cpm(i)*rho1(i)*chs_ice(i)
600 !flqc_ice(i)=rho1(i)*chs_ice(i)
601 !cqs2_ice(i)=chs2_ice(i)
602 end if !ice
603
604 if (wet(i)) then
605 wetc(i) = 1.0
606
607 !GJF: the lower boundary temperature passed in to MFLUX2 either follows GFS:
608 tstrc(i) = 0.5*(tskin_wat(i) + tsurf_wat(i)) !averaging tskin_wat and tsurf_wat as in GFS surface layer breaks ntsflg functionality
609 !GJF: or WRF module_sf_gfdl.F:
610 !tstrc(i) = tskin_wat(i)
611
612 ! DH* 20201009: these bounds on ocean roughness lengths are from Chunxi Zhang's module_sf_sfclayrev.f90 (in cm)
613 znt_wat(i)=min(2.85e-1,max(znt_wat(i),1.27e-5))
614
615 !GJF: from WRF's module_sf_gfdl.F
616 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then
617 wind10(i)=wspd(i)*alog(10.0/(0.01*znt_wat(i)))/alog(z1(i)/(0.01*znt_wat(i)))
618 end if
619 wind10(i)=wind10(i)*100.0 !! m/s to cm/s
620
621 !GJF: mflux2 expects negative roughness length for ocean points
622 znt_wat(i) = -znt_wat(i)
623
624 call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_wat(i:i), rib_wat(i:i), &
625 xxfh(i), znt_wat(i:i), mznt(i), tstrc(i), &
626 pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), &
627 scurx(i), scury(i), pert_cd, ens_random_seed, ens_cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), &
628 dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, &
629 errflg)
630 if (errflg /= 0) return
631
632 !GJF: this is broken when tstrc is set to an average of two variables
633 if (ntsflg==1) then
634 tskin_wat(i) = tstrc(i) ! gopal's doing
635 end if
636
637 znt_wat(i)= abs(znt_wat(i))
638 mznt(i)= abs(mznt(i))
639
640 !GJF: these bounds on ocean roughness lengths are from Chunxi Zhang's module_sf_sfclayrev.f90 (in cm)
641 znt_wat(i)=min(2.85e-1,max(znt_wat(i),1.27e-5))
642
643 if (diag_wind10m) then
644 u10_ocn(i) = u1(i)*(0.01*wind10(i)/wspd(i))
645 v10_ocn(i) = v1(i)*(0.01*wind10(i)/wspd(i))
646 end if
647
648 !GJF: these variables are not needed in a GFS-based suite, but are found in WRF's module_sf_gfdl.F and kept in comments for legacy
649 !gz1oz0(i) = alog(zkmax(i)/znt_wat(i))
650 !taux(i) = fxmx(i)/10. ! gopal's doing for Ocean coupling
651 !tauy(i) = fxmy(i)/10. ! gopal's doing for Ocean coupling
652
653 cdm_wat(i) = max(cdm_wat(i), cd_low_limit)
654 cdm_wat(i) = min(cdm_wat(i), cd_high_limit)
655 fm_wat(i) = karman/sqrt(cdm_wat(i))
656
657 !1) try fh_wat from MFLUX2
658 fh_wat(i) = karman*xxfh(i)
659
660 !2) calc ch_wat from fm_wat and fh_wat
661 ch_wat(i) = karman*karman/(fm_wat(i) * fh_wat(i))
662
663 !3) check if ch_lnd is out of bounds (if so, recalculate fh_lnd from bounded value)
664 ch_bound_excursion = .false.
665 if (ch_wat(i) < ch_low_limit) then
666 ch_bound_excursion = .true.
667 ch_wat(i) = ch_low_limit
668 else if (ch_wat(i) > ch_high_limit) then
669 ch_bound_excursion = .true.
670 ch_wat(i) = ch_high_limit
671 end if
672
673 fh2_wat(i) = karman*xxfh2(i)
674
675 if (ch_bound_excursion) then
676 fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0)
677 fh_wat(i) = karman*karman/(fm_wat(i)*ch_wat(i))
678 fh2_wat(i) = fh2_fh_ratio*fh_wat(i)
679 end if
680
681 !Other CCPP schemes (PBL) ask for fm/fh instead of psim/psih
682 !psim_ocn(i)=gz1oz0(i)-fm_wat(i)
683 !psih_ocn(i)=gz1oz0(i)-fh_wat(i)
684
685 ustar_wat(i) = 0.01*sqrt(cdm_wat(i)* &
686 (upc(i)*upc(i) + vpc(i)*vpc(i)))
687 !GJF: from Chunxi Zhang's module_sf_sfclayrev.f90 (I'm not sure it's necessary.)
688 ustar_wat(i) = amax1(ustar_wat(i),0.001)
689
690 stress_wat(i) = cdm_wat(i)*wspd(i)*wspd(i)
691
692 !GJF: from WRF's module_sf_gfdl.F
693 !!! convert cd, ch to values at 10m, for output
694 cd10 = cdm_wat(i)
695 if ( wind10(i) .ge. 0.1 ) then
696 cd10=cdm_wat(i)* (wspd(i)/(0.01*wind10(i)) )**2
697 !tmp9=0.01*abs(tzot(i))
698 !ch_out(i)=ch_wat(i)*(wspd(i)/(0.01*wind10(i)) ) * &
699 ! (alog(zkmax(i)/tmp9)/alog(10.0/tmp9))
700 end if
701 fm10_wat(i) = karman/sqrt(cd10)
702
703 !GJF: conductances aren't used in other CCPP schemes
704 !chs_ocn(i)=ch_wat(i)*wspd (i) !conductance
705 !chs2_ocn(i)=ustar_wat(i)*karman/fh2_wat(i) !2m conductance
706
707 if (diag_qss) then
708 esat = fpvs(tskin_wat(i))
709 qss_wat(i) = ep2*esat/(psfc(i)-esat)
710 end if
711 end if !wet
712
713 !flhc_ocn(i)=cpm(i)*rho1(i)*chs_ocn(i)
714 !flqc_ocn(i)=rho1(i)*chs_ocn(i)
715 !cqs2_ocn(i)=chs2_ocn(i)
716 end if !flag_iter
717 end do
718
719 !GJF: this code has not been updated since GFS suites don't require this; one would need to have different values of hfx, qfx, lh for each surface type
720 ! if (isfflx.eq.0) then
721 ! do i=its,ite
722 ! hfx(i)=0.
723 ! lh(i)=0.
724 ! qfx(i)=0.
725 ! enddo
726 ! else
727 ! do i=its,ite
728 ! if(islmsk == 0) then
729 ! !water
730 ! hfx(i)= -10.*cp*fxh(i)
731 ! else if (islmsk == 1) then
732 ! hfx(i)= -10.*cp*fxh(i)
733 ! hfx(i)=amax1(hfx(i),-250.)
734 ! end if
735 ! qfx(j)=-10.*fxe(i)
736 ! qfx(i)=amax1(qfx(i),0.)
737 ! lh(i)=xlv*qfx(i)
738 ! enddo
739 ! endif
740
741
742 end subroutine gfdl_sfc_layer_run
743
744!---------------------------------
745!GJF (2020/04/21): The starting point for the MFLUX2 subroutine here was module_sf_gfdl.F in WRF
746 SUBROUTINE mflux2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc KWON
747 pspc,pkmax,wetc,slwdc,z1, &
748 icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, &
749 pert_Cd, ens_random_seed, ens_Cdamp, &
750 upc,vpc,tpc,rpc,dt,wind10,xxfh2,ntsflg,sfenth, &
751 tzot, ep2, errmsg, errflg)
752
753!------------------------------------------------------------------------
754!
755! MFLUX2 computes surface fluxes of momentum, heat,and moisture
756! using monin-obukhov. the roughness length "z0" is prescribed
757! over land and over ocean "z0" is computed using charnocks formula.
758! the universal functions (from similarity theory approach) are
759! those of hicks. This is Bob's doing.
760!
761!------------------------------------------------------------------------
762
764 IMPLICIT NONE
765
766!-----------------------------------------------------------------------
767! user interface variables
768!-----------------------------------------------------------------------
769 !GJF: This subroutine was converted to expect data from a single point instead of a horizontal array to accommodate a fractional landmask
770 !integer,intent(in) :: ims,ime
771 !integer,intent(in) :: its,ite
772 integer, parameter :: ims = 1
773 integer, parameter :: ime = 1
774 integer, parameter :: its = 1
775 integer, parameter :: ite = 1
776 integer,intent(in) :: ntsflg
777 integer,intent(in) :: icoef_sf
778 integer,intent(in) :: iwavecpl
779 logical,intent(in) :: lcurr_sf
780 logical,intent(in) :: pert_Cd
781 integer,intent(in) :: ens_random_seed
782 real(kind=kind_phys),intent(in) :: ens_cdamp
783
784 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxh
785 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxe
786 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxmx
787 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxmy
788 real(kind=kind_phys), intent (inout), dimension (ims :ime ) :: cdm
789! real, intent (out), dimension (ims :ime ) :: cdm2
790 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: rib
791 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: xxfh
792 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: xxfh2
793 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: wind10
794
795 real(kind=kind_phys), intent ( inout), dimension (ims :ime ) :: zoc,mzoc !KWON
796 real(kind=kind_phys), intent ( inout), dimension (ims :ime ) :: tzot !WANG
797 real(kind=kind_phys), intent ( inout), dimension (ims :ime ) :: tstrc
798
799 real(kind=kind_phys), intent ( in) :: dt
800 real(kind=kind_phys), intent ( in) :: sfenth
801 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: pspc
802 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: pkmax
803 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: wetc
804 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: slwdc
805 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: alpha, gamma
806 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: xcur, ycur
807 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: z1
808
809 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: upc
810 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: vpc
811 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: tpc
812 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: rpc
813
814 real(kind=kind_phys), intent ( in) :: ep2
815
816 character(len=*), intent(out) :: errmsg
817 integer, intent(out) :: errflg
818
819!-----------------------------------------------------------------------
820! internal variables
821!-----------------------------------------------------------------------
822
823 integer, parameter :: icntx = 30
824
825 integer, dimension(1 :ime) :: ifz
826 integer, dimension(1 :ime) :: indx
827 integer, dimension(1 :ime) :: istb
828 integer, dimension(1 :ime) :: it
829 integer, dimension(1 :ime) :: iutb
830
831 real(kind=kind_phys), dimension(1 :ime) :: aap
832 real(kind=kind_phys), dimension(1 :ime) :: bq1
833 real(kind=kind_phys), dimension(1 :ime) :: bq1p
834 real(kind=kind_phys), dimension(1 :ime) :: delsrad
835 real(kind=kind_phys), dimension(1 :ime) :: ecof
836 real(kind=kind_phys), dimension(1 :ime) :: ecofp
837 real(kind=kind_phys), dimension(1 :ime) :: estso
838 real(kind=kind_phys), dimension(1 :ime) :: estsop
839 real(kind=kind_phys), dimension(1 :ime) :: fmz1
840 real(kind=kind_phys), dimension(1 :ime) :: fmz10
841 real(kind=kind_phys), dimension(1 :ime) :: fmz2
842 real(kind=kind_phys), dimension(1 :ime) :: fmzo1
843 real(kind=kind_phys), dimension(1 :ime) :: foft
844 real(kind=kind_phys), dimension(1 :ime) :: foftm
845 real(kind=kind_phys), dimension(1 :ime) :: frac
846 real(kind=kind_phys), dimension(1 :ime) :: land
847 real(kind=kind_phys), dimension(1 :ime) :: pssp
848 real(kind=kind_phys), dimension(1 :ime) :: qf
849 real(kind=kind_phys), dimension(1 :ime) :: rdiff
850 real(kind=kind_phys), dimension(1 :ime) :: rho
851 real(kind=kind_phys), dimension(1 :ime) :: rkmaxp
852 real(kind=kind_phys), dimension(1 :ime) :: rstso
853 real(kind=kind_phys), dimension(1 :ime) :: rstsop
854 real(kind=kind_phys), dimension(1 :ime) :: sf10
855 real(kind=kind_phys), dimension(1 :ime) :: sf2
856 real(kind=kind_phys), dimension(1 :ime) :: sfm
857 real(kind=kind_phys), dimension(1 :ime) :: sfzo
858 real(kind=kind_phys), dimension(1 :ime) :: sgzm
859 real(kind=kind_phys), dimension(1 :ime) :: slwa
860 real(kind=kind_phys), dimension(1 :ime) :: szeta
861 real(kind=kind_phys), dimension(1 :ime) :: szetam
862 real(kind=kind_phys), dimension(1 :ime) :: t1
863 real(kind=kind_phys), dimension(1 :ime) :: t2
864 real(kind=kind_phys), dimension(1 :ime) :: tab1
865 real(kind=kind_phys), dimension(1 :ime) :: tab2
866 real(kind=kind_phys), dimension(1 :ime) :: tempa1
867 real(kind=kind_phys), dimension(1 :ime) :: tempa2
868 real(kind=kind_phys), dimension(1 :ime) :: theta
869 real(kind=kind_phys), dimension(1 :ime) :: thetap
870 real(kind=kind_phys), dimension(1 :ime) :: tsg
871 real(kind=kind_phys), dimension(1 :ime) :: tsm
872 real(kind=kind_phys), dimension(1 :ime) :: tsp
873 real(kind=kind_phys), dimension(1 :ime) :: tss
874 real(kind=kind_phys), dimension(1 :ime) :: ucom
875 real(kind=kind_phys), dimension(1 :ime) :: uf10
876 real(kind=kind_phys), dimension(1 :ime) :: uf2
877 real(kind=kind_phys), dimension(1 :ime) :: ufh
878 real(kind=kind_phys), dimension(1 :ime) :: ufm
879 real(kind=kind_phys), dimension(1 :ime) :: ufzo
880 real(kind=kind_phys), dimension(1 :ime) :: ugzm
881 real(kind=kind_phys), dimension(1 :ime) :: uzeta
882 real(kind=kind_phys), dimension(1 :ime) :: uzetam
883 real(kind=kind_phys), dimension(1 :ime) :: vcom
884 real(kind=kind_phys), dimension(1 :ime) :: vrtkx
885 real(kind=kind_phys), dimension(1 :ime) :: vrts
886 real(kind=kind_phys), dimension(1 :ime) :: wind
887 real(kind=kind_phys), dimension(1 :ime) :: windp
888 real(kind=kind_phys), dimension(1 :ime) :: wind10p !WANG, 10m wind previous step
889 real(kind=kind_phys), dimension(1 :ime) :: uvs1
890! real(kind=kind_phys), dimension(1 :ime) :: xxfh
891 real(kind=kind_phys), dimension(1 :ime) :: xxfm
892 real(kind=kind_phys), dimension(1 :ime) :: xxsh
893 real(kind=kind_phys), dimension(1 :ime) :: z10
894 real(kind=kind_phys), dimension(1 :ime) :: z2
895 real(kind=kind_phys), dimension(1 :ime) :: zeta
896 real(kind=kind_phys), dimension(1 :ime) :: zkmax
897
898 real(kind=kind_phys), dimension(1 :ime) :: pss
899 real(kind=kind_phys), dimension(1 :ime) :: tstar
900 real(kind=kind_phys), dimension(1 :ime) :: ukmax
901 real(kind=kind_phys), dimension(1 :ime) :: vkmax
902 real(kind=kind_phys), dimension(1 :ime) :: tkmax
903 real(kind=kind_phys), dimension(1 :ime) :: rkmax
904 real(kind=kind_phys), dimension(1 :ime) :: zot
905 real(kind=kind_phys), dimension(1 :ime) :: fhzo1
906 real(kind=kind_phys), dimension(1 :ime) :: sfh
907
908 real(kind=kind_phys) :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
909 real(kind=kind_phys) :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
910 real(kind=kind_phys) :: szet2, zal2,ugz2
911 real(kind=kind_phys) :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
912 real(kind=kind_phys) :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
913 real(kind=kind_phys) :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
914 real(kind=kind_phys) :: shfx,sigt4,reflect
915 real(kind=kind_phys) :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
916 real(kind=kind_phys) :: wndm,ckg
917 real(kind=kind_phys) :: windmks,znott,znotm
918 real(kind=kind_phys) :: ubot, vbot
919 integer:: i,j,ii,iq,nnest,icnt,ngd,ip
920
921!-----------------------------------------------------------------------
922! internal variables
923!-----------------------------------------------------------------------
924
925 real(kind=kind_phys), dimension (223) :: tab
926 real(kind=kind_phys), dimension (223) :: table
927 real(kind=kind_phys), dimension (101) :: tab11
928 real(kind=kind_phys), dimension (41) :: table4
929 real(kind=kind_phys), dimension (42) :: tab3
930 real(kind=kind_phys), dimension (54) :: table2
931 real(kind=kind_phys), dimension (54) :: table3
932 real(kind=kind_phys), dimension (74) :: table1
933 real(kind=kind_phys), dimension (80) :: tab22
934
935 character(len=255) :: message
936
937 equivalence(tab(1),tab11(1))
938 equivalence(tab(102),tab22(1))
939 equivalence(tab(182),tab3(1))
940 equivalence(table(1),table1(1))
941 equivalence(table(75),table2(1))
942 equivalence(table(129),table3(1))
943 equivalence(table(183),table4(1))
944
945 data amask/ -98.0/
946!-----------------------------------------------------------------------
947! tables used to obtain the vapor pressures or saturated vapor
948! pressure
949!-----------------------------------------------------------------------
950
951 data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
952 &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
953 &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
954 &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
955 &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
956 &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
957 &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
958 &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
959
960 data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
961 &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
962 &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
963 &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
964 &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
965 &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
966 &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
967 &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
968 &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
969 &190220.,199260./
970
971 data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
972 &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
973 &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
974 &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
975 &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
976 &1013250.,1049940.,1087740./
977
978 data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
979 & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
980 &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
981 &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
982 &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
983 &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
984 &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
985 &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
986 &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
987 &.7220e+01/
988
989 data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
990 &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
991 &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
992 &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
993 &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
994 &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
995 &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
996 &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
997 &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
998 &.7090e+03/
999
1000 data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
1001 &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
1002 &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
1003 &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
1004 &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
1005 &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
1006 &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
1007 &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
1008 &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
1009 &.9780e+04/
1010
1011 data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
1012 &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
1013 &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
1014 &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
1015 &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
1016 &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
1017 &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
1018!
1019! spcify constants needed by MFLUX2
1020!
1021!GJF: should send through argument list, but these have nonstandard units
1022 real,parameter :: cp = 1.00464e7
1023 real,parameter :: g = 980.6
1024 real,parameter :: rgas = 2.87e6
1025 real,parameter :: og = 1./g
1026 integer :: ntstep = 0
1027
1028 ! Initialize CCPP error handling variables
1029 errmsg = ''
1030 errflg = 0
1031!
1032#if HWRF==1
1033 real*8 :: gasdev,ran1 !zhang
1034 real :: rr !zhang
1035 logical,save :: pert_Cd_local !zhang
1036 CHARACTER(len=3) :: env_memb,env_pp
1037 integer,save :: ens_random_seed_local,env_pp_local !zhang
1038 integer :: ensda_physics_pert !zhang
1039 real,save :: ens_Cdamp_local !zhang
1040 data ens_random_seed_local/0/
1041 data env_pp_local/0/
1042 if ( ens_random_seed_local .eq. 0 ) then
1043 CALL nl_get_ensda_physics_pert(1,ensda_physics_pert)
1044 ens_random_seed_local=ens_random_seed
1045 env_pp_local=ensda_physics_pert
1046 pert_cd_local=.false.
1047 ens_cdamp_local=0.0
1048! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99
1049 if ( env_pp_local .eq. 1 ) then
1050 if ( ens_random_seed .ne. 99 ) then
1051 pert_cd_local=.true.
1052 ens_cdamp_local=ens_cdamp
1053 else
1054! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero
1055 ens_random_seed_local=ens_random_seed
1056 pert_cd_local=pert_cd
1057 ens_cdamp_local=ens_cdamp
1058 endif
1059 else
1060 ens_random_seed_local=ens_random_seed
1061 pert_cd_local=pert_cd
1062 ens_cdamp_local=ens_cdamp
1063 endif
1064 print*, "Cd ===", ens_random_seed_local,pert_cd_local,ens_cdamp_local,ensda_physics_pert
1065 endif
1066#endif
1067
1068! character*10 routine
1069! routine = 'mflux2'
1070!
1071!------------------------------------------------------------------------
1072! set water availability constant "ecof" and land mask "land".
1073! limit minimum wind speed to 100 cm/s
1074!------------------------------------------------------------------------
1075! constants for 10 m winds (correction for knots
1076!
1077 cor1 = .120
1078 cor2 = 720.
1079! KWON : remove the artificial increase of 10m wind speed over 60kts
1080! which comes from GFDL hurricane model
1081 cor1 = 0.
1082 cor2 = 0.
1083!
1084
1085 do i = its,ite
1086 z10(i) = 1000.
1087 z2(i) = 200.
1088 pss(i) = pspc(i)
1089 tstar(i) = tstrc(i)
1090
1091 if ( lcurr_sf .and. zoc(i) .le. 0.0 ) then
1092 ubot = upc(i) - xcur(i) * 100.0
1093 vbot = vpc(i) - ycur(i) * 100.0
1094! ubot = upc(i)
1095! vbot = vpc(i)
1096 else
1097 ubot = upc(i)
1098 vbot = vpc(i)
1099 endif
1100 uvs1(i)= amax1( sqrt(ubot*ubot + &
1101 vbot*vbot), 100.0)
1102 if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then
1103 ukmax(i) = ( ubot * cos(gamma(i)) - &
1104 vbot * sin(gamma(i)) ) &
1105 * cos(gamma(i))
1106 vkmax(i) = ( vbot * cos(gamma(i)) - &
1107 ubot * sin(gamma(i)) ) &
1108 * cos(gamma(i))
1109
1110 else
1111 ukmax(i) = ubot
1112 vkmax(i) = vbot
1113 endif
1114
1115! ukmax(i) = upc(i)
1116! vkmax(i) = vpc(i)
1117 tkmax(i) = tpc(i)
1118 rkmax(i) = rpc(i)
1119 enddo
1120
1121 do i = its,ite
1122 windp(i) = sqrt(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
1123 wind(i) = amax1(windp(i),100.)
1124
1125!! use wind10 previous step
1126 wind10p(i) = wind10(i) !! cm/s
1127 wind10p(i) = amax1(wind10p(i),100.)
1128!!
1129
1130 if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind10p(i)*wind10p(i)*og
1131 if (zoc(i) .GT. 0.0) then
1132 ecof(i) = wetc(i)
1133 land(i) = 1.0
1134 zot(i) = zoc(i)
1135 else
1136 ecof(i) = wetc(i)
1137 land(i) = 0.0
1138 windmks=wind10p(i)*.01
1139 if ( iwavecpl .eq. 1 ) then
1140 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg)
1141 !Check if Charnock parameter ratio is received in a proper range.
1142 if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. ) then
1143 znotm = znotm*alpha(i)
1144 endif
1145 zoc(i) = -100.*znotm
1146 zot(i) = -100* znott
1147 else
1148 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg)
1149 zoc(i) = -100.*znotm
1150 zot(i) = -100* znott
1151 endif
1152 endif
1153!------------------------------------------------------------------------
1154! where necessary modify zo values over ocean.
1155!------------------------------------------------------------------------
1156!
1157 mzoc(i) = zoc(i) !FOR SAVE MOMENTUM Zo
1158 tzot(i) = zot(i) !output wang
1159 enddo
1160
1161!------------------------------------------------------------------------
1162! define constants:
1163! a and c = constants used in evaluating universal function for
1164! stable case
1165! ca = karmen constant
1166! cm01 = constant part of vertical integral of universal
1167! function; stable case ( 0.5 < zeta < or = 10.0)
1168! cm02 = constant part of vertical integral of universal
1169! function; stable case ( zeta > 10.0)
1170!------------------------------------------------------------------------
1171
1172 en = 2.
1173 c = .76
1174 a = 5.
1175 ca = .4
1176 cmo1 = .5*a - 1.648
1177 cmo2 = 17.193 + .5*a - 10.*c
1178 boycon = .61
1179 rovcp=rgas/cp
1180
1181 do i = its,ite
1182 theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
1183 vrtkx(i) = 1.0 + boycon*rkmax(i)
1184 !zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og
1185 zkmax(i) = z1(i) !use precalculated height of first model layer center
1186 enddo
1187
1188!------------------------------------------------------------------------
1189! get saturation mixing ratios at surface
1190!------------------------------------------------------------------------
1191
1192 do i = its,ite
1193 tsg(i) = tstar(i)
1194 tab1(i) = tstar(i) - 153.16
1195 it(i) = ifix(tab1(i))
1196 tab2(i) = tab1(i) - float(it(i))
1197 t1(i) = tab(min(223,max(1,it(i) + 1)))
1198 t2(i) = table(min(223,max(1,it(i) + 1)))
1199 estso(i) = t1(i) + tab2(i)*t2(i)
1200 psps1 = (pss(i) - estso(i))
1201 if(psps1 .EQ. 0.0)then
1202 psps1 = .1
1203 endif
1204 rstso(i) = ep2*estso(i)/psps1
1205 vrts(i) = 1. + boycon*ecof(i)*rstso(i)
1206 enddo
1207
1208!------------------------------------------------------------------------
1209! check if consideration of virtual temperature changes stability.
1210! if so, set "dthetav" to near neutral value (1.0e-4). also check
1211! for very small lapse rates; if ABS(tempa1) <1.0e-4 then
1212! tempa1=1.0e-4
1213!------------------------------------------------------------------------
1214
1215 do i = its,ite
1216 tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
1217 tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
1218 if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
1219 tab1(i) = abs(tempa1(i))
1220 if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
1221!------------------------------------------------------------------------
1222! compute bulk richardson number "rib" at each point. if "rib"
1223! exceeds 95% of critical richardson number "tab1" then "rib = tab1"
1224!------------------------------------------------------------------------
1225
1226 rib(i) = g*zkmax(i)*tempa1(i)/ &
1227 (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
1228 tab2(i) = abs(zoc(i))
1229 tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
1230 if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
1231 enddo
1232
1233 do i = its,ite
1234 zeta(i) = ca*rib(i)/0.03
1235 enddo
1236
1237!------------------------------------------------------------------------
1238! begin looping through points on line, solving wegsteins iteration
1239! for zeta at each point, and using hicks functions
1240!------------------------------------------------------------------------
1241
1242!------------------------------------------------------------------------
1243! set initial guess of zeta=non - dimensional height "szeta" for
1244! stable points
1245!------------------------------------------------------------------------
1246
1247 rca = 1./ca
1248 enrca = en*rca
1249! turn off interfacial layer by zeroing out enrca
1250 enrca = 0.0
1251 zog = .0185*og
1252
1253!------------------------------------------------------------------------
1254! stable points
1255!------------------------------------------------------------------------
1256
1257 ip = 0
1258 do i = its,ite
1259 if (zeta(i) .GE. 0.0) then
1260 ip = ip + 1
1261 istb(ip) = i
1262 endif
1263 enddo
1264
1265 if (ip .EQ. 0) go to 170
1266 do i = 1,ip
1267 szetam(i) = 1.0e+30
1268 sgzm(i) = 0.0e+00
1269 szeta(i) = zeta(istb(i))
1270 ifz(i) = 1
1271 enddo
1272
1273!------------------------------------------------------------------------
1274! begin wegstein iteration for "zeta" at stable points using
1275! hicks(1976)
1276!------------------------------------------------------------------------
1277
1278 do icnt = 1,icntx
1279 do i = 1,ip
1280 if (ifz(i) .EQ. 0) go to 80
1281 zal1g = alog(szeta(i))
1282 if (szeta(i) .LE. 0.5) then
1283 fmz1(i) = (zal1g + a*szeta(i))*rca
1284 else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
1285 rzeta1 = 1./szeta(i)
1286 fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
1287 0.5*rzeta1*rzeta1 + cmo1)*rca
1288 else if (szeta(i) .GT. 10.) then
1289 fmz1(i) = (c*szeta(i) + cmo2)*rca
1290 endif
1291 szetao = abs(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
1292 zal2g = alog(szetao)
1293 if (szetao .LE. 0.5) then
1294 fmzo1(i) = (zal2g + a*szetao)*rca
1295 sfzo(i) = 1. + a*szetao
1296 else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
1297 rzeta2 = 1./szetao
1298 fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
1299 0.5*rzeta2*rzeta2 + cmo1)*rca
1300 sfzo(i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1301 else if (szetao .GT. 10.) then
1302 fmzo1(i) = (c*szetao + cmo2)*rca
1303 sfzo(i) = c*szetao
1304 endif
1305
1306
1307! compute heat & moisture parts of zot.. for calculation of sfh
1308
1309 szetho = abs(zot(istb(i)))/zkmax(istb(i))*szeta(i)
1310 zal2gh = alog(szetho)
1311 if (szetho .LE. 0.5) then
1312 fhzo1(i) = (zal2gh + a*szetho)*rca
1313 sfzo(i) = 1. + a*szetho
1314 else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
1315 rzeta2 = 1./szetho
1316 fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
1317 0.5*rzeta2*rzeta2 + cmo1)*rca
1318 sfzo(i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1319 else if (szetho .GT. 10.) then
1320 fhzo1(i) = (c*szetho + cmo2)*rca
1321 sfzo(i) = c*szetho
1322 endif
1323
1324!------------------------------------------------------------------------
1325! compute universal function at 10 meters for diagnostic purposes
1326!------------------------------------------------------------------------
1327
1328 szet10 = abs(z10(istb(i)))/zkmax(istb(i))*szeta(i)
1329 zal10 = alog(szet10)
1330 if (szet10 .LE. 0.5) then
1331 fmz10(i) = (zal10 + a*szet10)*rca
1332 else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
1333 rzeta2 = 1./szet10
1334 fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
1335 0.5*rzeta2*rzeta2 + cmo1)*rca
1336 else if (szet10 .GT. 10.) then
1337 fmz10(i) = (c*szet10 + cmo2)*rca
1338 endif
1339 sf10(i) = fmz10(i) - fmzo1(i)
1340! compute 2m values for diagnostics in HWRF
1341 szet2 = abs(z2(istb(i)))/zkmax(istb(i))*szeta(i)
1342 zal2 = alog(szet2 )
1343 if (szet2 .LE. 0.5) then
1344 fmz2(i) = (zal2 + a*szet2 )*rca
1345 else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
1346 rzeta2 = 1./szet2
1347 fmz2(i) = (8.*zal2 + 4.25*rzeta2 - &
1348 0.5*rzeta2*rzeta2 + cmo1)*rca
1349 else if (szet2 .GT. 2.) then
1350 fmz2(i) = (c*szet2 + cmo2)*rca
1351 endif
1352 sf2(i) = fmz2(i) - fmzo1(i)
1353
1354 sfm(i) = fmz1(i) - fmzo1(i)
1355 sfh(i) = fmz1(i) - fhzo1(i)
1356 sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
1357 (sfh(i) + enrca*sfzo(i))
1358 fmz = (sgz - szeta(i))/szeta(i)
1359 fmzo = abs(fmz)
1360 if (fmzo .GE. 5.0e-5) then
1361 sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
1362 if(sq .EQ. 1) then
1363 write(errmsg,'(*(a))') 'NCO ERROR DIVIDE BY ZERO IN gfdl_sfc_layer.F90/MFLUX2 (STABLE CASE)'// &
1364 'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
1365 errflg = 1
1366 return
1367 endif
1368 szetam(i) = szeta(i)
1369 szeta(i) = (sgz - szeta(i)*sq)/(1.0 - sq)
1370 sgzm(i) = sgz
1371 else
1372 ifz(i) = 0
1373 endif
137480 continue
1375 enddo
1376 enddo
1377
1378 do i = 1,ip
1379 if (ifz(i) .GE. 1) go to 110
1380 enddo
1381
1382 go to 130
1383
1384110 continue
1385
1386 write(errmsg,'(*(a))') 'NON-CONVERGENCE FOR STABLE ZETA IN gfdl_sfc_layer.F90/MFLUX2'
1387 errflg = 1
1388 return
1389! call MPI_CLOSE(1,routine)
1390
1391!------------------------------------------------------------------------
1392! update "zo" for ocean points. "zo"cannot be updated within the
1393! wegsteins iteration as the scheme (for the near neutral case)
1394! can become unstable
1395!------------------------------------------------------------------------
1396
1397130 continue
1398 do i = 1,ip
1399 szo = zoc(istb(i))
1400 if (szo .LT. 0.0) then
1401 wndm=wind(istb(i))*0.01
1402 if(wndm.lt.15.0) then
1403 ckg=0.0185*og
1404 else
1405 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1406 endif
1407
1408 szo = - ckg*wind(istb(i))*wind(istb(i))/ &
1409 (sfm(i)*sfm(i))
1410 cons_p000001 = .000001
1411 cons_7 = 7.
1412 vis = 1.4e-1
1413
1414 ustar = sqrt( -szo / zog)
1415 restar = -ustar * szo / vis
1416 restar = max(restar,cons_p000001)
1417! Rat taken from Zeng, Zhao and Dickinson 1997
1418 rat = 2.67 * restar ** .25 - 2.57
1419 rat = min(rat ,cons_7) !constant
1420 rat=0.
1421 zot(istb(i)) = szo * exp(-rat)
1422 else
1423 zot(istb(i)) = zoc(istb(i))
1424 endif
1425
1426! in hwrf thermal znot is loaded back into the zoc array for next step
1427 zoc(istb(i)) = szo
1428 enddo
1429
1430 do i = 1,ip
1431 xxfm(istb(i)) = sfm(i)
1432 xxfh(istb(i)) = sfh(i)
1433 xxfh2(istb(i)) = sf2(i)
1434 xxsh(istb(i)) = sfzo(i)
1435 enddo
1436
1437!------------------------------------------------------------------------
1438! obtain wind at 10 meters for diagnostic purposes
1439!------------------------------------------------------------------------
1440
1441 do i = 1,ip
1442 wind10(istb(i)) = sf10(i)*uvs1(istb(i))/sfm(i)
1443 wind10(istb(i)) = wind10(istb(i)) * 1.944
1444 if(wind10(istb(i)) .GT. 6000.0) then
1445 wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
1446 - cor2
1447 endif
1448! the above correction done by GFDL in centi-kts!!!-change back
1449 wind10(istb(i)) = wind10(istb(i)) / 1.944
1450 enddo
1451
1452!------------------------------------------------------------------------
1453! unstable points
1454!------------------------------------------------------------------------
1455
1456170 continue
1457
1458 iq = 0
1459 do i = its,ite
1460 if (zeta(i) .LT. 0.0) then
1461 iq = iq + 1
1462 iutb(iq) = i
1463 endif
1464 enddo
1465
1466 if (iq .EQ. 0) go to 290
1467 do i = 1,iq
1468 uzeta(i) = zeta(iutb(i))
1469 ifz(i) = 1
1470 uzetam(i) = 1.0e+30
1471 ugzm(i) = 0.0e+00
1472 enddo
1473
1474!------------------------------------------------------------------------
1475! begin wegstein iteration for "zeta" at unstable points using
1476! hicks functions
1477!------------------------------------------------------------------------
1478
1479 do icnt = 1,icntx
1480 do i = 1,iq
1481 if (ifz(i) .EQ. 0) go to 200
1482 ugzzo = alog(zkmax(iutb(i))/abs(zot(iutb(i))))
1483 uzetao = abs(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1484 ux11 = 1. - 16.*uzeta(i)
1485 ux12 = 1. - 16.*uzetao
1486 y = sqrt(ux11)
1487 yo = sqrt(ux12)
1488 ufzo(i) = 1./yo
1489 ux13 = (1. + y)/(1. + yo)
1490 ux21 = alog(ux13)
1491 ufh(i) = (ugzzo - 2.*ux21)*rca
1492! recompute scalers for ufm in terms of mom znot... zoc
1493 ugzzo = alog(zkmax(iutb(i))/abs(zoc(iutb(i))))
1494 uzetao = abs(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1495 ux11 = 1. - 16.*uzeta(i)
1496 ux12 = 1. - 16.*uzetao
1497 y = sqrt(ux11)
1498 yo = sqrt(ux12)
1499 ux13 = (1. + y)/(1. + yo)
1500 ux21 = alog(ux13)
1501! ufzo(i) = 1./yo
1502 x = sqrt(y)
1503 xo = sqrt(yo)
1504 xnum = (x**2 + 1.)*((x + 1.)**2)
1505 xden = (xo**2 + 1.)*((xo + 1.)**2)
1506 xtan = atan(x) - atan(xo)
1507 ux3 = alog(xnum/xden)
1508 ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca
1509
1510!------------------------------------------------------------------------
1511! obtain ten meter winds for diagnostic purposes
1512!------------------------------------------------------------------------
1513
1514 ugz10 = alog(z10(iutb(i))/abs(zoc(iutb(i))))
1515 uzet1o = abs(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1516 uzetao = abs(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1517 ux11 = 1. - 16.*uzet1o
1518 ux12 = 1. - 16.*uzetao
1519 y = sqrt(ux11)
1520 y10 = sqrt(ux12)
1521 ux13 = (1. + y)/(1. + y10)
1522 ux21 = alog(ux13)
1523 x = sqrt(y)
1524 x10 = sqrt(y10)
1525 xnum = (x**2 + 1.)*((x + 1.)**2)
1526 xden = (x10**2 + 1.)*((x10 + 1.)**2)
1527 xtan = atan(x) - atan(x10)
1528 ux3 = alog(xnum/xden)
1529 uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca
1530
1531! obtain 2m values for diagnostics...
1532
1533
1534 ugz2 = alog(z2(iutb(i))/abs(zoc(iutb(i))))
1535 uzet1o = abs(z2(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1536 uzetao = abs(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1537 ux11 = 1. - 16.*uzet1o
1538 ux12 = 1. - 16.*uzetao
1539 y = sqrt(ux11)
1540 yo = sqrt(ux12)
1541 ux13 = (1. + y)/(1. + yo)
1542 ux21 = alog(ux13)
1543 uf2(i) = (ugzzo - 2.*ux21)*rca
1544
1545
1546 ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
1547 ux1 = (ugz - uzeta(i))/uzeta(i)
1548 ux2 = abs(ux1)
1549 if (ux2 .GE. 5.0e-5) then
1550 uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
1551 uzetam(i) = uzeta(i)
1552 if(uq .EQ. 1) then
1553 write(errmsg,'(*(a))') 'NCO ERROR DIVIDE BY ZERO IN gfdl_sfc_layer.F90/MFLUX2 (UNSTABLE CASE)'// &
1554 'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1555 errflg = 1
1556 return
1557 endif
1558 uzeta(i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
1559 ugzm(i) = ugz
1560 else
1561 ifz(i) = 0
1562 endif
1563200 continue
1564 enddo
1565 enddo
1566
1567
1568 do i = 1,iq
1569 if (ifz(i) .GE. 1) go to 230
1570 enddo
1571
1572 go to 250
1573
1574230 continue
1575 write(errmsg,'(*(a))') 'NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW'// &
1576 'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1577 errflg = 1
1578 return
1579
1580! call MPI_CLOSE(1,routine)
1581
1582!------------------------------------------------------------------------
1583! gather unstable values
1584!------------------------------------------------------------------------
1585
1586250 continue
1587
1588!------------------------------------------------------------------------
1589! update "zo" for ocean points. zo cannot be updated within the
1590! wegsteins iteration as the scheme (for the near neutral case)
1591! can become unstable.
1592!------------------------------------------------------------------------
1593
1594 do i = 1,iq
1595 uzo = zoc(iutb(i))
1596 if (zoc(iutb(i)) .LT. 0.0) then
1597 wndm=wind(iutb(i))*0.01
1598 if(wndm.lt.15.0) then
1599 ckg=0.0185*og
1600 else
1601 ckg=(4*0.000308*wndm)*og
1602 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1603 endif
1604 uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
1605 cons_p000001 = .000001
1606 cons_7 = 7.
1607 vis = 1.4e-1
1608
1609 ustar = sqrt( -uzo / zog)
1610 restar = -ustar * uzo / vis
1611 restar = max(restar,cons_p000001)
1612! Rat taken from Zeng, Zhao and Dickinson 1997
1613 rat = 2.67 * restar ** .25 - 2.57
1614 rat = min(rat ,cons_7) !constant
1615 rat=0.0
1616 zot(iutb(i)) = uzo * exp(-rat)
1617 else
1618 zot(iutb(i)) = zoc(iutb(i))
1619 endif
1620! in hwrf thermal znot is loaded back into the zoc array for next step
1621 zoc(iutb(i)) = uzo
1622 enddo
1623
1624!------------------------------------------------------------------------
1625! obtain wind at ten meters for diagnostic purposes
1626!------------------------------------------------------------------------
1627 do i = 1,iq
1628 wind10(iutb(i)) = uf10(i)*uvs1(iutb(i))/ufm(i)
1629 wind10(iutb(i)) = wind10(iutb(i)) * 1.944
1630 if(wind10(iutb(i)) .GT. 6000.0) then
1631 wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
1632 - cor2
1633 endif
1634! the above correction done by GFDL in centi-kts!!!-change back
1635 wind10(iutb(i)) = wind10(iutb(i)) / 1.944
1636 enddo
1637
1638 do i = 1,iq
1639 xxfm(iutb(i)) = ufm(i)
1640 xxfh(iutb(i)) = ufh(i)
1641 xxfh2(iutb(i)) = uf2(i)
1642 xxsh(iutb(i)) = ufzo(i)
1643 enddo
1644
1645290 continue
1646
1647 do i = its,ite
1648 ucom(i) = ukmax(i)
1649 vcom(i) = vkmax(i)
1650 if (windp(i) .EQ. 0.0) then
1651 windp(i) = 100.0
1652 ucom(i) = 100.0/sqrt(2.0)
1653 vcom(i) = 100.0/sqrt(2.0)
1654 endif
1655 rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - &
1656 tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i))))
1657 bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i)))
1658 enddo
1659
1660! do land sfc temperature prediction if ntsflg=1
1661! ntsflg = 1 ! gopal's doing
1662
1663 if (ntsflg .EQ. 0) go to 370
1664 alll = 600.
1665 xks = 0.01
1666 hcap = .5/2.39e-8
1667 pith = sqrt(4.*atan(1.0))
1668 alfus = alll/2.39e-8
1669 teps = 0.1
1670! slwdc... in units of cal/min ????
1671! slwa... in units of ergs/sec/cm*2
1672! 1 erg=2.39e-8 cal
1673!------------------------------------------------------------------------
1674! pack land and sea ice points
1675!------------------------------------------------------------------------
1676
1677 ip = 0
1678 do i = its,ite
1679 if (land(i) .EQ. 1) then
1680 ip = ip + 1
1681 indx(ip) = i
1682! slwa is defined as positive down....
1683 slwa(ip) = slwdc(i)/(2.39e-8*60.)
1684 tss(ip) = tstar(i)
1685 thetap(ip) = theta(i)
1686 rkmaxp(ip) = rkmax(i)
1687 aap(ip) = 5.673e-5
1688 pssp(ip) = pss(i)
1689 ecofp(ip) = ecof(i)
1690 estsop(ip) = estso(i)
1691 rstsop(ip) = rstso(i)
1692 bq1p(ip) = bq1(i)
1693 bq1p(ip) = amax1(bq1p(ip),0.1e-3)
1694 delsrad(ip) = dt *pith/(hcap*sqrt(3600.*24.*xks))
1695 endif
1696 enddo
1697
1698!------------------------------------------------------------------------
1699! initialize variables for first pass of iteration
1700!------------------------------------------------------------------------
1701
1702 do i = 1,ip
1703 ifz(i) = 1
1704 tsm(i) = tss(i)
1705 rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1706
1707300 format(2x, ' SURFACE EQUILIBRIUM CALCULATION ')
1708
1709 foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - &
1710 cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1711 rdiff(i))
1712 tsp(i) = foftm(i)
1713 enddo
1714
1715!------------------------------------------------------------------------
1716! do iteration to determine "tstar" at new time level
1717!------------------------------------------------------------------------
1718
1719 do icnt = 1,icntx
1720 do i = 1,ip
1721 if (ifz(i) .EQ. 0) go to 330
1722 tab1(i) = tsp(i) - 153.16
1723 it(i) = ifix(tab1(i))
1724 tab2(i) = tab1(i) - float(it(i))
1725 t1(i) = tab(min(223,max(1,it(i) + 1)))
1726 t2(i) = table(min(223,max(1,it(i) + 1)))
1727 estsop(i) = t1(i) + tab2(i)*t2(i)
1728 psps2 = (pssp(i) - estsop(i))
1729 if(psps2 .EQ. 0.0)then
1730 psps2 = .1
1731 endif
1732 rstsop(i) = ep2*estsop(i)/psps2
1733 rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1734
1735 foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - &
1736 cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1737 rdiff(i))
1738
1739 frac(i) = abs((foft(i) - tsp(i))/tsp(i))
1740
1741!------------------------------------------------------------------------
1742! check for convergence of all points use wegstein iteration
1743!------------------------------------------------------------------------
1744
1745 if (frac(i) .GE. teps) then
1746 qf(i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
1747 tsm(i) = tsp(i)
1748 tsp(i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
1749 foftm(i) = foft(i)
1750 else
1751 ifz(i) = 0
1752 endif
1753330 continue
1754 enddo
1755 enddo
1756
1757!------------------------------------------------------------------------
1758! check for convergence of "t star" prediction
1759!------------------------------------------------------------------------
1760
1761 do i = 1,ip
1762 if (ifz(i) .EQ. 1) then
1763 write(errmsg,'(*(a))') 'NON-CONVERGENCE OF T* PREDICTED (T*,I) = ', &
1764 tsp(i), i
1765 errflg = 1
1766 return
1767! call MPI_CLOSE(1,routine)
1768 endif
1769 enddo
1770
1771 do i = 1,ip
1772 ii = indx(i)
1773 tstrc(ii) = tsp(i)
1774 enddo
1775
1776!------------------------------------------------------------------------
1777! compute fluxes and momentum drag coef
1778!------------------------------------------------------------------------
1779
1780370 continue
1781 do i = its,ite
1782!!!
1783 if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then
1784 windmks = wind10(i) * 0.01
1785 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg)
1786 !Check if Charnock parameter ratio is received in a proper range.
1787 if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. ) then
1788 znotm = znotm*alpha(i)
1789 endif
1790 zoc(i) = -100.*znotm
1791 zot(i) = -100* znott
1792 endif
1793!!!!
1794 fxh(i) = bq1(i)*(theta(i) - tsg(i))
1795 fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i))
1796 if (fxe(i) .GT. 0.0) fxe(i) = 0.0
1797 fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ &
1798 windp(i)
1799 fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
1800 windp(i)
1801 cdm(i) = 1./(xxfm(i)*xxfm(i))
1802#if HWRF==1
1803! randomly perturb the Cd
1804!zzz if( pert_Cd_local .and. ens_random_seed_local .gt. 0 ) then
1805 if( pert_cd_local ) then
1806 ens_random_seed_local=ran1(-ens_random_seed_local)*1000
1807 rr=2.0*ens_cdamp_local*ran1(-ens_random_seed_local)-ens_cdamp_local
1808 cdm(i) = cdm(i) *(1.0+rr)
1809 endif
1810#endif
1811
1812 enddo
1813 ntstep = ntstep + 1
1814 return
1815 end subroutine mflux2
1816
1817 end module gfdl_sfc_layer
subroutine esat(t, esw, esi, desw, desi)
use polynomials to calculate saturation vapor pressure and derivative with respect to temperature: ov...
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains the CCPP-compliant GFDL surface layer scheme.
This module contains the CCPP-compliant Noah land surface scheme driver.
Definition lsm_noah.f:5
Data from MPTABLE.TBL, SOILPARM.TBL, GENPARM.TBL for NoahMP.