CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
fv_sat_adj.F90
1
5!***********************************************************************
6!* GNU Lesser General Public License
7!*
8!* This file is part of the GFDL Cloud Microphysics.
9!*
10!* The GFDL Cloud Microphysics is free software: you can
11!8 redistribute it and/or modify it under the terms of the
12!* GNU Lesser General Public License as published by the
13!* Free Software Foundation, either version 3 of the License, or
14!* (at your option) any later version.
15!*
16!* The GFDL Cloud Microphysics is distributed in the hope it will be
17!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
18!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
19!* See the GNU General Public License for more details.
20!*
21!* You should have received a copy of the GNU Lesser General Public
22!* License along with the GFDL Cloud Microphysics.
23!* If not, see <http://www.gnu.org/licenses/>.
24!***********************************************************************
25
28module fv_sat_adj
29! Modules Included:
30! <table>
31! <tr>
32! <th>Module Name</th>
33! <th>Functions Included</th>
34! </tr>
35! <tr>
36! <td>constants_mod</td>
37! <td>rvgas, rdgas, grav, hlv, hlf, cp_air</td>
38! </tr>
39! <tr>
40! <td>fv_arrays_mod</td>
41! <td> r_grid</td>
42! </tr>
43! <tr>
44! <tr>
45! <td>fv_mp_mod</td>
46! <td>is_master</td>
47! </tr>
48! <tr>
49! <td>gfdl_cloud_microphys_mod</td>
50! <td>ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt,
51! tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r,
52! rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs</td>
53! </tr>
54! </table>
55 ! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH
56 use physcons, only : rdgas => con_rd_dyn, &
57 rvgas => con_rv_dyn, &
58 grav => con_g_dyn, &
59 hlv => con_hvap_dyn, &
60 hlf => con_hfus_dyn, &
61 cp_air => con_cp_dyn
62 ! *DH
63 use machine, only: kind_grid, kind_dyn
64 use gfdl_cloud_microphys_mod, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt
65 use gfdl_cloud_microphys_mod, only: icloud_f, sat_adj0, t_sub, cld_min
66 use gfdl_cloud_microphys_mod, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r
67 use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
68#ifdef MULTI_GASES
69 use ccpp_multi_gases_mod, only: multi_gases_init, &
70 multi_gases_finalize, &
71 virq_qpz, vicpqd_qpz, &
72 vicvqd_qpz, num_gas
73#endif
74
75 implicit none
76
77 private
78
79 public fv_sat_adj_init, fv_sat_adj_run, fv_sat_adj_finalize
80
81 logical :: is_initialized = .false.
82
83 real(kind=kind_dyn), parameter :: rrg = -rdgas/grav
84 ! real, parameter :: cp_air = cp_air ! 1004.6, heat capacity of dry air at constant pressure, come from constants_mod
85 real(kind=kind_dyn), parameter :: cp_vap = 4.0 * rvgas
86 real(kind=kind_dyn), parameter :: cv_air = cp_air - rdgas
87 real(kind=kind_dyn), parameter :: cv_vap = 3.0 * rvgas
88 ! http: / / www.engineeringtoolbox.com / ice - thermal - properties - d_576.html
89 ! c_ice = 2050.0 at 0 deg c
90 ! c_ice = 1972.0 at - 15 deg c
91 ! c_ice = 1818.0 at - 40 deg c
92 ! http: / / www.engineeringtoolbox.com / water - thermal - properties - d_162.html
93 ! c_liq = 4205.0 at 4 deg c
94 ! c_liq = 4185.5 at 15 deg c
95 ! c_liq = 4178.0 at 30 deg c
96 ! real, parameter :: c_ice = 2106.0 ! ifs: heat capacity of ice at 0 deg c
97 ! real, parameter :: c_liq = 4218.0 ! ifs: heat capacity of liquid at 0 deg c
98 real(kind=kind_dyn), parameter :: c_ice = 1972.0
99 real(kind=kind_dyn), parameter :: c_liq = 4185.5
100 real(kind=kind_dyn), parameter :: dc_vap = cp_vap - c_liq
101 real(kind=kind_dyn), parameter :: dc_ice = c_liq - c_ice
102 real(kind=kind_dyn), parameter :: tice = 273.16
103 real(kind=kind_dyn), parameter :: t_wfr = tice - 40.
104 real(kind=kind_dyn), parameter :: lv0 = hlv - dc_vap * tice
105 real(kind=kind_dyn), parameter :: li00 = hlf - dc_ice * tice
106 ! real (kind_grid), parameter :: e00 = 610.71 ! gfdl: saturation vapor pressure at 0 deg c
107 real (kind_grid), parameter :: e00 = 611.21
108 real (kind_grid), parameter :: d2ice = dc_vap + dc_ice
109 real (kind_grid), parameter :: li2 = lv0 + li00
110 real(kind=kind_dyn), parameter :: lat2 = (hlv + hlf) ** 2
111 real(kind=kind_dyn) :: d0_vap
112 real(kind=kind_dyn) :: lv00
113 real(kind=kind_dyn), allocatable :: table(:), table2(:), tablew(:), des2(:), desw(:)
114
115contains
116
121subroutine fv_sat_adj_init(do_sat_adj, kmp, nwat, ngas, rilist, cpilist, &
122 mpirank, mpiroot, errmsg, errflg)
123
124 implicit none
125
126 ! Interface variables
127 logical, intent(in ) :: do_sat_adj
128 integer, intent(in ) :: kmp
129 integer, intent(in ) :: nwat
130 integer, intent(in ) :: ngas
131 real(kind_dyn), intent(in ) :: rilist(0:ngas)
132 real(kind_dyn), intent(in ) :: cpilist(0:ngas)
133 integer, intent(in ) :: mpirank
134 integer, intent(in ) :: mpiroot
135 character(len=*), intent( out) :: errmsg
136 integer, intent( out) :: errflg
137
138 ! Local variables
139 integer, parameter :: length = 2621
140 integer :: i
141
142 ! Initialize the CCPP error handling variables
143 errmsg = ''
144 errflg = 0
145
146 ! If saturation adjustment is not used, return immediately
147 if (.not.do_sat_adj) then
148 write(errmsg,'(a)') 'Logic error: fv_sat_adj_init is called but do_sat_adj is set to false'
149 errflg = 1
150 return
151 end if
152
153 if (.not.nwat==6) then
154 write(errmsg,'(a)') 'Logic error: fv_sat_adj requires six water species (nwat=6)'
155 errflg = 1
156 return
157 end if
158
159 if (is_initialized) return
160
161 ! generate es table (dt = 0.1 deg c)
162
163 allocate (table(length))
164 allocate (table2(length))
165 allocate (tablew(length))
166 allocate (des2(length))
167 allocate (desw(length))
168
169 call qs_table (length)
170 call qs_table2 (length)
171 call qs_tablew (length)
172
173 do i = 1, length - 1
174 des2(i) = max(0., table2(i + 1) - table2(i))
175 desw(i) = max(0., tablew(i + 1) - tablew(i))
176 enddo
177 des2(length) = des2(length - 1)
178 desw(length) = desw(length - 1)
179
180#ifdef MULTI_GASES
181 call multi_gases_init(ngas,nwat,rilist,cpilist,mpirank==mpiroot)
182#endif
183
184 is_initialized = .true.
185
186end subroutine fv_sat_adj_init
187
188!\ingroup fast_sat_adj
193subroutine fv_sat_adj_finalize (errmsg, errflg)
194
195 implicit none
196
197 character(len=*), intent(out) :: errmsg
198 integer, intent(out) :: errflg
199
200 ! Initialize the CCPP error handling variables
201 errmsg = ''
202 errflg = 0
203
204 if (.not.is_initialized) return
205
206 if (allocated(table )) deallocate(table )
207 if (allocated(table2)) deallocate(table2)
208 if (allocated(tablew)) deallocate(tablew)
209 if (allocated(des2 )) deallocate(des2 )
210 if (allocated(desw )) deallocate(desw )
211
212#ifdef MULTI_GASES
213 call multi_gases_finalize()
214#endif
215
216 is_initialized = .false.
217
218end subroutine fv_sat_adj_finalize
219
233subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, isc1, iec1, isc2, iec2, kmp, km, kmdelz, js, je, jsd, jed, jsc1, jec1, jsc2, jec2, &
234 ng, hydrostatic, fast_mp_consv, te0_2d, te0, ngas, qvi, qv, ql, qi, qr, &
235 qs, qg, hs, peln, delz, delp, pt, pkz, q_con, akap, cappa, area, dtdt, &
236 out_dt, last_step, do_qa, qa, &
237 nthreads, errmsg, errflg)
238
239 implicit none
240
241 ! Interface variables
242 real(kind=kind_dyn), intent(in) :: mdt
243 real(kind=kind_dyn), intent(in) :: zvir
244 integer, intent(in) :: is
245 integer, intent(in) :: ie
246 integer, intent(in) :: isd
247 integer, intent(in) :: ied
248 integer, intent(in) :: isc1
249 integer, intent(in) :: iec1
250 integer, intent(in) :: isc2
251 integer, intent(in) :: iec2
252 integer, intent(in) :: kmp
253 integer, intent(in) :: km
254 integer, intent(in) :: kmdelz
255 integer, intent(in) :: js
256 integer, intent(in) :: je
257 integer, intent(in) :: jsd
258 integer, intent(in) :: jed
259 integer, intent(in) :: jsc1
260 integer, intent(in) :: jec1
261 integer, intent(in) :: jsc2
262 integer, intent(in) :: jec2
263 integer, intent(in) :: ng
264 logical, intent(in) :: hydrostatic
265 logical, intent(in) :: fast_mp_consv
266 real(kind=kind_dyn), intent(inout) :: te0_2d(is:ie, js:je)
267 real(kind=kind_dyn), intent( out) :: te0(isd:ied, jsd:jed, 1:km)
268 ! If multi-gases physics are not used, ngas is one and qvi identical to qv
269 integer, intent(in) :: ngas
270#ifdef MULTI_GASES
271 real(kind=kind_dyn), intent(inout), optional :: qvi(isd:ied, jsd:jed, 1:km, 1:ngas)
272#else
273 real(kind=kind_dyn), intent(inout), optional :: qvi(:,:,:,:)
274#endif
275 real(kind=kind_dyn), intent(inout) :: qv(isd:ied, jsd:jed, 1:km)
276 real(kind=kind_dyn), intent(inout) :: ql(isd:ied, jsd:jed, 1:km)
277 real(kind=kind_dyn), intent(inout) :: qi(isd:ied, jsd:jed, 1:km)
278 real(kind=kind_dyn), intent(inout) :: qr(isd:ied, jsd:jed, 1:km)
279 real(kind=kind_dyn), intent(inout) :: qs(isd:ied, jsd:jed, 1:km)
280 real(kind=kind_dyn), intent(inout) :: qg(isd:ied, jsd:jed, 1:km)
281 real(kind=kind_dyn), intent(in) :: hs(isd:ied, jsd:jed)
282 real(kind=kind_dyn), intent(in) :: peln(is:ie, 1:km+1, js:je)
283 ! For hydrostatic build, kmdelz=1, otherwise kmdelz=km (see fv_arrays.F90)
284 real(kind=kind_dyn), intent(in) :: delz(is:ie, js:je, 1:kmdelz)
285 real(kind=kind_dyn), intent(in) :: delp(isd:ied, jsd:jed, 1:km)
286 real(kind=kind_dyn), intent(inout) :: pt(isd:ied, jsd:jed, 1:km)
287 real(kind=kind_dyn), intent(inout) :: pkz(is:ie, js:je, 1:km)
288#ifdef USE_COND
289 real(kind=kind_dyn), intent(inout) :: q_con(isd:ied, jsd:jed, 1:km)
290#else
291 real(kind=kind_dyn), intent(inout) :: q_con(isd:isd, jsd:jsd, 1)
292#endif
293 real(kind=kind_dyn), intent(in) :: akap
294#ifdef MOIST_CAPPA
295 real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1:km)
296#else
297 real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1)
298#endif
299 ! DH* WARNING, allocation in fv_arrays.F90 is area(isd_2d:ied_2d, jsd_2d:jed_2d),
300 ! where normally isd_2d = isd etc, but for memory optimization, these can be set
301 ! to isd_2d=0, ied_2d=-1 etc. I don't believe this optimization is actually used,
302 ! as it would break a whole lot of code (including the one below)!
303 ! Assume thus that isd_2d = isd etc.
304 real(kind_grid), intent(in) :: area(isd:ied, jsd:jed)
305 real(kind=kind_dyn), intent(inout) :: dtdt(is:ie, js:je, 1:km)
306 logical, intent(in) :: out_dt
307 logical, intent(in) :: last_step
308 logical, intent(in) :: do_qa
309 real(kind=kind_dyn), intent( out) :: qa(isd:ied, jsd:jed, 1:km)
310 integer, intent(in) :: nthreads
311 character(len=*), intent( out) :: errmsg
312 integer, intent( out) :: errflg
313
314 ! Local variables
315 real(kind=kind_dyn), dimension(is:ie,js:je) :: dpln
316 integer :: kdelz
317 integer :: k, j, i
318
319 ! Initialize the CCPP error handling variables
320 errmsg = ''
321 errflg = 0
322
323#ifndef FV3
324! Open parallel region if not already opened by host model
325!$OMP parallel num_threads(nthreads) default(none) &
326!$OMP shared(kmp,km,js,je,is,ie,peln,mdt, &
327!$OMP isd,jsd,delz,q_con,cappa,qa, &
328!$OMP do_qa,last_step,out_dt,dtdt, &
329!$OMP area,delp,pt,hs,qg,qs,qr,qi, &
330!$OMP ql,qv,te0,fast_mp_consv, &
331!$OMP hydrostatic,ng,zvir,pkz, &
332!$OMP akap,te0_2d,ngas,qvi) &
333!$OMP private(k,j,i,kdelz,dpln)
334#endif
335
336!$OMP do
337 do k=kmp,km
338 do j=js,je
339 do i=is,ie
340 dpln(i,j) = peln(i,k+1,j) - peln(i,k,j)
341 enddo
342 enddo
343 if (hydrostatic) then
344 kdelz = 1
345 else
346 kdelz = k
347 end if
348 call fv_sat_adj_work(abs(mdt), zvir, is, ie, js, je, ng, hydrostatic, fast_mp_consv, &
349 te0(isd,jsd,k), &
350#ifdef MULTI_GASES
351 qvi(isd,jsd,k,1:ngas), &
352#else
353 qv(isd,jsd,k), &
354#endif
355 ql(isd,jsd,k), qi(isd,jsd,k), &
356 qr(isd,jsd,k), qs(isd,jsd,k), qg(isd,jsd,k), &
357 hs, dpln, delz(is:,js:,kdelz), pt(isd,jsd,k), delp(isd,jsd,k),&
358 q_con(isd:,jsd:,k), cappa(isd:,jsd:,k), area, dtdt(is,js,k), &
359 out_dt, last_step, do_qa, qa(isd,jsd,k))
360 if ( .not. hydrostatic ) then
361 do j=js,je
362 do i=is,ie
363#ifdef MOIST_CAPPA
364 pkz(i,j,k) = exp(cappa(i,j,k)*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
365#else
366#ifdef MULTI_GASES
367 pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
368#else
369 pkz(i,j,k) = exp(akap*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
370#endif
371#endif
372 enddo
373 enddo
374 endif
375 enddo
376!$OMP end do
377
378 if ( fast_mp_consv ) then
379!$OMP do
380 do j=js,je
381 do i=is,ie
382 do k=kmp,km
383 te0_2d(i,j) = te0_2d(i,j) + te0(i,j,k)
384 enddo
385 enddo
386 enddo
387!$OMP end do
388 endif
389
390#ifndef FV3
391!$OMP end parallel
392#endif
393
394 return
395
396end subroutine fv_sat_adj_run
397
402subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0, &
403#ifdef MULTI_GASES
404 qvi, &
405#else
406 qv, &
407#endif
408 ql, qi, qr, qs, qg, hs, dpln, delz, pt, dp, q_con, cappa, &
409 area, dtdt, out_dt, last_step, do_qa, qa)
410
411 implicit none
412
413 ! Interface variables
414 integer, intent (in) :: is, ie, js, je, ng
415 logical, intent (in) :: hydrostatic, consv_te, out_dt, last_step, do_qa
416 real(kind=kind_dyn), intent (in) :: zvir, mdt ! remapping time step
417 real(kind=kind_dyn), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: dp, hs
418 real(kind=kind_dyn), intent (in), dimension (is:ie, js:je) :: dpln, delz
419 real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: pt
420#ifdef MULTI_GASES
421 real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng, 1:1, 1:num_gas) :: qvi
422#else
423 real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: qv
424#endif
425 real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: ql, qi, qr, qs, qg
426 real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: q_con, cappa
427 real(kind=kind_dyn), intent (inout), dimension (is:ie, js:je) :: dtdt
428 real(kind=kind_dyn), intent (out), dimension (is - ng:ie + ng, js - ng:je + ng) :: qa, te0
429 real (kind_grid), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: area
430
431 ! Local variables
432#ifdef MULTI_GASES
433 real, dimension (is - ng:ie + ng, js - ng:je + ng) :: qv
434#endif
435 real(kind=kind_dyn), dimension (is:ie) :: wqsat, dq2dt, qpz, cvm, t0, pt1, qstar
436 real(kind=kind_dyn), dimension (is:ie) :: icp2, lcp2, tcp2, tcp3
437 real(kind=kind_dyn), dimension (is:ie) :: den, q_liq, q_sol, q_cond, src, sink, hvar
438 real(kind=kind_dyn), dimension (is:ie) :: mc_air, lhl, lhi
439 real(kind=kind_dyn) :: qsw, rh
440 real(kind=kind_dyn) :: tc, qsi, dqsdt, dq, dq0, pidep, qi_crt, tmp, dtmp
441 real(kind=kind_dyn) :: tin, rqi, q_plus, q_minus
442 real(kind=kind_dyn) :: sdt, dt_bigg, adj_fac
443 real(kind=kind_dyn) :: fac_smlt, fac_r2g, fac_i2s, fac_imlt, fac_l2r, fac_v2l, fac_l2v
444 real(kind=kind_dyn) :: factor, qim, tice0, c_air, c_vap, dw
445 integer :: i, j
446
447#ifdef MULTI_GASES
448 qv(:,:) = qvi(:,:,1,1)
449#endif
450 sdt = 0.5 * mdt ! half remapping time step
451 dt_bigg = mdt ! bigg mechinism time step
452
453 tice0 = tice - 0.01 ! 273.15, standard freezing temperature
454
455 ! -----------------------------------------------------------------------
457 ! -----------------------------------------------------------------------
458
459 fac_i2s = 1. - exp(- mdt / tau_i2s)
460 fac_v2l = 1. - exp(- sdt / tau_v2l)
461 fac_r2g = 1. - exp(- mdt / tau_r2g)
462 fac_l2r = 1. - exp(- mdt / tau_l2r)
463
464 fac_l2v = 1. - exp(- sdt / tau_l2v)
465 fac_l2v = min(sat_adj0, fac_l2v)
466
467 fac_imlt = 1. - exp(- sdt / tau_imlt)
468 fac_smlt = 1. - exp(- mdt / tau_smlt)
469
470 ! -----------------------------------------------------------------------
472 ! -----------------------------------------------------------------------
473
474 if (hydrostatic) then
475 c_air = cp_air
476 c_vap = cp_vap
477 else
478 c_air = cv_air
479 c_vap = cv_vap
480 endif
481 d0_vap = c_vap - c_liq
482 lv00 = hlv - d0_vap * tice
483 ! dc_vap = cp_vap - c_liq ! - 2339.5
484 ! d0_vap = cv_vap - c_liq ! - 2801.0
485
486 do j = js, je ! start j loop
487
488 do i = is, ie
489 q_liq(i) = ql(i, j) + qr(i, j)
490 q_sol(i) = qi(i, j) + qs(i, j) + qg(i, j)
491 qpz(i) = q_liq(i) + q_sol(i)
492#ifdef MULTI_GASES
493 pt1(i) = pt(i, j) / virq_qpz(qvi(i,j,1,1:num_gas),qpz(i))
494#else
495#ifdef USE_COND
496 pt1(i) = pt(i, j) / ((1 + zvir * qv(i, j)) * (1 - qpz(i)))
497#else
498 pt1(i) = pt(i, j) / (1 + zvir * qv(i, j))
499#endif
500#endif
501 t0(i) = pt1(i) ! true temperature
502 qpz(i) = qpz(i) + qv(i, j) ! total_wat conserved in this routine
503 enddo
504
505 ! -----------------------------------------------------------------------
507 ! -----------------------------------------------------------------------
508
509 if (hydrostatic) then
510 do i = is, ie
511 den(i) = dp(i, j) / (dpln(i, j) * rdgas * pt(i, j))
512 enddo
513 else
514 do i = is, ie
515 den(i) = - dp(i, j) / (grav * delz(i, j)) ! moist_air density
516 enddo
517 endif
518
519 ! -----------------------------------------------------------------------
521 ! -----------------------------------------------------------------------
522
523 do i = is, ie
524#ifdef MULTI_GASES
525 if (hydrostatic) then
526 c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
527 else
528 c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
529 endif
530#endif
531 mc_air(i) = (1. - qpz(i)) * c_air ! constant
532 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
533 lhi(i) = li00 + dc_ice * pt1(i)
534 icp2(i) = lhi(i) / cvm(i)
535 enddo
536
537 ! -----------------------------------------------------------------------
539 ! -----------------------------------------------------------------------
540
541 if (consv_te) then
542 if (hydrostatic) then
543 do i = is, ie
544#ifdef MULTI_GASES
545 c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
546#endif
547 te0(i, j) = - c_air * t0(i)
548 enddo
549 else
550 do i = is, ie
551#ifdef USE_COND
552 te0(i, j) = - cvm(i) * t0(i)
553#else
554#ifdef MULTI_GASES
555 c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
556#endif
557 te0(i, j) = - c_air * t0(i)
558#endif
559 enddo
560 endif
561 endif
562
563 ! -----------------------------------------------------------------------
565 ! -----------------------------------------------------------------------
566
567 do i = is, ie
568 if (qi(i, j) < 0.) then
569 qs(i, j) = qs(i, j) + qi(i, j)
570 qi(i, j) = 0.
571 endif
572 enddo
573
574 ! -----------------------------------------------------------------------
576 ! -----------------------------------------------------------------------
577
578 do i = is, ie
579 if (qi(i, j) > 1.e-8 .and. pt1(i) > tice) then
580 sink(i) = min(qi(i, j), fac_imlt * (pt1(i) - tice) / icp2(i))
581 qi(i, j) = qi(i, j) - sink(i)
582 ! sjl, may 17, 2017
583 ! tmp = min (sink (i), dim (ql_mlt, ql (i, j))) ! max ql amount
584 ! ql (i, j) = ql (i, j) + tmp
585 ! qr (i, j) = qr (i, j) + sink (i) - tmp
586 ! sjl, may 17, 2017
587 ql(i, j) = ql(i, j) + sink(i)
588 q_liq(i) = q_liq(i) + sink(i)
589 q_sol(i) = q_sol(i) - sink(i)
590 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
591 pt1(i) = pt1(i) - sink(i) * lhi(i) / cvm(i)
592 endif
593 enddo
594
595 ! -----------------------------------------------------------------------
597 ! -----------------------------------------------------------------------
598
599 do i = is, ie
600 lhi(i) = li00 + dc_ice * pt1(i)
601 icp2(i) = lhi(i) / cvm(i)
602 enddo
603
604 ! -----------------------------------------------------------------------
606 ! -----------------------------------------------------------------------
607
608 do i = is, ie
609 if (qs(i, j) < 0.) then
610 qg(i, j) = qg(i, j) + qs(i, j)
611 qs(i, j) = 0.
612 elseif (qg(i, j) < 0.) then
613 tmp = min(- qg(i, j), max(0., qs(i, j)))
614 qg(i, j) = qg(i, j) + tmp
615 qs(i, j) = qs(i, j) - tmp
616 endif
617 enddo
618
619 ! after this point cloud ice & snow are positive definite
620
621 ! -----------------------------------------------------------------------
623 ! -----------------------------------------------------------------------
624
625 do i = is, ie
626 if (ql(i, j) < 0.) then
627 tmp = min(- ql(i, j), max(0., qr(i, j)))
628 ql(i, j) = ql(i, j) + tmp
629 qr(i, j) = qr(i, j) - tmp
630 elseif (qr(i, j) < 0.) then
631 tmp = min(- qr(i, j), max(0., ql(i, j)))
632 ql(i, j) = ql(i, j) - tmp
633 qr(i, j) = qr(i, j) + tmp
634 endif
635 enddo
636
637 ! -----------------------------------------------------------------------
639 ! -----------------------------------------------------------------------
640
641 do i = is, ie
642 dtmp = tice - 48. - pt1(i)
643 if (ql(i, j) > 0. .and. dtmp > 0.) then
644 sink(i) = min(ql(i, j), dtmp / icp2(i))
645 ql(i, j) = ql(i, j) - sink(i)
646 qi(i, j) = qi(i, j) + sink(i)
647 q_liq(i) = q_liq(i) - sink(i)
648 q_sol(i) = q_sol(i) + sink(i)
649 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
650 pt1(i) = pt1(i) + sink(i) * lhi(i) / cvm(i)
651 endif
652 enddo
653
654 ! -----------------------------------------------------------------------
656 ! -----------------------------------------------------------------------
657
658 do i = is, ie
659 lhl(i) = lv00 + d0_vap * pt1(i)
660 lhi(i) = li00 + dc_ice * pt1(i)
661 lcp2(i) = lhl(i) / cvm(i)
662 icp2(i) = lhi(i) / cvm(i)
663 tcp3(i) = lcp2(i) + icp2(i) * min(1., dim(tice, pt1(i)) /48.)
664 enddo
665
666 ! -----------------------------------------------------------------------
668 ! -----------------------------------------------------------------------
669
670 call wqs2_vect (is, ie, pt1, den, wqsat, dq2dt)
671
672 adj_fac = sat_adj0
673 do i = is, ie
674 dq0 = (qv(i, j) - wqsat(i)) / (1. + tcp3(i) * dq2dt(i))
675 if (dq0 > 0.) then ! whole grid - box saturated
676 src(i) = min(adj_fac * dq0, max(ql_gen - ql(i, j), fac_v2l * dq0))
677 else ! evaporation of ql
678 ! sjl 20170703 added ql factor to prevent the situation of high ql and rh<1
679 ! factor = - min (1., fac_l2v * sqrt (max (0., ql (i, j)) / 1.e-5) * 10. * (1. - qv (i, j) / wqsat (i)))
680 ! factor = - fac_l2v
681 ! factor = - 1
682 factor = - min(1., fac_l2v * 10. * (1. - qv(i, j) / wqsat(i))) ! the rh dependent factor = 1 at 90%
683 src(i) = - min(ql(i, j), factor * dq0)
684 endif
685 qv(i, j) = qv(i, j) - src(i)
686#ifdef MULTI_GASES
687 qvi(i,j,1,1) = qv(i, j)
688#endif
689 ql(i, j) = ql(i, j) + src(i)
690 q_liq(i) = q_liq(i) + src(i)
691 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
692 pt1(i) = pt1(i) + src(i) * lhl(i) / cvm(i)
693 enddo
694
695 ! -----------------------------------------------------------------------
697 ! -----------------------------------------------------------------------
698
699 do i = is, ie
700 lhl(i) = lv00 + d0_vap * pt1(i)
701 lhi(i) = li00 + dc_ice * pt1(i)
702 lcp2(i) = lhl(i) / cvm(i)
703 icp2(i) = lhi(i) / cvm(i)
704 tcp3(i) = lcp2(i) + icp2(i) * min(1., dim(tice, pt1(i)) / 48.)
705 enddo
706
707 if (last_step) then
708
709 ! -----------------------------------------------------------------------
712 ! final iteration:
713 ! -----------------------------------------------------------------------
714
715 call wqs2_vect (is, ie, pt1, den, wqsat, dq2dt)
716
717 do i = is, ie
718 dq0 = (qv(i, j) - wqsat(i)) / (1. + tcp3(i) * dq2dt(i))
719 if (dq0 > 0.) then ! remove super - saturation, prevent super saturation over water
720 src(i) = dq0
721 else ! evaporation of ql
722 ! factor = - min (1., fac_l2v * sqrt (max (0., ql (i, j)) / 1.e-5) * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90%
723 ! factor = - fac_l2v
724 ! factor = - 1
725 factor = - min(1., fac_l2v * 10. * (1. - qv(i, j) / wqsat(i))) ! the rh dependent factor = 1 at 90%
726 src(i) = - min(ql(i, j), factor * dq0)
727 endif
728 adj_fac = 1.
729 qv(i, j) = qv(i, j) - src(i)
730#ifdef MULTI_GASES
731 qvi(i,j,1,1) = qv(i,j)
732#endif
733 ql(i, j) = ql(i, j) + src(i)
734 q_liq(i) = q_liq(i) + src(i)
735 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
736 pt1(i) = pt1(i) + src(i) * lhl(i) / cvm(i)
737 enddo
738
739 ! -----------------------------------------------------------------------
741 ! -----------------------------------------------------------------------
742
743 do i = is, ie
744 lhl(i) = lv00 + d0_vap * pt1(i)
745 lhi(i) = li00 + dc_ice * pt1(i)
746 lcp2(i) = lhl(i) / cvm(i)
747 icp2(i) = lhi(i) / cvm(i)
748 enddo
749
750 endif
751
752 ! -----------------------------------------------------------------------
754 ! -----------------------------------------------------------------------
755
756 do i = is, ie
757 dtmp = t_wfr - pt1(i) ! [ - 40, - 48]
758 if (ql(i, j) > 0. .and. dtmp > 0.) then
759 sink(i) = min(ql(i, j), ql(i, j) * dtmp * 0.125, dtmp / icp2(i))
760 ql(i, j) = ql(i, j) - sink(i)
761 qi(i, j) = qi(i, j) + sink(i)
762 q_liq(i) = q_liq(i) - sink(i)
763 q_sol(i) = q_sol(i) + sink(i)
764 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
765 pt1(i) = pt1(i) + sink(i) * lhi(i) / cvm(i)
766 endif
767 enddo
768
769 ! -----------------------------------------------------------------------
771 ! -----------------------------------------------------------------------
772
773 do i = is, ie
774 lhi(i) = li00 + dc_ice * pt1(i)
775 icp2(i) = lhi(i) / cvm(i)
776 enddo
777
778 ! -----------------------------------------------------------------------
780 ! -----------------------------------------------------------------------
781
782 do i = is, ie
783 tc = tice0 - pt1(i)
784 if (ql(i, j) > 0.0 .and. tc > 0.) then
785 sink(i) = 3.3333e-10 * dt_bigg * (exp(0.66 * tc) - 1.) * den(i) * ql(i, j) ** 2
786 sink(i) = min(ql(i, j), tc / icp2(i), sink(i))
787 ql(i, j) = ql(i, j) - sink(i)
788 qi(i, j) = qi(i, j) + sink(i)
789 q_liq(i) = q_liq(i) - sink(i)
790 q_sol(i) = q_sol(i) + sink(i)
791 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
792 pt1(i) = pt1(i) + sink(i) * lhi(i) / cvm(i)
793 endif
794 enddo
795
796 ! -----------------------------------------------------------------------
798 ! -----------------------------------------------------------------------
799
800 do i = is, ie
801 lhi(i) = li00 + dc_ice * pt1(i)
802 icp2(i) = lhi(i) / cvm(i)
803 enddo
804
805 ! -----------------------------------------------------------------------
807 ! -----------------------------------------------------------------------
808
809 do i = is, ie
810 dtmp = (tice - 0.1) - pt1(i)
811 if (qr(i, j) > 1.e-7 .and. dtmp > 0.) then
812 tmp = min(1., (dtmp * 0.025) ** 2) * qr(i, j) ! no limit on freezing below - 40 deg c
813 sink(i) = min(tmp, fac_r2g * dtmp / icp2(i))
814 qr(i, j) = qr(i, j) - sink(i)
815 qg(i, j) = qg(i, j) + sink(i)
816 q_liq(i) = q_liq(i) - sink(i)
817 q_sol(i) = q_sol(i) + sink(i)
818 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
819 pt1(i) = pt1(i) + sink(i) * lhi(i) / cvm(i)
820 endif
821 enddo
822
823 ! -----------------------------------------------------------------------
825 ! -----------------------------------------------------------------------
826
827 do i = is, ie
828 lhi(i) = li00 + dc_ice * pt1(i)
829 icp2(i) = lhi(i) / cvm(i)
830 enddo
831
832 ! -----------------------------------------------------------------------
834 ! -----------------------------------------------------------------------
835
836 do i = is, ie
837 dtmp = pt1(i) - (tice + 0.1)
838 if (qs(i, j) > 1.e-7 .and. dtmp > 0.) then
839 tmp = min(1., (dtmp * 0.1) ** 2) * qs(i, j) ! no limter on melting above 10 deg c
840 sink(i) = min(tmp, fac_smlt * dtmp / icp2(i))
841 tmp = min(sink(i), dim(qs_mlt, ql(i, j))) ! max ql due to snow melt
842 qs(i, j) = qs(i, j) - sink(i)
843 ql(i, j) = ql(i, j) + tmp
844 qr(i, j) = qr(i, j) + sink(i) - tmp
845 ! qr (i, j) = qr (i, j) + sink (i)
846 q_liq(i) = q_liq(i) + sink(i)
847 q_sol(i) = q_sol(i) - sink(i)
848 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
849 pt1(i) = pt1(i) - sink(i) * lhi(i) / cvm(i)
850 endif
851 enddo
852
853 ! -----------------------------------------------------------------------
855 ! -----------------------------------------------------------------------
856
857 do i = is, ie
858 if (ql(i, j) > ql0_max) then
859 sink(i) = fac_l2r * (ql(i, j) - ql0_max)
860 qr(i, j) = qr(i, j) + sink(i)
861 ql(i, j) = ql(i, j) - sink(i)
862 endif
863 enddo
864
865 ! -----------------------------------------------------------------------
867 ! -----------------------------------------------------------------------
868
869 do i = is, ie
870 lhi(i) = li00 + dc_ice * pt1(i)
871 lhl(i) = lv00 + d0_vap * pt1(i)
872 lcp2(i) = lhl(i) / cvm(i)
873 icp2(i) = lhi(i) / cvm(i)
874 tcp2(i) = lcp2(i) + icp2(i)
875 enddo
876
877 ! -----------------------------------------------------------------------
879 ! -----------------------------------------------------------------------
880
881 do i = is, ie
882 src(i) = 0.
883 if (pt1(i) < t_sub) then ! too cold to be accurate; freeze qv as a fix
884 src(i) = dim(qv(i, j), 1.e-6)
885 elseif (pt1(i) < tice0) then
886 qsi = iqs2(pt1(i), den(i), dqsdt)
887 dq = qv(i, j) - qsi
888 sink(i) = adj_fac * dq / (1. + tcp2(i) * dqsdt)
889 if (qi(i, j) > 1.e-8) then
890 pidep = sdt * dq * 349138.78 * exp(0.875 * log(qi(i, j) * den(i))) &
891 / (qsi * den(i) * lat2 / (0.0243 * rvgas * pt1(i) ** 2) + 4.42478e4)
892 else
893 pidep = 0.
894 endif
895 if (dq > 0.) then ! vapor - > ice
896 tmp = tice - pt1(i)
897 qi_crt = qi_gen * min(qi_lim, 0.1 * tmp) / den(i)
898 src(i) = min(sink(i), max(qi_crt - qi(i, j), pidep), tmp / tcp2(i))
899 else
900 pidep = pidep * min(1., dim(pt1(i), t_sub) * 0.2)
901 src(i) = max(pidep, sink(i), - qi(i, j))
902 endif
903 endif
904 qv(i, j) = qv(i, j) - src(i)
905#ifdef MULTI_GASES
906 qvi(i,j,1,1) = qv(i,j)
907#endif
908 qi(i, j) = qi(i, j) + src(i)
909 q_sol(i) = q_sol(i) + src(i)
910 cvm(i) = mc_air(i) + qv(i, j) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
911 pt1(i) = pt1(i) + src(i) * (lhl(i) + lhi(i)) / cvm(i)
912 enddo
913
914 ! -----------------------------------------------------------------------
916 ! -----------------------------------------------------------------------
917
918 do i = is, ie
919#ifdef USE_COND
920 q_con(i, j) = q_liq(i) + q_sol(i)
921#ifdef MULTI_GASES
922 pt(i, j) = pt1(i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j))
923#else
924 tmp = 1. + zvir * qv(i, j)
925 pt(i, j) = pt1(i) * tmp * (1. - q_con(i, j))
926#endif
927 tmp = rdgas * tmp
928 cappa(i, j) = tmp / (tmp + cvm(i))
929#else
930#ifdef MULTI_GASES
931 q_con(i, j) = q_liq(i) + q_sol(i)
932 pt(i, j) = pt1(i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j)) * (1. - q_con(i,j))
933#else
934 pt(i, j) = pt1(i) * (1. + zvir * qv(i, j))
935#endif
936#endif
937 enddo
938
939 ! -----------------------------------------------------------------------
941 ! -----------------------------------------------------------------------
942
943 do i = is, ie
944 if (qg(i, j) < 0.) then
945 tmp = min(- qg(i, j), max(0., qi(i, j)))
946 qg(i, j) = qg(i, j) + tmp
947 qi(i, j) = qi(i, j) - tmp
948 endif
949 enddo
950
951 ! -----------------------------------------------------------------------
953 ! -----------------------------------------------------------------------
954
955 do i = is, ie
956 qim = qi0_max / den(i)
957 if (qi(i, j) > qim) then
958 sink(i) = fac_i2s * (qi(i, j) - qim)
959 qi(i, j) = qi(i, j) - sink(i)
960 qs(i, j) = qs(i, j) + sink(i)
961 endif
962 enddo
963
964 if (out_dt) then
965 do i = is, ie
966 dtdt(i, j) = dtdt(i, j) + pt1(i) - t0(i)
967 enddo
968 endif
969
970 ! -----------------------------------------------------------------------
972 ! -----------------------------------------------------------------------
973
974 if (consv_te) then
975 do i = is, ie
976 if (hydrostatic) then
977#ifdef MULTI_GASES
978 c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
979#endif
980 te0(i, j) = dp(i, j) * (te0(i, j) + c_air * pt1(i))
981 else
982#ifdef USE_COND
983 te0(i, j) = dp(i, j) * (te0(i, j) + cvm(i) * pt1(i))
984#else
985#ifdef MULTI_GASES
986 c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
987#endif
988 te0(i, j) = dp(i, j) * (te0(i, j) + c_air * pt1(i))
989#endif
990 endif
991 enddo
992 endif
993
994 ! -----------------------------------------------------------------------
996 ! -----------------------------------------------------------------------
997
998 do i = is, ie
999 lhi(i) = li00 + dc_ice * pt1(i)
1000 lhl(i) = lv00 + d0_vap * pt1(i)
1001 cvm(i) = mc_air(i) + (qv(i, j) + q_liq(i) + q_sol(i)) * c_vap
1002 lcp2(i) = lhl(i) / cvm(i)
1003 icp2(i) = lhi(i) / cvm(i)
1004 enddo
1005
1006 ! -----------------------------------------------------------------------
1008 ! -----------------------------------------------------------------------
1009
1010 if (do_qa .and. last_step) then
1011
1012 ! -----------------------------------------------------------------------
1014 ! -----------------------------------------------------------------------
1015
1016 if (rad_snow) then
1017 if (rad_graupel) then
1018 do i = is, ie
1019 q_sol(i) = qi(i, j) + qs(i, j) + qg(i, j)
1020 enddo
1021 else
1022 do i = is, ie
1023 q_sol(i) = qi(i, j) + qs(i, j)
1024 enddo
1025 endif
1026 else
1027 do i = is, ie
1028 q_sol(i) = qi(i, j)
1029 enddo
1030 endif
1031 if (rad_rain) then
1032 do i = is, ie
1033 q_liq(i) = ql(i, j) + qr(i, j)
1034 enddo
1035 else
1036 do i = is, ie
1037 q_liq(i) = ql(i, j)
1038 enddo
1039 endif
1040 do i = is, ie
1041 q_cond(i) = q_sol(i) + q_liq(i)
1042 enddo
1043
1044 ! -----------------------------------------------------------------------
1046 ! -----------------------------------------------------------------------
1047
1048 do i = is, ie
1049
1050 if(tintqs) then
1051 tin = pt1(i)
1052 else
1053 tin = pt1(i) - (lcp2(i) * q_cond(i) + icp2(i) * q_sol(i)) ! minimum temperature
1054 ! tin = pt1 (i) - ((lv00 + d0_vap * pt1 (i)) * q_cond (i) + &
1055 ! (li00 + dc_ice * pt1 (i)) * q_sol (i)) / (mc_air (i) + qpz (i) * c_vap)
1056 endif
1057
1058 ! -----------------------------------------------------------------------
1059 ! determine saturated specific humidity
1060 ! -----------------------------------------------------------------------
1061
1062 if (tin <= t_wfr) then
1063 ! ice phase:
1064 qstar(i) = iqs1(tin, den(i))
1065 elseif (tin >= tice) then
1066 ! liquid phase:
1067 qstar(i) = wqs1(tin, den(i))
1068 else
1069 ! mixed phase:
1070 qsi = iqs1(tin, den(i))
1071 qsw = wqs1(tin, den(i))
1072 if (q_cond(i) > 1.e-6) then
1073 rqi = q_sol(i) / q_cond(i)
1074 else
1075 ! mostly liquid water clouds at initial cloud development stage
1076 rqi = ((tice - tin) / (tice - t_wfr))
1077 endif
1078 qstar(i) = rqi * qsi + (1. - rqi) * qsw
1079 endif
1081 dw = dw_ocean + (dw_land - dw_ocean) * min(1., abs(hs(i, j)) / (10. * grav))
1083 hvar(i) = min(0.2, max(0.01, dw * sqrt(sqrt(area(i, j)) / 100.e3)))
1084
1085 ! -----------------------------------------------------------------------
1089 ! -----------------------------------------------------------------------
1090
1091 rh = qpz(i) / qstar(i)
1092
1093 ! -----------------------------------------------------------------------
1094 ! icloud_f = 0: bug - fixed
1095 ! icloud_f = 1: old fvgfs gfdl) mp implementation
1096 ! icloud_f = 2: binary cloud scheme (0 / 1)
1097 ! -----------------------------------------------------------------------
1098
1099 if (rh > 0.75 .and. qpz(i) > 1.e-8) then
1100 dq = hvar(i) * qpz(i)
1101 q_plus = qpz(i) + dq
1102 q_minus = qpz(i) - dq
1103 if (icloud_f == 2) then
1104 if (qpz(i) > qstar(i)) then
1105 qa(i, j) = 1.
1106 elseif (qstar(i) < q_plus .and. q_cond(i) > 1.e-8) then
1107 qa(i, j) = ((q_plus - qstar(i)) / dq) ** 2
1108 qa(i, j) = min(1., qa(i, j))
1109 else
1110 qa(i, j) = 0.
1111 endif
1112 else
1113 if (qstar(i) < q_minus) then
1114 qa(i, j) = 1.
1115 else
1116 if (qstar(i) < q_plus) then
1117 if (icloud_f == 0) then
1118 qa(i, j) = (q_plus - qstar(i)) / (dq + dq)
1119 else
1120 qa(i, j) = (q_plus - qstar(i)) / (2. * dq * (1. - q_cond(i)))
1121 endif
1122 else
1123 qa(i, j) = 0.
1124 endif
1125 ! impose minimum cloudiness if substantial q_cond (i) exist
1126 if (q_cond(i) > 1.e-8) then
1127 qa(i, j) = max(cld_min, qa(i, j))
1128 endif
1129 qa(i, j) = min(1., qa(i, j))
1130 endif
1131 endif
1132 else
1133 qa(i, j) = 0.
1134 endif
1135
1136 enddo
1137
1138 endif
1139
1140 enddo ! end j loop
1141
1142end subroutine fv_sat_adj_work
1144
1145! =======================================================================
1149! =======================================================================
1150real(kind=kind_dyn) function wqs1 (ta, den)
1151
1152 implicit none
1153
1154 ! pure water phase; universal dry / moist formular using air density
1155 ! input "den" can be either dry or moist air density
1156
1157 real(kind=kind_dyn), intent (in) :: ta, den
1158
1159 real(kind=kind_dyn) :: es, ap1, tmin
1160
1161 integer :: it
1162
1163 tmin = tice - 160.
1164 ap1 = 10. * dim(ta, tmin) + 1.
1165 ap1 = min(2621., ap1)
1166 it = ap1
1167 es = tablew(it) + (ap1 - it) * desw(it)
1168 wqs1 = es / (rvgas * ta * den)
1169
1170end function wqs1
1171
1172! =======================================================================
1176! =======================================================================
1177real(kind=kind_dyn) function iqs1 (ta, den)
1178
1179 implicit none
1180
1181 ! water - ice phase; universal dry / moist formular using air density
1182 ! input "den" can be either dry or moist air density
1183
1184 real(kind=kind_dyn), intent (in) :: ta, den
1185
1186 real(kind=kind_dyn) :: es, ap1, tmin
1187
1188 integer :: it
1189
1190 tmin = tice - 160.
1191 ap1 = 10. * dim(ta, tmin) + 1.
1192 ap1 = min(2621., ap1)
1193 it = ap1
1194 es = table2(it) + (ap1 - it) * des2(it)
1195 iqs1 = es / (rvgas * ta * den)
1196
1197end function iqs1
1198
1199! =======================================================================
1203! =======================================================================
1204real(kind=kind_dyn) function wqs2 (ta, den, dqdt)
1205
1206 implicit none
1207
1208 ! pure water phase; universal dry / moist formular using air density
1209 ! input "den" can be either dry or moist air density
1210
1211 real(kind=kind_dyn), intent (in) :: ta, den
1212
1213 real(kind=kind_dyn), intent (out) :: dqdt
1214
1215 real(kind=kind_dyn) :: es, ap1, tmin
1216
1217 integer :: it
1218
1219 tmin = tice - 160.
1220 ap1 = 10. * dim(ta, tmin) + 1.
1221 ap1 = min(2621., ap1)
1222 it = ap1
1223 es = tablew(it) + (ap1 - it) * desw(it)
1224 wqs2 = es / (rvgas * ta * den)
1225 it = ap1 - 0.5
1226 ! finite diff, del_t = 0.1:
1227 dqdt = 10. * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) / (rvgas * ta * den)
1228
1229end function wqs2
1230
1231! =======================================================================
1236! =======================================================================
1237subroutine wqs2_vect (is, ie, ta, den, wqsat, dqdt)
1238
1239 implicit none
1240
1241 ! pure water phase; universal dry / moist formular using air density
1242 ! input "den" can be either dry or moist air density
1243
1244 integer, intent (in) :: is, ie
1245
1246 real(kind=kind_dyn), intent (in), dimension (is:ie) :: ta, den
1247
1248 real(kind=kind_dyn), intent (out), dimension (is:ie) :: wqsat, dqdt
1249
1250 real(kind=kind_dyn) :: es, ap1, tmin
1251
1252 integer :: i, it
1253
1254 tmin = tice - 160.
1255
1256 do i = is, ie
1257 ap1 = 10. * dim(ta(i), tmin) + 1.
1258 ap1 = min(2621., ap1)
1259 it = ap1
1260 es = tablew(it) + (ap1 - it) * desw(it)
1261 wqsat(i) = es / (rvgas * ta(i) * den(i))
1262 it = ap1 - 0.5
1263 ! finite diff, del_t = 0.1:
1264 dqdt(i) = 10. * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) / (rvgas * ta(i) * den(i))
1265 enddo
1266
1267end subroutine wqs2_vect
1268
1269! =======================================================================
1273! =======================================================================
1274real(kind=kind_dyn) function iqs2 (ta, den, dqdt)
1275
1276 implicit none
1277
1278 ! water - ice phase; universal dry / moist formular using air density
1279 ! input "den" can be either dry or moist air density
1280
1281 real(kind=kind_dyn), intent (in) :: ta, den
1282
1283 real(kind=kind_dyn), intent (out) :: dqdt
1284
1285 real(kind=kind_dyn) :: es, ap1, tmin
1286
1287 integer :: it
1288
1289 tmin = tice - 160.
1290 ap1 = 10. * dim(ta, tmin) + 1.
1291 ap1 = min(2621., ap1)
1292 it = ap1
1293 es = table2(it) + (ap1 - it) * des2(it)
1294 iqs2 = es / (rvgas * ta * den)
1295 it = ap1 - 0.5
1296 ! finite diff, del_t = 0.1:
1297 dqdt = 10. * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) / (rvgas * ta * den)
1298
1299end function iqs2
1300
1301! =======================================================================
1304! 3 - phase table
1305! =======================================================================
1306
1307subroutine qs_table (n)
1308
1309 implicit none
1310
1311 integer, intent (in) :: n
1312 real (kind_grid) :: delt = 0.1
1313 real (kind_grid) :: tmin, tem, esh20
1314 real (kind_grid) :: wice, wh2o, fac0, fac1, fac2
1315 real (kind_grid) :: esupc (200)
1316 integer :: i
1317
1318 tmin = tice - 160.
1319
1320 ! -----------------------------------------------------------------------
1321 ! compute es over ice between - 160 deg c and 0 deg c.
1322 ! -----------------------------------------------------------------------
1323
1324 do i = 1, 1600
1325 tem = tmin + delt * real(i - 1)
1326 fac0 = (tem - tice) / (tem * tice)
1327 fac1 = fac0 * li2
1328 fac2 = (d2ice * log(tem / tice) + fac1) / rvgas
1329 table(i) = e00 * exp(fac2)
1330 enddo
1331
1332 ! -----------------------------------------------------------------------
1333 ! compute es over water between - 20 deg c and 102 deg c.
1334 ! -----------------------------------------------------------------------
1335
1336 do i = 1, 1221
1337 tem = 253.16 + delt * real(i - 1)
1338 fac0 = (tem - tice) / (tem * tice)
1339 fac1 = fac0 * lv0
1340 fac2 = (dc_vap * log(tem / tice) + fac1) / rvgas
1341 esh20 = e00 * exp(fac2)
1342 if (i <= 200) then
1343 esupc(i) = esh20
1344 else
1345 table(i + 1400) = esh20
1346 endif
1347 enddo
1348
1349 ! -----------------------------------------------------------------------
1350 ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c
1351 ! -----------------------------------------------------------------------
1352
1353 do i = 1, 200
1354 tem = 253.16 + delt * real(i - 1)
1355 wice = 0.05 * (tice - tem)
1356 wh2o = 0.05 * (tem - 253.16)
1357 table(i + 1400) = wice * table(i + 1400) + wh2o * esupc(i)
1358 enddo
1359
1360end subroutine qs_table
1361
1362! =======================================================================
1365! 1 - phase table
1366! =======================================================================
1367
1368subroutine qs_tablew (n)
1369
1370 implicit none
1371
1372 integer, intent (in) :: n
1373 real (kind_grid) :: delt = 0.1
1374 real (kind_grid) :: tmin, tem, fac0, fac1, fac2
1375 integer :: i
1376
1377 tmin = tice - 160.
1378
1379 ! -----------------------------------------------------------------------
1380 ! compute es over water
1381 ! -----------------------------------------------------------------------
1382
1383 do i = 1, n
1384 tem = tmin + delt * real(i - 1)
1385 fac0 = (tem - tice) / (tem * tice)
1386 fac1 = fac0 * lv0
1387 fac2 = (dc_vap * log(tem / tice) + fac1) / rvgas
1388 tablew(i) = e00 * exp(fac2)
1389 enddo
1390
1391end subroutine qs_tablew
1392
1393! =======================================================================
1396! 2 - phase table
1397! =======================================================================
1398
1399subroutine qs_table2 (n)
1400
1401 implicit none
1402
1403 integer, intent (in) :: n
1404 real (kind_grid) :: delt = 0.1
1405 real (kind_grid) :: tmin, tem0, tem1, fac0, fac1, fac2
1406 integer :: i, i0, i1
1407
1408 tmin = tice - 160.
1409
1410 do i = 1, n
1411 tem0 = tmin + delt * real(i - 1)
1412 fac0 = (tem0 - tice) / (tem0 * tice)
1413 if (i <= 1600) then
1414 ! -----------------------------------------------------------------------
1415 ! compute es over ice between - 160 deg c and 0 deg c.
1416 ! -----------------------------------------------------------------------
1417 fac1 = fac0 * li2
1418 fac2 = (d2ice * log(tem0 / tice) + fac1) / rvgas
1419 else
1420 ! -----------------------------------------------------------------------
1421 ! compute es over water between 0 deg c and 102 deg c.
1422 ! -----------------------------------------------------------------------
1423 fac1 = fac0 * lv0
1424 fac2 = (dc_vap * log(tem0 / tice) + fac1) / rvgas
1425 endif
1426 table2(i) = e00 * exp(fac2)
1427 enddo
1428
1429 ! -----------------------------------------------------------------------
1430 ! smoother around 0 deg c
1431 ! -----------------------------------------------------------------------
1432
1433 i0 = 1600
1434 i1 = 1601
1435 tem0 = 0.25 * (table2(i0 - 1) + 2. * table(i0) + table2(i0 + 1))
1436 tem1 = 0.25 * (table2(i1 - 1) + 2. * table(i1) + table2(i1 + 1))
1437 table2(i0) = tem0
1438 table2(i1) = tem1
1439
1440end subroutine qs_table2
1441
1442end module fv_sat_adj
fv_sat_adj::qs_table
subroutine qs_table(n)
saturation water vapor pressure table i
Definition fv_sat_adj.F90:1308
fv_sat_adj::qs_tablew
subroutine qs_tablew(n)
saturation water vapor pressure table ii.
Definition fv_sat_adj.F90:1369
fv_sat_adj::qs_table2
subroutine qs_table2(n)
saturation water vapor pressure table iii.
Definition fv_sat_adj.F90:1400
fv_sat_adj::wqs2
real(kind=kind_dyn) function wqs2(ta, den, dqdt)
The function 'wqs2'computes the gradient of saturated specific humidity for table ii.
Definition fv_sat_adj.F90:1205
fv_sat_adj::wqs1
real(kind=kind_dyn) function wqs1(ta, den)
the function 'wqs1' computes the saturated specific humidity for table ii.
Definition fv_sat_adj.F90:1151
fv_sat_adj::iqs2
real(kind=kind_dyn) function iqs2(ta, den, dqdt)
The function 'iqs2' computes the gradient of saturated specific humidity for table iii.
Definition fv_sat_adj.F90:1275
fv_sat_adj::fv_sat_adj_work
subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0, ifdef multi_gases
This subroutine includes the entity of the fast saturation adjustment processes.
Definition fv_sat_adj.F90:404
fv_sat_adj::iqs1
real(kind=kind_dyn) function iqs1(ta, den)
the function 'wqs1' computes the saturated specific humidity for table iii
Definition fv_sat_adj.F90:1178
fv_sat_adj::wqs2_vect
subroutine wqs2_vect(is, ie, ta, den, wqsat, dqdt)
The function wqs2_vect computes the gradient of saturated specific humidity for table ii....
Definition fv_sat_adj.F90:1238
fv_sat_adj::fv_sat_adj_run
subroutine, public fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, isc1, iec1, isc2, iec2, kmp, km, kmdelz, js, je, jsd, jed, jsc1, jec1, jsc2, jec2, ng, hydrostatic, fast_mp_consv, te0_2d, te0, ngas, qvi, qv, ql, qi, qr, qs, qg, hs, peln, delz, delp, pt, pkz, q_con, akap, cappa, area, dtdt, out_dt, last_step, do_qa, qa, nthreads, errmsg, errflg)
Definition fv_sat_adj.F90:238
ccpp_multi_gases_mod
The module 'multi_gases' peforms multi constitutents computations.
Definition multi_gases.F90:25
fv_sat_adj
This module contains the GFDL in-core fast saturation adjustment called in FV3 dynamics solver.
Definition fv_sat_adj.F90:28
gfdl_cloud_microphys_mod
This module contains the column GFDL Cloud microphysics scheme.
Definition module_gfdl_cloud_microphys.F90:32