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