CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
GFS_surface_composites_post.F90
1
3
5
6 use machine, only: kind_phys
7
8 ! For consistent calculations of composite surface properties
9 use sfc_diff, only: stability
10
11 implicit none
12
13 private
14
15 public gfs_surface_composites_post_run
16
17 real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys, &
18 half = 0.5_kind_phys, qmin = 1.0e-8_kind_phys
19
20contains
21
25 subroutine gfs_surface_composites_post_run ( &
26 im, kice, km, rd, rvrdm1, cplflx, cplwav2atm, frac_grid, flag_cice, thsfc_loc, islmsk, dry, wet, icy, wind, t1, q1, prsl1, &
27 landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli, garea, frac_ice, &
28 cd, cd_wat, cd_lnd, cd_ice, cdq, cdq_wat, cdq_lnd, cdq_ice, rb, rb_wat, rb_lnd, rb_ice, stress, stress_wat, stress_lnd, &
29 stress_ice, ffmm, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh, ffhh_wat, ffhh_lnd, ffhh_ice, uustar, uustar_wat, uustar_lnd, &
30 uustar_ice, fm10, fm10_wat, fm10_lnd, fm10_ice, fh2, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice, &
31 cmm, cmm_wat, cmm_lnd, cmm_ice, chh, chh_wat, chh_lnd, chh_ice, gflx, gflx_wat, gflx_lnd, gflx_ice, ep1d, ep1d_wat, &
32 ep1d_lnd, ep1d_ice, weasd, weasd_lnd, weasd_ice, snowd, snowd_lnd, snowd_ice, tprcp, tprcp_wat, &
33 tprcp_lnd, tprcp_ice, evap, evap_wat, evap_lnd, evap_ice, hflx, hflx_wat, hflx_lnd, hflx_ice, qss, qss_wat, qss_lnd, &
34 qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tisfc, hice, cice, tiice, &
35 sigmaf, zvfun, lheatstrg, h0facu, h0facs, hflxq, hffac, stc, lkm, iopt_lake, iopt_lake_clm, use_lake_model, &
36 grav, prsik1, prslk1, prslki, z1, ztmax_wat, ztmax_lnd, ztmax_ice, huge, errmsg, errflg)
37
38 implicit none
39
40 integer, intent(in) :: im, kice, km, lkm, iopt_lake, iopt_lake_clm
41 logical, intent(in) :: cplflx, frac_grid, cplwav2atm, frac_ice
42 logical, intent(in) :: lheatstrg
43 logical, dimension(:), intent(in) :: flag_cice, dry, icy
44 logical, dimension(:), intent(in) :: wet
45 integer, dimension(:), intent(in) :: islmsk, use_lake_model
46 real(kind=kind_phys), dimension(:), intent(in) :: wind, t1, q1, prsl1, landfrac, lakefrac, oceanfrac, &
47 cd_wat, cd_lnd, cd_ice, cdq_wat, cdq_lnd, cdq_ice, rb_wat, rb_lnd, rb_ice, stress_wat, &
48 stress_lnd, stress_ice, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh_wat, ffhh_lnd, ffhh_ice, uustar_wat, uustar_lnd, uustar_ice, &
49 fm10_wat, fm10_lnd, fm10_ice, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice, cmm_wat, cmm_lnd, cmm_ice, &
50 chh_wat, chh_lnd, chh_ice, gflx_wat, gflx_lnd, gflx_ice, ep1d_wat, ep1d_lnd, ep1d_ice, weasd_lnd, weasd_ice, &
51 snowd_lnd, snowd_ice, tprcp_wat, tprcp_lnd, tprcp_ice, evap_wat, evap_lnd, evap_ice, hflx_wat, hflx_lnd, &
52 hflx_ice, qss_wat, qss_lnd, qss_ice, tsfc_wat, zorlo, zorll, zorli, garea
53
54 real(kind=kind_phys), dimension(:), intent(inout) :: zorl, cd, cdq, rb, stress, ffmm, ffhh, uustar, fm10, &
55 fh2, cmm, chh, gflx, ep1d, weasd, snowd, tprcp, evap, hflx, qss, tsfc, tsfco, tsfcl, tisfc
56
57 real(kind=kind_phys), dimension(:), intent(inout) :: hice, cice
58 real(kind=kind_phys), dimension(:), intent(inout) :: sigmaf, zvfun, hflxq, hffac
59 real(kind=kind_phys), intent(in ) :: h0facu, h0facs
60! real(kind=kind_phys), intent(in ) :: min_seaice
61 real(kind=kind_phys), intent(in ) :: rd, rvrdm1, huge
62
63 real(kind=kind_phys), dimension(:,:), intent(in ) :: tiice
64 real(kind=kind_phys), dimension(:,:), intent(inout) :: stc
65
66 ! Additional data needed for calling "stability"
67 logical, intent(in ) :: thsfc_loc
68 real(kind=kind_phys), intent(in ) :: grav
69 real(kind=kind_phys), dimension(:), intent(in ) :: prsik1, prslk1, prslki, z1
70 real(kind=kind_phys), dimension(:), intent(in ) :: ztmax_wat, ztmax_lnd, ztmax_ice
71
72 character(len=*), intent(out) :: errmsg
73 integer, intent(out) :: errflg
74
75 ! Local variables
76 integer :: i, k
77 real(kind=kind_phys) :: txl, txi, txo, wfrac, q0, rho
78 ! For calling "stability"
79 real(kind=kind_phys) :: tsurf, virtfac, tv1, thv1, tvs, z0max, ztmax
80 real(kind=kind_phys) :: lnzorll, lnzorli, lnzorlo
81!
82 real(kind=kind_phys) :: tem1, tem2, gdx
83 real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0
84!
85
86 ! Initialize CCPP error handling variables
87 errmsg = ''
88 errflg = 0
89
90 ! --- generate ocean/land/ice composites
91
92 fractional_grid: if (frac_grid) then
93
94 do i=1, im
95
96 ! Three-way composites (fields from sfc_diff)
97 txl = landfrac(i) ! land fraction
98 wfrac = one - txl ! ocean fraction
99 txi = cice(i) * wfrac ! txi = ice fraction wrt whole cell
100 txo = max(zero, wfrac-txi) ! txo = open water fraction
101
102 !gflx(i) = txl*gflx_lnd(i) + txi*gflx_ice(i) + txo*gflx_wat(i)
103 ep1d(i) = txl*ep1d_lnd(i) + txi*ep1d_ice(i) + txo*ep1d_wat(i)
104 weasd(i) = txl*weasd_lnd(i) + txi*weasd_ice(i)
105 snowd(i) = txl*snowd_lnd(i) + txi*snowd_ice(i)
106 !tprcp(i) = txl*tprcp_lnd(i) + txi*tprcp_ice(i) + txo*tprcp_wat(i)
107!
108 sigmaf(i) = txl*sigmaf(i)
109
110 evap(i) = txl*evap_lnd(i) + txi*evap_ice(i) + txo*evap_wat(i)
111 hflx(i) = txl*hflx_lnd(i) + txi*hflx_ice(i) + txo*hflx_wat(i)
112 qss(i) = txl*qss_lnd(i) + txi*qss_ice(i) + txo*qss_wat(i)
113 gflx(i) = txl*gflx_lnd(i) + txi*gflx_ice(i) + txo*gflx_wat(i)
114
115!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116! Call stability for consistent surface properties. Currently this comes from !
117! the GFS surface layere scheme (sfc_diff), regardless of the actual surface !
118! layer parameterization being used - to be extended in the future !
119!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
120
121 tsfc(i) = ( txl * cdq_lnd(i) * tsfcl(i) &
122 + txi * cdq_ice(i) * tisfc(i) & ! DH* Ben had tsurf_ice(i), but GFS_surface_composites_post_run uses tice instead
123 + txo * cdq_wat(i) * tsfc_wat(i)) &
124 / (txl * cdq_lnd(i) + txi * cdq_ice(i) + txo * cdq_wat(i) )
125 tsurf = ( txl * cdq_lnd(i) * tsurf_lnd(i) &
126 + txi * cdq_ice(i) * tsurf_ice(i) &
127 + txo * cdq_wat(i) * tsurf_wat(i)) &
128 / (txl * cdq_lnd(i) + txi * cdq_ice(i) + txo * cdq_wat(i) )
129
130 q0 = max( q1(i), qmin )
131 virtfac = one + rvrdm1 * q0
132
133 tv1 = t1(i) * virtfac ! Virtual temperature in middle of lowest layer
134 if(thsfc_loc) then ! Use local potential temperature
135 thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level
136 tvs = half * (tsfc(i)+tsurf) * virtfac
137 else ! Use potential temperature referenced to 1000 hPa
138 thv1 = t1(i) / prslk1(i) * virtfac ! Theta-v at lowest level
139 tvs = half * (tsfc(i)+tsurf)/prsik1(i) * virtfac
140 endif
141
142 lnzorll = zero ; lnzorli = zero ; lnzorlo = zero
143 if (zorll(i) /= huge) then
144 lnzorll = log(zorll(i))
145 endif
146 if (zorli(i) /= huge) then
147 lnzorli = log(zorli(i))
148 endif
149 if (zorlo(i) /= huge) then
150 lnzorlo = log(zorlo(i))
151 endif
152 zorl(i) = exp(txl*lnzorll + txi*lnzorli + txo*lnzorlo)
153 ! zorl(i) = exp(txl*log(zorll(i)) + txi*log(zorli(i)) + txo*log(zorlo(i)))
154 z0max = 0.01_kind_phys * zorl(i)
155 ztmax = exp(txl*log(ztmax_lnd(i)) + txi*log(ztmax_ice(i)) + txo*log(ztmax_wat(i)))
156
157 ! Only actually need to call "stability" if multiple surface types exist...
158 if(txl == one) then ! 100% land
159 rb(i) = rb_lnd(i)
160 ffmm(i) = ffmm_lnd(i)
161 ffhh(i) = ffhh_lnd(i)
162 fm10(i) = fm10_lnd(i)
163 fh2(i) = fh2_lnd(i)
164 cd(i) = cd_lnd(i)
165 cdq(i) = cdq_lnd(i)
166 stress(i) = stress_lnd(i)
167 uustar(i) = uustar_lnd(i)
168 elseif(txo == one) then ! 100% open water
169 rb(i) = rb_wat(i)
170 ffmm(i) = ffmm_wat(i)
171 ffhh(i) = ffhh_wat(i)
172 fm10(i) = fm10_wat(i)
173 fh2(i) = fh2_wat(i)
174 cd(i) = cd_wat(i)
175 cdq(i) = cdq_wat(i)
176 stress(i) = stress_wat(i)
177 uustar(i) = uustar_wat(i)
178 elseif(txi == one) then ! 100% ice
179 rb(i) = rb_ice(i)
180 ffmm(i) = ffmm_ice(i)
181 ffhh(i) = ffhh_ice(i)
182 fm10(i) = fm10_ice(i)
183 fh2(i) = fh2_ice(i)
184 cd(i) = cd_ice(i)
185 cdq(i) = cdq_ice(i)
186 stress(i) = stress_ice(i)
187 uustar(i) = uustar_ice(i)
188 else ! Mix of multiple surface types (land, water, and/or ice)
189!
190! re-compute zvfun with composite surface roughness & green vegetation fraction
191!
192 tem1 = (z0max - z0lo) / (z0up - z0lo)
193 tem1 = min(max(tem1, zero), one)
194 tem2 = max(sigmaf(i), 0.1)
195 zvfun(i) = sqrt(tem1 * tem2)
196 gdx = sqrt(garea(i))
197!
198! re-compute variables for canopy heat storage parameterization with the updated zvfun
199! in the fractional grid
200!
201 hflxq(i) = hflx(i)
202 hffac(i) = 1.0
203 if (lheatstrg) then
204 if(hflx(i) > 0.) then
205 hffac(i) = h0facu * zvfun(i)
206 else
207 hffac(i) = h0facs * zvfun(i)
208 endif
209 hffac(i) = 1. + hffac(i)
210 hflxq(i) = hflx(i) / hffac(i)
211 endif
212!
213 call stability(z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & ! inputs
214 z0max, ztmax, tvs, grav, thsfc_loc, & ! inputs
215 rb(i), ffmm(i), ffhh(i), fm10(i), fh2(i), cd(i), cdq(i), & ! outputs
216 stress(i), uustar(i))
217 endif ! Checking to see if point is one or multiple surface types
218
219 ! BWG, 2021/02/25: cmm=cd*wind, chh=cdq*wind, so use composite cd, cdq
220 rho = prsl1(i) / (rd*tv1)
221 cmm(i) = cd(i)*wind(i) !txl*cmm_lnd(i) + txi*cmm_ice(i) + txo*cmm_wat(i)
222 chh(i) = rho*cdq(i)*wind(i) !txl*chh_lnd(i) + txi*chh_ice(i) + txo*chh_wat(i)
223
224!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
225
226 if (dry(i)) then
227 elseif (wet(i)) then
228 tsfcl(i) = tsfc_wat(i) ! over water
229 else
230 tsfcl(i) = tisfc(i) ! over ice
231 endif
232 if (wet(i)) then
233 tsfco(i) = tsfc_wat(i) ! over lake or ocean when uncoupled
234 elseif (icy(i)) then
235 tsfco(i) = tisfc(i) ! over lake or ocean ice when uncoupled
236 else
237 tsfco(i) = tsfcl(i) ! over land
238 endif
239 if (icy(i)) then
240 !tisfc(i) = tisfc(i) ! over lake or ocean ice when uncoupled
241 elseif (wet(i)) then
242 tisfc(i) = tsfc_wat(i) ! over lake or ocean when uncoupled
243 else
244 tisfc(i) = tsfcl(i) ! over land
245 endif
246 ! for coupled model ocean will replace this
247! endif
248
249! if (.not. flag_cice(i)) then
250! if (islmsk(i) == 2) then ! return updated lake ice thickness & concentration to global array
251! tisfc(i) = tice(i)
252! else ! this would be over open ocean or land (no ice fraction)
253! hice(i) = zero
254! cice(i) = zero
255! tisfc(i) = tsfc(i)
256! endif
257! endif
258 if (.not. icy(i)) then
259 hice(i) = zero
260 cice(i) = zero
261 endif
262 enddo
263
264 else ! not fractional grid
265
266 do i=1,im
267
268 if (use_lake_model(i)>0) then
269 if(frac_ice .and. icy(i)) then
270 call composite_icy(iopt_lake==iopt_lake_clm)
271 call composite_wet_and_icy
272 else
273 call composite_wet
274 endif
275 else if (islmsk(i) == 1) then
276 !-- land
277 call composite_land
278 elseif (islmsk(i) == 0) then
279 !-- water
280 call composite_wet
281 else ! islmsk(i) == 2
282 !-- ice
283 call composite_icy(.false.)
284 call composite_wet_and_icy
285 endif
286 enddo
287
288 endif fractional_grid
289
290 ! --- compositing done
291
292 contains
293
294 subroutine composite_land
295 implicit none
296 zorl(i) = zorll(i)
297 cd(i) = cd_lnd(i)
298 cdq(i) = cdq_lnd(i)
299 rb(i) = rb_lnd(i)
300 stress(i) = stress_lnd(i)
301 ffmm(i) = ffmm_lnd(i)
302 ffhh(i) = ffhh_lnd(i)
303 uustar(i) = uustar_lnd(i)
304 fm10(i) = fm10_lnd(i)
305 fh2(i) = fh2_lnd(i)
306 tsfc(i) = tsfcl(i)
307 tsfco(i) = tsfc(i)
308 tisfc(i) = tsfc(i)
309 cmm(i) = cmm_lnd(i)
310 chh(i) = chh_lnd(i)
311 gflx(i) = gflx_lnd(i)
312 ep1d(i) = ep1d_lnd(i)
313 weasd(i) = weasd_lnd(i)
314 snowd(i) = snowd_lnd(i)
315 evap(i) = evap_lnd(i)
316 hflx(i) = hflx_lnd(i)
317 qss(i) = qss_lnd(i)
318 hice(i) = zero
319 cice(i) = zero
320 end subroutine composite_land
321
322 subroutine composite_wet
323 implicit none
324 zorl(i) = zorlo(i)
325 cd(i) = cd_wat(i)
326 cdq(i) = cdq_wat(i)
327 rb(i) = rb_wat(i)
328 stress(i) = stress_wat(i)
329 ffmm(i) = ffmm_wat(i)
330 ffhh(i) = ffhh_wat(i)
331 uustar(i) = uustar_wat(i)
332 fm10(i) = fm10_wat(i)
333 fh2(i) = fh2_wat(i)
334 tsfco(i) = tsfc_wat(i) ! over lake (and ocean when uncoupled)
335 tsfc(i) = tsfco(i)
336 tsfcl(i) = tsfc(i)
337 tisfc(i) = tsfc(i)
338 cmm(i) = cmm_wat(i)
339 chh(i) = chh_wat(i)
340 gflx(i) = gflx_wat(i)
341 ep1d(i) = ep1d_wat(i)
342 weasd(i) = zero
343 snowd(i) = zero
344 evap(i) = evap_wat(i)
345 hflx(i) = hflx_wat(i)
346 qss(i) = qss_wat(i)
347 hice(i) = zero
348 cice(i) = zero
349 end subroutine composite_wet
350
351 subroutine composite_icy(is_clm)
352 implicit none
353 logical, intent(in) :: is_clm
354 zorl(i) = zorli(i)
355 cd(i) = cd_ice(i)
356 cdq(i) = cdq_ice(i)
357 rb(i) = rb_ice(i)
358 ffmm(i) = ffmm_ice(i)
359 ffhh(i) = ffhh_ice(i)
360 uustar(i) = uustar_ice(i)
361 fm10(i) = fm10_ice(i)
362 fh2(i) = fh2_ice(i)
363 stress(i) = stress_ice(i)
364 cmm(i) = cmm_ice(i)
365 chh(i) = chh_ice(i)
366 gflx(i) = gflx_ice(i)
367 ep1d(i) = ep1d_ice(i)
368 if(is_clm) then
369 weasd(i) = weasd_ice(i)
370 snowd(i) = snowd_ice(i)
371 else
372 weasd(i) = weasd_ice(i) * cice(i)
373 snowd(i) = snowd_ice(i) * cice(i)
374 endif
375 qss(i) = qss_ice(i)
376 evap(i) = evap_ice(i)
377 hflx(i) = hflx_ice(i)
378 end subroutine composite_icy
379
380 subroutine composite_wet_and_icy
381 implicit none
382 txi = cice(i)
383 txo = one - txi
384 evap(i) = txi * evap_ice(i) + txo * evap_wat(i)
385 hflx(i) = txi * hflx_ice(i) + txo * hflx_wat(i)
386 tsfc(i) = txi * tisfc(i) + txo * tsfc_wat(i)
387 stress(i) = txi * stress_ice(i) + txo * stress_wat(i)
388 qss(i) = txi * qss_ice(i) + txo * qss_wat(i)
389 ep1d(i) = txi * ep1d_ice(i) + txo * ep1d_wat(i)
390
391 lnzorli = zero ; lnzorlo = zero
392 if (zorli(i) /= huge) then
393 lnzorli = log(zorli(i))
394 endif
395 if (zorlo(i) /= huge) then
396 lnzorlo = log(zorlo(i))
397 endif
398 zorl(i) = exp(txi*lnzorli + txo*lnzorlo)
399! zorl(i) = exp(txi*log(zorli(i)) + txo*log(zorlo(i)))
400!
401 if (wet(i)) then
402 tsfco(i) = tsfc_wat(i)
403 else
404 tsfco(i) = tsfc(i)
405 endif
406 tsfcl(i) = tsfc(i)
407 do k=1,min(kice,km) ! store tiice in stc to reduce output in the nonfrac grid case
408 stc(i,k) = tiice(i,k)
409 enddo
410 end subroutine composite_wet_and_icy
411
412 end subroutine gfs_surface_composites_post_run
413
subroutine, public stability(z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, thsfc_loc, rb, fm, fh, fm10, fh2, cm, ch, stress, ustar)
Definition sfc_diff.f:481
This module contains the CCPP-compliant GFS surface layer scheme.
Definition sfc_diff.f:7