109 integer,
public,
parameter ::
qa_mp = 11
127 'Ratio of Water Vapor mass to total mass (Specific humidity)', &
128 'Ratio of Cloud Water mass to total mass ', &
129 'Ratio of Rain Water mass to total mass ', &
130 'Ratio of Cloud Ice mass ratio to total mass ', &
131 'Ratio of Snow miass ratio to total mass ', &
132 'Ratio of Graupel mass ratio to total mass ', &
133 'Cloud Water Number Density ', &
134 'Rain Water Number Density ', &
135 'Cloud Ice Number Density ', &
136 'Snow Number Density ', &
137 'Graupel Number Density '/)
155 private :: mp_sn14_init
162 integer,
private,
parameter :: i_qv = 1
163 integer,
private,
parameter :: i_qc = 2
164 integer,
private,
parameter :: i_qr = 3
165 integer,
private,
parameter :: i_qi = 4
166 integer,
private,
parameter :: i_qs = 5
167 integer,
private,
parameter :: i_qg = 6
168 integer,
private,
parameter :: i_nc = 7
169 integer,
private,
parameter :: i_nr = 8
170 integer,
private,
parameter :: i_ni = 9
171 integer,
private,
parameter :: i_ns = 10
172 integer,
private,
parameter :: i_ng = 11
174 integer,
private,
parameter :: hydro_max = 5
176 integer,
private,
parameter :: i_mp_qc = 1
177 integer,
private,
parameter :: i_mp_qr = 2
178 integer,
private,
parameter :: i_mp_qi = 3
179 integer,
private,
parameter :: i_mp_qs = 4
180 integer,
private,
parameter :: i_mp_qg = 5
181 integer,
private,
parameter :: i_mp_nc = 6
182 integer,
private,
parameter :: i_mp_nr = 7
183 integer,
private,
parameter :: i_mp_ni = 8
184 integer,
private,
parameter :: i_mp_ns = 9
185 integer,
private,
parameter :: i_mp_ng = 10
189 integer,
parameter :: i_lcccn = 1
190 integer,
parameter :: i_ncccn = 2
191 integer,
parameter :: i_liccn = 3
192 integer,
parameter :: i_niccn = 4
194 integer,
parameter :: i_lchom = 5
195 integer,
parameter :: i_nchom = 6
196 integer,
parameter :: i_lchet = 7
197 integer,
parameter :: i_nchet = 8
198 integer,
parameter :: i_lrhet = 9
199 integer,
parameter :: i_nrhet = 10
201 integer,
parameter :: i_limlt = 11
202 integer,
parameter :: i_nimlt = 12
203 integer,
parameter :: i_lsmlt = 13
204 integer,
parameter :: i_nsmlt = 14
205 integer,
parameter :: i_lgmlt = 15
206 integer,
parameter :: i_ngmlt = 16
208 integer,
parameter :: i_lrdep = 17
209 integer,
parameter :: i_nrdep = 18
210 integer,
parameter :: i_lidep = 19
211 integer,
parameter :: i_nidep = 20
212 integer,
parameter :: i_lsdep = 21
213 integer,
parameter :: i_nsdep = 22
214 integer,
parameter :: i_lgdep = 23
215 integer,
parameter :: i_ngdep = 24
216 integer,
parameter :: i_lcdep = 25
220 integer,
parameter :: i_lcaut = 26
221 integer,
parameter :: i_ncaut = 27
222 integer,
parameter :: i_nraut = 28
224 integer,
parameter :: i_lcacc = 29
225 integer,
parameter :: i_ncacc = 30
227 integer,
parameter :: i_nrslc = 31
228 integer,
parameter :: i_nrbrk = 32
231 integer,
parameter :: i_licon = 33
232 integer,
parameter :: i_nicon = 34
233 integer,
parameter :: i_lscon = 35
234 integer,
parameter :: i_nscon = 36
236 integer,
parameter :: i_liacm = 37
237 integer,
parameter :: i_niacm = 38
238 integer,
parameter :: i_liarm = 39
239 integer,
parameter :: i_niarm = 40
240 integer,
parameter :: i_lsacm = 41
241 integer,
parameter :: i_nsacm = 42
242 integer,
parameter :: i_lsarm = 43
243 integer,
parameter :: i_nsarm = 44
244 integer,
parameter :: i_lgacm = 45
245 integer,
parameter :: i_ngacm = 46
246 integer,
parameter :: i_lgarm = 47
247 integer,
parameter :: i_ngarm = 48
249 integer,
parameter :: i_lgspl = 49
250 integer,
parameter :: i_lsspl = 50
251 integer,
parameter :: i_nispl = 51
253 integer,
parameter :: pq_max = 51
257 integer,
parameter :: i_liaclc2li = 1
258 integer,
parameter :: i_niacnc2ni = 2
259 integer,
parameter :: i_lsaclc2ls = 3
260 integer,
parameter :: i_nsacnc2ns = 4
261 integer,
parameter :: i_lgaclc2lg = 5
262 integer,
parameter :: i_ngacnc2ng = 6
263 integer,
parameter :: i_lracli2lg_i = 7
264 integer,
parameter :: i_nracni2ng_i = 8
265 integer,
parameter :: i_lracli2lg_r = 9
266 integer,
parameter :: i_nracni2ng_r = 10
267 integer,
parameter :: i_lracls2lg_s = 11
268 integer,
parameter :: i_nracns2ng_s = 12
269 integer,
parameter :: i_lracls2lg_r = 13
270 integer,
parameter :: i_nracns2ng_r = 14
271 integer,
parameter :: i_lraclg2lg = 15
272 integer,
parameter :: i_nracng2ng = 16
273 integer,
parameter :: i_liacli2ls = 17
274 integer,
parameter :: i_niacni2ns = 18
275 integer,
parameter :: i_liacls2ls = 19
276 integer,
parameter :: i_niacns2ns = 20
277 integer,
parameter :: i_nsacns2ns = 21
278 integer,
parameter :: i_ngacng2ng = 22
279 integer,
parameter :: i_lgacls2lg = 23
280 integer,
parameter :: i_ngacns2ng = 24
282 integer,
parameter :: pac_max = 24
285 integer,
parameter :: i_cgngacns2ng = 25
286 integer,
parameter :: i_cgngacni2ng = 26
287 integer,
parameter :: i_ngspl = 49
288 integer,
parameter :: i_nsspl = 50
289 integer,
parameter :: pcrg_max = 26
291 character(len=H_SHORT),
save :: wlabel(hydro_max)
294 real(
rp),
private,
parameter :: nqmin(hydro_max) = (/ 1.e+4_rp, 1.0_rp, 1.0_rp, 1.e-4_rp, 1.e-4_rp /)
297 real(
rp),
private,
parameter :: xqmax(hydro_max) = (/ 2.6e-10_rp, 5.0e-6_rp, 1.377e-6_rp, 7.519e-6_rp, 4.90e-5_rp /)
300 real(
rp),
private,
parameter :: xqmin(hydro_max) = (/ 4.20e-15_rp, 2.60e-10_rp, 3.382e-13_rp, 1.847e-12_rp, 1.230e-10_rp /)
306 real(
rp),
private,
parameter :: xc_min = 4.20e-15_rp
307 real(
rp),
private,
parameter :: xr_min = 2.60e-10_rp
308 real(
rp),
private,
parameter :: xi_min = 3.382e-13_rp
309 real(
rp),
private,
parameter :: xs_min = 1.847e-12_rp
310 real(
rp),
private,
parameter :: xg_min = 1.230e-10_rp
312 real(
rp),
private,
parameter :: xc_max = 2.6e-10_rp
313 real(
rp),
private,
parameter :: xr_max = 5.00e-6_rp
314 real(
rp),
private,
parameter :: xi_max = 1.377e-6_rp
315 real(
rp),
private,
parameter :: xs_max = 7.519e-6_rp
316 real(
rp),
private,
parameter :: xg_max = 4.900e-5_rp
318 real(
rp),
private,
parameter :: xmin_filter= xc_min
320 real(
rp),
private,
parameter :: rmin_re= 1.e-6_rp
323 real(
rp),
private,
parameter :: n0r_min= 2.5e+5_rp
324 real(
rp),
private,
parameter :: n0r_max= 2.0e+7_rp
325 real(
rp),
private,
parameter :: lambdar_min= 1.e+3_rp
326 real(
rp),
private,
parameter :: lambdar_max= 1.e+4_rp
328 real(
rp),
private,
parameter :: nc_min = 1.e+4_rp
329 real(
rp),
private,
parameter :: nr_min = 1.0_rp
330 real(
rp),
private,
parameter :: ni_min = 1.0_rp
331 real(
rp),
private,
parameter :: ns_min = 1.e-4_rp
332 real(
rp),
private,
parameter :: ng_min = 1.e-4_rp
334 real(
rp),
private,
parameter :: lc_min = xc_min*nc_min
335 real(
rp),
private,
parameter :: lr_min = xr_min*nr_min
336 real(
rp),
private,
parameter :: li_min = xi_min*ni_min
337 real(
rp),
private,
parameter :: ls_min = xs_min*ns_min
338 real(
rp),
private,
parameter :: lg_min = xg_min*ng_min
340 real(
rp),
private,
parameter :: x_sep = 2.6e-10_rp
342 real(
rp),
private,
parameter :: tem_min=100.0_rp
343 real(
rp),
private,
parameter :: rho_min=1.e-5_rp
344 real(
rp),
private,
parameter :: rhoi = 916.70_rp
346 integer,
private,
save :: ntmax_phase_change = 1
347 integer,
private,
save :: ntmax_collection = 1
350 real(
rp),
private,
parameter :: rho_0 = 1.280_rp
352 real(
rp),
allocatable,
private,
save :: nc_uplim_d(:,:,:)
355 real(
rp),
private,
parameter :: ka0 = 2.428e-2_rp
357 real(
rp),
private,
parameter :: dka_dt = 7.47e-5_rp
362 real(
rp),
private,
parameter :: mua0 = 1.718e-5_rp
364 real(
rp),
private,
parameter :: dmua_dt = 5.28e-8_rp
368 real(
rp),
private,
save :: xc_ccn = 1.e-12_rp
369 real(
rp),
private,
save :: xi_ccn = 1.e-12_rp
373 real(
rp),
private,
save :: cap(hydro_max)
377 real(
rp),
private,
save :: a_m(hydro_max)
378 real(
rp),
private,
save :: b_m(hydro_max)
381 real(
rp),
private,
save :: alpha_v(hydro_max,2)
382 real(
rp),
private,
save :: beta_v(hydro_max,2)
383 real(
rp),
private,
save :: alpha_vn(hydro_max,2)
384 real(
rp),
private,
save :: beta_vn(hydro_max,2)
385 real(
rp),
private,
save :: gamma_v(hydro_max)
388 real(
rp),
private,
parameter :: pre0_vt = 300.e+2_rp
389 real(
rp),
private,
parameter :: tem0_vt = 233.0_rp
390 real(
rp),
private,
parameter :: a_pre0_vt = -0.1780_rp
391 real(
rp),
private,
parameter :: a_tem0_vt = -0.3940_rp
399 real(
rp),
private,
save :: nu(hydro_max)
400 real(
rp),
private,
save :: mu(hydro_max)
406 real(
rp),
private,
save :: a_area(hydro_max)
407 real(
rp),
private,
save :: b_area(hydro_max)
408 real(
rp),
private,
save :: ax_area(hydro_max)
409 real(
rp),
private,
save :: bx_area(hydro_max)
412 real(
rp),
private,
save :: a_rea(hydro_max)
413 real(
rp),
private,
save :: b_rea(hydro_max)
414 real(
rp),
private,
save :: a_rea2(hydro_max)
415 real(
rp),
private,
save :: b_rea2(hydro_max)
416 real(
rp),
private,
save :: a_rea3(hydro_max)
417 real(
rp),
private,
save :: b_rea3(hydro_max)
419 real(
rp),
private,
save :: a_d2vt(hydro_max)
420 real(
rp),
private,
save :: b_d2vt(hydro_max)
424 real(
rp),
private,
save :: coef_m2(hydro_max)
426 real(
rp),
private,
save :: coef_d6(hydro_max)
428 real(
rp),
private,
save :: coef_d3(hydro_max)
430 real(
rp),
private,
save :: coef_d(hydro_max)
432 real(
rp),
private,
save :: coef_d2v(hydro_max)
434 real(
rp),
private,
save :: coef_md2v(hydro_max)
437 real(
rp),
private,
save :: coef_r2(hydro_max)
438 real(
rp),
private,
save :: coef_r3(hydro_max)
439 real(
rp),
private,
save :: coef_re(hydro_max)
441 real(
rp),
private,
save :: coef_rea2(hydro_max)
442 real(
rp),
private,
save :: coef_rea3(hydro_max)
443 logical,
private,
save :: opt_m96_ice=.true.
444 logical,
private,
save :: opt_m96_column_ice=.false.
449 real(
rp),
private,
save :: coef_vt0(hydro_max,2)
450 real(
rp),
private,
save :: coef_vt1(hydro_max,2)
451 real(
rp),
private,
save :: coef_deplc
452 real(
rp),
private,
save :: coef_dave_n(hydro_max)
453 real(
rp),
private,
save :: coef_dave_l(hydro_max)
456 real(
rp),
private,
save :: d0_ni=261.76e-6_rp
457 real(
rp),
private,
save :: d0_li=398.54e-6_rp
458 real(
rp),
private,
parameter :: d0_ns=270.03e-6_rp
459 real(
rp),
private,
parameter :: d0_ls=397.47e-6_rp
460 real(
rp),
private,
parameter :: d0_ng=269.08e-6_rp
461 real(
rp),
private,
parameter :: d0_lg=376.36e-6_rp
463 real(
rp),
private,
parameter :: coef_vtr_ar1=9.65_rp
465 real(
rp),
private,
parameter :: coef_vtr_br1=10.43_rp
466 real(
rp),
private,
parameter :: coef_vtr_cr1=600.0_rp
467 real(
rp),
private,
parameter :: coef_vtr_ar2=4.e+3_rp
468 real(
rp),
private,
parameter :: coef_vtr_br2=12.e+3_rp
469 real(
rp),
private,
parameter :: d_vtr_branch=0.745e-3_rp
471 real(
rp),
private,
parameter :: dr_eq = 1.10e-3_rp
476 real(
rp),
private,
save :: coef_a(hydro_max)
477 real(
rp),
private,
save :: coef_lambda(hydro_max)
481 real(
rp),
private,
save :: ah_vent (hydro_max,2)
482 real(
rp),
private,
save :: bh_vent (hydro_max,2)
483 real(
rp),
private,
save :: ah_vent0 (hydro_max,2)
484 real(
rp),
private,
save :: bh_vent0 (hydro_max,2)
485 real(
rp),
private,
save :: ah_vent1 (hydro_max,2)
486 real(
rp),
private,
save :: bh_vent1 (hydro_max,2)
488 real(
rp),
private,
save :: delta_b0 (hydro_max)
489 real(
rp),
private,
save :: delta_b1 (hydro_max)
490 real(
rp),
private,
save :: delta_ab0(hydro_max,hydro_max)
491 real(
rp),
private,
save :: delta_ab1(hydro_max,hydro_max)
493 real(
rp),
private,
save :: theta_b0 (hydro_max)
494 real(
rp),
private,
save :: theta_b1 (hydro_max)
495 real(
rp),
private,
save :: theta_ab0(hydro_max,hydro_max)
496 real(
rp),
private,
save :: theta_ab1(hydro_max,hydro_max)
498 logical,
private,
save :: opt_debug=.false.
500 logical,
private,
save :: opt_debug_inc=.true.
501 logical,
private,
save :: opt_debug_act=.true.
502 logical,
private,
save :: opt_debug_ree=.true.
503 logical,
private,
save :: opt_debug_bcs=.true.
505 logical,
private,
save :: mp_doautoconversion = .true.
506 logical,
private,
save :: mp_couple_aerosol = .false.
507 real(
rp),
private,
save :: mp_ssw_lim = 1.e+1_rp
513 real(
rp),
private,
parameter :: c_ccn_ocean= 1.00e+8_rp
514 real(
rp),
private,
parameter :: c_ccn_land = 1.26e+9_rp
515 real(
rp),
private,
save :: c_ccn = 1.00e+8_rp
517 real(
rp),
private,
parameter :: kappa_ocean= 0.462_rp
518 real(
rp),
private,
parameter :: kappa_land = 0.308_rp
519 real(
rp),
private,
save :: kappa = 0.462_rp
520 real(
rp),
private,
save :: c_in = 1.0_rp
522 real(
rp),
private,
save :: nm_m92 = 1.e+3_rp
523 real(
rp),
private,
save :: am_m92 = -0.639_rp
524 real(
rp),
private,
save :: bm_m92 = 12.96_rp
526 real(
rp),
private,
save :: in_max = 1000.e+3_rp
527 real(
rp),
private,
save :: ssi_max= 0.60_rp
528 real(
rp),
private,
save :: ssw_max= 1.1_rp
530 real(
rp),
private,
save :: qke_min = 0.03_rp
531 real(
rp),
private,
save :: tem_ccn_low=233.150_rp
532 real(
rp),
private,
save :: tem_in_low =173.150_rp
533 logical,
private,
save :: nucl_twomey = .false.
534 logical,
private,
save :: inucl_w = .false.
538 real(
rp),
private,
parameter :: rc_cr= 12.e-6_rp
539 real(
rp),
private,
save :: xc_cr
540 real(
rp),
private,
save :: alpha
541 real(
rp),
private,
save :: gm, lgm
546 real(
rp),
private,
save :: dc0 = 15.0e-6_rp
547 real(
rp),
private,
save :: dc1 = 40.0e-6_rp
548 real(
rp),
private,
save :: di0 = 150.0e-6_rp
549 real(
rp),
private,
save :: ds0 = 150.0e-6_rp
550 real(
rp),
private,
save :: dg0 = 150.0e-6_rp
552 real(
rp),
private,
save :: sigma_c=0.00_rp
553 real(
rp),
private,
save :: sigma_r=0.00_rp
554 real(
rp),
private,
save :: sigma_i=0.2_rp
555 real(
rp),
private,
save :: sigma_s=0.2_rp
556 real(
rp),
private,
save :: sigma_g=0.00_rp
558 real(
rp),
private,
save :: e_im = 0.80_rp
559 real(
rp),
private,
save :: e_sm = 0.80_rp
560 real(
rp),
private,
save :: e_gm = 1.00_rp
562 real(
rp),
private,
save :: e_ir=1.0_rp
563 real(
rp),
private,
save :: e_sr=1.0_rp
564 real(
rp),
private,
save :: e_gr=1.0_rp
565 real(
rp),
private,
save :: e_ii=1.0_rp
566 real(
rp),
private,
save :: e_si=1.0_rp
567 real(
rp),
private,
save :: e_gi=1.0_rp
568 real(
rp),
private,
save :: e_ss=1.0_rp
569 real(
rp),
private,
save :: e_gs=1.0_rp
570 real(
rp),
private,
save :: e_gg=1.0_rp
573 integer,
private,
save :: i_iconv2g=1
574 integer,
private,
save :: i_sconv2g=1
576 real(
rp),
private,
save :: rho_g = 900.0_rp
578 real(
rp),
private,
save :: cfill_i = 0.68_rp
579 real(
rp),
private,
save :: cfill_s = 0.01_rp
581 real(
rp),
private,
save :: di_cri = 500.e-6_rp
582 logical,
private,
save :: opt_stick_ks96=.false.
583 logical,
private,
save :: opt_stick_co86=.false.
586 real(
rp),
private,
save :: fac_cndc = 1.0_rp
587 logical,
private,
save :: opt_fix_taucnd_c=.false.
601 integer,
intent(in) :: ka
602 integer,
intent(in) :: ia
603 integer,
intent(in) :: ja
605 namelist / param_atmos_phy_mp_sn14 / &
606 mp_doautoconversion, &
614 log_info(
"ATMOS_PHY_MP_sn14_setup",*)
'Setup'
615 log_info(
"ATMOS_PHY_MP_sn14_setup",*)
'Seiki and Nakajima (2014) 2-moment bulk 6 category'
619 read(
io_fid_conf,nml=param_atmos_phy_mp_sn14,iostat=ierr)
621 log_info(
"ATMOS_PHY_MP_sn14_setup",*)
'Not found namelist. Default used.'
622 elseif( ierr > 0 )
then
623 log_error(
"ATMOS_PHY_MP_sn14_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14. Check!'
626 log_nml(param_atmos_phy_mp_sn14)
632 wlabel(5) =
"GRAUPEL"
636 allocate(nc_uplim_d(1,ia,ja))
637 nc_uplim_d(:,:,:) = 150.e6_rp
647 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
674 integer,
intent(in) :: ka, ks, ke
675 integer,
intent(in) :: ia, is, ie
676 integer,
intent(in) :: ja, js, je
678 real(
rp),
intent(in) :: dens (ka,ia,ja)
679 real(
rp),
intent(in) :: w (ka,ia,ja)
680 real(
rp),
intent(in) :: qtrc (ka,ia,ja,
qa_mp)
681 real(
rp),
intent(in) :: pres(ka,ia,ja)
682 real(
rp),
intent(in) :: temp(ka,ia,ja)
683 real(
rp),
intent(in) :: qdry(ka,ia,ja)
684 real(
rp),
intent(in) :: cptot(ka,ia,ja)
685 real(
rp),
intent(in) :: cvtot(ka,ia,ja)
686 real(
rp),
intent(in) :: ccn (ka,ia,ja)
687 real(
dp),
intent(in) :: dt
688 real(
rp),
intent(in) :: cz( ka,ia,ja)
689 real(
rp),
intent(in) :: fz(0:ka,ia,ja)
691 real(
rp),
intent(out) :: rhoq_t (ka,ia,ja,
qa_mp)
692 real(
rp),
intent(out) :: rhoe_t (ka,ia,ja)
693 real(
rp),
intent(out) :: cptot_t(ka,ia,ja)
694 real(
rp),
intent(out) :: cvtot_t(ka,ia,ja)
695 real(
rp),
intent(out) :: evaporate(ka,ia,ja)
698 logical,
intent(in),
optional :: flg_lt
699 real(
rp),
intent(in),
optional :: d0_crg, v0_crg
700 real(
rp),
intent(in),
optional :: dqcrg(ka,ia,ja)
701 real(
rp),
intent(in),
optional :: beta_crg(ka,ia,ja)
702 real(
rp),
intent(in),
optional :: qtrc_crg(ka,ia,ja,hydro_max)
703 real(
rp),
intent(out),
optional :: qsplt_in(ka,ia,ja,3)
704 real(
rp),
intent(out),
optional :: sarea(ka,ia,ja,hydro_max)
705 real(
rp),
intent(out),
optional :: rhoqcrg_t(ka,ia,ja,hydro_max)
708 log_progress(*)
'atmosphere / physics / microphysics / SN14'
716 ka, ks, ke, ia, is, ie, ja, js, je, &
717 dens(:,:,:), w(:,:,:), qtrc(:,:,:,:), pres(:,:,:), temp(:,:,:), &
718 qdry(:,:,:), cptot(:,:,:), cvtot(:,:,:), ccn(:,:,:), &
719 real(dt,
rp), cz(:,:,:), fz(:,:,:), &
720 rhoq_t(:,:,:,:), rhoe_t(:,:,:), cptot_t(:,:,:), cvtot_t(:,:,:), &
722 flg_lt, d0_crg, v0_crg, dqcrg(:,:,:), beta_crg(:,:,:), &
724 qsplt_in(:,:,:,:), sarea(:,:,:,:), rhoqcrg_t(:,:,:,:) )
738 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
743 integer,
intent(in) :: ka, ks, ke
744 integer,
intent(in) :: ia, is, ie
745 integer,
intent(in) :: ja, js, je
747 real(
rp),
intent(in) :: qtrc (ka,ia,ja,
qa_mp-1)
748 real(
rp),
intent(in) :: mask_criterion
750 real(
rp),
intent(out) :: cldfrac(ka,ia,ja)
753 integer :: k, i, j, iq
763 qhydro = qhydro + qtrc(k,i,j,iq)
765 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
777 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
778 DENS0, TEMP0, QTRC0, &
789 integer,
intent(in) :: ka, ks, ke
790 integer,
intent(in) :: ia, is, ie
791 integer,
intent(in) :: ja, js, je
793 real(
rp),
intent(in) :: dens0(ka,ia,ja)
794 real(
rp),
intent(in) :: temp0(ka,ia,ja)
795 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,i_qc:i_ng)
797 real(
rp),
intent(out) :: re (ka,ia,ja,
n_hyd)
806 real(
rp) :: dc_ave(ka)
807 real(
rp) :: dr_ave(ka)
811 real(
rp) :: ri2m(ka), ri3m(ka)
812 real(
rp) :: rs2m(ka), rs3m(ka)
813 real(
rp) :: rg2m(ka), rg3m(ka)
815 real(
rp),
parameter :: coef_fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
818 real(
rp),
parameter :: r2m_min=1.e-12_rp
819 real(
rp),
parameter :: um2cm = 100.0_rp
821 real(
rp) :: limitsw, zerosw
833 xc(k) = min( xc_max, max( xc_min, dens0(k,i,j)*qtrc0(k,i,j,i_qc)/(qtrc0(k,i,j,i_nc)+nc_min) ) )
834 xr(k) = min( xr_max, max( xr_min, dens0(k,i,j)*qtrc0(k,i,j,i_qr)/(qtrc0(k,i,j,i_nr)+nr_min) ) )
835 xi(k) = min( xi_max, max( xi_min, dens0(k,i,j)*qtrc0(k,i,j,i_qi)/(qtrc0(k,i,j,i_ni)+ni_min) ) )
836 xs(k) = min( xs_max, max( xs_min, dens0(k,i,j)*qtrc0(k,i,j,i_qs)/(qtrc0(k,i,j,i_ns)+ns_min) ) )
837 xg(k) = min( xg_max, max( xg_min, dens0(k,i,j)*qtrc0(k,i,j,i_qg)/(qtrc0(k,i,j,i_ng)+ng_min) ) )
842 dc_ave(k) = a_m(i_mp_qc) * xc(k)**b_m(i_mp_qc)
843 dr_ave(k) = a_m(i_mp_qr) * xr(k)**b_m(i_mp_qr)
848 rc = 0.5_rp * dc_ave(k)
849 limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
850 re(k,i,j,
i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
855 rr = 0.5_rp * dr_ave(k)
856 limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
857 re(k,i,j,
i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
861 ri2m(k) = pi * coef_rea2(i_mp_qi) * qtrc0(k,i,j,i_ni) * a_rea2(i_mp_qi) * xi(k)**b_rea2(i_mp_qi)
862 rs2m(k) = pi * coef_rea2(i_mp_qs) * qtrc0(k,i,j,i_ns) * a_rea2(i_mp_qs) * xs(k)**b_rea2(i_mp_qs)
863 rg2m(k) = pi * coef_rea2(i_mp_qg) * qtrc0(k,i,j,i_ng) * a_rea2(i_mp_qg) * xg(k)**b_rea2(i_mp_qg)
868 ri3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k)
869 rs3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k)
870 rg3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k)
875 zerosw = 0.5_rp - sign(0.5_rp, ri2m(k) - r2m_min )
876 re(k,i,j,
i_hi) = ri3m(k) / ( ri2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
881 zerosw = 0.5_rp - sign(0.5_rp, rs2m(k) - r2m_min )
882 re(k,i,j,
i_hs) = rs3m(k) / ( rs2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
887 zerosw = 0.5_rp - sign(0.5_rp, rg2m(k) - r2m_min )
888 re(k,i,j,
i_hg) = rg3m(k) / ( rg2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
892 re(k,i,j,
i_hh) = 0.0_rp
906 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
918 integer,
intent(in) :: ka, ks, ke
919 integer,
intent(in) :: ia, is, ie
920 integer,
intent(in) :: ja, js, je
922 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,
qa_mp-1)
924 real(
rp),
intent(out) :: qe (ka,ia,ja,
n_hyd)
934 qe(k,i,j,
i_hc) = qtrc0(k,i,j,i_mp_qc)
935 qe(k,i,j,
i_hr) = qtrc0(k,i,j,i_mp_qr)
936 qe(k,i,j,
i_hi) = qtrc0(k,i,j,i_mp_qi)
937 qe(k,i,j,
i_hs) = qtrc0(k,i,j,i_mp_qs)
938 qe(k,i,j,
i_hg) = qtrc0(k,i,j,i_mp_qg)
939 qe(k,i,j,
i_hh) = 0.0_rp
949 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
961 integer,
intent(in) :: ka, ks, ke
962 integer,
intent(in) :: ia, is, ie
963 integer,
intent(in) :: ja, js, je
965 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,
qa_mp-1)
967 real(
rp),
intent(out) :: ne (ka,ia,ja,
n_hyd)
977 ne(k,i,j,
i_hc) = qtrc0(k,i,j,i_mp_nc)
978 ne(k,i,j,
i_hr) = qtrc0(k,i,j,i_mp_nr)
979 ne(k,i,j,
i_hi) = qtrc0(k,i,j,i_mp_ni)
980 ne(k,i,j,
i_hs) = qtrc0(k,i,j,i_mp_ns)
981 ne(k,i,j,
i_hg) = qtrc0(k,i,j,i_mp_ng)
982 ne(k,i,j,
i_hh) = 0.0_rp
991 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1007 integer,
intent(in) :: ka, ks, ke
1008 integer,
intent(in) :: ia, is, ie
1009 integer,
intent(in) :: ja, js, je
1011 real(
rp),
intent(in) :: qe(ka,ia,ja,
n_hyd)
1013 real(
rp),
intent(out) :: qtrc(ka,ia,ja,
qa_mp-1)
1015 real(
rp),
intent(in),
optional :: qnum(ka,ia,ja,
n_hyd)
1017 real(
rp),
parameter :: dc = 20.e-6_rp
1018 real(
rp),
parameter :: dr = 200.e-6_rp
1019 real(
rp),
parameter :: di = 80.e-6_rp
1020 real(
rp),
parameter :: ds = 80.e-6_rp
1021 real(
rp),
parameter :: dg = 200.e-6_rp
1022 real(
rp),
parameter :: rhow = 1000.0_rp
1023 real(
rp),
parameter :: rhof = 100.0_rp
1024 real(
rp),
parameter :: rhog = 400.0_rp
1025 real(
rp),
parameter :: b = 3.0_rp
1038 qtrc(k,i,j,i_mp_qc) = qe(k,i,j,
i_hc)
1048 qtrc(k,i,j,i_mp_qr) = qe(k,i,j,
i_hr)
1058 qtrc(k,i,j,i_mp_qi) = qe(k,i,j,
i_hi)
1068 qtrc(k,i,j,i_mp_qs) = qe(k,i,j,
i_hs)
1078 qtrc(k,i,j,i_mp_qg) = qe(k,i,j,
i_hg) + qe(k,i,j,
i_hh)
1086 if (
present(qnum) )
then
1093 if ( qnum(k,i,j,
i_hc) .ne. undef )
then
1094 qtrc(k,i,j,i_mp_nc) = qnum(k,i,j,
i_hc)
1096 qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1107 if ( qnum(k,i,j,
i_hr) .ne. undef )
then
1108 qtrc(k,i,j,i_mp_nr) = qnum(k,i,j,
i_hr)
1110 qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1121 if ( qnum(k,i,j,
i_hi) .ne. undef )
then
1122 qtrc(k,i,j,i_mp_ni) = qnum(k,i,j,
i_hi)
1124 qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1135 if ( qnum(k,i,j,
i_hs) .ne. undef )
then
1136 qtrc(k,i,j,i_mp_ns) = qnum(k,i,j,
i_hs)
1138 qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1149 if ( qnum(k,i,j,
i_hg) .ne. undef )
then
1150 if ( qnum(k,i,j,
i_hh) .ne. undef )
then
1151 qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,
i_hg) + qnum(k,i,j,
i_hh)
1153 qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,
i_hg)
1156 qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1169 qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1179 qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1189 qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1199 qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1209 qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1235 integer,
intent(in) :: ka, ks, ke
1237 real(
rp),
intent(in) :: rhoq(ka,i_qc:i_ng)
1238 real(
rp),
intent(in) :: dens(ka)
1239 real(
rp),
intent(in) :: temp(ka)
1240 real(
rp),
intent(in) :: pres(ka)
1242 real(
rp),
intent(out) :: vterm(ka,
qa_mp-1)
1247 real(
rp) :: rhofac_q(hydro_max)
1249 real(
rp) :: rlambdar
1259 integer :: k, i, j, iq
1262 mud_r = 3.0_rp * nu(i_mp_qr) + 2.0_rp
1266 rhofac = rho_0 / max( dens(k), rho_min )
1269 rhofac_q(i_mp_qc) = rhofac ** gamma_v(i_mp_qc)
1270 xq = max( xqmin(i_mp_qc), min( xqmax(i_mp_qc), rhoq(k,i_qc) / ( rhoq(k,i_nc) + nqmin(i_mp_qc) ) ) )
1272 vterm(k,i_mp_qc) = -rhofac_q(i_mp_qc) * coef_vt1(i_mp_qc,1) * xq**beta_v(i_mp_qc,1)
1274 vterm(k,i_mp_nc) = -rhofac_q(i_mp_qc) * coef_vt0(i_mp_qc,1) * xq**beta_vn(i_mp_qc,1)
1277 rhofac_q(i_mp_qr) = rhofac ** gamma_v(i_mp_qr)
1278 xq = max( xqmin(i_mp_qr), min( xqmax(i_mp_qr), rhoq(k,i_qr) / ( rhoq(k,i_nr) + nqmin(i_mp_qr) ) ) )
1280 rlambdar = a_m(i_mp_qr) * xq**b_m(i_mp_qr) &
1281 * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
1282 dq = ( 4.0_rp + mud_r ) * rlambdar
1284 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
1286 velq_s = coef_vtr_ar2 * dq &
1287 * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-5.0_rp-mud_r) )
1288 velq_l = coef_vtr_ar1 - coef_vtr_br1 &
1289 * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-4.0_rp-mud_r)
1290 vterm(k,i_mp_qr) = -rhofac_q(i_mp_qr) &
1291 * ( velq_l * ( weight ) &
1292 + velq_s * ( 1.0_rp - weight ) )
1294 dq = ( 1.0_rp + mud_r ) * rlambdar
1295 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
1297 velq_s = coef_vtr_ar2 * dql &
1298 * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-2.0_rp-mud_r) )
1299 velq_l = coef_vtr_ar1 - coef_vtr_br1 &
1300 * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-1.0_rp-mud_r)
1301 vterm(k,i_mp_nr) = -rhofac_q(i_mp_qr) &
1302 * ( velq_l * ( weight ) &
1303 + velq_s * ( 1.0_rp - weight ) )
1306 rhofac_q(i_mp_qi) = ( pres(k)/pre0_vt )**a_pre0_vt * ( temp(k)/tem0_vt )**a_tem0_vt
1307 xq = max( xqmin(i_mp_qi), min( xqmax(i_mp_qi), rhoq(k,i_qi) / ( rhoq(k,i_ni) + nqmin(i_mp_qi) ) ) )
1309 tmp = a_m(i_mp_qi) * xq**b_m(i_mp_qi)
1310 dq = coef_dave_l(i_mp_qi) * tmp
1311 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_li ) ) ) )
1313 velq_s = coef_vt1(i_mp_qi,1) * xq**beta_v(i_mp_qi,1)
1314 velq_l = coef_vt1(i_mp_qi,2) * xq**beta_v(i_mp_qi,2)
1315 vterm(k,i_mp_qi) = -rhofac_q(i_mp_qi) &
1316 * ( velq_l * ( weight ) &
1317 + velq_s * ( 1.0_rp - weight ) )
1319 dq = coef_dave_n(i_mp_qi) * tmp
1320 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ni ) ) ) )
1322 velq_s = coef_vt0(i_mp_qi,1) * xq**beta_vn(i_mp_qi,1)
1323 velq_l = coef_vt0(i_mp_qi,2) * xq**beta_vn(i_mp_qi,2)
1324 vterm(k,i_mp_ni) = -rhofac_q(i_mp_qi) &
1325 * ( velq_l * ( weight ) &
1326 + velq_s * ( 1.0_rp - weight ) )
1329 rhofac_q(i_mp_qs) = rhofac_q(i_mp_qi)
1330 xq = max( xqmin(i_mp_qs), min( xqmax(i_mp_qs), rhoq(k,i_qs) / ( rhoq(k,i_ns) + nqmin(i_mp_qs) ) ) )
1332 tmp = a_m(i_mp_qs) * xq**b_m(i_mp_qs)
1333 dq = coef_dave_l(i_mp_qs) * tmp
1334 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ls ) ) ) )
1336 velq_s = coef_vt1(i_mp_qs,1) * xq**beta_v(i_mp_qs,1)
1337 velq_l = coef_vt1(i_mp_qs,2) * xq**beta_v(i_mp_qs,2)
1338 vterm(k,i_mp_qs) = -rhofac_q(i_mp_qs) &
1339 * ( velq_l * ( weight ) &
1340 + velq_s * ( 1.0_rp - weight ) )
1342 dq = coef_dave_n(i_mp_qs) * tmp
1343 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ns ) ) ) )
1345 velq_s = coef_vt0(i_mp_qs,1) * xq**beta_vn(i_mp_qs,1)
1346 velq_l = coef_vt0(i_mp_qs,2) * xq**beta_vn(i_mp_qs,2)
1347 vterm(k,i_mp_ns) = -rhofac_q(i_mp_qs) &
1348 * ( velq_l * ( weight ) &
1349 + velq_s * ( 1.0_rp - weight ) )
1352 rhofac_q(i_mp_qg) = rhofac_q(i_mp_qi)
1353 xq = max( xqmin(i_mp_qg), min( xqmax(i_mp_qg), rhoq(k,i_qg) / ( rhoq(k,i_ng) + nqmin(i_mp_qg) ) ) )
1355 tmp = a_m(i_mp_qg) * xq**b_m(i_mp_qg)
1356 dq = coef_dave_l(i_mp_qg) * tmp
1357 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_lg ) ) ) )
1359 velq_s = coef_vt1(i_mp_qg,1) * xq**beta_v(i_mp_qg,1)
1360 velq_l = coef_vt1(i_mp_qg,2) * xq**beta_v(i_mp_qg,2)
1361 vterm(k,i_mp_qg) = -rhofac_q(i_mp_qg) &
1362 * ( velq_l * ( weight ) &
1363 + velq_s * ( 1.0_rp - weight ) )
1365 dq = coef_dave_n(i_mp_qg) * tmp
1366 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ng ) ) ) )
1368 velq_s = coef_vt0(i_mp_qg,1) * xq**beta_vn(i_mp_qg,1)
1369 velq_l = coef_vt0(i_mp_qg,2) * xq**beta_vn(i_mp_qg,2)
1370 vterm(k,i_mp_ng) = -rhofac_q(i_mp_qg) &
1371 * ( velq_l * ( weight ) &
1372 + velq_s * ( 1.0_rp - weight ) )
1377 vterm(ks-1,iq) = vterm(ks,iq)
1387 subroutine mp_sn14_init
1394 real(
rp),
parameter :: eps_gamma=1.e-30_rp
1396 real(
rp) :: w1(hydro_max)
1397 real(
rp) :: w2(hydro_max)
1398 real(
rp) :: w3(hydro_max)
1399 real(
rp) :: w4(hydro_max)
1400 real(
rp) :: w5(hydro_max)
1401 real(
rp) :: w6(hydro_max)
1402 real(
rp) :: w7(hydro_max)
1403 real(
rp) :: w8(hydro_max)
1406 real(
rp) :: ar_ice_fix = 0.7_rp
1407 real(
rp) :: wcap1, wcap2
1409 logical :: flag_vent0(hydro_max), flag_vent1(hydro_max)
1411 integer :: iw, ia, ib
1414 namelist / param_atmos_phy_mp_sn14_init / &
1420 ntmax_phase_change, &
1423 namelist / param_atmos_phy_mp_sn14_particles / &
1424 a_m, b_m, alpha_v, beta_v, gamma_v, &
1425 alpha_vn, beta_vn, &
1426 a_area, b_area, cap, &
1428 opt_m96_column_ice, &
1432 namelist / param_atmos_phy_mp_sn14_nucleation / &
1435 nm_m92, am_m92, bm_m92, &
1440 nucl_twomey, inucl_w
1442 namelist / param_atmos_phy_mp_sn14_collection / &
1443 dc0, dc1, di0, ds0, dg0, &
1444 sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
1448 e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
1449 i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
1452 namelist / param_atmos_phy_mp_sn14_condensation / &
1453 opt_fix_taucnd_c, fac_cndc
1458 alpha_v(:,:) = undef8
1459 beta_v(:,:) = undef8
1460 alpha_vn(:,:) = undef8
1461 beta_vn(:,:) = undef8
1479 coef_dave_n(:) = undef8
1480 coef_dave_l(:) = undef8
1484 coef_d2v(:) = undef8
1485 coef_md2v(:) = undef8
1489 coef_rea2(:) = undef8
1490 coef_rea3(:) = undef8
1493 coef_lambda(:) = undef8
1494 coef_vt0(:,:) = undef8
1495 coef_vt1(:,:) = undef8
1496 delta_b0(:) = undef8
1497 delta_b1(:) = undef8
1498 delta_ab0(:,:) = undef8
1499 delta_ab1(:,:) = undef8
1500 theta_b0(:) = undef8
1501 theta_b1(:) = undef8
1502 theta_ab0(:,:) = undef8
1503 theta_ab1(:,:) = undef8
1505 ah_vent(:,:) = undef8
1506 ah_vent0(:,:) = undef8
1507 ah_vent1(:,:) = undef8
1508 bh_vent(:,:) = undef8
1509 bh_vent0(:,:) = undef8
1510 bh_vent1(:,:) = undef8
1514 read(
io_fid_conf,nml=param_atmos_phy_mp_sn14_init,iostat=ierr)
1517 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.'
1518 elseif( ierr > 0 )
then
1519 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_init. Check!'
1522 log_nml(param_atmos_phy_mp_sn14_init)
1528 a_area(i_mp_qc) = pi/4.0_rp
1529 a_area(i_mp_qr) = pi/4.0_rp
1530 a_area(i_mp_qi) = 0.65_rp*1.e-4_rp*100.0_rp**(2.00_rp)
1531 a_area(i_mp_qs) = 0.2285_rp*1.e-4_rp*100.0_rp**(1.88_rp)
1532 a_area(i_mp_qg) = 0.50_rp*1.e-4_rp*100.0_rp**(2.0_rp)
1533 b_area(i_mp_qc) = 2.0_rp
1534 b_area(i_mp_qr) = 2.0_rp
1535 b_area(i_mp_qi) = 2.0_rp
1536 b_area(i_mp_qs) = 1.88_rp
1537 b_area(i_mp_qg) = 2.0_rp
1543 a_m(i_mp_qc) = 0.124_rp
1544 a_m(i_mp_qr) = 0.124_rp
1545 a_m(i_mp_qi) = 0.217_rp
1546 a_m(i_mp_qs) = 8.156_rp
1547 a_m(i_mp_qg) = 0.190_rp
1548 b_m(i_mp_qc) = 1.0_rp/3.0_rp
1549 b_m(i_mp_qr) = 1.0_rp/3.0_rp
1550 b_m(i_mp_qi) = 0.302_rp
1551 b_m(i_mp_qs) = 0.526_rp
1552 b_m(i_mp_qg) = 0.323_rp
1556 alpha_v(i_mp_qc,:)= 3.75e+5_rp
1557 alpha_v(i_mp_qr,:)= 159.0_rp
1558 alpha_v(i_mp_qi,:)= 317.0_rp
1559 alpha_v(i_mp_qs,:)= 27.70_rp
1560 alpha_v(i_mp_qg,:)= 40.0_rp
1561 beta_v(i_mp_qc,:) = 2.0_rp/3.0_rp
1562 beta_v(i_mp_qr,:) = 0.266_rp
1563 beta_v(i_mp_qi,:) = 0.363_rp
1564 beta_v(i_mp_qs,:) = 0.216_rp
1565 beta_v(i_mp_qg,:) = 0.230_rp
1566 gamma_v(i_mp_qc) = 1.0_rp
1568 gamma_v(i_mp_qr) = 1.0_rp/2.0_rp
1569 gamma_v(i_mp_qi) = 1.0_rp/2.0_rp
1570 gamma_v(i_mp_qs) = 1.0_rp/2.0_rp
1571 gamma_v(i_mp_qg) = 1.0_rp/2.0_rp
1579 nu(i_mp_qc) = 1.0_rp
1580 nu(i_mp_qr) = -1.0_rp/3.0_rp
1581 nu(i_mp_qi) = 1.0_rp
1582 nu(i_mp_qs) = 1.0_rp
1583 nu(i_mp_qg) = 1.0_rp
1585 mu(i_mp_qc) = 1.0_rp
1586 mu(i_mp_qr) = 1.0_rp/3.0_rp
1587 mu(i_mp_qi) = 1.0_rp/3.0_rp
1588 mu(i_mp_qs) = 1.0_rp/3.0_rp
1589 mu(i_mp_qg) = 1.0_rp/3.0_rp
1598 cap(i_mp_qc) = 2.0_rp
1599 cap(i_mp_qr) = 2.0_rp
1601 cap(i_mp_qs) = 2.0_rp
1602 cap(i_mp_qg) = 2.0_rp
1604 alpha_vn(:,:) = alpha_v(:,:)
1605 beta_vn(:,:) = beta_v(:,:)
1613 read(
io_fid_conf,nml=param_atmos_phy_mp_sn14_particles,iostat=ierr)
1615 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.'
1616 elseif( ierr > 0 )
then
1617 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_particles. Check!'
1620 log_nml(param_atmos_phy_mp_sn14_particles)
1624 if( opt_m96_ice )
then
1628 a_area(i_mp_qi) = 0.120284936_rp
1629 a_area(i_mp_qs) = 0.131488_rp
1630 a_area(i_mp_qg) = 0.5_rp
1631 b_area(i_mp_qi) = 1.850000_rp
1632 b_area(i_mp_qs) = 1.880000_rp
1633 b_area(i_mp_qg) = 2.0_rp
1634 a_m(i_mp_qi) = 1.23655360084766_rp
1635 a_m(i_mp_qs) = a_m(i_mp_qi)
1636 a_m(i_mp_qg) = 0.346111225718402_rp
1637 b_m(i_mp_qi) = 0.408329930583912_rp
1638 b_m(i_mp_qs) = b_m(i_mp_qi)
1639 b_m(i_mp_qg) = 0.357142857142857_rp
1641 if( opt_m96_column_ice )
then
1644 a_area(i_mp_qi)= (0.684_rp*1.e-4_rp)*10.0_rp**(2.0_rp*2.00_rp)
1645 b_area(i_mp_qi)= 2.0_rp
1646 a_m(i_mp_qi) = 0.19834046116844_rp
1647 b_m(i_mp_qi) = 0.343642611683849_rp
1650 wcap1 = sqrt(1.0_rp-ar_ice_fix**2)
1651 wcap2 = log( (1.0_rp+wcap1)/ar_ice_fix )
1652 cap(i_mp_qi) = 2.0_rp*wcap2/wcap1
1661 alpha_v(i_mp_qi,:) =(/ 5798.60107421875_rp, 167.347076416016_rp/)
1662 alpha_vn(i_mp_qi,:) =(/ 12408.177734375_rp, 421.799865722656_rp/)
1663 if( opt_m96_column_ice )
then
1664 alpha_v(i_mp_qi,:) = (/2901.0_rp, 32.20_rp/)
1665 alpha_vn(i_mp_qi,:) = (/9675.2_rp, 64.16_rp/)
1667 alpha_v(i_mp_qs,:) =(/ 15173.3916015625_rp, 305.678619384766_rp/)
1668 alpha_vn(i_mp_qs,:) =(/ 29257.1601562500_rp, 817.985717773438_rp/)
1669 alpha_v(i_mp_qg,:) =(/ 15481.6904296875_rp, 311.642242431641_rp/)
1670 alpha_vn(i_mp_qg,:) =(/ 27574.6562500000_rp, 697.536132812500_rp/)
1672 beta_v(i_mp_qi,:) =(/ 0.504873454570770_rp, 0.324817866086960_rp/)
1673 beta_vn(i_mp_qi,:) =(/ 0.548495233058929_rp, 0.385287821292877_rp/)
1674 if( opt_m96_column_ice )
then
1675 beta_v(i_mp_qi,:) =(/ 0.465552181005478_rp, 0.223826110363007_rp/)
1676 beta_vn(i_mp_qi,:) =(/ 0.530453503131866_rp, 0.273761242628098_rp/)
1678 beta_v(i_mp_qs,:) =(/ 0.528109610080719_rp, 0.329863965511322_rp/)
1679 beta_vn(i_mp_qs,:) =(/ 0.567154467105865_rp, 0.393876969814301_rp/)
1680 beta_v(i_mp_qg,:) =(/ 0.534656763076782_rp, 0.330253750085831_rp/)
1681 beta_vn(i_mp_qg,:) =(/ 0.570551633834839_rp, 0.387124240398407_rp/)
1685 ax_area(:) = a_area(:)*a_m(:)**b_area(:)
1686 bx_area(:) = b_area(:)*b_m(:)
1690 a_rea(:) = sqrt(ax_area(:)/pi)
1691 b_rea(:) = bx_area(:)/2.0_rp
1692 a_rea2(:) = a_rea(:)**2
1693 b_rea2(:) = b_rea(:)*2.0_rp
1694 a_rea3(:) = a_rea(:)**3
1695 b_rea3(:) = b_rea(:)*3.0_rp
1697 a_d2vt(:)=alpha_v(:,2)*(1.0_rp/alpha_v(:,2))**(beta_v(:,2)/b_m(:))
1698 b_d2vt(:)=(beta_v(:,2)/b_m(:))
1718 w1(iw) = gammafunc( (n+nu(iw)+1.0_rp)/mu(iw) )
1719 w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1720 w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1721 coef_m2(iw) = w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**n
1723 w4(iw) = gammafunc( (b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1724 coef_d(iw) = a_m(iw) * w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**b_m(iw)
1725 w5(iw) = gammafunc( (2.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
1726 w6(iw) = gammafunc( (3.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
1727 coef_d2v(iw) = a_m(iw) * w6(iw)/w5(iw)* ( w2(iw)/w3(iw) )**b_m(iw)
1728 coef_md2v(iw)= w5(iw)/w2(iw)* ( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw)+beta_v(iw,2))
1730 w7(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1731 coef_d3(iw) = a_m(iw)**3 * w7(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
1732 w8(iw) = gammafunc( (6.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1733 coef_d6(iw) = a_m(iw)**6 * w8(iw)/w2(iw)*( w2(iw)/w3(iw) )**(6.0_rp*b_m(iw))
1736 coef_deplc = coef_d(i_mp_qc)/a_m(i_mp_qc)
1742 w1(iw) = gammafunc( (2.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1743 w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1744 w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1746 w4(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1748 coef_r2(iw)=w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw))
1749 coef_r3(iw)=w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
1750 coef_re(iw)=coef_r3(iw)/coef_r2(iw)
1757 w1(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1758 w2(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1759 w3(iw) = gammafunc( (b_rea2(iw)+nu(iw)+1.0_rp)/mu(iw) )
1760 w4(iw) = gammafunc( (b_rea3(iw)+nu(iw)+1.0_rp)/mu(iw) )
1762 coef_rea2(iw) = w3(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea2(iw)
1763 coef_rea3(iw) = w4(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea3(iw)
1769 w1(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1770 w2(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1771 coef_a(iw) = mu(iw)/w1(iw)
1773 coef_lambda(iw) = (w1(iw)/w2(iw))**(-mu(iw))
1781 w1(iw) = gammafunc( (beta_vn(iw,ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
1782 w2(iw) = gammafunc( ( nu(iw) + 1.0_rp + n)/mu(iw) )
1783 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1784 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1786 coef_vt0(iw,ia)=alpha_vn(iw,ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_vn(iw,ia)
1788 w1(iw) = gammafunc( (beta_v(iw,ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
1789 w2(iw) = gammafunc( ( nu(iw) + 1.0_rp + n)/mu(iw) )
1791 coef_vt1(iw,ia)=alpha_v(iw,ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_v(iw,ia)
1796 w1(iw) = gammafunc( ( b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
1797 w2(iw) = gammafunc( (1.0_rp + b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
1798 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1799 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1800 coef_dave_n(iw) = (w1(iw)/w3(iw))*(w3(iw)/w4(iw))** b_m(iw)
1801 coef_dave_l(iw) = (w2(iw)/w3(iw))*(w3(iw)/w4(iw))**(1.0_rp+b_m(iw))
1805 ah_vent(i_mp_qc,1:2) = (/1.0000_rp,1.0000_rp/)
1806 ah_vent(i_mp_qr,1:2) = (/1.0000_rp,0.780_rp/)
1807 ah_vent(i_mp_qi,1:2) = (/1.0000_rp,0.860_rp/)
1808 ah_vent(i_mp_qs,1:2) = (/1.0000_rp,0.780_rp/)
1809 ah_vent(i_mp_qg,1:2) = (/1.0000_rp,0.780_rp/)
1810 bh_vent(i_mp_qc,1:2) = (/0.0000_rp,0.0000_rp/)
1811 bh_vent(i_mp_qr,1:2) = (/0.108_rp,0.308_rp/)
1812 bh_vent(i_mp_qi,1:2) = (/0.140_rp,0.280_rp/)
1813 bh_vent(i_mp_qs,1:2) = (/0.108_rp,0.308_rp/)
1814 bh_vent(i_mp_qg,1:2) = (/0.108_rp,0.308_rp/)
1818 if( (nu(iw) + b_m(iw) + n) > eps_gamma )
then
1819 w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
1820 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1821 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1822 ah_vent0(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1823 ah_vent0(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1824 flag_vent0(iw)=.true.
1826 ah_vent0(iw,1)= 1.0_rp
1827 ah_vent0(iw,2)= 1.0_rp
1828 flag_vent0(iw)=.false.
1831 if( (nu(iw) + b_m(iw) + n) > eps_gamma )
then
1832 w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
1833 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1834 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1835 ah_vent1(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1836 ah_vent1(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1837 flag_vent1(iw)=.true.
1839 ah_vent1(iw,1)= 1.0_rp
1840 ah_vent1(iw,2)= 1.0_rp
1841 flag_vent1(iw)=.true.
1846 if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )
then
1847 flag_vent0(iw)=.false.
1849 if(flag_vent0(iw))
then
1850 w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
1851 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1852 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1854 w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
1855 bh_vent0(iw,1)=bh_vent(iw,1)*(w4(iw)/w2(iw))*(w2(iw)/w3(iw))**(2.00_rp*b_m(iw)+beta_v(iw,1)+n-1.0_rp)
1856 w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
1857 bh_vent0(iw,2)=bh_vent(iw,2)*(w5(iw)/w2(iw))*(w2(iw)/w3(iw))**(1.5_rp*b_m(iw)+0.5_rp*beta_v(iw,2)+n-1.0_rp)
1859 bh_vent0(iw,1) = 0.0_rp
1860 bh_vent0(iw,2) = 0.0_rp
1864 if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )
then
1865 flag_vent1(iw)=.false.
1867 if(flag_vent1(iw))
then
1868 w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
1869 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1870 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1872 w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
1873 bh_vent1(iw,1)=bh_vent(iw,1)*(w4(iw)/w2(iw))*(w2(iw)/w3(iw))**(2.00_rp*b_m(iw)+beta_v(iw,1)+n-1.0_rp)
1875 w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
1876 bh_vent1(iw,2)=bh_vent(iw,2)*(w5(iw)/w2(iw))*(w2(iw)/w3(iw))**(1.5_rp*b_m(iw)+0.5_rp*beta_v(iw,2)+n-1.0_rp)
1878 bh_vent1(iw,1) = 0.0_rp
1879 bh_vent1(iw,2) = 0.0_rp
1888 w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1889 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1890 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1891 delta_b0(iw) = w1(iw)/w2(iw) &
1892 *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
1894 w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1895 delta_b1(iw) = w1(iw)/w2(iw) &
1896 *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
1902 w1(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1903 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1904 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1905 w4(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1907 w5(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1915 delta_ab0(ia,ib) = 2.0_rp*(w1(ib)/w2(ib))*(w4(ia)/w2(ia)) &
1916 * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
1917 * ( w2(ia)/w3(ia) )**(b_rea(ia) )
1919 delta_ab1(ia,ib) = 2.0_rp*(w5(ib)/w2(ib))*(w4(ia)/w2(ia)) &
1920 * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
1921 * ( w2(ia)/w3(ia) )**(b_rea(ia) )
1929 w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1930 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1931 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1932 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1933 theta_b0(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
1935 w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1936 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1937 theta_b1(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
1944 w1(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1945 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1946 w3(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1947 w4(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1949 w5(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1950 w6(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1952 w7(iw) = gammafunc( (beta_v(iw,2) + b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1953 w8(iw) = gammafunc( ( b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1958 theta_ab0(ia,ib) = 2.0_rp * (w1(ib)/w2(ib))*(w3(ia)/w4(ia)) &
1959 * (w5(ia)/w6(ia))**beta_v(ia,2) &
1960 * (w5(ib)/w6(ib))**beta_v(ib,2)
1961 theta_ab1(ia,ib) = 2.0_rp * (w7(ib)/w8(ib))*(w3(ia)/w4(ia)) &
1962 * (w5(ia)/w6(ia))**beta_v(ia,2) &
1963 * (w5(ib)/w6(ib))**beta_v(ib,2)
1968 read(
io_fid_conf, nml=param_atmos_phy_mp_sn14_nucleation, iostat=ierr)
1970 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.'
1971 elseif( ierr > 0 )
then
1972 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_nucleation. Check!'
1975 log_nml(param_atmos_phy_mp_sn14_nucleation)
1976 if ( mp_couple_aerosol .AND. nucl_twomey )
then
1977 log_error(
"ATMOS_PHY_MP_SN14_nucleation_kij",*)
"nucl_twomey should be false when MP_couple_aerosol is true, stop"
1983 read(
io_fid_conf, nml=param_atmos_phy_mp_sn14_collection, iostat=ierr )
1985 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.'
1986 elseif( ierr > 0 )
then
1987 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_collection. Check!'
1990 log_nml(param_atmos_phy_mp_sn14_collection)
1995 read (
io_fid_conf,nml=param_atmos_phy_mp_sn14_condensation, iostat=ierr )
1997 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.'
1998 elseif( ierr > 0 )
then
1999 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_condensation. Check!'
2002 log_nml(param_atmos_phy_mp_sn14_condensation)
2007 xc_cr = (2.0_rp*rc_cr/a_m(i_mp_qc))**(1.0_rp/b_m(i_mp_qc))
2008 alpha = (nu(i_mp_qc)+1.0_rp)/mu(i_mp_qc)
2009 gm = gammafunc(alpha)
2014 log_info(
"ATMOS_PHY_MP_sn14_init",
'(100a16)')
"LABEL ",wlabel(:)
2015 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"capacity ",cap(:)
2016 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_m2 ",coef_m2(:)
2017 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d ",coef_d(:)
2019 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d3 ",coef_d3(:)
2020 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d6 ",coef_d6(:)
2021 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d2v ",coef_d2v(:)
2022 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_md2v ",coef_md2v(:)
2023 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_d2vt ",a_d2vt(:)
2024 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_d2vt ",b_d2vt(:)
2026 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_r2 ",coef_r2(:)
2027 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_r3 ",coef_r3(:)
2028 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_re ",coef_re(:)
2030 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_area ",a_area(:)
2031 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_area ",b_area(:)
2032 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ax_area ",ax_area(:)
2033 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bx_area ",bx_area(:)
2034 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_rea ",a_rea(:)
2035 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_rea ",b_rea(:)
2036 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_rea3 ",a_rea3(:)
2037 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_rea3 ",b_rea3(:)
2039 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_rea2 ",coef_rea2(:)
2040 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_rea3 ",coef_rea3(:)
2041 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_vt0 ",coef_vt0(:,1)
2042 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_vt1 ",coef_vt1(:,1)
2043 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_A ",coef_a(:)
2044 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_lambda ",coef_lambda(:)
2046 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent0 sml",ah_vent0(:,1)
2047 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent0 lrg",ah_vent0(:,2)
2048 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent1 sml",ah_vent1(:,1)
2049 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent1 lrg",ah_vent1(:,2)
2050 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent0 sml",bh_vent0(:,1)
2051 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent0 lrg",bh_vent0(:,2)
2052 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent1 sml",bh_vent1(:,1)
2053 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent1 lrg",bh_vent1(:,2)
2055 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"delta_b0 ",delta_b0(:)
2056 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"delta_b1 ",delta_b1(:)
2057 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"theta_b0 ",theta_b0(:)
2058 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"theta_b1 ",theta_b1(:)
2061 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,a10,a,100ES16.6)')
"delta0(a,b)=(",trim(wlabel(ia)),
",b)=",(delta_ab0(ia,ib),ib=1,hydro_max)
2064 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,a10,a,100ES16.6)')
"delta1(a,b)=(",trim(wlabel(ia)),
",b)=",(delta_ab1(ia,ib),ib=1,hydro_max)
2067 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,a10,a,100ES16.6)')
"theta0(a,b)=(",trim(wlabel(ia)),
",b)=",(theta_ab0(ia,ib),ib=1,hydro_max)
2070 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,a10,a,100ES16.6)')
"theta1(a,b)=(",trim(wlabel(ia)),
",b)=",(theta_ab1(ia,ib),ib=1,hydro_max)
2074 end subroutine mp_sn14_init
2076 subroutine mp_sn14 ( &
2077 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2113 moist_psat_liq => atmos_saturation_psat_liq, &
2114 moist_psat_ice => atmos_saturation_psat_ice
2117 integer,
intent(in) :: ka, ks, ke
2118 integer,
intent(in) :: ia, is, ie
2119 integer,
intent(in) :: ja, js, je
2121 real(
rp),
intent(in) :: dens (ka,ia,ja)
2122 real(
rp),
intent(in) :: w (ka,ia,ja)
2123 real(
rp),
intent(in) :: qtrc (ka,ia,ja,
qa_mp)
2124 real(
rp),
intent(in) :: pres0 (ka,ia,ja)
2125 real(
rp),
intent(in) :: temp0 (ka,ia,ja)
2126 real(
rp),
intent(in) :: qdry (ka,ia,ja)
2127 real(
rp),
intent(in) :: cptot0(ka, ia, ja)
2128 real(
rp),
intent(in) :: cvtot0(ka, ia, ja)
2129 real(
rp),
intent(in) :: ccn (ka,ia,ja)
2130 real(
rp),
intent(in) :: dt
2131 real(
rp),
intent(in) :: cz( ka,ia,ja)
2132 real(
rp),
intent(in) :: fz(0:ka,ia,ja)
2134 real(
rp),
intent(out) :: rhoq_t(ka,ia,ja,
qa_mp)
2135 real(
rp),
intent(out) :: rhoe_t(ka,ia,ja)
2136 real(
rp),
intent(out) :: cptot_t(ka,ia,ja)
2137 real(
rp),
intent(out) :: cvtot_t(ka,ia,ja)
2139 real(
rp),
intent(out) :: evaporate(ka,ia,ja)
2142 logical,
intent(in),
optional :: flg_lt
2143 real(
rp),
intent(in),
optional :: d0_crg, v0_crg
2144 real(
rp),
intent(in),
optional :: dqcrg(ka,ia,ja)
2145 real(
rp),
intent(in),
optional :: beta_crg(ka,ia,ja)
2146 real(
rp),
intent(in),
optional :: qtrc_crg(ka,ia,ja,hydro_max)
2147 real(
rp),
intent(out),
optional :: qsplt_in(ka,ia,ja,3)
2148 real(
rp),
intent(out),
optional :: sarea(ka,ia,ja,hydro_max)
2149 real(
rp),
intent(out),
optional :: rhoqcrg_t_mp(ka,ia,ja,hydro_max)
2151 real(
rp) :: pres (ka)
2152 real(
rp) :: temp (ka)
2153 real(
rp) :: cva (ka)
2154 real(
rp) :: cpa (ka)
2155 real(
rp) :: rrho (ka)
2156 real(
rp) :: rhoe (ka)
2157 real(
rp) :: rhoq (ka,i_qv:i_ng)
2158 real(
rp) :: rhoq2(ka,i_qv:i_ng)
2160 real(
rp) :: rhoq0_t (ka,
qa_mp)
2161 real(
rp) :: rhoe0_t (ka)
2162 real(
rp) :: cptot0_t(ka)
2163 real(
rp) :: cvtot0_t(ka)
2165 real(
rp) :: xq(ka,hydro_max)
2167 real(
rp) :: dq_xa(ka,hydro_max)
2168 real(
rp) :: vt_xa(ka,hydro_max,2)
2170 real(
rp) :: wtemp(ka)
2171 real(
rp) :: esw (ka)
2172 real(
rp) :: esi (ka)
2175 real(
rp) :: rho_fac_q(ka,hydro_max)
2178 real(
rp) :: drhoqc, drhonc
2179 real(
rp) :: drhoqr, drhonr
2180 real(
rp) :: drhoqi, drhoni
2181 real(
rp) :: drhoqs, drhons
2182 real(
rp) :: drhoqg, drhong
2185 real(
rp) :: pq(ka,pq_max)
2186 real(
rp) :: wrm_dqc, wrm_dnc
2187 real(
rp) :: wrm_dqr, wrm_dnr
2190 real(
rp) :: pac(ka,pac_max)
2191 real(
rp) :: gc_dqc, gc_dnc
2192 real(
rp) :: sc_dqc, sc_dnc
2193 real(
rp) :: ic_dqc, ic_dnc
2194 real(
rp) :: rg_dqg, rg_dng
2195 real(
rp) :: rg_dqr, rg_dnr
2196 real(
rp) :: rs_dqr, rs_dnr, rs_dqs, rs_dns
2197 real(
rp) :: ri_dqr, ri_dnr
2198 real(
rp) :: ri_dqi, ri_dni
2199 real(
rp) :: ii_dqi, ii_dni
2200 real(
rp) :: is_dqi, is_dni, ss_dns
2201 real(
rp) :: gs_dqs, gs_dns, gg_dng
2204 real(
rp) :: clp_dqc, clp_dnc, clm_dqc, clm_dnc
2205 real(
rp) :: clp_dqr, clp_dnr, clm_dqr, clm_dnr
2206 real(
rp) :: clp_dqi, clp_dni, clm_dqi, clm_dni
2207 real(
rp) :: clp_dqs, clp_dns, clm_dqs, clm_dns
2208 real(
rp) :: clp_dqg, clp_dng, clm_dqg, clm_dng
2209 real(
rp) :: fac1, fac3, fac4, fac6, fac7, fac9, fac10
2212 real(
rp) :: pco_dqi, pco_dni
2213 real(
rp) :: pco_dqs, pco_dns
2214 real(
rp) :: pco_dqg, pco_dng
2217 real(
rp) :: eml_dqc, eml_dnc
2218 real(
rp) :: eml_dqr, eml_dnr
2219 real(
rp) :: eml_dqi, eml_dni
2220 real(
rp) :: eml_dqs, eml_dns
2221 real(
rp) :: eml_dqg, eml_dng
2224 real(
rp) :: spl_dqi, spl_dni
2225 real(
rp) :: spl_dqg, spl_dqs
2230 real(
rp) :: dtdt_equiv_d(ka)
2237 real(
rp) :: sl_plcdep
2238 real(
rp) :: sl_plrdep, sl_pnrdep
2240 real(
rp) :: qke_d(ka)
2242 real(
rp),
parameter :: eps = 1.e-19_rp
2243 real(
rp),
parameter :: eps_qv = 1.e-19_rp
2244 real(
rp),
parameter :: eps_rhoe = 1.e-19_rp
2245 real(
rp),
parameter :: eps_rho = 1.e-19_rp
2248 real(
rp) :: di2l, dtem
2253 integer :: k, i, j, iq
2255 real(
rp) :: dqv, dql, dqi
2256 real(
rp) :: dcv, dcp
2262 real(
rp) :: v0_crg_l=0.0_rp, d0_crg_l=0.0_rp
2263 real(
rp) :: facq(i_qc:i_qg), f_crg
2264 integer :: grid(2), pp, qq
2265 real(
rp) :: drhoqcrg_c, drhoqcrg_r
2266 real(
rp) :: drhoqcrg_i, drhoqcrg_s, drhoqcrg_g
2269 real(
rp) :: pcrg1(pq_max,ka)
2270 real(
rp) :: pcrg2(pcrg_max,ka)
2271 real(
rp) :: rhoq_crg(i_qc:i_qg,ka)
2272 real(
rp) :: rhoq2_crg(i_qc:i_qg,ka)
2273 real(
rp) :: qtrc0(ka,ia,ja,
qa_mp)
2274 real(
rp) :: crs(ka,ia,ja,hydro_max)
2275 real(
rp),
allocatable :: rhoqcrg0_t(:,:,:,:)
2277 real(
rp) :: crg_split_s
2278 real(
rp) :: crg_split_g
2279 real(
rp) :: crg_split_i
2280 real(
rp) :: wrm_dnc_crg
2281 real(
rp) :: wrm_dnr_crg
2282 real(
rp) :: gc_dnc_crg
2283 real(
rp) :: sc_dnc_crg
2284 real(
rp) :: ic_dnc_crg
2285 real(
rp) :: rg_dng_crg
2286 real(
rp) :: rg_dnr_crg
2287 real(
rp) :: rs_dnr_crg
2288 real(
rp) :: rs_dns_crg
2289 real(
rp) :: ri_dnr_crg
2290 real(
rp) :: ri_dni_crg
2291 real(
rp) :: ii_dni_crg
2292 real(
rp) :: is_dni_crg
2293 real(
rp) :: ss_dns_crg
2294 real(
rp) :: gs_dns_crg
2295 real(
rp) :: gi_dni_crg
2296 real(
rp) :: gg_dng_crg
2298 real(
rp) :: clp_dnc_crg, clm_dnc_crg
2299 real(
rp) :: clp_dnr_crg, clm_dnr_crg
2300 real(
rp) :: clp_dni_crg, clm_dni_crg
2301 real(
rp) :: clp_dns_crg, clm_dns_crg
2302 real(
rp) :: clp_dng_crg, clm_dng_crg
2304 real(
rp) :: pco_dni_crg
2305 real(
rp) :: pco_dns_crg
2306 real(
rp) :: pco_dng_crg
2308 real(
rp) :: eml_dnc_crg
2309 real(
rp) :: eml_dnr_crg
2310 real(
rp) :: eml_dni_crg
2311 real(
rp) :: eml_dns_crg
2312 real(
rp) :: eml_dng_crg
2314 real(
rp) :: spl_dni_crg
2315 real(
rp) :: spl_dns_crg
2316 real(
rp) :: spl_dng_crg
2320 real(
rp) :: sw1, sw2
2326 if (
present(flg_lt) )
then
2336 crg_split_i = 0.0_rp
2337 crg_split_s = 0.0_rp
2338 crg_split_g = 0.0_rp
2341 qsplt_in(:,:,:,:) = 0.0_rp
2342 allocate(rhoqcrg0_t(ka,ia,ja,hydro_max))
2343 rhoqcrg_t_mp(:,:,:,:) = 0.0_rp
2381 rhoq_t(k,i,j,:) = 0.0_rp
2382 rhoe_t(k,i,j) = 0.0_rp
2383 cptot_t(k,i,j) = 0.0_rp
2384 cvtot_t(k,i,j) = 0.0_rp
2389 cpa(k) = cptot0(k,i,j)
2390 cva(k) = cvtot0(k,i,j)
2391 pres(k) = pres0(k,i,j)
2392 temp(k) = temp0(k,i,j)
2411 rhoq(k,iq) = dens(k,i,j) * qtrc(k,i,j,iq)
2413 rhoq2(k,i_qv) = dens(k,i,j)*qtrc(k,i,j,i_qv)
2414 rhoq2(k,i_ni) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,i_ni) )
2415 rhoq2(k,i_nc) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,i_nc) )
2419 rrho(k) = 1.0_rp / dens(k,i,j)
2420 rhoe(k) = dens(k,i,j) * temp(k) * cva(k)
2421 wtemp(k) = max(temp(k), tem_min)
2426 1, i, j, temp(:), dens(:,i,j), pres(:), qtrc(:,i,j,i_qv) )
2430 rho_fac = rho_0 / max(dens(k,i,j),rho_min)
2431 rho_fac_q(k,i_mp_qc) = rho_fac**gamma_v(i_mp_qc)
2432 rho_fac_q(k,i_mp_qr) = rho_fac**gamma_v(i_mp_qr)
2433 rho_fac_q(k,i_mp_qi) = (pres(k)/pre0_vt)**a_pre0_vt * (temp(k)/tem0_vt)**a_tem0_vt
2434 rho_fac_q(k,i_mp_qs) = rho_fac_q(k,i_mp_qi)
2435 rho_fac_q(k,i_mp_qg) = rho_fac_q(k,i_mp_qi)
2450 dtdt_equiv_d(k) = 0.0_rp
2454 nc_uplim_d(1,i,j) = c_ccn*1.5_rp
2458 cz(:,i,j), fz(:,i,j), &
2459 w(:,i,j), dens(:,i,j), &
2460 wtemp(:), pres(:), qdry(:,i,j), &
2461 rhoq2(:,:), cpa(:), &
2464 ccn(:,i,j), nc_uplim_d(1,i,j), &
2471 drhoqc = dt * pq(k,i_lcccn)
2472 drhonc = dt * pq(k,i_ncccn)
2473 drhoqi = dt * pq(k,i_liccn)
2474 drhoni = dt * pq(k,i_niccn)
2475 drhoqv = max( -rhoq(k,i_qv), -drhoqc-drhoqi )
2478 fac1 = drhoqv / min( -drhoqc-drhoqi, -eps )
2480 drhoqc = drhoqc * fac1
2481 drhoqi = drhoqi * fac1
2483 rhoq0_t(k,i_qv) = drhoqv / dt
2484 rhoq0_t(k,i_qc) = drhoqc / dt
2485 rhoq0_t(k,i_qi) = drhoqi / dt
2486 rhoq0_t(k,i_nc) = drhonc / dt
2487 rhoq0_t(k,i_ni) = drhoni / dt
2489 rhoe0_t(k) = ( - lhv * drhoqv + lhf * drhoqi ) / dt
2491 dqv = rrho(k) * drhoqv
2492 dql = rrho(k) * drhoqc
2493 dqi = rrho(k) * drhoqi
2498 cvtot0_t(k) = dcv / dt
2499 cptot0_t(k) = dcp / dt
2504 rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe0_t(k)
2505 cvtot_t(k,i,j) = cvtot_t(k,i,j) + cvtot0_t(k)
2506 cptot_t(k,i,j) = cptot_t(k,i,j) + cptot0_t(k)
2511 rhoq(k,i_qv) = rhoq(k,i_qv) + rhoq0_t(k,i_qv)*dt
2512 rhoq(k,i_qc) = max(0.0_rp, rhoq(k,i_qc) + rhoq0_t(k,i_qc)*dt )
2513 rhoq(k,i_qi) = max(0.0_rp, rhoq(k,i_qi) + rhoq0_t(k,i_qi)*dt )
2514 rhoq(k,i_nc) = max(0.0_rp, rhoq(k,i_nc) + rhoq0_t(k,i_nc)*dt )
2515 rhoq(k,i_ni) = max(0.0_rp, rhoq(k,i_ni) + rhoq0_t(k,i_ni)*dt )
2518 rhoq(k,i_nc) = min( rhoq(k,i_nc), nc_uplim_d(1,i,j) )
2520 rhoe(k) = rhoe(k) + rhoe0_t(k)*dt
2521 cva(k) = cva(k) + cvtot0_t(k)*dt
2522 cpa(k) = cpa(k) + cptot0_t(k)*dt
2524 temp(k) = rhoe(k) / ( dens(k,i,j) * cva(k) )
2525 pres(k) = dens(k,i,j) * (cpa(k)-cva(k)) * temp(k)
2526 wtemp(k) = max( temp(k), tem_min )
2532 2, i, j, temp(:), dens(:,i,j), pres(:), qtrc(:,i,j,i_qv) )
2542 rhoq2(k,i_qr) = rhoq(k,i_qr)
2543 rhoq2(k,i_nr) = rhoq(k,i_nr)
2544 xq(k,i_mp_qr) = max(xr_min, min(xr_max, rhoq2(k,i_qr)/(rhoq2(k,i_nr)+nr_min) ))
2546 dq_xa(k,i_mp_qr) = a_m(i_mp_qr)*xq(k,i_mp_qr)**b_m(i_mp_qr)
2547 vt_xa(k,i_mp_qr,1) = alpha_v(i_mp_qr,1)*(xq(k,i_mp_qr)**beta_v(i_mp_qr,1))*rho_fac_q(k,i_mp_qr)
2548 vt_xa(k,i_mp_qr,2) = vt_xa(k,i_mp_qr,1)
2552 rhoq2(k,i_qv) = rhoq(k,i_qv)
2553 rhoq2(k,i_qc) = rhoq(k,i_qc)
2554 rhoq2(k,i_qi) = rhoq(k,i_qi)
2555 rhoq2(k,i_qs) = rhoq(k,i_qs)
2556 rhoq2(k,i_qg) = rhoq(k,i_qg)
2558 rhoq2(k,i_nc) = rhoq(k,i_nc)
2559 rhoq2(k,i_ni) = rhoq(k,i_ni)
2560 rhoq2(k,i_ns) = rhoq(k,i_ns)
2561 rhoq2(k,i_ng) = rhoq(k,i_ng)
2566 xq(k,i_mp_qc) = min(xc_max, max(xc_min, rhoq2(k,i_qc)/(rhoq2(k,i_nc)+nc_min) ))
2567 xq(k,i_mp_qi) = min(xi_max, max(xi_min, rhoq2(k,i_qi)/(rhoq2(k,i_ni)+ni_min) ))
2568 xq(k,i_mp_qs) = min(xs_max, max(xs_min, rhoq2(k,i_qs)/(rhoq2(k,i_ns)+ns_min) ))
2569 xq(k,i_mp_qg) = min(xg_max, max(xg_min, rhoq2(k,i_qg)/(rhoq2(k,i_ng)+ng_min) ))
2572 dq_xa(k,i_mp_qc) = a_m(i_mp_qc)*xq(k,i_mp_qc)**b_m(i_mp_qc)
2573 dq_xa(k,i_mp_qi) = a_m(i_mp_qi)*xq(k,i_mp_qi)**b_m(i_mp_qi)
2574 dq_xa(k,i_mp_qs) = a_m(i_mp_qs)*xq(k,i_mp_qs)**b_m(i_mp_qs)
2575 dq_xa(k,i_mp_qg) = a_m(i_mp_qg)*xq(k,i_mp_qg)**b_m(i_mp_qg)
2578 vt_xa(k,i_mp_qc,1) = alpha_v(i_mp_qc,1)*(xq(k,i_mp_qc)**beta_v(i_mp_qc,1))*rho_fac_q(k,i_mp_qc)
2579 vt_xa(k,i_mp_qi,1) = alpha_v(i_mp_qi,1)*(xq(k,i_mp_qi)**beta_v(i_mp_qi,1))*rho_fac_q(k,i_mp_qi)
2580 vt_xa(k,i_mp_qs,1) = alpha_v(i_mp_qs,1)*(xq(k,i_mp_qs)**beta_v(i_mp_qs,1))*rho_fac_q(k,i_mp_qs)
2581 vt_xa(k,i_mp_qg,1) = alpha_v(i_mp_qg,1)*(xq(k,i_mp_qg)**beta_v(i_mp_qg,1))*rho_fac_q(k,i_mp_qg)
2582 vt_xa(k,i_mp_qc,2) = alpha_v(i_mp_qc,2)*(xq(k,i_mp_qc)**beta_v(i_mp_qc,2))*rho_fac_q(k,i_mp_qc)
2583 vt_xa(k,i_mp_qi,2) = alpha_v(i_mp_qi,2)*(xq(k,i_mp_qi)**beta_v(i_mp_qi,2))*rho_fac_q(k,i_mp_qi)
2584 vt_xa(k,i_mp_qs,2) = alpha_v(i_mp_qs,2)*(xq(k,i_mp_qs)**beta_v(i_mp_qs,2))*rho_fac_q(k,i_mp_qs)
2585 vt_xa(k,i_mp_qg,2) = alpha_v(i_mp_qg,2)*(xq(k,i_mp_qg)**beta_v(i_mp_qg,2))*rho_fac_q(k,i_mp_qg)
2588 rhoq_crg(:,:) = 0.0_rp
2589 rhoq2_crg(:,:) = 0.0_rp
2593 rhoq_crg(iq,k) = dens(k,i,j) * qtrc_crg(k,i,j,iq-1)
2594 rhoq2_crg(iq,k) = rhoq_crg(iq,k)
2599 call moist_psat_liq( ka, ks, ke, &
2601 call moist_psat_ice( ka, ks, ke, &
2607 rhoq2(:,:), xq(:,:), temp(:), &
2610 call dep_vapor_melt_ice( &
2612 dens(:,i,j), wtemp(:), pres(:), qdry(:,i,j), rhoq2(:,:), &
2613 esw(:), esi(:), xq(:,:), vt_xa(:,:,:), dq_xa(:,:), &
2620 call update_by_phase_change( &
2622 ntmax_phase_change, dt, &
2623 cz(:,i,j), fz(:,i,j), &
2626 dens(:,i,j), qdry(:,i,j), &
2628 rhoq2(:,:), pres(:), temp(:), &
2632 sl_plcdep, sl_plrdep, sl_pnrdep, &
2633 rhoq0_t(:,:), rhoe0_t(:), &
2634 cptot0_t(:), cvtot0_t(:), &
2636 rhoq2_crg(i_qc:i_qg,:), &
2637 rhoqcrg0_t(:,i,j,:) )
2639 call update_by_phase_change( &
2641 ntmax_phase_change, dt, &
2642 cz(:,i,j), fz(:,i,j), &
2645 dens(:,i,j), qdry(:,i,j), &
2647 rhoq2(:,:), pres(:), temp(:), &
2651 sl_plcdep, sl_plrdep, sl_pnrdep, &
2652 rhoq0_t(:,:), rhoe0_t(:), &
2653 cptot0_t(:), cvtot0_t(:), &
2659 rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe0_t(k)
2660 cvtot_t(k,i,j) = cvtot_t(k,i,j) + cvtot0_t(k)
2661 cptot_t(k,i,j) = cptot_t(k,i,j) + cptot0_t(k)
2667 rhoq(k,iq) = max(0.0_rp, rhoq(k,iq) + rhoq0_t(k,iq)*dt )
2672 rhoe(k) = rhoe(k) + rhoe0_t(k)*dt
2673 cva(k) = cva(k) + cvtot0_t(k)*dt
2674 cpa(k) = cpa(k) + cptot0_t(k)*dt
2675 temp(k) = rhoe(k) / ( dens(k,i,j) * cva(k) )
2676 pres(k) = dens(k,i,j) * ( cpa(k) - cva(k) ) * temp(k)
2682 rhoq_crg(iq,k) = rhoq_crg(iq,k) + rhoqcrg0_t(k,i,j,iq-1)*dt
2689 3, i, j, temp(:), dens(:,i,j), pres(:), qtrc(:,i,j,i_qv) )
2702 rhoq2(k,i_qv) = rhoq(k,i_qv)
2703 rhoq2(k,i_qc) = rhoq(k,i_qc)
2704 rhoq2(k,i_qr) = rhoq(k,i_qr)
2705 rhoq2(k,i_qi) = rhoq(k,i_qi)
2706 rhoq2(k,i_qs) = rhoq(k,i_qs)
2707 rhoq2(k,i_qg) = rhoq(k,i_qg)
2709 rhoq2(k,i_nc) = rhoq(k,i_nc)
2710 rhoq2(k,i_nr) = rhoq(k,i_nr)
2711 rhoq2(k,i_ni) = rhoq(k,i_ni)
2712 rhoq2(k,i_ns) = rhoq(k,i_ns)
2713 rhoq2(k,i_ng) = rhoq(k,i_ng)
2716 xq(k,i_mp_qc) = min(xc_max, max(xc_min, rhoq2(k,i_qc)/(rhoq2(k,i_nc)+nc_min) ) )
2717 xq(k,i_mp_qr) = min(xr_max, max(xr_min, rhoq2(k,i_qr)/(rhoq2(k,i_nr)+nr_min) ) )
2718 xq(k,i_mp_qi) = min(xi_max, max(xi_min, rhoq2(k,i_qi)/(rhoq2(k,i_ni)+ni_min) ) )
2719 xq(k,i_mp_qs) = min(xs_max, max(xs_min, rhoq2(k,i_qs)/(rhoq2(k,i_ns)+ns_min) ) )
2720 xq(k,i_mp_qg) = min(xg_max, max(xg_min, rhoq2(k,i_qg)/(rhoq2(k,i_ng)+ng_min) ) )
2723 dq_xa(k,i_mp_qc) = 2.0_rp*a_rea(i_mp_qc)*xq(k,i_mp_qc)**b_rea(i_mp_qc)
2724 dq_xa(k,i_mp_qr) = 2.0_rp*a_rea(i_mp_qr)*xq(k,i_mp_qr)**b_rea(i_mp_qr)
2725 dq_xa(k,i_mp_qi) = 2.0_rp*a_rea(i_mp_qi)*xq(k,i_mp_qi)**b_rea(i_mp_qi)
2726 dq_xa(k,i_mp_qs) = 2.0_rp*a_rea(i_mp_qs)*xq(k,i_mp_qs)**b_rea(i_mp_qs)
2727 dq_xa(k,i_mp_qg) = 2.0_rp*a_rea(i_mp_qg)*xq(k,i_mp_qg)**b_rea(i_mp_qg)
2731 vt_xa(k,i_mp_qc,2) = alpha_v(i_mp_qc,2)*(xq(k,i_mp_qc)**beta_v(i_mp_qc,2))*rho_fac_q(k,i_mp_qc)
2732 vt_xa(k,i_mp_qr,2) = alpha_v(i_mp_qr,2)*(xq(k,i_mp_qr)**beta_v(i_mp_qr,2))*rho_fac_q(k,i_mp_qr)
2733 vt_xa(k,i_mp_qi,2) = alpha_v(i_mp_qi,2)*(xq(k,i_mp_qi)**beta_v(i_mp_qi,2))*rho_fac_q(k,i_mp_qi)
2734 vt_xa(k,i_mp_qs,2) = alpha_v(i_mp_qs,2)*(xq(k,i_mp_qs)**beta_v(i_mp_qs,2))*rho_fac_q(k,i_mp_qs)
2735 vt_xa(k,i_mp_qg,2) = alpha_v(i_mp_qg,2)*(xq(k,i_mp_qg)**beta_v(i_mp_qg,2))*rho_fac_q(k,i_mp_qg)
2743 rhoq2_crg(iq,k) = rhoq_crg(iq,k)
2749 if ( mp_doautoconversion )
then
2754 rhoq2_crg(i_qc:i_qg,:), &
2755 xq(:,:), dq_xa(:,:), &
2763 pq(k,i_lcaut) = 0.0_rp
2764 pq(k,i_ncaut) = 0.0_rp
2765 pq(k,i_nraut) = 0.0_rp
2766 pq(k,i_lcacc) = 0.0_rp
2767 pq(k,i_ncacc) = 0.0_rp
2768 pq(k,i_nrslc) = 0.0_rp
2769 pq(k,i_nrbrk) = 0.0_rp
2771 pcrg1(i_lcaut,k) = 0.0_rp
2772 pcrg1(i_ncaut,k) = 0.0_rp
2773 pcrg1(i_nraut,k) = 0.0_rp
2774 pcrg1(i_lcacc,k) = 0.0_rp
2775 pcrg1(i_ncacc,k) = 0.0_rp
2776 pcrg1(i_nrslc,k) = 0.0_rp
2777 pcrg1(i_nrbrk,k) = 0.0_rp
2785 d0_crg_l, v0_crg_l, &
2788 temp(:), rhoq2(:,:), &
2789 rhoq2_crg(i_qc:i_qg,:), &
2790 xq(:,:), dq_xa(:,:), vt_xa(:,:,:), &
2800 temp(:), rhoq2(:,:), &
2801 rhoq2_crg(i_qc:i_qg,:), &
2811 wrm_dqc = max( dt*( pq(k,i_lcaut)+pq(k,i_lcacc) ), -rhoq2(k,i_qc) )
2812 wrm_dnc = max( dt*( pq(k,i_ncaut)+pq(k,i_ncacc) ), -rhoq2(k,i_nc) )
2813 wrm_dnr = max( dt*( pq(k,i_nraut)+pq(k,i_nrslc)+pq(k,i_nrbrk) ), -rhoq2(k,i_nr) )
2817 wrm_dnc_crg = dt*( pcrg1(i_ncaut,k)+pcrg1(i_ncacc,k) )
2819 sw1 = min( abs(rhoq2_crg(i_qc,k)),abs(wrm_dnc_crg) )
2820 wrm_dnc_crg = sign( sw1,wrm_dnc_crg )
2821 wrm_dnr_crg = - wrm_dnc_crg
2830 gc_dqc = max( dt*pac(k,i_lgaclc2lg), min(0.0_rp, -rhoq2(k,i_qc)-wrm_dqc ))
2831 sc_dqc = max( dt*pac(k,i_lsaclc2ls), min(0.0_rp, -rhoq2(k,i_qc)-wrm_dqc-gc_dqc ))
2832 ic_dqc = max( dt*pac(k,i_liaclc2li), min(0.0_rp, -rhoq2(k,i_qc)-wrm_dqc-gc_dqc-sc_dqc ))
2834 gc_dnc = max( dt*pac(k,i_ngacnc2ng), min(0.0_rp, -rhoq2(k,i_nc)-wrm_dnc ))
2835 sc_dnc = max( dt*pac(k,i_nsacnc2ns), min(0.0_rp, -rhoq2(k,i_nc)-wrm_dnc-gc_dnc ))
2836 ic_dnc = max( dt*pac(k,i_niacnc2ni), min(0.0_rp, -rhoq2(k,i_nc)-wrm_dnc-gc_dnc-sc_dnc ))
2839 gc_dnc_crg = dt*pcrg2(i_ngacnc2ng,k)
2840 sc_dnc_crg = dt*pcrg2(i_nsacnc2ns,k)
2841 ic_dnc_crg = dt*pcrg2(i_niacnc2ni,k)
2843 sw1 = min( abs(rhoq2_crg(i_qc,k)+wrm_dnc_crg ),abs(gc_dnc_crg) )
2844 gc_dnc_crg = sign( sw1,gc_dnc_crg )
2845 sw1 = min( abs(rhoq2_crg(i_qc,k)+wrm_dnc_crg+gc_dnc_crg ),abs(sc_dnc_crg) )
2846 sc_dnc_crg = sign( sw1,sc_dnc_crg )
2847 sw1 = min( abs(rhoq2_crg(i_qc,k)+wrm_dnc_crg+gc_dnc_crg+sc_dnc_crg),abs(ic_dnc_crg) )
2848 ic_dnc_crg = sign( sw1,ic_dnc_crg )
2852 sw = sign(0.5_rp, t00-temp(k)) + 0.5_rp
2853 rg_dqr = max( dt*pac(k,i_lraclg2lg ), min(0.0_rp, -rhoq2(k,i_qr)-wrm_dqr )) * sw
2854 rg_dqg = max( dt*pac(k,i_lraclg2lg ), min(0.0_rp, -rhoq2(k,i_qg) )) * ( 1.0_rp - sw )
2855 rs_dqr = max( dt*pac(k,i_lracls2lg_r), min(0.0_rp, -rhoq2(k,i_qr)-wrm_dqr-rg_dqr )) * sw
2856 ri_dqr = max( dt*pac(k,i_lracli2lg_r), min(0.0_rp, -rhoq2(k,i_qr)-wrm_dqr-rg_dqr-rs_dqr )) * sw
2858 rg_dnr = max( dt*pac(k,i_nracng2ng ), min(0.0_rp, -rhoq2(k,i_nr)-wrm_dnr )) * sw
2859 rg_dng = max( dt*pac(k,i_nracng2ng ), min(0.0_rp, -rhoq2(k,i_ng) )) * ( 1.0_rp - sw )
2860 rs_dnr = max( dt*pac(k,i_nracns2ng_r), min(0.0_rp, -rhoq2(k,i_nr)-wrm_dnr-rg_dnr )) * sw
2861 ri_dnr = max( dt*pac(k,i_nracni2ng_r), min(0.0_rp, -rhoq2(k,i_nr)-wrm_dnr-rg_dnr-rs_dnr )) * sw
2864 rg_dnr_crg = dt*pcrg2(i_nracng2ng, k)* sw
2865 rg_dng_crg = dt*pcrg2(i_nracng2ng, k)* ( 1.0_rp - sw )
2866 rs_dnr_crg = dt*pcrg2(i_nracns2ng_r,k)* sw
2867 ri_dnr_crg = dt*pcrg2(i_nracni2ng_r,k)* sw
2869 sw1 = min( abs(rhoq2_crg(i_qr,k)+wrm_dnr_crg),abs(rg_dnr_crg) )
2870 rg_dnr_crg = sign( sw1,rg_dnr_crg )
2871 sw1 = min( abs(rhoq2_crg(i_qg,k) ),abs(rg_dng_crg) )
2872 rg_dng_crg = sign( sw1,rg_dng_crg )
2873 sw1 = min( abs(rhoq2_crg(i_qr,k)+wrm_dnr_crg),abs(rs_dnr_crg) )
2874 rs_dnr_crg = sign( sw1,rs_dnr_crg )
2875 sw1 = min( abs(rhoq2_crg(i_qr,k)+wrm_dnr_crg),abs(ri_dnr_crg) )
2876 ri_dnr_crg = sign( sw1,ri_dnr_crg )
2879 fac1 = (ri_dqr-eps)/ (dt*pac(k,i_lracli2lg_r)-eps)
2880 ri_dqi = max( dt*pac(k,i_lracli2lg_i)*fac1, min(0.0_rp, -rhoq2(k,i_qi)+ic_dqc ))
2881 ii_dqi = max( dt*pac(k,i_liacli2ls ) , min(0.0_rp, -rhoq2(k,i_qi)+ic_dqc-ri_dqi ))
2882 is_dqi = max( dt*pac(k,i_liacls2ls ) , min(0.0_rp, -rhoq2(k,i_qi)+ic_dqc-ri_dqi-ii_dqi ))
2886 fac4 = (ri_dnr-eps)/ (dt*pac(k,i_nracni2ng_r)-eps)
2887 ri_dni = max( dt*pac(k,i_nracni2ng_i)*fac4, min(0.0_rp, -rhoq2(k,i_ni) ))
2888 ii_dni = max( dt*pac(k,i_niacni2ns ) , min(0.0_rp, -rhoq2(k,i_ni)-ri_dni ))
2889 is_dni = max( dt*pac(k,i_niacns2ns ) , min(0.0_rp, -rhoq2(k,i_ni)-ri_dni-ii_dni ))
2894 ri_dni_crg = dt*pcrg2(i_nracni2ng_i,k)*fac4
2895 ii_dni_crg = dt*pcrg2(i_niacni2ns,k)
2896 is_dni_crg = dt*pcrg2(i_niacns2ns,k)
2900 sw1 = min( abs(rhoq2_crg(i_qi,k)-ic_dnc_crg) ,abs(ri_dni_crg) )
2901 ri_dni_crg = sign( sw1,ri_dni_crg )
2902 sw1 = min( abs(rhoq2_crg(i_qi,k)-ic_dnc_crg+ri_dni_crg) ,abs(ii_dni_crg) )
2903 ii_dni_crg = sign( sw1,ii_dni_crg )
2904 sw1 = min( abs(rhoq2_crg(i_qi,k)-ic_dnc_crg+ri_dni_crg+ii_dni_crg),abs(is_dni_crg) )
2905 is_dni_crg = sign( sw1,is_dni_crg )
2911 fac3 = (rs_dqr-eps)/(dt*pac(k,i_lracls2lg_r)-eps)
2912 rs_dqs = max( dt*pac(k,i_lracls2lg_s)*fac3, min(0.0_rp, -rhoq2(k,i_qs)+sc_dqc+ii_dqi+is_dqi ))
2913 gs_dqs = max( dt*pac(k,i_lgacls2lg ) , min(0.0_rp, -rhoq2(k,i_qs)+sc_dqc+ii_dqi+is_dqi-rs_dqs ))
2915 fac6 = (rs_dnr-eps)/(dt*pac(k,i_nracns2ng_r)-eps)
2917 rs_dns = max( dt*pac(k,i_nracns2ng_s)*fac6, min(0.0_rp, -rhoq2(k,i_ns)+0.50_rp*ii_dni+is_dni ))
2918 gs_dns = max( dt*pac(k,i_ngacns2ng ) , min(0.0_rp, -rhoq2(k,i_ns)+0.50_rp*ii_dni+is_dni-rs_dns ))
2919 ss_dns = max( dt*pac(k,i_nsacns2ns ) , min(0.0_rp, -rhoq2(k,i_ns)+0.50_rp*ii_dni+is_dni-rs_dns-gs_dns ))
2921 gg_dng = max( dt*pac(k,i_ngacng2ng) , min(0.0_rp, -rhoq2(k,i_ng) ))
2924 rs_dns_crg = dt*pcrg2(i_nracns2ng_s,k)*fac6
2925 gs_dns_crg = dt*pcrg2(i_ngacns2ng,k)
2929 sw1 = min( abs(rhoq2_crg(i_qs,k)-sc_dnc_crg-ii_dni-is_dni_crg), abs(rs_dns_crg) )
2930 rs_dns_crg = sign( sw1,rs_dns_crg )
2931 sw1 = min( abs(rhoq2_crg(i_qs,k)-sc_dnc_crg-ii_dni-is_dni_crg+rs_dns_crg),abs(gs_dns_crg) )
2932 gs_dns_crg = sign( sw1,gs_dns_crg )
2934 sw1 = sign(0.5_rp, abs( pcrg2(i_cgngacns2ng,k) )-eps ) + 0.5_rp
2935 sw2 = sign(0.5_rp, abs( pcrg2(i_cgngacni2ng,k) )-eps ) + 0.5_rp
2936 crg_split_g = dt*pcrg2(i_cgngacns2ng,k)*sw1 &
2937 + dt*pcrg2(i_cgngacni2ng,k)*sw2
2938 crg_split_s = -dt*pcrg2(i_cgngacns2ng,k)*sw1
2940 qsplt_in(k,i,j,1) = crg_split_g / dt
2941 qsplt_in(k,i,j,3) = crg_split_s / dt
2942 qsplt_in(k,i,j,2) = crg_split_i / dt
2949 clp_dqr = (-rg_dqg-rs_dqs-ri_dqi) * (1.0_rp-sw)
2951 clp_dqs = -sc_dqc-ii_dqi-is_dqi
2952 clp_dqg = -gc_dqc -gs_dqs + (-rg_dqr-rs_dqr-rs_dqs-ri_dqr-ri_dqi) * sw
2957 clp_dns = -ii_dni*0.5_rp
2958 clp_dng = (-rs_dnr-ri_dnr) * sw
2961 clp_dnc_crg = 0.0_rp
2962 clp_dnr_crg = -rg_dng_crg*(1.0_rp-sw)
2963 clp_dni_crg = -ic_dnc_crg
2965 clp_dns_crg = -sc_dnc_crg-ii_dni_crg-is_dni_crg-ss_dns_crg &
2967 clp_dng_crg = -gc_dnc_crg+(-rg_dnr_crg-rs_dnr_crg-ri_dnr_crg)*sw &
2968 -ri_dni_crg-rs_dns_crg-gs_dns_crg-gg_dng_crg &
2975 clm_dqc = gc_dqc+sc_dqc+ic_dqc
2976 clm_dqr = (rg_dqr+rs_dqr+ri_dqr) * sw
2977 clm_dqi = ri_dqi+ii_dqi+is_dqi
2978 clm_dqs = rs_dqs+gs_dqs
2979 clm_dqg = rg_dqg * (1.0_rp-sw)
2981 clm_dnc = gc_dnc+sc_dnc+ic_dnc
2982 clm_dnr = (rg_dnr+rs_dnr+ri_dnr) * sw
2983 clm_dni = ri_dni+ii_dni+is_dni
2984 clm_dns = rs_dns+ss_dns+gs_dns
2985 clm_dng = gg_dng + rg_dng * (1.0_rp-sw)
2988 clm_dnc_crg = gc_dnc_crg+sc_dnc_crg+ic_dnc_crg
2989 clm_dnr_crg = (rg_dnr_crg+rs_dnr_crg+ri_dnr_crg) * sw
2990 clm_dni_crg = ri_dni_crg+ii_dni_crg+is_dni_crg
2992 clm_dns_crg = rs_dns_crg+gs_dns_crg+ss_dns_crg
2993 clm_dng_crg = gg_dng_crg+rg_dng_crg*(1.0_rp-sw)
2998 pco_dqi = max( dt*pq(k,i_licon), -clp_dqi )
2999 pco_dqs = max( dt*pq(k,i_lscon), -clp_dqs )
3000 pco_dqg = -pco_dqi-pco_dqs
3002 pco_dni = max( dt*pq(k,i_nicon), -clp_dni )
3003 pco_dns = max( dt*pq(k,i_nscon), -clp_dns )
3004 pco_dng = -pco_dni-pco_dns
3007 pco_dni_crg = dt*pcrg1(i_nicon,k)
3008 pco_dns_crg = dt*pcrg1(i_nscon,k)
3010 sw1 = min( abs(rhoq2_crg(i_qi,k)+clp_dni_crg ),abs(pco_dni_crg) )
3011 pco_dni_crg = sign( sw1,pco_dni_crg )
3012 sw1 = min( abs(rhoq2_crg(i_qs,k)+clp_dns_crg ),abs(pco_dns_crg) )
3013 pco_dns_crg = sign( sw1,pco_dns_crg )
3014 pco_dng_crg = -pco_dni_crg-pco_dns_crg
3018 eml_dqi = max( dt*pq(k,i_liacm), min(0.0_rp, -rhoq2(k,i_qi)-(clp_dqi+clm_dqi)-pco_dqi ))
3019 eml_dqs = max( dt*pq(k,i_lsacm), min(0.0_rp, -rhoq2(k,i_qs)-(clp_dqs+clm_dqs)-pco_dqs ))
3020 eml_dqg = max( dt*(pq(k,i_lgacm)+pq(k,i_lgarm)+pq(k,i_lsarm)+pq(k,i_liarm)), &
3021 min(0.0_rp, -rhoq2(k,i_qg)-(clp_dqg+clm_dqg)-pco_dqg ))
3023 eml_dqr = -eml_dqs-eml_dqg
3025 eml_dni = max( dt*pq(k,i_niacm), min(0.0_rp, -rhoq2(k,i_ni)-(clp_dni+clm_dni)-pco_dni ))
3026 eml_dns = max( dt*pq(k,i_nsacm), min(0.0_rp, -rhoq2(k,i_ns)-(clp_dns+clm_dns)-pco_dns ))
3027 eml_dng = max( dt*(pq(k,i_ngacm)+pq(k,i_ngarm)+pq(k,i_nsarm)+pq(k,i_niarm)), &
3028 min(0.0_rp, -rhoq2(k,i_ng)-(clp_dng+clm_dng)-pco_dng ))
3030 eml_dnr = -eml_dns-eml_dng
3033 eml_dni_crg = dt*pcrg1(i_niacm,k)
3034 eml_dns_crg = dt*pcrg1(i_nsacm,k)
3035 eml_dng_crg = dt*(pcrg1(i_ngacm,k)+pcrg1(i_ngarm,k)+pcrg1(i_nsarm,k)+pcrg1(i_niarm,k))
3037 sw1 = min( abs(rhoq2_crg(i_qi,k)+clp_dni_crg+clm_dni_crg+pco_dni_crg ),abs(eml_dni_crg) )
3038 eml_dni_crg = sign( sw1,eml_dni_crg )
3039 sw1 = min( abs(rhoq2_crg(i_qs,k)+clp_dns_crg+clm_dns_crg+pco_dns_crg ),abs(eml_dns_crg) )
3040 eml_dns_crg = sign( sw1,eml_dns_crg )
3041 sw1 = min( abs(rhoq2_crg(i_qg,k)+clp_dng_crg+clm_dng_crg+pco_dng_crg ),abs(eml_dng_crg) )
3042 eml_dng_crg = sign( sw1,eml_dng_crg )
3044 eml_dnc_crg = -eml_dni_crg
3045 eml_dnr_crg = -eml_dns_crg-eml_dng_crg
3049 spl_dqg = max( dt*pq(k,i_lgspl), min(0.0_rp, -rhoq2(k,i_qg)-(clp_dqg+clm_dqg)-pco_dqg-eml_dqg ))
3050 spl_dqs = max( dt*pq(k,i_lsspl), min(0.0_rp, -rhoq2(k,i_qs)-(clp_dqs+clm_dqs)-pco_dqs-eml_dqs ))
3051 spl_dqi = -spl_dqg-spl_dqs
3052 fac9 = (spl_dqg-eps)/(dt*pq(k,i_lgspl)-eps)
3053 fac10 = (spl_dqs-eps)/(dt*pq(k,i_lsspl)-eps)
3054 spl_dni = dt*pq(k,i_nispl)*fac9*fac10
3057 spl_dns_crg = dt*pcrg1(i_nsspl,k)*fac9*fac10
3058 spl_dng_crg = dt*pcrg1(i_ngspl,k)*fac9*fac10
3060 sw1 = min( abs(rhoq2_crg(i_qs,k)+clp_dns_crg+pco_dns_crg+eml_dns_crg ),abs(spl_dns_crg) )
3061 spl_dns_crg = sign( sw1,spl_dns_crg )
3062 sw1 = min( abs(rhoq2_crg(i_qg,k)+clp_dng_crg+pco_dng_crg+eml_dng_crg ),abs(spl_dng_crg) )
3063 spl_dng_crg = sign( sw1,spl_dng_crg )
3064 spl_dni_crg = -spl_dns_crg-spl_dng_crg
3070 di2l = clp_dqc + clp_dqr + clm_dqc + clm_dqr + eml_dqc + eml_dqr
3071 dtem = - di2l * lhf0 / ( cva(k) * dens(k,i,j) )
3072 if ( abs(dtem) < eps )
then
3075 fact = min( 1.0_rp, max( 0.0_rp, ( t00 - temp(k) ) / dtem ) )
3080 drhoqc = wrm_dqc + ( clp_dqc + clm_dqc + eml_dqc ) * fact
3081 drhonc = wrm_dnc + ( clp_dnc + clm_dnc + eml_dnc ) * fact
3083 drhoqr = wrm_dqr + ( clp_dqr + clm_dqr + eml_dqr ) * fact
3084 drhonr = wrm_dnr + ( clp_dnr + clm_dnr + eml_dnr ) * fact
3086 drhoqi = ( clp_dqi + clm_dqi + eml_dqi ) * fact + pco_dqi + spl_dqi
3087 drhoni = ( clp_dni + clm_dni + eml_dni ) * fact + pco_dni + spl_dni
3089 drhoqs = ( clp_dqs + clm_dqs + eml_dqs ) * fact + pco_dqs + spl_dqs
3090 drhons = ( clp_dns + clm_dns + eml_dns ) * fact + pco_dns
3092 drhoqg = ( clp_dqg + clm_dqg + eml_dqg ) * fact + pco_dqg + spl_dqg
3093 drhong = ( clp_dng + clm_dng + eml_dng ) * fact + pco_dng
3096 drhoqcrg_c = wrm_dnc_crg + ( clp_dnc_crg + clm_dnc_crg + eml_dnc_crg ) * fact
3097 drhoqcrg_r = wrm_dnr_crg + ( clp_dnr_crg + clm_dnr_crg + eml_dnr_crg ) * fact
3098 drhoqcrg_i = ( clp_dni_crg + clm_dni_crg + eml_dni_crg ) * fact + pco_dni_crg + spl_dni_crg
3099 drhoqcrg_s = ( clp_dns_crg + clm_dns_crg + eml_dns_crg ) * fact + pco_dns_crg + spl_dns_crg
3100 drhoqcrg_g = ( clp_dng_crg + clm_dng_crg + eml_dng_crg ) * fact + pco_dng_crg + spl_dng_crg
3105 rhoq0_t(k,i_qc) = drhoqc / dt
3106 rhoq0_t(k,i_nc) = drhonc / dt
3107 rhoq0_t(k,i_qr) = drhoqr / dt
3108 rhoq0_t(k,i_nr) = drhonr / dt
3109 rhoq0_t(k,i_qi) = drhoqi / dt
3110 rhoq0_t(k,i_ni) = drhoni / dt
3111 rhoq0_t(k,i_qs) = drhoqs / dt
3112 rhoq0_t(k,i_ns) = drhons / dt
3113 rhoq0_t(k,i_qg) = drhoqg / dt
3114 rhoq0_t(k,i_ng) = drhong / dt
3117 rhoqcrg0_t(k,i,j,i_qc-1) = drhoqcrg_c / dt
3118 rhoqcrg0_t(k,i,j,i_qr-1) = drhoqcrg_r / dt
3119 rhoqcrg0_t(k,i,j,i_qi-1) = drhoqcrg_i / dt
3120 rhoqcrg0_t(k,i,j,i_qs-1) = drhoqcrg_s / dt
3121 rhoqcrg0_t(k,i,j,i_qg-1) = drhoqcrg_g / dt
3124 rhoe0_t(k) = lhf * ( drhoqi + drhoqs + drhoqg ) / dt
3126 dql = rrho(k) * ( drhoqc + drhoqr )
3127 dqi = rrho(k) * ( drhoqi + drhoqs + drhoqg )
3132 cvtot0_t(k) = dcv / dt
3133 cptot0_t(k) = dcp / dt
3139 rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe0_t(k)
3140 cvtot_t(k,i,j) = cvtot_t(k,i,j) + cvtot0_t(k)
3141 cptot_t(k,i,j) = cptot_t(k,i,j) + cptot0_t(k)
3147 rhoq(k,iq) = max(0.0_rp, rhoq(k,iq) + rhoq0_t(k,iq)*dt )
3155 rhoq_crg(iq,k) = rhoq_crg(iq,k) + rhoqcrg0_t(k,i,j,iq-1)*dt
3158 qtrc0(k,i,j,iq) = rhoq(k,iq) / dens(k,i,j)
3168 sarea(k,i,j,i_mp_qc) = crs(k,i,j,i_mp_qc)
3169 sarea(k,i,j,i_mp_qr) = crs(k,i,j,i_mp_qr)
3170 sarea(k,i,j,i_mp_qi) = crs(k,i,j,i_mp_qi)
3171 sarea(k,i,j,i_mp_qs) = crs(k,i,j,i_mp_qs)
3172 sarea(k,i,j,i_mp_qg) = crs(k,i,j,i_mp_qg)
3179 rhoq_t(k,i,j,iq) = ( rhoq(k,iq) - dens(k,i,j)*qtrc(k,i,j,iq) )/dt
3186 rhoqcrg_t_mp(k,i,j,iq-1) = ( rhoq_crg(iq,k) - dens(k,i,j)*qtrc_crg(k,i,j,iq-1) )/dt
3193 4, i, j, temp(:), dens(:,i,j), pres(:), qtrc(:,i,j,i_qv) )
3200 end subroutine mp_sn14
3211 integer,
intent(in) :: KA, KS, KE
3213 integer,
intent(in) :: point
3214 integer,
intent(in) :: i, j
3215 real(RP),
intent(in) :: tem(KA)
3216 real(RP),
intent(in) :: rho(KA)
3217 real(RP),
intent(in) :: pre(KA)
3218 real(RP),
intent(in) :: qv (KA)
3224 if ( tem(k) < tem_min &
3225 .OR. rho(k) < rho_min &
3226 .OR. pre(k) < 1.0_rp )
then
3228 log_info(
"ATMOS_PHY_MP_SN14_debug_tem_kij",
'(A,I3,A,4(F16.5),3(I6))') &
3229 "point: ", point,
" low tem,rho,pre:", tem(k), rho(k), pre(k), qv(k), k, i, j,
prc_myrank
3240 rho, tem, pre, qdry, &
3245 CCN, nc_uplim_d, & ! in
3251 moist_psat_liq => atmos_saturation_psat_liq, &
3252 moist_psat_ice => atmos_saturation_psat_ice, &
3253 moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
3254 moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
3255 moist_dqsi_dtem_dens => atmos_saturation_dqs_dtem_dens_liq
3258 integer,
intent(in) :: KA, KS, KE
3260 real(RP),
intent(in) :: cz( KA)
3261 real(RP),
intent(in) :: fz(0:KA)
3262 real(RP),
intent(in) :: w (KA)
3263 real(RP),
intent(in) :: rho (KA)
3264 real(RP),
intent(in) :: tem (KA)
3265 real(RP),
intent(in) :: pre (KA)
3266 real(RP),
intent(in) :: qdry(KA)
3268 real(RP),
intent(in) :: rhoq(KA,I_QV:I_NG)
3270 real(RP),
intent(in) :: cpa(KA)
3271 real(RP),
intent(in) :: dTdt_rad(KA)
3272 real(RP),
intent(in) :: qke(KA)
3273 real(RP),
intent(in) :: dt
3274 real(RP),
intent(in) :: CCN(KA)
3275 real(RP),
intent(in) :: nc_uplim_d
3277 real(RP),
intent(out) :: PQ(KA,PQ_MAX)
3288 real(RP) :: w_dssidz(KA)
3290 real(RP) :: ssi_below(KA)
3291 real(RP) :: z_below(KA)
3297 real(RP) :: dqsidtem_rho(KA)
3298 real(RP) :: dssidt_rad(KA)
3299 real(RP) :: wssi, wdssi
3305 real(RP) :: sigma_w(KA)
3306 real(RP) :: weff(KA)
3307 real(RP) :: weff_max(KA)
3308 real(RP) :: velz(KA)
3310 real(RP) :: coef_ccn
3311 real(RP) :: slope_ccn
3312 real(RP) :: nc_new(KA)
3313 real(RP) :: nc_new_below(KA)
3315 real(RP) :: nc_new_max
3318 logical :: flag_nucleation(KA)
3320 real(RP) :: r_gravity
3321 real(RP),
parameter :: r_sqrt3=0.577350269_rp
3322 real(RP),
parameter :: eps=1.e-30_rp
3325 real(RP) :: dlcdt_max, dli_max
3326 real(RP) :: dncdt_max, dni_max
3340 r_gravity = 1.0_rp/grav
3342 call moist_psat_liq ( ka, ks, ke, &
3344 call moist_psat_ice ( ka, ks, ke, &
3346 call moist_pres2qsat_liq ( ka, ks, ke, &
3347 tem(:), pre(:), qdry(:), &
3349 call moist_pres2qsat_ice ( ka, ks, ke, &
3350 tem(:), pre(:), qdry(:), &
3352 call moist_dqsi_dtem_dens( ka, ks, ke, &
3357 a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
3363 pv = rhoq(k,i_qv)*rvap*tem(k)
3364 ssw(k) = min( mp_ssw_lim, ( pv/esw(k)-1.0_rp ) )*100.0_rp
3365 ssi(k) = ( pv/esi(k) - 1.00_rp )
3367 ssi_below(k+1) = ssi(k)
3368 z_below(k+1) = cz(k)
3371 ssi_below(ks) = ssi(ks)
3372 z_below(ks) = cz(ks-1)
3377 coef_ccn = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
3379 * (70.0_rp)**(kappa/(kappa + 2.0_rp))
3381 slope_ccn = 1.5_rp*kappa/(kappa + 2.0_rp)
3384 sigma_w(k) = r_sqrt3*sqrt(max(qke(k),qke_min))
3386 sigma_w(ks-1) = sigma_w(ks)
3387 sigma_w(ke+1) = sigma_w(ke)
3390 weff(k) = w(k) - cpa(k)*r_gravity*dtdt_rad(k)
3393 if( mp_couple_aerosol )
then
3396 if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp )
then
3397 nc_new(k) = max( ccn(k), c_ccn )
3405 if( nucl_twomey )
then
3409 weff_max(k) = weff(k) + sigma_w(k)
3411 if( (weff(k) > 1.e-8_rp) .AND. (ssw(k) > 1.e-10_rp) .AND. pre(k) > 300.e+2_rp )
then
3413 nc_new_max = coef_ccn*weff_max(k)**slope_ccn
3414 nc_new(k) = a_max*nc_new_max**b_max
3423 if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp )
then
3424 nc_new(k) = c_ccn*ssw(k)**kappa
3435 if( nc_new(k) > nc_uplim_d )
then
3436 flag_nucleation(k) = .false.
3437 nc_new_below(k+1) = 1.e+30_rp
3438 else if( nc_new(k) > eps )
then
3439 flag_nucleation(k) = .true.
3440 nc_new_below(k+1) = nc_new(k)
3442 flag_nucleation(k) = .false.
3443 nc_new_below(k+1) = 0.0_rp
3446 nc_new_below(ks) = 0.0_rp
3456 if( mp_couple_aerosol )
then
3461 if ( flag_nucleation(k) .AND. &
3462 tem(k) > tem_ccn_low )
then
3463 dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3464 dlcdt_max = max( dlcdt_max, 0.0_rp )
3465 dncdt_max = dlcdt_max/xc_min
3468 pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3469 pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3471 pq(k,i_ncccn) = 0.0_rp
3472 pq(k,i_lcccn) = 0.0_rp
3477 if( nucl_twomey )
then
3482 if ( flag_nucleation(k) .AND. &
3483 tem(k) > tem_ccn_low .AND. &
3484 nc_new(k) > rhoq(k,i_nc) )
then
3485 dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3486 dlcdt_max = max( dlcdt_max, 0.0_rp )
3487 dncdt_max = dlcdt_max/xc_min
3488 dnc_new = nc_new(k)-rhoq(k,i_nc)
3489 pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3490 pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3492 pq(k,i_ncccn) = 0.0_rp
3493 pq(k,i_lcccn) = 0.0_rp
3499 if( tem(k) > tem_ccn_low .AND. &
3500 nc_new(k) > rhoq(k,i_nc) )
then
3501 dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3502 dlcdt_max = max( dlcdt_max, 0.0_rp )
3503 dncdt_max = dlcdt_max/xc_min
3504 dnc_new = nc_new(k)-rhoq(k,i_nc)
3505 pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3506 pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3508 pq(k,i_ncccn) = 0.0_rp
3509 pq(k,i_lcccn) = 0.0_rp
3525 velz(k) = ( w(k) * ( cz(k+1) - fz(k) ) + w(k+1) * ( fz(k) - cz(k) ) ) / ( cz(k+1) - cz(k) )
3529 dzh = cz(k) - z_below(k)
3530 w_dssidz(k) = velz(k) * (ssi(k) - ssi_below(k))/dzh
3531 dssidt_rad(k) = -rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))*dqsidtem_rho(k)*dtdt_rad(k)
3532 dli_max = ( rhoq(k,i_qv) - esi(k) / ( rvap * tem(k) ) ) * rdt
3533 dli_max = max( dli_max, 0.0_rp )
3534 dni_max = min( dli_max/xi_ccn, (in_max-rhoq(k,i_ni))*rdt )
3535 wdssi = min( w_dssidz(k)+dssidt_rad(k), 0.01_rp)
3536 wssi = min( ssi(k), ssi_max)
3538 if( ( wdssi > eps ) .AND. &
3539 (tem(k) < 273.15_rp ) .AND. &
3540 (rhoq(k,i_ni) < in_max ) .AND. &
3541 (wssi >= eps ) )
then
3542 tmp = c_in * nm_m92 * exp( 0.3_rp * bm_m92 * ( wssi - 0.1_rp ) )
3544 tmp = bm_m92 * 0.3_rp * tmp * wdssi
3546 tmp = max( tmp - rhoq(k,i_ni), 0.0_rp ) * rdt
3548 pq(k,i_niccn) = min(dni_max, tmp)
3549 pq(k,i_liccn) = min(dli_max, pq(k,i_niccn)*xi_ccn )
3551 pq(k,i_niccn) = 0.0_rp
3552 pq(k,i_liccn) = 0.0_rp
3565 rhoq_crg, xq, & ! in
3574 integer,
intent(in) :: KA, KS, KE
3576 real(RP),
intent(in) :: Pac(KA,Pac_MAX)
3577 real(RP),
intent(in) :: tem(KA)
3578 real(RP),
intent(in) :: rhoq(KA,I_QV:I_NG)
3579 real(RP),
intent(in) :: xq(KA,HYDRO_MAX)
3581 real(RP),
intent(inout):: PQ(KA,PQ_MAX)
3583 logical,
intent(in) :: flg_lt
3584 real(RP),
intent(in) :: rhoq_crg(I_QC:I_QG,KA)
3585 real(RP),
intent(inout):: Pcrg1(PQ_MAX,KA)
3588 real(RP),
parameter :: pice = 350.0e+6_rp
3590 real(RP),
parameter :: pnc = 250.0_rp
3597 real(RP) :: a0,a1,a2,a3,a4,a5
3598 real(RP) :: a6,a7,a8,a9,a10
3599 real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
3600 real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
3601 real(RP),
parameter :: eps=1.0e-30_rp
3605 real(RP) :: wn, wni, wns, wng
3614 if (tem(k) > 270.16_rp)
then
3616 else if(tem(k) >= 268.16_rp)
then
3617 fp = (270.16_rp-tem(k))*0.5_rp
3618 else if(tem(k) >= 265.16_rp)
then
3619 fp = (tem(k)-265.16_rp)*0.333333333_rp
3627 x = coef_lambda(i_mp_qc)*(xc_cr/xq(k,i_mp_qc))**mu(i_mp_qc)
3629 if(x<1.e-2_rp*alpha)
then
3631 else if(x<alpha+1.0_rp)
then
3634 a1 = a0*x/(alpha+1.0_rp)
3635 a2 = a1*x/(alpha+2.0_rp)
3636 a3 = a2*x/(alpha+3.0_rp)
3637 a4 = a3*x/(alpha+4.0_rp)
3638 a5 = a4*x/(alpha+5.0_rp)
3639 a6 = a5*x/(alpha+6.0_rp)
3640 a7 = a6*x/(alpha+7.0_rp)
3641 a8 = a7*x/(alpha+8.0_rp)
3642 a9 = a8*x/(alpha+9.0_rp)
3643 a10 = a9*x/(alpha+10.0_rp)
3644 igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
3645 else if(x<alpha*1.d2)
then
3653 an1 = -(1.0_rp-alpha)
3655 d1 = 1.0_rp/(an1*d0+b1)
3660 an2 = -2.0_rp*(2.0_rp-alpha)
3662 d2 = 1.0_rp/(an2*d1+b2)
3667 igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
3672 n12 = rhoq(k,i_nc)*(1.0_rp-igm)
3674 wn = (pice + n12/((rhoq(k,i_qc)+xc_min)*pnc) )*fp
3675 wni = wn*(-pac(k,i_liaclc2li) )
3676 wns = wn*(-pac(k,i_lsaclc2ls) )
3677 wng = wn*(-pac(k,i_lgaclc2lg) )
3678 pq(k,i_nispl) = wni+wns+wng
3680 pq(k,i_lsspl) = - wns*xq(k,i_mp_qi)
3681 pq(k,i_lgspl) = - wng*xq(k,i_mp_qi)
3683 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
3684 pcrg1(i_nsspl,k) = - wns*(1.0_rp-sw1) &
3685 / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
3686 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
3687 pcrg1(i_ngspl,k) = - wng*(1.0_rp-sw1) &
3688 / (rhoq(k,i_ng)+sw1)*rhoq_crg(i_qg,k)
3689 pcrg1(i_nispl,k) = - ( pcrg1(i_nsspl,k) + pcrg1(i_ngspl,k) )
3699 ! collection process
3702 d0_crg, v0_crg, & ! in
3703 beta_crg, dqcrg, & ! in
3704 wtem, rhoq, rhoq_crg, & ! in
3705 xq, dq_xave, vt_xave, & ! in
3706 ! rho ! [Add] 11/08/30
3708 Pcrg1, Pcrg2, & ! inout
3711 moist_psat_ice => atmos_saturation_psat_ice
3714 integer,
intent(in) :: KA, KS, KE
3719 real(RP),
intent(in) :: wtem(KA)
3721 real(RP),
intent(in) :: rhoq(KA,I_QV:I_NG)
3723 real(RP),
intent(in) :: xq(KA,HYDRO_MAX)
3725 real(RP),
intent(in) :: dq_xave(KA,HYDRO_MAX)
3727 real(RP),
intent(in) :: vt_xave(KA,HYDRO_MAX,2)
3731 real(RP),
intent(inout):: PQ(KA,PQ_MAX)
3733 real(RP),
intent(out):: Pac(KA,Pac_MAX)
3735 logical,
intent(in) :: flg_lt
3736 real(RP),
intent(in) :: beta_crg(KA)
3737 real(RP),
intent(in) :: dqcrg(KA)
3738 real(RP),
intent(in) :: d0_crg, v0_crg
3739 real(RP),
intent(in) :: rhoq_crg(I_QC:I_QG,KA)
3740 real(RP),
intent(inout):: Pcrg1(PQ_MAX,KA)
3741 real(RP),
intent(inout):: Pcrg2(Pcrg_MAX,KA)
3743 real(RP),
parameter :: a_dec = 0.883_rp
3744 real(RP),
parameter :: b_dec = 0.093_rp
3745 real(RP),
parameter :: c_dec = 0.00348_rp
3746 real(RP),
parameter :: d_dec = 4.5185e-5_rp
3752 real(RP) :: E_c, E_r, E_i, E_s, E_g
3753 real(RP) :: E_ic, E_sc, E_gc
3755 real(RP) :: E_stick(KA)
3757 real(RP) :: temc, temc2, temc3
3762 real(RP) :: temc_p, temc_m
3773 real(RP) :: coef_acc_LCI, coef_acc_NCI
3774 real(RP) :: coef_acc_LCS, coef_acc_NCS
3776 real(RP) :: coef_acc_LCG, coef_acc_NCG
3777 real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I
3778 real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R
3779 real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S
3780 real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R
3781 real(RP) :: coef_acc_LRG, coef_acc_NRG
3782 real(RP) :: coef_acc_LII, coef_acc_NII
3783 real(RP) :: coef_acc_LIS, coef_acc_NIS
3784 real(RP) :: coef_acc_NSS
3785 real(RP) :: coef_acc_NGG
3786 real(RP) :: coef_acc_LSG, coef_acc_NSG
3788 real(RP) :: dcdc, dcdi, dcds, dcdg
3789 real(RP) :: drdr, drdi, drds, drdg
3790 real(RP) :: didi, dids, didg
3791 real(RP) :: dsds, dsdg
3794 real(RP) :: vcvc, vcvi, vcvs, vcvg
3795 real(RP) :: vrvr, vrvi, vrvs, vrvg
3796 real(RP) :: vivi, vivs, vivg
3797 real(RP) :: vsvs, vsvg
3800 real(RP) :: wx_cri, wx_crs
3801 real(RP) :: coef_emelt
3804 real(RP) :: sw, sw1, sw2
3805 real(RP) :: alpha_lt
3811 tem(k) = max( wtem(k), tem_min )
3814 call moist_psat_ice( ka, ks, ke, &
3817 if( opt_stick_ks96 )
then
3823 e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
3824 esi_rat = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
3825 e_stick(k) = min(1.0_rp, e_dec*esi_rat)
3827 else if( opt_stick_co86 )
then
3830 temc = min(tem(k) - t00,0.0_rp)
3831 w1 = 0.035_rp*temc-0.7_rp
3832 e_stick(k) = 10._rp**w1
3837 temc_m = min(tem(k) - t00,0.0_rp)
3838 e_stick(k) = exp(0.09_rp*temc_m)
3846 temc_p = max(tem(k) - t00,0.0_rp)
3848 ave_dc = coef_d(i_mp_qc)*xq(k,i_mp_qc)**b_m(i_mp_qc)
3849 ave_di = coef_d(i_mp_qi)*xq(k,i_mp_qi)**b_m(i_mp_qi)
3850 ave_ds = coef_d(i_mp_qs)*xq(k,i_mp_qs)**b_m(i_mp_qs)
3851 ave_dg = coef_d(i_mp_qg)*xq(k,i_mp_qg)**b_m(i_mp_qg)
3854 e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
3855 sw = 0.5_rp - sign(0.5_rp, di0-ave_di)
3857 sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds)
3859 sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg)
3866 dcdc = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qc)
3867 drdr = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qr)
3868 didi = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qi)
3869 dsds = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qs)
3870 dgdg = dq_xave(k,i_mp_qg) * dq_xave(k,i_mp_qg)
3871 dcdi = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qi)
3872 dcds = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qs)
3873 dcdg = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qg)
3874 drdi = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qi)
3875 drds = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qs)
3876 drdg = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qg)
3877 dids = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qs)
3879 dsdg = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qg)
3881 vcvc = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qc,2)
3882 vrvr = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qr,2)
3883 vivi = vt_xave(k,i_mp_qi,2)* vt_xave(k,i_mp_qi,2)
3884 vsvs = vt_xave(k,i_mp_qs,2)* vt_xave(k,i_mp_qs,2)
3885 vgvg = vt_xave(k,i_mp_qg,2)* vt_xave(k,i_mp_qg,2)
3886 vcvi = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qi,2)
3887 vcvs = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qs,2)
3888 vcvg = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qg,2)
3889 vrvi = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qi,2)
3890 vrvs = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qs,2)
3891 vrvg = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qg,2)
3892 vivs = vt_xave(k,i_mp_qi,2)* vt_xave(k,i_mp_qs,2)
3894 vsvg = vt_xave(k,i_mp_qs,2)* vt_xave(k,i_mp_qg,2)
3903 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3904 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3905 + sigma_i + sigma_c )**0.5_rp
3907 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3908 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3909 + sigma_i + sigma_c )**0.5_rp
3910 pac(k,i_liaclc2li)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_qc)*coef_acc_lci
3911 pac(k,i_niacnc2ni)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_nc)*coef_acc_nci
3914 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
3915 pcrg2(i_niacnc2ni,k) = pac(k,i_niacnc2ni)*(1.0_rp-sw1) &
3916 / (rhoq(k,i_nc)+sw1)*rhoq_crg(i_qc,k)
3921 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3922 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3923 + sigma_s + sigma_c )**0.5_rp
3925 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3926 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3927 + sigma_s + sigma_c )**0.5_rp
3928 pac(k,i_lsaclc2ls)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_qc)*coef_acc_lcs
3929 pac(k,i_nsacnc2ns)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_nc)*coef_acc_ncs
3932 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
3933 pcrg2(i_nsacnc2ns,k) = pac(k,i_nsacnc2ns)*(1.0_rp-sw1) &
3934 / (rhoq(k,i_nc)+sw1)*rhoq_crg(i_qc,k)
3939 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3940 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3941 + sigma_g + sigma_c )**0.5_rp
3943 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3944 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3945 + sigma_g + sigma_c )**0.5_rp
3946 pac(k,i_lgaclc2lg)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_qc)*coef_acc_lcg
3947 pac(k,i_ngacnc2ng)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_nc)*coef_acc_ncg
3950 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
3951 pcrg2(i_ngacnc2ng,k) = pac(k,i_ngacnc2ng)*(1.0_rp-sw1) &
3952 / (rhoq(k,i_nc)+sw1)*rhoq_crg(i_qc,k)
3956 ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3957 * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3958 + sigma_g + sigma_s )**0.5_rp
3960 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3962 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3963 + sigma_g + sigma_s )**0.5_rp
3964 pac(k,i_lgacls2lg)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_qs)*coef_acc_lsg
3965 pac(k,i_ngacns2ng)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg
3968 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
3969 pcrg2(i_ngacns2ng,k) = pac(k,i_ngacns2ng)*(1.0_rp-sw1) &
3970 / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
3972 alpha_lt = 5.0_rp * ( dq_xave(k,i_mp_qs) / d0_crg )**2 * vt_xave(k,i_mp_qg,2) / v0_crg
3973 alpha_lt = min( alpha_lt, 10.0_rp )
3974 pcrg2(i_cgngacns2ng,k)= 0.25_rp*pi*( 1.0_rp - e_stick(k) )*e_gs &
3975 * rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg &
3976 * ( dqcrg(k)*alpha_lt ) &
4012 ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
4013 * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
4014 + sigma_i + sigma_s )**0.5_rp
4016 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
4017 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
4018 + sigma_i + sigma_s )**0.5_rp
4019 pac(k,i_liacls2ls)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_qi)*coef_acc_lis
4020 pac(k,i_niacns2ns)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_ni)*coef_acc_nis
4024 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
4025 pcrg2(i_niacns2ns,k) = pac(k,i_niacns2ns)*(1.0_rp-sw1) &
4026 / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4028 sw = sign(0.5_rp, t00-tem(k)) + 0.5_rp
4038 ( ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) * sw &
4039 + ( delta_b1(i_mp_qg)*dgdg + delta_ab1(i_mp_qr,i_mp_qg)*drdg + delta_b0(i_mp_qr)*drdr ) * (1.0_rp-sw) ) &
4040 * sqrt( ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg ) * sw &
4041 + ( theta_b1(i_mp_qg)*vgvg - theta_ab1(i_mp_qr,i_mp_qg)*vrvg + theta_b0(i_mp_qr)*vrvr ) * (1.0_rp-sw) &
4042 + sigma_r + sigma_g )
4043 pac(k,i_lraclg2lg) = -0.25_rp*pi*e_gr*coef_acc_lrg &
4044 * ( rhoq(k,i_ng)*rhoq(k,i_qr) * sw &
4045 + rhoq(k,i_nr)*rhoq(k,i_qg) * (1.0_rp-sw) )
4047 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) &
4048 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg &
4049 + sigma_r + sigma_g )**0.5_rp
4050 pac(k,i_nracng2ng) = -0.25_rp*pi*e_gr*rhoq(k,i_ng)*rhoq(k,i_nr)*coef_acc_nrg
4053 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
4054 sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
4055 pcrg2(i_nracng2ng,k) = pac(k,i_nracng2ng)*(1.0_rp-sw1)/(rhoq(k,i_nr)+sw1) &
4056 *rhoq_crg(i_qr,k) * sw &
4057 + pac(k,i_nracng2ng)*(1.0_rp-sw2)/(rhoq(k,i_ng)+sw2) &
4058 *rhoq_crg(i_qg,k) * (1.0_rp-sw)
4069 ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
4070 * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
4071 + sigma_r + sigma_i )**0.5_rp
4073 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
4074 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
4075 + sigma_r + sigma_i )**0.5_rp
4076 pac(k,i_lracli2lg_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_qi)*coef_acc_lri_i
4077 pac(k,i_nracni2ng_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_ni)*coef_acc_nri_i
4080 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
4081 pcrg2(i_nracni2ng_i,k) = pac(k,i_nracni2ng_i)*(1.0_rp-sw1) &
4082 / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4086 ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
4087 * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
4088 + sigma_r + sigma_i )**0.5_rp
4090 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
4091 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
4092 + sigma_r + sigma_i )**0.5_rp
4093 pac(k,i_lracli2lg_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_qr)*coef_acc_lri_r
4094 pac(k,i_nracni2ng_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_nr)*coef_acc_nri_r
4097 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
4098 pcrg2(i_nracni2ng_r,k) = pac(k,i_nracni2ng_r)*(1.0_rp-sw1)&
4099 / (rhoq(k,i_nr)+sw1)*rhoq_crg(i_qr,k)
4104 ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
4105 * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
4106 + sigma_r + sigma_s )**0.5_rp
4108 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
4109 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
4110 + sigma_r + sigma_s )**0.5_rp
4111 pac(k,i_lracls2lg_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_qs)*coef_acc_lrs_s
4112 pac(k,i_nracns2ng_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_ns)*coef_acc_nrs_s
4115 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
4116 pcrg2(i_nracns2ng_s,k) = pac(k,i_nracns2ng_s)*(1.0_rp-sw1)&
4117 / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4121 ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
4122 * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
4123 + sigma_r + sigma_s )**0.5_rp
4125 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
4126 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
4127 + sigma_r + sigma_s )**0.5_rp
4128 pac(k,i_lracls2lg_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_qr)*coef_acc_lrs_r
4129 pac(k,i_nracns2ng_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_nr)*coef_acc_nrs_r
4132 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
4133 pcrg2(i_nracns2ng_r,k) = pac(k,i_nracns2ng_r)*(1.0_rp-sw1) &
4134 / (rhoq(k,i_nr)+sw1)*rhoq_crg(i_qr,k)
4143 ( delta_b0(i_mp_qi)*didi + delta_ab1(i_mp_qi,i_mp_qi)*didi + delta_b1(i_mp_qi)*didi ) &
4144 * ( theta_b0(i_mp_qi)*vivi - theta_ab1(i_mp_qi,i_mp_qi)*vivi + theta_b1(i_mp_qi)*vivi &
4145 + sigma_i + sigma_i )**0.5_rp
4147 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qi,i_mp_qi)*didi + delta_b0(i_mp_qi)*didi ) &
4148 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qi,i_mp_qi)*vivi + theta_b0(i_mp_qi)*vivi &
4149 + sigma_i + sigma_i )**0.5_rp
4150 pac(k,i_liacli2ls)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_qi)*coef_acc_lii
4151 pac(k,i_niacni2ns)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_ni)*coef_acc_nii
4154 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
4155 pcrg2(i_niacni2ns,k) = pac(k,i_niacni2ns)*(1.0_rp-sw1) &
4156 / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4169 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qs,i_mp_qs)*dsds + delta_b0(i_mp_qs)*dsds ) &
4170 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qs,i_mp_qs)*vsvs + theta_b0(i_mp_qs)*vsvs &
4171 + sigma_s + sigma_s )**0.5_rp
4172 pac(k,i_nsacns2ns)= -0.125_rp*pi*e_stick(k)*e_ss*rhoq(k,i_ns)*rhoq(k,i_ns)*coef_acc_nss
4176 ( delta_b0(i_mp_qg)*dgdg + delta_ab0(i_mp_qg,i_mp_qg)*dgdg + delta_b0(i_mp_qg)*dgdg ) &
4177 * ( theta_b0(i_mp_qg)*vgvg - theta_ab0(i_mp_qg,i_mp_qg)*vgvg + theta_b0(i_mp_qg)*vgvg &
4178 + sigma_g + sigma_g )**0.5_rp
4179 pac(k,i_ngacng2ng)= -0.125_rp*pi*e_stick(k)*e_gg*rhoq(k,i_ng)*rhoq(k,i_ng)*coef_acc_ngg
4186 sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di)
4187 wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di*ave_di*ave_di/xq(k,i_mp_qi) - 1.0_rp ) * sw
4188 pq(k,i_licon) = i_iconv2g * pac(k,i_liaclc2li)/max(1.0_rp, wx_cri) * sw
4189 pq(k,i_nicon) = i_iconv2g * pq(k,i_licon)/xq(k,i_mp_qi) * sw
4191 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
4192 pcrg1(i_nicon,k) = i_iconv2g * pq(k,i_nicon)*(1.0_rp-sw) &
4193 / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4197 wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds*ave_ds*ave_ds/xq(k,i_mp_qs) - 1.0_rp )
4198 pq(k,i_lscon) = i_sconv2g * (pac(k,i_lsaclc2ls))/max(1.0_rp, wx_crs)
4199 pq(k,i_nscon) = i_sconv2g * pq(k,i_lscon)/xq(k,i_mp_qs)
4201 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
4202 pcrg1(i_nscon,k) = i_sconv2g * pq(k,i_nscon)*(1.0_rp-sw) &
4203 / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4213 coef_emelt = cl/lhf0*temc_p
4215 pq(k,i_lgacm) = coef_emelt*pac(k,i_lgaclc2lg)
4216 pq(k,i_ngacm) = pq(k,i_lgacm)/xq(k,i_mp_qg)
4218 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
4219 pcrg1(i_ngacm,k) = pq(k,i_ngacm)*(1.0_rp-sw1) &
4220 / (rhoq(k,i_ng)+sw1)*rhoq_crg(i_qg,k)
4223 pq(k,i_lgarm) = coef_emelt*pac(k,i_lraclg2lg)
4224 pq(k,i_ngarm) = pq(k,i_lgarm)/xq(k,i_mp_qg)
4226 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
4227 pcrg1(i_ngarm,k) = pq(k,i_ngarm)*(1.0_rp-sw1) &
4228 / (rhoq(k,i_ng)+sw1)*rhoq_crg(i_qg,k)
4231 pq(k,i_lsacm) = coef_emelt*(pac(k,i_lsaclc2ls))
4232 pq(k,i_nsacm) = pq(k,i_lsacm)/xq(k,i_mp_qs)
4234 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
4235 pcrg1(i_nsacm,k) = pq(k,i_nsacm)*(1.0_rp-sw1) &
4236 / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4239 pq(k,i_lsarm) = coef_emelt*(pac(k,i_lracls2lg_r)+pac(k,i_lracls2lg_s))
4240 pq(k,i_nsarm) = pq(k,i_lsarm)/xq(k,i_mp_qg)
4242 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
4243 pcrg1(i_nsarm,k) = pq(k,i_nsarm)*(1.0_rp-sw1) &
4244 / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4247 pq(k,i_liacm) = coef_emelt*pac(k,i_liaclc2li)
4248 pq(k,i_niacm) = pq(k,i_liacm)/xq(k,i_mp_qi)
4250 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
4251 pcrg1(i_niacm,k) = pq(k,i_niacm)*(1.0_rp-sw1) &
4252 / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4255 pq(k,i_liarm) = coef_emelt*(pac(k,i_lracli2lg_r)+pac(k,i_lracli2lg_i))
4256 pq(k,i_niarm) = pq(k,i_liarm)/xq(k,i_mp_qg)
4258 sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
4259 pcrg1(i_niarm,k) = pq(k,i_niarm)*(1.0_rp-sw1) &
4260 / (rhoq(k,i_nr)+sw1)*rhoq_crg(i_qi,k)
4278 integer,
intent(in) :: KA, KS, KE
4280 real(RP),
intent(in) :: rhoq(KA,I_QV:I_NG)
4281 real(RP),
intent(in) :: rhoq_crg(I_QC:I_QG,KA)
4282 logical,
intent(in) :: flg_lt
4283 real(RP),
intent(in) :: xq(KA,HYDRO_MAX)
4284 real(RP),
intent(in) :: dq_xave(KA,HYDRO_MAX)
4285 real(RP),
intent(in) :: rho(KA)
4287 real(RP),
intent(inout) :: PQ(KA,PQ_MAX)
4288 real(RP),
intent(inout) :: Pcrg(PQ_MAX,KA)
4291 real(RP),
parameter :: kcc = 4.44e+9_rp
4292 real(RP),
parameter :: tau_min = 1.e-20_rp
4293 real(RP),
parameter :: rx_sep = 1.0_rp/x_sep
4296 real(RP),
parameter :: kcr = 5.8_rp
4297 real(RP),
parameter :: thr_acc = 5.e-5_rp
4300 real(RP),
parameter :: krr = 4.33_rp
4301 real(RP),
parameter :: kaprr = 60.7_rp
4302 real(RP),
parameter :: kbr = 1000._rp
4303 real(RP),
parameter :: kapbr = 2.3e+3_rp
4304 real(RP),
parameter :: dr_min = 0.35e-3_rp
4307 real(RP) :: coef_nuc0
4308 real(RP) :: coef_nuc1
4309 real(RP) :: coef_aut0
4310 real(RP) :: coef_aut1
4322 coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
4323 coef_nuc1 = (nu(i_mp_qc)+2.0_rp)*(nu(i_mp_qc)+4.0_rp)/(nu(i_mp_qc)+1.0_rp)/(nu(i_mp_qc)+1.0_rp)
4324 coef_aut0 = -kcc*coef_nuc0
4325 coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
4328 lwc = rhoq(k,i_qr) + rhoq(k,i_qc)
4329 if( lwc > xc_min )
then
4330 tau = max(tau_min, rhoq(k,i_qr)/lwc)
4334 rho_fac = sqrt(rho_0/max(rho(k),rho_min))
4337 psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3
4338 pq(k,i_ncaut) = coef_aut0*rhoq(k,i_qc)*rhoq(k,i_qc)*rho_fac*rho_fac
4340 pq(k,i_lcaut) = coef_aut1*lwc*lwc*xq(k,i_mp_qc)*xq(k,i_mp_qc) &
4341 *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac
4342 pq(k,i_nraut) = -rx_sep*pq(k,i_lcaut)
4345 sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
4346 pcrg(i_ncaut,k) = pq(k,i_ncaut)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq_crg(i_qc,k)
4347 pcrg(i_nraut,k) = -pcrg(i_ncaut,k)
4351 psi_acc =(tau/(tau+thr_acc))**4
4352 pq(k,i_lcacc) = -kcr*rhoq(k,i_qc)*rhoq(k,i_qr)*rho_fac*psi_acc
4353 pq(k,i_ncacc) = -kcr*rhoq(k,i_nc)*rhoq(k,i_qr)*rho_fac*psi_acc
4356 sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
4357 pcrg(i_ncacc,k) = pq(k,i_ncacc)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq(k,i_qc)
4361 pq(k,i_nrslc) = -krr*rhoq(k,i_nr)*rhoq(k,i_qr)*rho_fac
4364 ddr = min(1.e-3_rp, dq_xave(k,i_mp_qr) - dr_eq )
4365 if ( dq_xave(k,i_mp_qr) < dr_min )
then
4367 else if ( dq_xave(k,i_mp_qr) <= dr_eq )
then
4370 psi_brk = exp(kapbr*ddr) - 1.0_rp
4372 pq(k,i_nrbrk) = - (psi_brk + 1.0_rp)*pq(k,i_nrslc)
4380 subroutine dep_vapor_melt_ice( &
4382 rho, tem, pre, qd, & ! in
4385 xq, vt_xave, dq_xave, & ! in
4391 integer,
intent(in) :: KA, KS, KE
4394 real(RP),
intent(inout) :: PQ(KA,PQ_MAX)
4396 real(RP),
intent(in) :: rho(KA)
4397 real(RP),
intent(in) :: tem(KA)
4398 real(RP),
intent(in) :: pre(KA)
4399 real(RP),
intent(in) :: qd (KA)
4400 real(RP),
intent(in) :: esw(KA)
4401 real(RP),
intent(in) :: esi(KA)
4402 real(RP),
intent(in) :: rhoq(KA,I_QV:I_NG)
4403 real(RP),
intent(in) :: xq(KA,HYDRO_MAX)
4407 real(RP),
intent(in) :: vt_xave(KA,HYDRO_MAX,2)
4409 real(RP),
intent(in) :: dq_xave(KA,HYDRO_MAX)
4412 real(RP) :: temc_lim
4419 real(RP) :: nua, r_nua
4425 real(RP) :: Gwr, Gii, Gis, Gig
4430 real(RP) :: Nrers_r2, Nreis_r2
4431 real(RP) :: Nress_r2, Nregs_r2
4433 real(RP) :: Nrerl_r2, Nreil_r2
4434 real(RP) :: Nresl_r2, Nregl_r2
4435 real(RP) :: NscNrer_s, NscNrer_l
4436 real(RP) :: NscNrei_s, NscNrei_l
4437 real(RP) :: NscNres_s, NscNres_l
4438 real(RP) :: NscNreg_s, NscNreg_l
4439 real(RP) :: ventLR_s, ventLR_l
4440 real(RP) :: ventNI_s, ventNI_l, ventLI_s, ventLI_l
4441 real(RP) :: ventNS_s, ventNS_l, ventLS_s, ventLS_l
4442 real(RP) :: ventNG_s, ventNG_l, ventLG_s, ventLG_l
4444 real(RP) :: wtr, wti, wts, wtg
4445 real(RP),
parameter :: r_14=1.0_rp/1.4_rp
4446 real(RP),
parameter :: r_15=1.0_rp/1.5_rp
4449 real(RP) :: ventNI, ventLI
4450 real(RP) :: ventNS, ventLS
4451 real(RP) :: ventNG, ventLG
4453 real(RP),
parameter :: Re_max=1.e+3_rp
4454 real(RP),
parameter :: Re_min=1.e-4_rp
4468 temc_lim= max(temc, -40._rp )
4469 rho_lim = max(rho(k),rho_min)
4470 qv = rhoq(k,i_qv)/rho_lim
4471 pre_lim = rho_lim*(qd(k)*rdry + qv*rvap)*(temc_lim+t00)
4479 dw = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(p00/pre_lim)
4480 kalfa = ka0 + temc_lim*dka_dt
4481 mua = mua0 + temc_lim*dmua_dt
4484 gw = (lhv0/kalfa/tem(k))*(lhv0/rvap/tem(k)-1.0_rp)+(rvap*tem(k)/dw/esw(k))
4485 gi = (lhs0/kalfa/tem(k))*(lhs0/rvap/tem(k)-1.0_rp)+(rvap*tem(k)/dw/esi(k))
4487 gwr = 4.0_rp*pi/cap(i_mp_qr)/gw
4488 gii = 4.0_rp*pi/cap(i_mp_qi)/gi
4489 gis = 4.0_rp*pi/cap(i_mp_qs)/gi
4490 gig = 4.0_rp*pi/cap(i_mp_qg)/gi
4493 nsc_r3 = (nua/dw)**(0.33333333_rp)
4496 nrers_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qr,1)*dq_xave(k,i_mp_qr)*r_nua)))
4497 nreis_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qi,1)*dq_xave(k,i_mp_qi)*r_nua)))
4498 nress_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qs,1)*dq_xave(k,i_mp_qs)*r_nua)))
4499 nregs_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qg,1)*dq_xave(k,i_mp_qg)*r_nua)))
4502 nrerl_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qr,2)*dq_xave(k,i_mp_qr)*r_nua)))
4503 nreil_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qi,2)*dq_xave(k,i_mp_qi)*r_nua)))
4504 nresl_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qs,2)*dq_xave(k,i_mp_qs)*r_nua)))
4505 nregl_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qg,2)*dq_xave(k,i_mp_qg)*r_nua)))
4506 nscnrer_s=nsc_r3*nrers_r2
4507 nscnrer_l=nsc_r3*nrerl_r2
4509 nscnrei_s=nsc_r3*nreis_r2
4510 nscnrei_l=nsc_r3*nreil_r2
4512 nscnres_s=nsc_r3*nress_r2
4513 nscnres_l=nsc_r3*nresl_r2
4515 nscnreg_s=nsc_r3*nregs_r2
4516 nscnreg_l=nsc_r3*nregl_r2
4518 ventlr_s = ah_vent1(i_mp_qr,1) + bh_vent1(i_mp_qr,1)*nscnrer_s
4519 ventlr_l = ah_vent1(i_mp_qr,2) + bh_vent1(i_mp_qr,2)*nscnrer_l
4521 ventni_s = ah_vent0(i_mp_qi,1) + bh_vent0(i_mp_qi,1)*nscnrei_s
4522 ventni_l = ah_vent0(i_mp_qi,2) + bh_vent0(i_mp_qi,2)*nscnrei_l
4523 ventli_s = ah_vent1(i_mp_qi,1) + bh_vent1(i_mp_qi,1)*nscnrei_s
4524 ventli_l = ah_vent1(i_mp_qi,2) + bh_vent1(i_mp_qi,2)*nscnrei_l
4526 ventns_s = ah_vent0(i_mp_qs,1) + bh_vent0(i_mp_qs,1)*nscnres_s
4527 ventns_l = ah_vent0(i_mp_qs,2) + bh_vent0(i_mp_qs,2)*nscnres_l
4528 ventls_s = ah_vent1(i_mp_qs,1) + bh_vent1(i_mp_qs,1)*nscnres_s
4529 ventls_l = ah_vent1(i_mp_qs,2) + bh_vent1(i_mp_qs,2)*nscnres_l
4531 ventng_s = ah_vent0(i_mp_qg,1) + bh_vent0(i_mp_qg,1)*nscnreg_s
4532 ventng_l = ah_vent0(i_mp_qg,2) + bh_vent0(i_mp_qg,2)*nscnreg_l
4533 ventlg_s = ah_vent1(i_mp_qg,1) + bh_vent1(i_mp_qg,1)*nscnreg_s
4534 ventlg_l = ah_vent1(i_mp_qg,2) + bh_vent1(i_mp_qg,2)*nscnreg_l
4538 wtr = ( min(max( nscnrer_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
4539 wti = ( min(max( nscnrei_s , 0.5_rp), 2.0_rp) -0.5_rp )*r_15
4540 wts = ( min(max( nscnres_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
4541 wtg = ( min(max( nscnreg_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
4543 ventni = (1.0_rp-wti)*ventni_s + wti*ventni_l
4544 ventns = (1.0_rp-wts)*ventns_s + wts*ventns_l
4545 ventng = (1.0_rp-wtg)*ventng_s + wtg*ventng_l
4547 ventlr = (1.0_rp-wtr)*ventlr_s + wtr*ventlr_l
4548 ventli = (1.0_rp-wti)*ventli_s + wti*ventli_l
4549 ventls = (1.0_rp-wts)*ventls_s + wts*ventls_l
4550 ventlg = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
4569 pq(k,i_lcdep) = gwr*rhoq(k,i_nc)*dq_xave(k,i_mp_qc)*coef_deplc
4570 pq(k,i_lrdep) = gwr*rhoq(k,i_nr)*dq_xave(k,i_mp_qr)*ventlr
4571 pq(k,i_lidep) = gii*rhoq(k,i_ni)*dq_xave(k,i_mp_qi)*ventli
4572 pq(k,i_lsdep) = gis*rhoq(k,i_ns)*dq_xave(k,i_mp_qs)*ventls
4573 pq(k,i_lgdep) = gig*rhoq(k,i_ng)*dq_xave(k,i_mp_qg)*ventlg
4574 pq(k,i_nrdep) = pq(k,i_lrdep)/xq(k,i_mp_qr)
4575 pq(k,i_nidep) = 0.0_rp
4576 pq(k,i_nsdep) = pq(k,i_lsdep)/xq(k,i_mp_qs)
4577 pq(k,i_ngdep) = pq(k,i_lgdep)/xq(k,i_mp_qg)
4585 dt = kalfa/(cpvap*rho_0)
4591 gm = 2.0_rp*pi/emelt&
4592 * ( (kalfa*dt/dw)*(temc) + (dw*lhs0/rvap)*(esi(k)/tem(k)-psat0/t00) )
4597 sw = ( sign(0.5_rp,temc) + 0.5_rp ) * ( sign(0.5_rp,gm-eps) + 0.5_rp )
4601 pq(k,i_limlt) = - gm * rhoq(k,i_qi)*dq_xave(k,i_mp_qi)*ventli/xq(k,i_mp_qi) * sw
4602 pq(k,i_nimlt) = - gm * rhoq(k,i_ni)*dq_xave(k,i_mp_qi)*ventni/xq(k,i_mp_qi) * sw
4603 pq(k,i_lsmlt) = - gm * rhoq(k,i_qs)*dq_xave(k,i_mp_qs)*ventls/xq(k,i_mp_qs) * sw
4604 pq(k,i_nsmlt) = - gm * rhoq(k,i_ns)*dq_xave(k,i_mp_qs)*ventns/xq(k,i_mp_qs) * sw
4605 pq(k,i_lgmlt) = - gm * rhoq(k,i_qg)*dq_xave(k,i_mp_qg)*ventlg/xq(k,i_mp_qg) * sw
4606 pq(k,i_ngmlt) = - gm * rhoq(k,i_ng)*dq_xave(k,i_mp_qg)*ventng/xq(k,i_mp_qg) * sw
4611 end subroutine dep_vapor_melt_ice
4624 integer,
intent(in) :: KA, KS, KE
4626 real(RP),
intent(in) :: dt
4628 real(RP),
intent(in) :: tem(KA)
4630 real(RP),
intent(in) :: rhoq(KA,I_QV:I_NG)
4631 real(RP),
intent(in) :: xq(KA,HYDRO_MAX)
4633 real(RP),
intent(inout):: PQ(KA,PQ_MAX)
4635 real(RP),
parameter :: temc_min = -65.0_rp
4636 real(RP),
parameter :: a_het = 0.2_rp
4637 real(RP),
parameter :: b_het = 0.65_rp
4639 real(RP) :: coef_m2_c
4640 real(RP) :: coef_m2_r
4642 real(RP) :: temc, temc2, temc3, temc4
4644 real(RP) :: Jhom, Jhet
4652 coef_m2_c = coef_m2(i_mp_qc)
4653 coef_m2_r = coef_m2(i_mp_qr)
4657 temc = max( tem(k) - t00, temc_min )
4660 jhet = a_het*(exp( -b_het*temc) - 1.0_rp )
4663 if( temc < -65.0_rp )
then
4664 jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
4665 jhet = a_het*(exp( 65.0_rp*b_het) - 1.0_rp )
4666 else if( temc < -30.0_rp )
then
4671 - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
4672 - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
4673 else if( temc < 0.0_rp)
then
4674 jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
4686 pq(k,i_lchom) = 0.0_rp
4687 pq(k,i_nchom) = 0.0_rp
4689 #if defined(PGI) || defined(SX)
4690 tmp = min( xq(k,i_mp_qc)*(jhet+jhom)*dt, 1.e+3_rp)
4691 pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
4692 pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( - tmp ) )
4694 tmp = min( xq(k,i_mp_qr)*(jhet+jhom)*dt, 1.e+3_rp)
4695 pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
4696 pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( - tmp ) )
4698 pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*xq(k,i_mp_qc)*(jhet+jhom)*dt ) )
4699 pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( - xq(k,i_mp_qc)*(jhet+jhom)*dt ) )
4700 pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*xq(k,i_mp_qr)*(jhet+jhom)*dt ) )
4701 pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( - xq(k,i_mp_qr)*(jhet+jhom)*dt ) )
4711 subroutine update_by_phase_change( &
4721 esw, esi, rhoq2, & ! in
4726 sl_PLCdep, & ! inout
4727 sl_PLRdep, sl_PNRdep, & ! inout
4732 qc_evaporate, & ! out
4733 rhoq2_crg, & ! in:optional
4744 moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
4745 moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
4746 moist_dqs_dtem_dens_liq => atmos_saturation_dqs_dtem_dens_liq, &
4747 moist_dqs_dtem_dens_ice => atmos_saturation_dqs_dtem_dens_ice, &
4748 moist_dqs_dtem_dpre_liq => atmos_saturation_dqs_dtem_dpre_liq, &
4749 moist_dqs_dtem_dpre_ice => atmos_saturation_dqs_dtem_dpre_ice
4752 integer,
intent(in) :: KA, KS, KE
4754 integer,
intent(in) :: ntmax
4756 real(RP),
intent(in) :: dt
4757 real(RP),
intent(in) :: cz(KA)
4758 real(RP),
intent(in) :: fz(0:KA)
4759 real(RP),
intent(in) :: w (KA)
4760 real(RP),
intent(in) :: dTdt_rad(KA)
4761 real(RP),
intent(in) :: rho (KA)
4762 real(RP),
intent(in) :: qdry(KA)
4763 real(RP),
intent(in) :: esw (KA)
4764 real(RP),
intent(in) :: esi (KA)
4765 real(RP),
intent(in) :: rhoq2(KA,I_QV:I_NG)
4767 real(RP),
intent(in) :: tem(KA)
4768 real(RP),
intent(in) :: pre(KA)
4769 real(RP),
intent(in) :: cpa(KA)
4770 real(RP),
intent(in) :: cva(KA)
4773 real(RP),
intent(inout) :: PQ(KA,PQ_MAX)
4775 real(RP),
intent(inout) :: sl_PLCdep
4776 real(RP),
intent(inout) :: sl_PLRdep, sl_PNRdep
4778 real(RP),
intent(out) :: RHOQ_t(KA,QA_MP)
4779 real(RP),
intent(out) :: RHOE_t(KA)
4780 real(RP),
intent(out) :: CPtot_t(KA)
4781 real(RP),
intent(out) :: CVtot_t(KA)
4784 real(RP),
intent(out) :: qc_evaporate(KA)
4787 logical,
intent(in) :: flg_lt
4788 real(RP),
intent(in),
optional :: rhoq2_crg(I_QC:I_QG,KA)
4789 real(RP),
intent(inout),
optional :: RHOQcrg_t(KA,HYDRO_MAX)
4793 real(RP) :: wtem(KA)
4800 real(RP) :: dqswdtem_rho(KA)
4801 real(RP) :: dqsidtem_rho(KA)
4802 real(RP) :: dqswdtem_pre(KA)
4803 real(RP) :: dqsidtem_pre(KA)
4804 real(RP) :: dqswdpre_tem(KA)
4805 real(RP) :: dqsidpre_tem(KA)
4810 real(RP) :: aliqliq, asolliq
4811 real(RP) :: aliqsol, asolsol
4816 real(RP) :: taucnd, r_taucnd
4817 real(RP) :: taudep, r_taudep
4818 real(RP) :: taucnd_c(KA), r_taucnd_c
4819 real(RP) :: taucnd_r(KA), r_taucnd_r
4820 real(RP) :: taudep_i(KA), r_taudep_i
4821 real(RP) :: taudep_s(KA), r_taudep_s
4822 real(RP) :: taudep_g(KA), r_taudep_g
4825 real(RP) :: PLR2NR, PLI2NI, PLS2NS, PLG2NG
4826 real(RP) :: coef_a_cnd, coef_b_cnd
4827 real(RP) :: coef_a_dep, coef_b_dep
4849 real(RP) :: r_xc_ccn, r_xi_ccn
4852 real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
4853 real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
4855 real(RP) :: drhoqcrg_c, drhoqcrg_r
4856 real(RP) :: drhoqcrg_i, drhoqcrg_s, drhoqcrg_g
4857 real(RP) :: frz_dnc_crg
4858 real(RP) :: frz_dnr_crg
4859 real(RP) :: mlt_dni_crg
4860 real(RP) :: mlt_dns_crg
4861 real(RP) :: mlt_dng_crg
4862 real(RP) :: dep_dni_crg
4863 real(RP) :: dep_dns_crg
4864 real(RP) :: dep_dng_crg
4865 real(RP) :: dep_dnr_crg
4866 real(RP) :: dep_dnc_crg
4868 real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
4869 real(RP) :: r_rvaptem
4871 real(RP) :: lvsw, lvsi
4872 real(RP) :: dlvsw, dlvsi
4874 real(RP) :: dcnd, ddep
4875 real(RP) :: uplim_cnd
4876 real(RP) :: lowlim_cnd
4878 real(RP) :: uplim_dep
4879 real(RP) :: lowlim_dep
4880 real(RP) :: ssw, ssi
4881 real(RP) :: r_esw, r_esi
4882 real(RP) :: r_lvsw, r_lvsi
4884 real(RP) :: ssw_o, ssi_o
4891 real(RP) :: fac_cndc_wrk
4893 real(RP),
parameter :: tau100day = 1.e+7_rp
4894 real(RP),
parameter :: r_tau100day = 1.e-7_rp
4895 real(RP),
parameter :: eps=1.e-30_rp
4901 real(RP) :: dqv, dql, dqi
4902 real(RP) :: dcv, dcp
4903 real(RP) :: dqc_crg, dqr_crg, dqi_crg, dqs_crg, dqg_crg
4913 r_xc_ccn=1.0_rp/xc_ccn
4916 if( opt_fix_taucnd_c )
then
4917 fac_cndc_wrk = fac_cndc**(1.0_rp-b_m(i_mp_qc))
4919 pq(k,i_lcdep) = pq(k,i_lcdep)*fac_cndc_wrk
4921 log_info(
"ATMOS_PHY_MP_SN14_update_by_phase_change",*)
"taucnd:fac_cndc_wrk=",fac_cndc_wrk
4928 wtem(k) = max( tem(k), tem_min )
4931 call moist_pres2qsat_liq( ka, ks, ke, &
4932 wtem(:), pre(:), qdry(:), &
4934 call moist_pres2qsat_ice( ka, ks, ke, &
4935 wtem(:), pre(:), qdry(:), &
4937 call moist_dqs_dtem_dens_liq( ka, ks, ke, &
4940 call moist_dqs_dtem_dens_ice( ka, ks, ke, &
4943 call moist_dqs_dtem_dpre_liq( ka, ks, ke, &
4944 wtem(:), pre(:), qdry(:), &
4945 dqswdtem_pre(:), dqswdpre_tem(:) )
4946 call moist_dqs_dtem_dpre_ice( ka, ks, ke, &
4947 wtem(:), pre(:), qdry(:), &
4948 dqsidtem_pre(:), dqsidpre_tem(:) )
4951 if( cz(k) <= 25000.0_rp )
then
4956 if( pre(k) < esw(k)+1.e-10_rp )
then
4958 dqswdtem_rho(k) = 0.0_rp
4959 dqswdtem_pre(k) = 0.0_rp
4960 dqswdpre_tem(k) = 0.0_rp
4962 if( pre(k) < esi(k)+1.e-10_rp )
then
4964 dqsidtem_rho(k) = 0.0_rp
4965 dqsidtem_pre(k) = 0.0_rp
4966 dqsidpre_tem(k) = 0.0_rp
4969 r_rvaptem = 1.0_rp/(rvap*wtem(k))
4970 lvsw = esw(k)*r_rvaptem
4971 lvsi = esi(k)*r_rvaptem
4972 pv = rhoq2(k,i_qv)*rvap*tem(k)
4973 r_esw = 1.0_rp/esw(k)
4974 r_esi = 1.0_rp/esi(k)
4975 ssw = min( mp_ssw_lim, ( pv*r_esw-1.0_rp ) )
4976 ssi = pv*r_esi - 1.0_rp
4977 r_lvsw = 1.0_rp/lvsw
4978 r_lvsi = 1.0_rp/lvsi
4979 r_taucnd_c = pq(k,i_lcdep)*r_lvsw
4980 r_taucnd_r = pq(k,i_lrdep)*r_lvsw
4981 r_taudep_i = pq(k,i_lidep)*r_lvsi
4982 r_taudep_s = pq(k,i_lsdep)*r_lvsi
4983 r_taudep_g = pq(k,i_lgdep)*r_lvsi
4990 r_cva = 1.0_rp / cva(k)
4991 r_cpa = 1.0_rp / cpa(k)
4995 + r_cva*( lhv00 + (cvvap-cl)*tem(k) )*dqswdtem_rho(k)
4998 + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k) )*dqswdtem_rho(k)
5001 + r_cva*( lhv00 + (cvvap-cl)*tem(k) )*dqsidtem_rho(k)
5004 + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k) )*dqsidtem_rho(k)
5005 pdynliq = w2 * grav * ( r_cpa*dqswdtem_pre(k) + rho(k)*dqswdpre_tem(k) )
5006 pdynsol = w2 * grav * ( r_cpa*dqsidtem_pre(k) + rho(k)*dqsidpre_tem(k) )
5007 pradliq = -dtdt_rad(k) * dqswdtem_rho(k)
5008 pradsol = -dtdt_rad(k) * dqsidtem_rho(k)
5015 acnd = pdynliq + pradliq &
5016 - ( r_taudep_i+r_taudep_s+r_taudep_g ) * ( qsw(k) - qsi(k) )
5017 adep = pdynsol + pradsol &
5018 + ( r_taucnd_c+r_taucnd_r ) * ( qsw(k) - qsi(k) )
5020 + aliqliq*( r_taucnd_c+r_taucnd_r ) &
5021 + asolliq*( r_taudep_i+r_taudep_s+r_taudep_g )
5023 + aliqsol*( r_taucnd_c+r_taucnd_r )&
5024 + asolsol*( r_taudep_i+r_taudep_s+r_taudep_g )
5026 uplim_cnd = max( rho(k)*ssw_o*qsw(k)*r_dt, 0.0_rp )
5027 lowlim_cnd = min( rho(k)*ssw_o*qsw(k)*r_dt, 0.0_rp )
5028 if( r_taucnd < r_tau100day )
then
5030 pq(k,i_lcdep) = max(lowlim_cnd, min(uplim_cnd, pq(k,i_lcdep)*ssw_o ))
5031 pq(k,i_lrdep) = max(lowlim_cnd, min(uplim_cnd, pq(k,i_lrdep)*ssw_o ))
5032 pq(k,i_nrdep) = min(0.0_rp, pq(k,i_nrdep)*ssw_o )
5035 taucnd = 1.0_rp/r_taucnd
5037 coef_a_cnd = rho(k)*acnd*taucnd
5038 coef_b_cnd = rho(k)*taucnd*r_dt*(ssw_o*qsw(k)-acnd*taucnd) * ( exp(-dt*r_taucnd) - 1.0_rp )
5039 pq(k,i_lcdep) = coef_a_cnd*r_taucnd_c - coef_b_cnd*r_taucnd_c
5040 plr2nr = pq(k,i_nrdep)/(pq(k,i_lrdep)+1.e-30_rp)
5041 pq(k,i_lrdep) = coef_a_cnd*r_taucnd_r - coef_b_cnd*r_taucnd_r
5042 pq(k,i_nrdep) = min(0.0_rp, pq(k,i_lrdep)*plr2nr )
5045 uplim_dep = max( rho(k)*ssi_o*qsi(k)*r_dt, 0.0_rp )
5046 lowlim_dep = min( rho(k)*ssi_o*qsi(k)*r_dt, 0.0_rp )
5047 if( r_taudep < r_tau100day )
then
5049 pq(k,i_lidep) = max(lowlim_dep, min(uplim_dep, pq(k,i_lidep)*ssi_o ))
5050 pq(k,i_lsdep) = max(lowlim_dep, min(uplim_dep, pq(k,i_lsdep)*ssi_o ))
5051 pq(k,i_lgdep) = max(lowlim_dep, min(uplim_dep, pq(k,i_lgdep)*ssi_o ))
5052 pq(k,i_nidep) = min(0.0_rp, pq(k,i_nidep)*ssi_o )
5053 pq(k,i_nsdep) = min(0.0_rp, pq(k,i_nsdep)*ssi_o )
5054 pq(k,i_ngdep) = min(0.0_rp, pq(k,i_ngdep)*ssi_o )
5056 taudep = 1.0_rp/r_taudep
5058 coef_a_dep = rho(k)*adep*taudep
5059 coef_b_dep = rho(k)*taudep*r_dt*(ssi_o*qsi(k)-adep*taudep) * ( exp(-dt*r_taudep) - 1.0_rp )
5060 pli2ni = pq(k,i_nidep)/max(pq(k,i_lidep),1.e-30_rp)
5061 pls2ns = pq(k,i_nsdep)/max(pq(k,i_lsdep),1.e-30_rp)
5062 plg2ng = pq(k,i_ngdep)/max(pq(k,i_lgdep),1.e-30_rp)
5063 pq(k,i_lidep) = coef_a_dep*r_taudep_i - coef_b_dep*r_taudep_i
5064 pq(k,i_lsdep) = coef_a_dep*r_taudep_s - coef_b_dep*r_taudep_s
5065 pq(k,i_lgdep) = coef_a_dep*r_taudep_g - coef_b_dep*r_taudep_g
5066 pq(k,i_nidep) = min(0.0_rp, pq(k,i_lidep)*pli2ni )
5067 pq(k,i_nsdep) = min(0.0_rp, pq(k,i_lsdep)*pls2ns )
5068 pq(k,i_ngdep) = min(0.0_rp, pq(k,i_lgdep)*plg2ng )
5071 sw = 0.5_rp - sign(0.5_rp, pq(k,i_lcdep)+eps)
5072 pncdep = min(0.0_rp, ((rhoq2(k,i_qc)+pq(k,i_lcdep)*dt)*r_xc_ccn - rhoq2(k,i_nc))*r_dt ) * sw
5085 r_rvaptem = 1.0_rp/(rvap*wtem(k))
5086 lvsw = esw(k)*r_rvaptem
5087 dlvsw = rhoq2(k,i_qv)-lvsw
5088 dcnd = dt*(pq(k,i_lcdep)+pq(k,i_lrdep))
5090 sw = ( sign(0.5_rp,dcnd) + sign(0.5_rp,dlvsw) ) &
5091 * ( 0.5_rp + sign(0.5_rp,abs(dcnd)-eps) )
5095 fac1 = min(dlvsw*sw,dcnd*sw)*sw / (abs(sw)-1.0_rp+dcnd) &
5097 dep_dqc = max( dt*pq(k,i_lcdep)*fac1, &
5098 -rhoq2(k,i_qc) - 1e30_rp*(sw+1.0_rp) )*abs(sw)
5099 dep_dqr = max( dt*pq(k,i_lrdep)*fac1, &
5100 -rhoq2(k,i_qr) - 1e30_rp*(sw+1.0_rp) )*abs(sw)
5119 dep_dnc = max( dt*pncdep*fac1, -rhoq2(k,i_nc) )
5120 dep_dnr = max( dt*pq(k,i_nrdep)*fac1, -rhoq2(k,i_nr) )
5122 qc_evaporate(k) = - dep_dnc
5126 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_nc)-small )
5127 dep_dnc_crg = dep_dnc*( 1.0_rp-sw )/( rhoq2(k,i_nc)+sw )*rhoq2_crg(i_qc,k)
5128 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_nr)-small )
5129 dep_dnr_crg = dep_dnr*( 1.0_rp-sw )/( rhoq2(k,i_nr)+sw )*rhoq2_crg(i_qr,k)
5131 sw = min( abs(rhoq2_crg(i_qc,k)),abs(dep_dnc_crg) )
5132 dep_dnc_crg = sign( sw,dep_dnc_crg )
5133 sw = min( abs(rhoq2_crg(i_qr,k)),abs(dep_dnr_crg) )
5134 dep_dnr_crg = sign( sw,dep_dnr_crg )
5138 lvsi = esi(k)*r_rvaptem
5139 ddep = dt*(pq(k,i_lidep)+pq(k,i_lsdep)+pq(k,i_lgdep))
5140 dlvsi = rhoq2(k,i_qv)-lvsi
5142 sw = ( sign(0.5_rp,ddep) + sign(0.5_rp,dlvsi) ) &
5143 * ( 0.5_rp + sign(0.5_rp,abs(ddep)-eps) )
5147 fac2 = min(dlvsi*sw,ddep*sw)*sw / (abs(sw)-1.0_rp+ddep) &
5149 dep_dqi = max( dt*pq(k,i_lidep) &
5150 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
5151 -rhoq2(k,i_qi) - 1e30_rp*(sw+1.0_rp) )
5152 dep_dqs = max( dt*pq(k,i_lsdep) &
5153 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
5154 -rhoq2(k,i_qs) - 1e30_rp*(sw+1.0_rp) )
5155 dep_dqg = max( dt*pq(k,i_lgdep) &
5156 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
5157 -rhoq2(k,i_qg) - 1e30_rp*(sw+1.0_rp) )
5179 dep_dni = max( dt*pq(k,i_nidep)*fac2, -rhoq2(k,i_ni) )
5180 dep_dns = max( dt*pq(k,i_nsdep)*fac2, -rhoq2(k,i_ns) )
5181 dep_dng = max( dt*pq(k,i_ngdep)*fac2, -rhoq2(k,i_ng) )
5184 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_ni)-small )
5185 dep_dni_crg = dep_dni*( 1.0_rp-sw )/( rhoq2(k,i_ni)+sw ) * rhoq2_crg(i_qi,k)
5186 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_ns)-small )
5187 dep_dns_crg = dep_dns*( 1.0_rp-sw )/( rhoq2(k,i_ns)+sw ) * rhoq2_crg(i_qs,k)
5188 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_ng)-small )
5189 dep_dng_crg = dep_dng*( 1.0_rp-sw )/( rhoq2(k,i_ng)+sw ) * rhoq2_crg(i_qg,k)
5191 sw = min( abs(rhoq2_crg(i_qi,k)),abs(dep_dni_crg) )
5192 dep_dni_crg = sign( sw,dep_dni_crg )
5193 sw = min( abs(rhoq2_crg(i_qs,k)),abs(dep_dns_crg) )
5194 dep_dns_crg = sign( sw,dep_dns_crg )
5195 sw = min( abs(rhoq2_crg(i_qg,k)),abs(dep_dng_crg) )
5196 dep_dng_crg = sign( sw,dep_dng_crg )
5200 frz_dqc = max( dt*(pq(k,i_lchom)+pq(k,i_lchet)), -rhoq2(k,i_qc)-dep_dqc )
5201 frz_dnc = max( dt*(pq(k,i_nchom)+pq(k,i_nchet)), -rhoq2(k,i_nc)-dep_dnc )
5202 fac3 = ( frz_dqc-eps )/( dt*(pq(k,i_lchom)+pq(k,i_lchet))-eps )
5203 fac4 = ( frz_dnc-eps )/( dt*(pq(k,i_nchom)+pq(k,i_nchet))-eps )
5204 pq(k,i_lchom) = fac3*pq(k,i_lchom)
5205 pq(k,i_lchet) = fac3*pq(k,i_lchet)
5206 pq(k,i_nchom) = fac4*pq(k,i_nchom)
5207 pq(k,i_nchet) = fac4*pq(k,i_nchet)
5209 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_nc)-small )
5210 frz_dnc_crg = frz_dnc*( 1.0_rp-sw )/( rhoq2(k,i_nc)+sw ) * rhoq2_crg(i_qc,k)
5212 sw = min( abs(rhoq2_crg(i_qc,k)+dep_dnc_crg),abs(frz_dnc_crg) )
5213 frz_dnc_crg = sign( sw,frz_dnc_crg )
5218 mlt_dqi = max( dt*pq(k,i_limlt), -rhoq2(k,i_qi)-dep_dqi )
5219 mlt_dni = max( dt*pq(k,i_nimlt), -rhoq2(k,i_ni)-dep_dni )
5221 mlt_dqs = max( dt*pq(k,i_lsmlt), -rhoq2(k,i_qs)-dep_dqs )
5222 mlt_dns = max( dt*pq(k,i_nsmlt), -rhoq2(k,i_ns)-dep_dns )
5224 mlt_dqg = max( dt*pq(k,i_lgmlt), -rhoq2(k,i_qg)-dep_dqg )
5225 mlt_dng = max( dt*pq(k,i_ngmlt), -rhoq2(k,i_ng)-dep_dng )
5228 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_ni)-small )
5229 mlt_dni_crg = mlt_dni*( 1.0_rp-sw ) / ( rhoq2(k,i_ni)+sw ) * rhoq2_crg(i_qi,k)
5230 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_ns)-small )
5231 mlt_dns_crg = mlt_dns*( 1.0_rp-sw ) / ( rhoq2(k,i_ns)+sw ) * rhoq2_crg(i_qs,k)
5232 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_ng)-small )
5233 mlt_dng_crg = mlt_dng*( 1.0_rp-sw ) / ( rhoq2(k,i_ng)+sw ) * rhoq2_crg(i_qg,k)
5236 sw = min( abs(rhoq2_crg(i_qi,k)+dep_dni_crg-frz_dnc_crg),abs(mlt_dni_crg) )
5237 mlt_dni_crg = sign( sw, mlt_dni_crg )
5238 sw = min( abs(rhoq2_crg(i_qs,k)+dep_dns_crg ),abs(mlt_dns_crg) )
5239 mlt_dns_crg = sign( sw, mlt_dns_crg )
5240 sw = min( abs(rhoq2_crg(i_qg,k)+dep_dng_crg ),abs(mlt_dng_crg) )
5241 mlt_dng_crg = sign( sw, mlt_dng_crg )
5245 frz_dqr = max( dt*(pq(k,i_lrhet)), min(0.0_rp, -rhoq2(k,i_qr)-dep_dqr) )
5246 frz_dnr = max( dt*(pq(k,i_nrhet)), min(0.0_rp, -rhoq2(k,i_nr)-dep_dnr) )
5249 sw = 0.5_rp - sign( 0.5_rp, rhoq2(k,i_nr)-small )
5250 frz_dnr_crg = frz_dnr*( 1.0_rp-sw ) /( rhoq2(k,i_nr)+sw ) * rhoq2_crg(i_qr,k)
5252 sw = min( abs(rhoq2_crg(i_qr,k)+dep_dnr_crg ),abs(frz_dnr_crg) )
5253 frz_dnr_crg = sign( sw,frz_dnr_crg )
5256 fac5 = ( frz_dqr-eps )/( dt*pq(k,i_lrhet)-eps )
5257 pq(k,i_lrhet) = fac5*pq(k,i_lrhet)
5258 fac6 = ( frz_dnr-eps )/( dt*pq(k,i_nrhet)-eps )
5259 pq(k,i_nrhet) = fac6*pq(k,i_nrhet)
5262 drhoqv = -(dep_dqc+dep_dqi+dep_dqs+dep_dqg+dep_dqr)
5265 sw = 0.5_rp - sign(0.5_rp, abs(drhoqv) - eps)
5266 fact = ( max( rhoq2(k,i_qv) + drhoqv * dt, 0.0_rp ) - rhoq2(k,i_qv) ) / dt / ( drhoqv + sw ) * ( 1.0_rp - sw ) &
5268 fact = min( 1.0_rp, max( 0.0_rp, fact ) )
5270 drhoqv = drhoqv * fact
5271 dep_dqc = dep_dqc * fact
5272 dep_dnc = dep_dnc * fact
5273 dep_dqr = dep_dqr * fact
5274 dep_dnr = dep_dnr * fact
5275 dep_dqi = dep_dqi * fact
5276 dep_dni = dep_dni * fact
5277 dep_dqs = dep_dqs * fact
5278 dep_dns = dep_dns * fact
5279 dep_dqg = dep_dqg * fact
5280 dep_dng = dep_dng * fact
5283 xi = min(xi_max, max(xi_min, rhoq2(k,i_qi)/(rhoq2(k,i_ni)+ni_min) ))
5284 sw = 0.5_rp + sign(0.5_rp,xi-x_sep)
5288 drhoqc = ( frz_dqc - mlt_dqi*(1.0_rp-sw) + dep_dqc )
5289 drhonc = ( frz_dnc - mlt_dni*(1.0_rp-sw) + dep_dnc )
5292 drhoqr = ( frz_dqr - mlt_dqg - mlt_dqs - mlt_dqi*sw + dep_dqr )
5293 drhonr = ( frz_dnr - mlt_dng - mlt_dns - mlt_dni*sw + dep_dnr )
5296 drhoqi = (-frz_dqc + mlt_dqi + dep_dqi )
5297 drhoni = (-frz_dnc + mlt_dni + dep_dni )
5300 drhoqs = ( mlt_dqs + dep_dqs )
5301 drhons = ( mlt_dns + dep_dns )
5304 drhoqg = (-frz_dqr + mlt_dqg + dep_dqg )
5305 drhong = (-frz_dnr + mlt_dng + dep_dng )
5309 drhoqcrg_c = ( frz_dnc_crg - mlt_dni_crg*(1.0_rp-sw) + dep_dnc_crg )
5310 drhoqcrg_r = ( frz_dnr_crg - mlt_dni_crg*sw - mlt_dng_crg - mlt_dns_crg + dep_dnr_crg )
5311 drhoqcrg_i = ( -frz_dnc_crg + mlt_dni_crg + dep_dni_crg )
5312 drhoqcrg_s = ( mlt_dns_crg + dep_dns_crg )
5313 drhoqcrg_g = ( -frz_dnr_crg + mlt_dng_crg + dep_dng_crg )
5317 rhoq_t(k,i_qv) = drhoqv / dt
5318 rhoq_t(k,i_qc) = drhoqc / dt
5319 rhoq_t(k,i_nc) = drhonc / dt
5320 rhoq_t(k,i_qr) = drhoqr / dt
5321 rhoq_t(k,i_nr) = drhonr / dt
5322 rhoq_t(k,i_qi) = drhoqi / dt
5323 rhoq_t(k,i_ni) = drhoni / dt
5324 rhoq_t(k,i_qs) = drhoqs / dt
5325 rhoq_t(k,i_ns) = drhons / dt
5326 rhoq_t(k,i_qg) = drhoqg / dt
5327 rhoq_t(k,i_ng) = drhong / dt
5329 rhoe_t(k) = ( - lhv * drhoqv + lhf * ( drhoqi + drhoqs + drhoqg ) ) / dt
5333 rhoqcrg_t(k,i_qc-1) = drhoqcrg_c / dt
5334 rhoqcrg_t(k,i_qr-1) = drhoqcrg_r / dt
5335 rhoqcrg_t(k,i_qi-1) = drhoqcrg_i / dt
5336 rhoqcrg_t(k,i_qs-1) = drhoqcrg_s / dt
5337 rhoqcrg_t(k,i_qg-1) = drhoqcrg_g / dt
5340 rrho = 1.0_rp/rho(k)
5342 dql = rrho * ( drhoqc + drhoqr )
5343 dqi = rrho * ( drhoqi + drhoqs + drhoqg )
5351 dz = fz(k) - fz(k-1)
5352 sl_plcdep = sl_plcdep + dep_dqc*dz
5353 sl_plrdep = sl_plrdep + dep_dqr*dz
5354 sl_pnrdep = sl_pnrdep + dep_dnr*dz
5358 end subroutine update_by_phase_change
5369 integer,
intent(in) :: KA, KS, KE
5370 integer,
intent(in) :: QA_MP
5371 real(RP),
intent(out) :: Crs (KA,HYDRO_MAX)
5372 real(RP),
intent(in) :: QTRC0(KA,QA_MP)
5373 real(RP),
intent(in) :: DENS0(KA)
5382 real(RP) :: dc_ave(KA)
5383 real(RP) :: dr_ave(KA)
5387 real(RP) :: coef_Fuetal1998
5389 real(RP),
parameter :: r2m_min=1.e-12_rp
5390 real(RP),
parameter :: um2cm = 100.0_rp
5392 real(RP) :: limitsw, zerosw
5398 xc(k) = min( xc_max, max( xc_min, dens0(k)*qtrc0(k,i_qc)/(qtrc0(k,i_nc)+nc_min) ) )
5399 xr(k) = min( xr_max, max( xr_min, dens0(k)*qtrc0(k,i_qr)/(qtrc0(k,i_nr)+nr_min) ) )
5400 xi(k) = min( xi_max, max( xi_min, dens0(k)*qtrc0(k,i_qi)/(qtrc0(k,i_ni)+ni_min) ) )
5401 xs(k) = min( xs_max, max( xs_min, dens0(k)*qtrc0(k,i_qs)/(qtrc0(k,i_ns)+ns_min) ) )
5402 xg(k) = min( xg_max, max( xg_min, dens0(k)*qtrc0(k,i_qg)/(qtrc0(k,i_ng)+ng_min) ) )
5407 dc_ave(k) = a_m(i_qc) * xc(k)**b_m(i_qc)
5408 dr_ave(k) = a_m(i_qr) * xr(k)**b_m(i_qr)
5412 crs(k,i_mp_qc) = pi * coef_r2(i_mp_qc) * qtrc0(k,i_nc) * a_rea2(i_mp_qc) * xc(k)**b_rea2(i_mp_qc)
5416 crs(k,i_mp_qr) = pi * coef_r2(i_mp_qr) * qtrc0(k,i_nr) * a_rea2(i_mp_qr) * xr(k)**b_rea2(i_mp_qr)
5420 crs(k,i_mp_qi) = pi * coef_rea2(i_mp_qi) * qtrc0(k,i_ni) * a_rea2(i_mp_qi) * xi(k)**b_rea2(i_mp_qi)
5421 crs(k,i_mp_qs) = pi * coef_rea2(i_mp_qs) * qtrc0(k,i_ns) * a_rea2(i_mp_qs) * xs(k)**b_rea2(i_mp_qs)
5422 crs(k,i_mp_qg) = pi * coef_rea2(i_mp_qg) * qtrc0(k,i_ng) * a_rea2(i_mp_qg) * xg(k)**b_rea2(i_mp_qg)