51 #define PROFILE_START(name) call fapp_start(name, 1, 1) 52 #define PROFILE_STOP(name) call fapp_stop (name, 1, 1) 53 #elif defined(PROFILE_FINEPA) 54 #define PROFILE_START(name) call start_collection(name) 55 #define PROFILE_STOP(name) call stop_collection (name) 57 #define PROFILE_START(name) 58 #define PROFILE_STOP(name) 119 integer,
public,
parameter ::
qa_mp = 11
137 'Ratio of Water Vapor mass to total mass (Specific humidity)', &
138 'Ratio of Cloud Water mass to total mass ', &
139 'Ratio of Rain Water mass to total mass ', &
140 'Ratio of Cloud Ice mass ratio to total mass ', &
141 'Ratio of Snow miass ratio to total mass ', &
142 'Ratio of Graupel mass ratio to total mass ', &
143 'Cloud Water Number Density ', &
144 'Rain Water Number Density ', &
145 'Cloud Ice Number Density ', &
146 'Snow Number Density ', &
147 'Graupel Number Density '/)
165 private :: mp_sn14_init
172 integer,
private,
parameter :: i_qv = 1
173 integer,
private,
parameter :: i_qc = 2
174 integer,
private,
parameter :: i_qr = 3
175 integer,
private,
parameter :: i_qi = 4
176 integer,
private,
parameter :: i_qs = 5
177 integer,
private,
parameter :: i_qg = 6
178 integer,
private,
parameter :: i_nc = 7
179 integer,
private,
parameter :: i_nr = 8
180 integer,
private,
parameter :: i_ni = 9
181 integer,
private,
parameter :: i_ns = 10
182 integer,
private,
parameter :: i_ng = 11
184 integer,
private,
parameter :: hydro_max = 5
186 integer,
private,
parameter :: i_mp_qc = 1
187 integer,
private,
parameter :: i_mp_qr = 2
188 integer,
private,
parameter :: i_mp_qi = 3
189 integer,
private,
parameter :: i_mp_qs = 4
190 integer,
private,
parameter :: i_mp_qg = 5
191 integer,
private,
parameter :: i_mp_nc = 6
192 integer,
private,
parameter :: i_mp_nr = 7
193 integer,
private,
parameter :: i_mp_ni = 8
194 integer,
private,
parameter :: i_mp_ns = 9
195 integer,
private,
parameter :: i_mp_ng = 10
199 integer,
parameter :: i_lcccn = 1
200 integer,
parameter :: i_ncccn = 2
201 integer,
parameter :: i_liccn = 3
202 integer,
parameter :: i_niccn = 4
204 integer,
parameter :: i_lchom = 5
205 integer,
parameter :: i_nchom = 6
206 integer,
parameter :: i_lchet = 7
207 integer,
parameter :: i_nchet = 8
208 integer,
parameter :: i_lrhet = 9
209 integer,
parameter :: i_nrhet = 10
211 integer,
parameter :: i_limlt = 11
212 integer,
parameter :: i_nimlt = 12
213 integer,
parameter :: i_lsmlt = 13
214 integer,
parameter :: i_nsmlt = 14
215 integer,
parameter :: i_lgmlt = 15
216 integer,
parameter :: i_ngmlt = 16
218 integer,
parameter :: i_lrdep = 17
219 integer,
parameter :: i_nrdep = 18
220 integer,
parameter :: i_lidep = 19
221 integer,
parameter :: i_nidep = 20
222 integer,
parameter :: i_lsdep = 21
223 integer,
parameter :: i_nsdep = 22
224 integer,
parameter :: i_lgdep = 23
225 integer,
parameter :: i_ngdep = 24
226 integer,
parameter :: i_lcdep = 25
230 integer,
parameter :: i_lcaut = 26
231 integer,
parameter :: i_ncaut = 27
232 integer,
parameter :: i_nraut = 28
234 integer,
parameter :: i_lcacc = 29
235 integer,
parameter :: i_ncacc = 30
237 integer,
parameter :: i_nrslc = 31
238 integer,
parameter :: i_nrbrk = 32
241 integer,
parameter :: i_licon = 33
242 integer,
parameter :: i_nicon = 34
243 integer,
parameter :: i_lscon = 35
244 integer,
parameter :: i_nscon = 36
246 integer,
parameter :: i_liacm = 37
247 integer,
parameter :: i_niacm = 38
248 integer,
parameter :: i_liarm = 39
249 integer,
parameter :: i_niarm = 40
250 integer,
parameter :: i_lsacm = 41
251 integer,
parameter :: i_nsacm = 42
252 integer,
parameter :: i_lsarm = 43
253 integer,
parameter :: i_nsarm = 44
254 integer,
parameter :: i_lgacm = 45
255 integer,
parameter :: i_ngacm = 46
256 integer,
parameter :: i_lgarm = 47
257 integer,
parameter :: i_ngarm = 48
259 integer,
parameter :: i_lgspl = 49
260 integer,
parameter :: i_lsspl = 50
261 integer,
parameter :: i_nispl = 51
263 integer,
parameter :: pq_max = 51
267 integer,
parameter :: i_liaclc2li = 1
268 integer,
parameter :: i_niacnc2ni = 2
269 integer,
parameter :: i_lsaclc2ls = 3
270 integer,
parameter :: i_nsacnc2ns = 4
271 integer,
parameter :: i_lgaclc2lg = 5
272 integer,
parameter :: i_ngacnc2ng = 6
273 integer,
parameter :: i_lracli2lg_i = 7
274 integer,
parameter :: i_nracni2ng_i = 8
275 integer,
parameter :: i_lracli2lg_r = 9
276 integer,
parameter :: i_nracni2ng_r = 10
277 integer,
parameter :: i_lracls2lg_s = 11
278 integer,
parameter :: i_nracns2ng_s = 12
279 integer,
parameter :: i_lracls2lg_r = 13
280 integer,
parameter :: i_nracns2ng_r = 14
281 integer,
parameter :: i_lraclg2lg = 15
282 integer,
parameter :: i_nracng2ng = 16
283 integer,
parameter :: i_liacli2ls = 17
284 integer,
parameter :: i_niacni2ns = 18
285 integer,
parameter :: i_liacls2ls = 19
286 integer,
parameter :: i_niacns2ns = 20
287 integer,
parameter :: i_nsacns2ns = 21
288 integer,
parameter :: i_ngacng2ng = 22
289 integer,
parameter :: i_lgacls2lg = 23
290 integer,
parameter :: i_ngacns2ng = 24
292 integer,
parameter :: pac_max = 24
296 character(len=H_SHORT),
save :: wlabel(hydro_max)
299 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 /)
302 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 /)
305 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 /)
311 real(RP),
private,
parameter :: xc_min = 4.20e-15_rp
312 real(RP),
private,
parameter :: xr_min = 2.60e-10_rp
313 real(RP),
private,
parameter :: xi_min = 3.382e-13_rp
314 real(RP),
private,
parameter :: xs_min = 1.847e-12_rp
315 real(RP),
private,
parameter :: xg_min = 1.230e-10_rp
317 real(RP),
private,
parameter :: xc_max = 2.6e-10_rp
318 real(RP),
private,
parameter :: xr_max = 5.00e-6_rp
319 real(RP),
private,
parameter :: xi_max = 1.377e-6_rp
320 real(RP),
private,
parameter :: xs_max = 7.519e-6_rp
321 real(RP),
private,
parameter :: xg_max = 4.900e-5_rp
323 real(RP),
private,
parameter :: xmin_filter= xc_min
325 real(RP),
private,
parameter :: rmin_re= 1.e-6_rp
328 real(RP),
private,
parameter :: n0r_min= 2.5e+5_rp
329 real(RP),
private,
parameter :: n0r_max= 2.0e+7_rp
330 real(RP),
private,
parameter :: lambdar_min= 1.e+3_rp
331 real(RP),
private,
parameter :: lambdar_max= 1.e+4_rp
333 real(RP),
private,
parameter :: nc_min = 1.e+4_rp
334 real(RP),
private,
parameter :: nr_min = 1.0_rp
335 real(RP),
private,
parameter :: ni_min = 1.0_rp
336 real(RP),
private,
parameter :: ns_min = 1.e-4_rp
337 real(RP),
private,
parameter :: ng_min = 1.e-4_rp
339 real(RP),
private,
parameter :: lc_min = xc_min*nc_min
340 real(RP),
private,
parameter :: lr_min = xr_min*nr_min
341 real(RP),
private,
parameter :: li_min = xi_min*ni_min
342 real(RP),
private,
parameter :: ls_min = xs_min*ns_min
343 real(RP),
private,
parameter :: lg_min = xg_min*ng_min
345 real(RP),
private,
parameter :: x_sep = 2.6e-10_rp
347 real(RP),
private,
parameter :: tem_min=100.0_rp
348 real(RP),
private,
parameter :: rho_min=1.e-5_rp
349 real(RP),
private,
parameter :: rhoi = 916.70_rp
351 integer,
private,
save :: ntmax_phase_change = 1
352 integer,
private,
save :: ntmax_collection = 1
355 real(RP),
private,
parameter :: rho_0 = 1.280_rp
357 real(RP),
allocatable,
private,
save :: nc_uplim_d(:,:,:)
360 real(RP),
private,
parameter :: ka0 = 2.428e-2_rp
362 real(RP),
private,
parameter :: dka_dt = 7.47e-5_rp
367 real(RP),
private,
parameter :: mua0 = 1.718e-5_rp
369 real(RP),
private,
parameter :: dmua_dt = 5.28e-8_rp
373 real(RP),
private,
save :: xc_ccn = 1.e-12_rp
374 real(RP),
private,
save :: xi_ccn = 1.e-12_rp
378 real(RP),
private,
save :: cap(hydro_max)
382 real(RP),
private,
save :: a_m(hydro_max)
383 real(RP),
private,
save :: b_m(hydro_max)
386 real(RP),
private,
save :: alpha_v(hydro_max,2)
387 real(RP),
private,
save :: beta_v(hydro_max,2)
388 real(RP),
private,
save :: alpha_vn(hydro_max,2)
389 real(RP),
private,
save :: beta_vn(hydro_max,2)
390 real(RP),
private,
save :: gamma_v(hydro_max)
393 real(RP),
private,
parameter :: pre0_vt = 300.e+2_rp
394 real(RP),
private,
parameter :: tem0_vt = 233.0_rp
395 real(RP),
private,
parameter :: a_pre0_vt = -0.1780_rp
396 real(RP),
private,
parameter :: a_tem0_vt = -0.3940_rp
404 real(RP),
private,
save :: nu(hydro_max)
405 real(RP),
private,
save :: mu(hydro_max)
411 real(RP),
private,
save :: a_area(hydro_max)
412 real(RP),
private,
save :: b_area(hydro_max)
413 real(RP),
private,
save :: ax_area(hydro_max)
414 real(RP),
private,
save :: bx_area(hydro_max)
417 real(RP),
private,
save :: a_rea(hydro_max)
418 real(RP),
private,
save :: b_rea(hydro_max)
419 real(RP),
private,
save :: a_rea2(hydro_max)
420 real(RP),
private,
save :: b_rea2(hydro_max)
421 real(RP),
private,
save :: a_rea3(hydro_max)
422 real(RP),
private,
save :: b_rea3(hydro_max)
424 real(RP),
private,
save :: a_d2vt(hydro_max)
425 real(RP),
private,
save :: b_d2vt(hydro_max)
429 real(RP),
private,
save :: coef_m2(hydro_max)
431 real(RP),
private,
save :: coef_d6(hydro_max)
433 real(RP),
private,
save :: coef_d3(hydro_max)
435 real(RP),
private,
save :: coef_d(hydro_max)
437 real(RP),
private,
save :: coef_d2v(hydro_max)
439 real(RP),
private,
save :: coef_md2v(hydro_max)
442 real(RP),
private,
save :: coef_r2(hydro_max)
443 real(RP),
private,
save :: coef_r3(hydro_max)
444 real(RP),
private,
save :: coef_re(hydro_max)
446 real(RP),
private,
save :: coef_rea2(hydro_max)
447 real(RP),
private,
save :: coef_rea3(hydro_max)
448 logical,
private,
save :: opt_m96_ice=.true.
449 logical,
private,
save :: opt_m96_column_ice=.false.
454 real(RP),
private,
save :: coef_vt0(hydro_max,2)
455 real(RP),
private,
save :: coef_vt1(hydro_max,2)
456 real(RP),
private,
save :: coef_deplc
457 real(RP),
private,
save :: coef_dave_n(hydro_max)
458 real(RP),
private,
save :: coef_dave_l(hydro_max)
461 real(RP),
private,
save :: d0_ni=261.76e-6_rp
462 real(RP),
private,
save :: d0_li=398.54e-6_rp
463 real(RP),
private,
parameter :: d0_ns=270.03e-6_rp
464 real(RP),
private,
parameter :: d0_ls=397.47e-6_rp
465 real(RP),
private,
parameter :: d0_ng=269.08e-6_rp
466 real(RP),
private,
parameter :: d0_lg=376.36e-6_rp
468 real(RP),
private,
parameter :: coef_vtr_ar1=9.65_rp
470 real(RP),
private,
parameter :: coef_vtr_br1=10.43_rp
471 real(RP),
private,
parameter :: coef_vtr_cr1=600.0_rp
472 real(RP),
private,
parameter :: coef_vtr_ar2=4.e+3_rp
473 real(RP),
private,
parameter :: coef_vtr_br2=12.e+3_rp
474 real(RP),
private,
parameter :: d_vtr_branch=0.745e-3_rp
476 real(RP),
private,
parameter :: dr_eq = 1.10e-3_rp
481 real(RP),
private,
save :: coef_a(hydro_max)
482 real(RP),
private,
save :: coef_lambda(hydro_max)
486 real(RP),
private,
save :: ah_vent (hydro_max,2)
487 real(RP),
private,
save :: bh_vent (hydro_max,2)
488 real(RP),
private,
save :: ah_vent0 (hydro_max,2)
489 real(RP),
private,
save :: bh_vent0 (hydro_max,2)
490 real(RP),
private,
save :: ah_vent1 (hydro_max,2)
491 real(RP),
private,
save :: bh_vent1 (hydro_max,2)
493 real(RP),
private,
save :: delta_b0 (hydro_max)
494 real(RP),
private,
save :: delta_b1 (hydro_max)
495 real(RP),
private,
save :: delta_ab0(hydro_max,hydro_max)
496 real(RP),
private,
save :: delta_ab1(hydro_max,hydro_max)
498 real(RP),
private,
save :: theta_b0 (hydro_max)
499 real(RP),
private,
save :: theta_b1 (hydro_max)
500 real(RP),
private,
save :: theta_ab0(hydro_max,hydro_max)
501 real(RP),
private,
save :: theta_ab1(hydro_max,hydro_max)
503 logical,
private,
save :: opt_debug=.false.
505 logical,
private,
save :: opt_debug_tem=.false.
506 logical,
private,
save :: opt_debug_inc=.true.
507 logical,
private,
save :: opt_debug_act=.true.
508 logical,
private,
save :: opt_debug_ree=.true.
509 logical,
private,
save :: opt_debug_bcs=.true.
511 logical,
private,
save :: mp_doautoconversion = .true.
512 logical,
private,
save :: mp_couple_aerosol = .false.
513 real(RP),
private,
save :: mp_ssw_lim = 1.e+1_rp
527 integer,
intent(in) :: KA
528 integer,
intent(in) :: IA
529 integer,
intent(in) :: JA
531 namelist / param_atmos_phy_mp_sn14 / &
532 mp_doautoconversion, &
540 log_info(
"ATMOS_PHY_MP_sn14_setup",*)
'Setup' 541 log_info(
"ATMOS_PHY_MP_sn14_setup",*)
'Seiki and Nakajima (2014) 2-moment bulk 6 category' 545 read(
io_fid_conf,nml=param_atmos_phy_mp_sn14,iostat=ierr)
547 log_info(
"ATMOS_PHY_MP_sn14_setup",*)
'Not found namelist. Default used.' 548 elseif( ierr > 0 )
then 549 log_error(
"ATMOS_PHY_MP_sn14_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14. Check!' 552 log_nml(param_atmos_phy_mp_sn14)
558 wlabel(5) =
"GRAUPEL" 562 allocate(nc_uplim_d(1,ia,ja))
563 nc_uplim_d(:,:,:) = 150.e6_rp
573 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
593 integer,
intent(in) :: KA, KS, KE
594 integer,
intent(in) :: IA, IS, IE
595 integer,
intent(in) :: JA, JS, JE
597 real(RP),
intent(in) :: DENS (ka,ia,ja)
598 real(RP),
intent(in) :: W (ka,ia,ja)
599 real(RP),
intent(in) :: QTRC (ka,ia,ja,
qa_mp)
600 real(RP),
intent(in) :: PRES(ka,ia,ja)
601 real(RP),
intent(in) :: TEMP(ka,ia,ja)
602 real(RP),
intent(in) :: Qdry(ka,ia,ja)
603 real(RP),
intent(in) :: CPtot(ka,ia,ja)
604 real(RP),
intent(in) :: CVtot(ka,ia,ja)
605 real(RP),
intent(in) :: CCN (ka,ia,ja)
606 real(DP),
intent(in) :: dt
607 real(RP),
intent(in) :: cz( ka,ia,ja)
608 real(RP),
intent(in) :: fz(0:ka,ia,ja)
610 real(RP),
intent(out) :: RHOQ_t (ka,ia,ja,
qa_mp)
611 real(RP),
intent(out) :: RHOE_t (ka,ia,ja)
612 real(RP),
intent(out) :: CPtot_t(ka,ia,ja)
613 real(RP),
intent(out) :: CVtot_t(ka,ia,ja)
614 real(RP),
intent(out) :: EVAPORATE(ka,ia,ja)
618 log_progress(*)
'atmosphere / physics / microphysics / SN14' 626 ka, ks, ke, ia, is, ie, ja, js, je, &
627 dens(:,:,:), w(:,:,:), qtrc(:,:,:,:), pres(:,:,:), temp(:,:,:), &
628 qdry(:,:,:), cptot(:,:,:), cvtot(:,:,:), ccn(:,:,:), &
629 real(dt,RP), cz(:,:,:), fz(:,:,:), &
630 RHOQ_t(:,:,:,:), RHOE_t(:,:,:), CPtot_t(:,:,:), CVtot_t(:,:,:), &
645 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
650 integer,
intent(in) :: KA, KS, KE
651 integer,
intent(in) :: IA, IS, IE
652 integer,
intent(in) :: JA, JS, JE
654 real(RP),
intent(in) :: QTRC (ka,ia,ja,
qa_mp-1)
655 real(RP),
intent(in) :: mask_criterion
657 real(RP),
intent(out) :: cldfrac(ka,ia,ja)
660 integer :: k, i, j, iq
668 qhydro = qhydro + qtrc(k,i,j,iq)
670 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
682 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
683 DENS0, TEMP0, QTRC0, &
693 integer,
intent(in) :: KA, KS, KE
694 integer,
intent(in) :: IA, IS, IE
695 integer,
intent(in) :: JA, JS, JE
697 real(RP),
intent(in) :: DENS0(ka,ia,ja)
698 real(RP),
intent(in) :: TEMP0(ka,ia,ja)
699 real(RP),
intent(in) :: QTRC0(ka,ia,ja,i_qc:i_ng)
701 real(RP),
intent(out) :: Re (ka,ia,ja,
n_hyd)
704 real(RP) :: xc(ka,ia,ja)
705 real(RP) :: xr(ka,ia,ja)
706 real(RP) :: xi(ka,ia,ja)
707 real(RP) :: xs(ka,ia,ja)
708 real(RP) :: xg(ka,ia,ja)
710 real(RP) :: dc_ave(ka,ia,ja)
711 real(RP) :: dr_ave(ka,ia,ja)
715 real(RP) :: ri2m(ka,ia,ja), ri3m(ka,ia,ja)
716 real(RP) :: rs2m(ka,ia,ja), rs3m(ka,ia,ja)
717 real(RP) :: rg2m(ka,ia,ja), rg3m(ka,ia,ja)
719 real(RP) :: coef_Fuetal1998
721 real(RP),
parameter :: r2m_min=1.e-12_rp
722 real(RP),
parameter :: um2cm = 100.0_rp
724 real(RP) :: limitsw, zerosw
732 xc(k,i,j) = min( xc_max, max( xc_min, dens0(k,i,j)*qtrc0(k,i,j,i_qc)/(qtrc0(k,i,j,i_nc)+nc_min) ) )
733 xr(k,i,j) = min( xr_max, max( xr_min, dens0(k,i,j)*qtrc0(k,i,j,i_qr)/(qtrc0(k,i,j,i_nr)+nr_min) ) )
734 xi(k,i,j) = min( xi_max, max( xi_min, dens0(k,i,j)*qtrc0(k,i,j,i_qi)/(qtrc0(k,i,j,i_ni)+ni_min) ) )
735 xs(k,i,j) = min( xs_max, max( xs_min, dens0(k,i,j)*qtrc0(k,i,j,i_qs)/(qtrc0(k,i,j,i_ns)+ns_min) ) )
736 xg(k,i,j) = min( xg_max, max( xg_min, dens0(k,i,j)*qtrc0(k,i,j,i_qg)/(qtrc0(k,i,j,i_ng)+ng_min) ) )
745 dc_ave(k,i,j) = a_m(i_mp_qc) * xc(k,i,j)**b_m(i_mp_qc)
746 dr_ave(k,i,j) = a_m(i_mp_qr) * xr(k,i,j)**b_m(i_mp_qr)
755 rc = 0.5_rp * dc_ave(k,i,j)
756 limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
757 re(k,i,j,
i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
766 rr = 0.5_rp * dr_ave(k,i,j)
767 limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
768 re(k,i,j,
i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
776 ri2m(k,i,j) = pi * coef_rea2(i_mp_qi) * qtrc0(k,i,j,i_ni) * a_rea2(i_mp_qi) * xi(k,i,j)**b_rea2(i_mp_qi)
777 rs2m(k,i,j) = pi * coef_rea2(i_mp_qs) * qtrc0(k,i,j,i_ns) * a_rea2(i_mp_qs) * xs(k,i,j)**b_rea2(i_mp_qs)
778 rg2m(k,i,j) = pi * coef_rea2(i_mp_qg) * qtrc0(k,i,j,i_ng) * a_rea2(i_mp_qg) * xg(k,i,j)**b_rea2(i_mp_qg)
784 coef_fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
788 ri3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k,i,j)
789 rs3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k,i,j)
790 rg3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k,i,j)
799 zerosw = 0.5_rp - sign(0.5_rp, ri2m(k,i,j) - r2m_min )
800 re(k,i,j,
i_hi) = ri3m(k,i,j) / ( ri2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
809 zerosw = 0.5_rp - sign(0.5_rp, rs2m(k,i,j) - r2m_min )
810 re(k,i,j,
i_hs) = rs3m(k,i,j) / ( rs2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
819 zerosw = 0.5_rp - sign(0.5_rp, rg2m(k,i,j) - r2m_min )
820 re(k,i,j,
i_hg) = rg3m(k,i,j) / ( rg2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
825 re(:,:,:,
i_hg+1:) = 0.0_rp
834 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
845 integer,
intent(in) :: KA, KS, KE
846 integer,
intent(in) :: IA, IS, IE
847 integer,
intent(in) :: JA, JS, JE
849 real(RP),
intent(in) :: QTRC0(ka,ia,ja,
qa_mp-1)
851 real(RP),
intent(out) :: Qe (ka,ia,ja,
n_hyd)
855 qe(:,:,:,
i_hc) = qtrc0(:,:,:,i_mp_qc)
857 qe(:,:,:,
i_hr) = qtrc0(:,:,:,i_mp_qr)
859 qe(:,:,:,
i_hi) = qtrc0(:,:,:,i_mp_qi)
861 qe(:,:,:,
i_hs) = qtrc0(:,:,:,i_mp_qs)
863 qe(:,:,:,
i_hg) = qtrc0(:,:,:,i_mp_qg)
865 qe(:,:,:,
i_hg+1:) = 0.0_rp
872 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
883 integer,
intent(in) :: KA, KS, KE
884 integer,
intent(in) :: IA, IS, IE
885 integer,
intent(in) :: JA, JS, JE
887 real(RP),
intent(in) :: QTRC0(ka,ia,ja,
qa_mp-1)
889 real(RP),
intent(out) :: Ne (ka,ia,ja,
n_hyd)
893 ne(:,:,:,
i_hc) = qtrc0(:,:,:,i_mp_nc)
895 ne(:,:,:,
i_hr) = qtrc0(:,:,:,i_mp_nr)
897 ne(:,:,:,
i_hi) = qtrc0(:,:,:,i_mp_ni)
899 ne(:,:,:,
i_hs) = qtrc0(:,:,:,i_mp_ns)
901 ne(:,:,:,
i_hg) = qtrc0(:,:,:,i_mp_ng)
903 ne(:,:,:,
i_hg+1:) = 0.0_rp
909 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
924 integer,
intent(in) :: KA, KS, KE
925 integer,
intent(in) :: IA, IS, IE
926 integer,
intent(in) :: JA, JS, JE
928 real(RP),
intent(in) :: Qe(ka,ia,ja,
n_hyd)
930 real(RP),
intent(out) :: QTRC(ka,ia,ja,
qa_mp-1)
932 real(RP),
intent(in),
optional :: QNUM(ka,ia,ja,
n_hyd)
934 real(RP),
parameter :: Dc = 20.e-6_rp
935 real(RP),
parameter :: Dr = 200.e-6_rp
936 real(RP),
parameter :: Di = 80.e-6_rp
937 real(RP),
parameter :: Ds = 80.e-6_rp
938 real(RP),
parameter :: Dg = 200.e-6_rp
939 real(RP),
parameter :: RHOw = 1000.0_rp
940 real(RP),
parameter :: RHOf = 100.0_rp
941 real(RP),
parameter :: RHOg = 400.0_rp
942 real(RP),
parameter :: b = 3.0_rp
955 qtrc(k,i,j,i_mp_qc) = qe(k,i,j,
i_hc)
965 qtrc(k,i,j,i_mp_qr) = qe(k,i,j,
i_hr)
975 qtrc(k,i,j,i_mp_qi) = qe(k,i,j,
i_hi)
985 qtrc(k,i,j,i_mp_qs) = qe(k,i,j,
i_hs)
995 qtrc(k,i,j,i_mp_qg) = qe(k,i,j,
i_hg) + qe(k,i,j,
i_hh)
1000 if (
present(qnum) )
then 1007 qtrc(k,i,j,i_mp_nc) = qnum(k,i,j,
i_hc)
1017 qtrc(k,i,j,i_mp_nr) = qnum(k,i,j,
i_hr)
1027 qtrc(k,i,j,i_mp_ni) = qnum(k,i,j,
i_hi)
1037 qtrc(k,i,j,i_mp_ns) = qnum(k,i,j,
i_hs)
1047 qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,
i_hg) + qnum(k,i,j,
i_hh)
1060 qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1070 qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1080 qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1090 qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1100 qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1113 subroutine mp_sn14_init
1120 real(RP) :: w1(hydro_max)
1121 real(RP) :: w2(hydro_max)
1122 real(RP) :: w3(hydro_max)
1123 real(RP) :: w4(hydro_max)
1124 real(RP) :: w5(hydro_max)
1125 real(RP) :: w6(hydro_max)
1126 real(RP) :: w7(hydro_max)
1127 real(RP) :: w8(hydro_max)
1130 real(RP) :: ar_ice_fix = 0.7_rp
1131 real(RP) :: wcap1, wcap2
1133 logical :: flag_vent0(hydro_max), flag_vent1(hydro_max)
1135 integer :: iw, ia, ib
1138 namelist / param_atmos_phy_mp_sn14_init / &
1145 ntmax_phase_change, &
1148 namelist / param_atmos_phy_mp_sn14_particles / &
1149 a_m, b_m, alpha_v, beta_v, gamma_v, &
1150 alpha_vn, beta_vn, &
1151 a_area, b_area, cap, &
1153 opt_m96_column_ice, &
1156 real(RP),
parameter :: eps_gamma=1.e-30_rp
1160 alpha_v(:,:) = undef8
1161 beta_v(:,:) = undef8
1162 alpha_vn(:,:) = undef8
1163 beta_vn(:,:) = undef8
1181 coef_dave_n(:) = undef8
1182 coef_dave_l(:) = undef8
1186 coef_d2v(:) = undef8
1187 coef_md2v(:) = undef8
1191 coef_rea2(:) = undef8
1192 coef_rea3(:) = undef8
1195 coef_lambda(:) = undef8
1196 coef_vt0(:,:) = undef8
1197 coef_vt1(:,:) = undef8
1198 delta_b0(:) = undef8
1199 delta_b1(:) = undef8
1200 delta_ab0(:,:) = undef8
1201 delta_ab1(:,:) = undef8
1202 theta_b0(:) = undef8
1203 theta_b1(:) = undef8
1204 theta_ab0(:,:) = undef8
1205 theta_ab1(:,:) = undef8
1207 ah_vent(:,:) = undef8
1208 ah_vent0(:,:) = undef8
1209 ah_vent1(:,:) = undef8
1210 bh_vent(:,:) = undef8
1211 bh_vent0(:,:) = undef8
1212 bh_vent1(:,:) = undef8
1216 read(
io_fid_conf,nml=param_atmos_phy_mp_sn14_init,iostat=ierr)
1219 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.' 1220 elseif( ierr > 0 )
then 1221 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_init. Check!' 1224 log_nml(param_atmos_phy_mp_sn14_init)
1230 a_area(i_mp_qc) = pi/4.0_rp
1231 a_area(i_mp_qr) = pi/4.0_rp
1232 a_area(i_mp_qi) = 0.65_rp*1.e-4_rp*100.0_rp**(2.00_rp)
1233 a_area(i_mp_qs) = 0.2285_rp*1.e-4_rp*100.0_rp**(1.88_rp)
1234 a_area(i_mp_qg) = 0.50_rp*1.e-4_rp*100.0_rp**(2.0_rp)
1235 b_area(i_mp_qc) = 2.0_rp
1236 b_area(i_mp_qr) = 2.0_rp
1237 b_area(i_mp_qi) = 2.0_rp
1238 b_area(i_mp_qs) = 1.88_rp
1239 b_area(i_mp_qg) = 2.0_rp
1245 a_m(i_mp_qc) = 0.124_rp
1246 a_m(i_mp_qr) = 0.124_rp
1247 a_m(i_mp_qi) = 0.217_rp
1248 a_m(i_mp_qs) = 8.156_rp
1249 a_m(i_mp_qg) = 0.190_rp
1250 b_m(i_mp_qc) = 1.0_rp/3.0_rp
1251 b_m(i_mp_qr) = 1.0_rp/3.0_rp
1252 b_m(i_mp_qi) = 0.302_rp
1253 b_m(i_mp_qs) = 0.526_rp
1254 b_m(i_mp_qg) = 0.323_rp
1258 alpha_v(i_mp_qc,:)= 3.75e+5_rp
1259 alpha_v(i_mp_qr,:)= 159.0_rp
1260 alpha_v(i_mp_qi,:)= 317.0_rp
1261 alpha_v(i_mp_qs,:)= 27.70_rp
1262 alpha_v(i_mp_qg,:)= 40.0_rp
1263 beta_v(i_mp_qc,:) = 2.0_rp/3.0_rp
1264 beta_v(i_mp_qr,:) = 0.266_rp
1265 beta_v(i_mp_qi,:) = 0.363_rp
1266 beta_v(i_mp_qs,:) = 0.216_rp
1267 beta_v(i_mp_qg,:) = 0.230_rp
1268 gamma_v(i_mp_qc) = 1.0_rp
1270 gamma_v(i_mp_qr) = 1.0_rp/2.0_rp
1271 gamma_v(i_mp_qi) = 1.0_rp/2.0_rp
1272 gamma_v(i_mp_qs) = 1.0_rp/2.0_rp
1273 gamma_v(i_mp_qg) = 1.0_rp/2.0_rp
1281 nu(i_mp_qc) = 1.0_rp
1282 nu(i_mp_qr) = -1.0_rp/3.0_rp
1283 nu(i_mp_qi) = 1.0_rp
1284 nu(i_mp_qs) = 1.0_rp
1285 nu(i_mp_qg) = 1.0_rp
1287 mu(i_mp_qc) = 1.0_rp
1288 mu(i_mp_qr) = 1.0_rp/3.0_rp
1289 mu(i_mp_qi) = 1.0_rp/3.0_rp
1290 mu(i_mp_qs) = 1.0_rp/3.0_rp
1291 mu(i_mp_qg) = 1.0_rp/3.0_rp
1300 cap(i_mp_qc) = 2.0_rp
1301 cap(i_mp_qr) = 2.0_rp
1303 cap(i_mp_qs) = 2.0_rp
1304 cap(i_mp_qg) = 2.0_rp
1306 alpha_vn(:,:) = alpha_v(:,:)
1307 beta_vn(:,:) = beta_v(:,:)
1315 read(
io_fid_conf,nml=param_atmos_phy_mp_sn14_particles,iostat=ierr)
1318 log_info(
"ATMOS_PHY_MP_sn14_init",*)
'Not found namelist. Default used.' 1319 elseif( ierr > 0 )
then 1320 log_error(
"ATMOS_PHY_MP_sn14_init",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14_particles. Check!' 1323 log_nml(param_atmos_phy_mp_sn14_particles)
1327 if( opt_m96_ice )
then 1331 a_area(i_mp_qi) = 0.120284936_rp
1332 a_area(i_mp_qs) = 0.131488_rp
1333 a_area(i_mp_qg) = 0.5_rp
1334 b_area(i_mp_qi) = 1.850000_rp
1335 b_area(i_mp_qs) = 1.880000_rp
1336 b_area(i_mp_qg) = 2.0_rp
1337 a_m(i_mp_qi) = 1.23655360084766_rp
1338 a_m(i_mp_qs) = a_m(i_mp_qi)
1339 a_m(i_mp_qg) = 0.346111225718402_rp
1340 b_m(i_mp_qi) = 0.408329930583912_rp
1341 b_m(i_mp_qs) = b_m(i_mp_qi)
1342 b_m(i_mp_qg) = 0.357142857142857_rp
1344 if( opt_m96_column_ice )
then 1347 a_area(i_mp_qi)= (0.684_rp*1.e-4_rp)*10.0_rp**(2.0_rp*2.00_rp)
1348 b_area(i_mp_qi)= 2.0_rp
1349 a_m(i_mp_qi) = 0.19834046116844_rp
1350 b_m(i_mp_qi) = 0.343642611683849_rp
1353 wcap1 = sqrt(1.0_rp-ar_ice_fix**2)
1354 wcap2 = log( (1.0_rp+wcap1)/ar_ice_fix )
1355 cap(i_mp_qi) = 2.0_rp*wcap2/wcap1
1364 alpha_v(i_mp_qi,:) =(/ 5798.60107421875_rp, 167.347076416016_rp/)
1365 alpha_vn(i_mp_qi,:) =(/ 12408.177734375_rp, 421.799865722656_rp/)
1366 if( opt_m96_column_ice )
then 1367 alpha_v(i_mp_qi,:) = (/2901.0_rp, 32.20_rp/)
1368 alpha_vn(i_mp_qi,:) = (/9675.2_rp, 64.16_rp/)
1370 alpha_v(i_mp_qs,:) =(/ 15173.3916015625_rp, 305.678619384766_rp/)
1371 alpha_vn(i_mp_qs,:) =(/ 29257.1601562500_rp, 817.985717773438_rp/)
1372 alpha_v(i_mp_qg,:) =(/ 15481.6904296875_rp, 311.642242431641_rp/)
1373 alpha_vn(i_mp_qg,:) =(/ 27574.6562500000_rp, 697.536132812500_rp/)
1375 beta_v(i_mp_qi,:) =(/ 0.504873454570770_rp, 0.324817866086960_rp/)
1376 beta_vn(i_mp_qi,:) =(/ 0.548495233058929_rp, 0.385287821292877_rp/)
1377 if( opt_m96_column_ice )
then 1378 beta_v(i_mp_qi,:) =(/ 0.465552181005478_rp, 0.223826110363007_rp/)
1379 beta_vn(i_mp_qi,:) =(/ 0.530453503131866_rp, 0.273761242628098_rp/)
1381 beta_v(i_mp_qs,:) =(/ 0.528109610080719_rp, 0.329863965511322_rp/)
1382 beta_vn(i_mp_qs,:) =(/ 0.567154467105865_rp, 0.393876969814301_rp/)
1383 beta_v(i_mp_qg,:) =(/ 0.534656763076782_rp, 0.330253750085831_rp/)
1384 beta_vn(i_mp_qg,:) =(/ 0.570551633834839_rp, 0.387124240398407_rp/)
1388 ax_area(:) = a_area(:)*a_m(:)**b_area(:)
1389 bx_area(:) = b_area(:)*b_m(:)
1393 a_rea(:) = sqrt(ax_area(:)/pi)
1394 b_rea(:) = bx_area(:)/2.0_rp
1395 a_rea2(:) = a_rea(:)**2
1396 b_rea2(:) = b_rea(:)*2.0_rp
1397 a_rea3(:) = a_rea(:)**3
1398 b_rea3(:) = b_rea(:)*3.0_rp
1400 a_d2vt(:)=alpha_v(:,2)*(1.0_rp/alpha_v(:,2))**(beta_v(:,2)/b_m(:))
1401 b_d2vt(:)=(beta_v(:,2)/b_m(:))
1421 w1(iw) = gammafunc( (n+nu(iw)+1.0_rp)/mu(iw) )
1422 w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1423 w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1424 coef_m2(iw) = w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**n
1426 w4(iw) = gammafunc( (b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1427 coef_d(iw) = a_m(iw) * w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**b_m(iw)
1428 w5(iw) = gammafunc( (2.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
1429 w6(iw) = gammafunc( (3.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
1430 coef_d2v(iw) = a_m(iw) * w6(iw)/w5(iw)* ( w2(iw)/w3(iw) )**b_m(iw)
1431 coef_md2v(iw)= w5(iw)/w2(iw)* ( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw)+beta_v(iw,2))
1433 w7(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1434 coef_d3(iw) = a_m(iw)**3 * w7(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
1435 w8(iw) = gammafunc( (6.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1436 coef_d6(iw) = a_m(iw)**6 * w8(iw)/w2(iw)*( w2(iw)/w3(iw) )**(6.0_rp*b_m(iw))
1439 coef_deplc = coef_d(i_mp_qc)/a_m(i_mp_qc)
1445 w1(iw) = gammafunc( (2.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1446 w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1447 w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1449 w4(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1451 coef_r2(iw)=w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw))
1452 coef_r3(iw)=w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
1453 coef_re(iw)=coef_r3(iw)/coef_r2(iw)
1460 w1(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1461 w2(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1462 w3(iw) = gammafunc( (b_rea2(iw)+nu(iw)+1.0_rp)/mu(iw) )
1463 w4(iw) = gammafunc( (b_rea3(iw)+nu(iw)+1.0_rp)/mu(iw) )
1465 coef_rea2(iw) = w3(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea2(iw)
1466 coef_rea3(iw) = w4(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea3(iw)
1472 w1(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1473 w2(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1474 coef_a(iw) = mu(iw)/w1(iw)
1476 coef_lambda(iw) = (w1(iw)/w2(iw))**(-mu(iw))
1484 w1(iw) = gammafunc( (beta_vn(iw,ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
1485 w2(iw) = gammafunc( ( nu(iw) + 1.0_rp + n)/mu(iw) )
1486 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1487 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1489 coef_vt0(iw,ia)=alpha_vn(iw,ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_vn(iw,ia)
1491 w1(iw) = gammafunc( (beta_v(iw,ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
1492 w2(iw) = gammafunc( ( nu(iw) + 1.0_rp + n)/mu(iw) )
1494 coef_vt1(iw,ia)=alpha_v(iw,ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_v(iw,ia)
1499 w1(iw) = gammafunc( ( b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
1500 w2(iw) = gammafunc( (1.0_rp + b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
1501 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1502 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1503 coef_dave_n(iw) = (w1(iw)/w3(iw))*(w3(iw)/w4(iw))** b_m(iw)
1504 coef_dave_l(iw) = (w2(iw)/w3(iw))*(w3(iw)/w4(iw))**(1.0_rp+b_m(iw))
1508 ah_vent(i_mp_qc,1:2) = (/1.0000_rp,1.0000_rp/)
1509 ah_vent(i_mp_qr,1:2) = (/1.0000_rp,0.780_rp/)
1510 ah_vent(i_mp_qi,1:2) = (/1.0000_rp,0.860_rp/)
1511 ah_vent(i_mp_qs,1:2) = (/1.0000_rp,0.780_rp/)
1512 ah_vent(i_mp_qg,1:2) = (/1.0000_rp,0.780_rp/)
1513 bh_vent(i_mp_qc,1:2) = (/0.0000_rp,0.0000_rp/)
1514 bh_vent(i_mp_qr,1:2) = (/0.108_rp,0.308_rp/)
1515 bh_vent(i_mp_qi,1:2) = (/0.140_rp,0.280_rp/)
1516 bh_vent(i_mp_qs,1:2) = (/0.108_rp,0.308_rp/)
1517 bh_vent(i_mp_qg,1:2) = (/0.108_rp,0.308_rp/)
1521 if( (nu(iw) + b_m(iw) + n) > eps_gamma )
then 1522 w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
1523 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1524 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1525 ah_vent0(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1526 ah_vent0(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1527 flag_vent0(iw)=.true.
1529 ah_vent0(iw,1)= 1.0_rp
1530 ah_vent0(iw,2)= 1.0_rp
1531 flag_vent0(iw)=.false.
1534 if( (nu(iw) + b_m(iw) + n) > eps_gamma )
then 1535 w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
1536 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1537 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1538 ah_vent1(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1539 ah_vent1(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1540 flag_vent1(iw)=.true.
1542 ah_vent1(iw,1)= 1.0_rp
1543 ah_vent1(iw,2)= 1.0_rp
1544 flag_vent1(iw)=.true.
1549 if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )
then 1550 flag_vent0(iw)=.false.
1552 if(flag_vent0(iw))
then 1553 w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
1554 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1555 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1557 w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
1558 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)
1559 w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
1560 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)
1562 bh_vent0(iw,1) = 0.0_rp
1563 bh_vent0(iw,2) = 0.0_rp
1567 if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )
then 1568 flag_vent1(iw)=.false.
1570 if(flag_vent1(iw))
then 1571 w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
1572 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1573 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1575 w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
1576 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)
1578 w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
1579 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)
1581 bh_vent1(iw,1) = 0.0_rp
1582 bh_vent1(iw,2) = 0.0_rp
1591 w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1592 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1593 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1594 delta_b0(iw) = w1(iw)/w2(iw) &
1595 *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
1597 w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1598 delta_b1(iw) = w1(iw)/w2(iw) &
1599 *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
1605 w1(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1606 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1607 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1608 w4(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1610 w5(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1618 delta_ab0(ia,ib) = 2.0_rp*(w1(ib)/w2(ib))*(w4(ia)/w2(ia)) &
1619 * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
1620 * ( w2(ia)/w3(ia) )**(b_rea(ia) )
1622 delta_ab1(ia,ib) = 2.0_rp*(w5(ib)/w2(ib))*(w4(ia)/w2(ia)) &
1623 * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
1624 * ( w2(ia)/w3(ia) )**(b_rea(ia) )
1632 w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1633 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1634 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1635 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1636 theta_b0(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
1638 w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1639 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1640 theta_b1(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
1647 w1(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1648 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1649 w3(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1650 w4(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1652 w5(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1653 w6(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1655 w7(iw) = gammafunc( (beta_v(iw,2) + b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1656 w8(iw) = gammafunc( ( b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1661 theta_ab0(ia,ib) = 2.0_rp * (w1(ib)/w2(ib))*(w3(ia)/w4(ia)) &
1662 * (w5(ia)/w6(ia))**beta_v(ia,2) &
1663 * (w5(ib)/w6(ib))**beta_v(ib,2)
1664 theta_ab1(ia,ib) = 2.0_rp * (w7(ib)/w8(ib))*(w3(ia)/w4(ia)) &
1665 * (w5(ia)/w6(ia))**beta_v(ia,2) &
1666 * (w5(ib)/w6(ib))**beta_v(ib,2)
1670 log_info(
"ATMOS_PHY_MP_sn14_init",
'(100a16)')
"LABEL ",wlabel(:)
1671 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"capacity ",cap(:)
1672 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_m2 ",coef_m2(:)
1673 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d ",coef_d(:)
1675 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d3 ",coef_d3(:)
1676 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d6 ",coef_d6(:)
1677 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_d2v ",coef_d2v(:)
1678 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_md2v ",coef_md2v(:)
1679 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_d2vt ",a_d2vt(:)
1680 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_d2vt ",b_d2vt(:)
1682 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_r2 ",coef_r2(:)
1683 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_r3 ",coef_r3(:)
1684 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_re ",coef_re(:)
1686 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_area ",a_area(:)
1687 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_area ",b_area(:)
1688 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ax_area ",ax_area(:)
1689 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bx_area ",bx_area(:)
1690 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_rea ",a_rea(:)
1691 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_rea ",b_rea(:)
1692 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"a_rea3 ",a_rea3(:)
1693 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"b_rea3 ",b_rea3(:)
1695 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_rea2 ",coef_rea2(:)
1696 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_rea3 ",coef_rea3(:)
1697 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_vt0 ",coef_vt0(:,1)
1698 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_vt1 ",coef_vt1(:,1)
1699 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_A ",coef_a(:)
1700 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"coef_lambda ",coef_lambda(:)
1702 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent0 sml",ah_vent0(:,1)
1703 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent0 lrg",ah_vent0(:,2)
1704 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent1 sml",ah_vent1(:,1)
1705 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"ah_vent1 lrg",ah_vent1(:,2)
1706 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent0 sml",bh_vent0(:,1)
1707 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent0 lrg",bh_vent0(:,2)
1708 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent1 sml",bh_vent1(:,1)
1709 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"bh_vent1 lrg",bh_vent1(:,2)
1711 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"delta_b0 ",delta_b0(:)
1712 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"delta_b1 ",delta_b1(:)
1713 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"theta_b0 ",theta_b0(:)
1714 log_info(
"ATMOS_PHY_MP_sn14_init",
'(a,100ES16.6)')
"theta_b1 ",theta_b1(:)
1717 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)
1720 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)
1723 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)
1726 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)
1730 end subroutine mp_sn14_init
1732 subroutine mp_sn14 ( &
1733 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1760 moist_psat_liq => atmos_saturation_psat_liq, &
1761 moist_psat_ice => atmos_saturation_psat_ice
1764 integer,
intent(in) :: KA, KS, KE
1765 integer,
intent(in) :: IA, IS, IE
1766 integer,
intent(in) :: JA, JS, JE
1768 real(RP),
intent(in) :: DENS (ka,ia,ja)
1769 real(RP),
intent(in) :: W (ka,ia,ja)
1770 real(RP),
intent(in) :: QTRC (ka,ia,ja,
qa_mp)
1771 real(RP),
intent(in) :: PRES0 (ka,ia,ja)
1772 real(RP),
intent(in) :: TEMP0 (ka,ia,ja)
1773 real(RP),
intent(in) :: Qdry (ka,ia,ja)
1774 real(RP),
intent(in) :: CPtot0(ka, ia, ja)
1775 real(RP),
intent(in) :: CVtot0(ka, ia, ja)
1776 real(RP),
intent(in) :: CCN (ka,ia,ja)
1777 real(RP),
intent(in) :: dt
1778 real(RP),
intent(in) :: cz( ka,ia,ja)
1779 real(RP),
intent(in) :: fz(0:ka,ia,ja)
1781 real(RP),
intent(out) :: RHOQ_t(ka, ia, ja,
qa_mp)
1782 real(RP),
intent(out) :: RHOE_t(ka, ia, ja)
1783 real(RP),
intent(out) :: CPtot_t(ka, ia, ja)
1784 real(RP),
intent(out) :: CVtot_t(ka, ia, ja)
1786 real(RP),
intent(out) :: EVAPORATE(ka,ia,ja)
1788 real(RP) :: pres(ka,ia,ja)
1789 real(RP) :: temp(ka,ia,ja)
1790 real(RP) :: cva(ka,ia,ja)
1791 real(RP) :: cpa(ka,ia,ja)
1792 real(RP) :: RHOQ0_t(ka, ia, ja,
qa_mp)
1793 real(RP) :: RHOE0_t(ka, ia, ja)
1794 real(RP) :: CPtot0_t(ka, ia, ja)
1795 real(RP) :: CVtot0_t(ka, ia, ja)
1799 real(RP) :: rhoe(ka,ia,ja)
1800 real(RP) :: rhoq(i_qv:i_ng,ka,ia,ja)
1801 real(RP) :: rhoq2(i_qv:i_ng,ka,ia,ja)
1803 real(RP) :: xq(hydro_max,ka,ia,ja)
1805 real(RP) :: dq_xa(hydro_max,ka,ia,ja)
1806 real(RP) :: vt_xa(hydro_max,2,ka,ia,ja)
1808 real(RP) :: wtemp(ka,ia,ja)
1809 real(RP) :: esw(ka,ia,ja)
1810 real(RP) :: esi(ka,ia,ja)
1813 real(RP) :: rho_fac_q(hydro_max,ka,ia,ja)
1816 real(RP) :: drhogqc, drhognc
1817 real(RP) :: drhogqr, drhognr
1818 real(RP) :: drhogqi, drhogni
1819 real(RP) :: drhogqs, drhogns
1820 real(RP) :: drhogqg, drhogng
1823 real(RP) :: PQ(pq_max,ka,ia,ja)
1825 real(RP) :: wrm_dqc, wrm_dnc
1826 real(RP) :: wrm_dqr, wrm_dnr
1829 real(RP) :: Pac(pac_max,ka,ia,ja)
1831 real(RP) :: gc_dqc, gc_dnc
1832 real(RP) :: sc_dqc, sc_dnc
1833 real(RP) :: ic_dqc, ic_dnc
1834 real(RP) :: rg_dqg, rg_dng
1835 real(RP) :: rg_dqr, rg_dnr
1836 real(RP) :: rs_dqr, rs_dnr, rs_dqs, rs_dns
1837 real(RP) :: ri_dqr, ri_dnr
1838 real(RP) :: ri_dqi, ri_dni
1839 real(RP) :: ii_dqi, ii_dni
1840 real(RP) :: is_dqi, is_dni, ss_dns
1841 real(RP) :: gs_dqs, gs_dns, gg_dng
1843 real(RP) :: clp_dqc, clp_dnc, clm_dqc, clm_dnc
1844 real(RP) :: clp_dqr, clp_dnr, clm_dqr, clm_dnr
1845 real(RP) :: clp_dqi, clp_dni, clm_dqi, clm_dni
1846 real(RP) :: clp_dqs, clp_dns, clm_dqs, clm_dns
1847 real(RP) :: clp_dqg, clp_dng, clm_dqg, clm_dng
1848 real(RP) :: fac1, fac3, fac4, fac6, fac7, fac9, fac10
1850 real(RP) :: pco_dqi, pco_dni
1851 real(RP) :: pco_dqs, pco_dns
1852 real(RP) :: pco_dqg, pco_dng
1854 real(RP) :: eml_dqc, eml_dnc
1855 real(RP) :: eml_dqr, eml_dnr
1856 real(RP) :: eml_dqi, eml_dni
1857 real(RP) :: eml_dqs, eml_dns
1858 real(RP) :: eml_dqg, eml_dng
1860 real(RP) :: spl_dqi, spl_dni
1861 real(RP) :: spl_dqg, spl_dqs
1863 real(RP) :: rrho(ka,ia,ja)
1868 real(RP) :: dTdt_equiv_d(ka,ia,ja)
1875 real(RP) :: sl_PLCdep(ia,ja)
1876 real(RP) :: sl_PLRdep(ia,ja), sl_PNRdep(ia,ja)
1878 real(RP) :: qke_d(ka,ia,ja)
1880 real(RP),
parameter :: eps = 1.e-19_rp
1881 real(RP),
parameter :: eps_qv = 1.e-19_rp
1882 real(RP),
parameter :: eps_rhoge = 1.e-19_rp
1883 real(RP),
parameter :: eps_rhog = 1.e-19_rp
1888 integer :: k, i, j, iq
1890 real(RP) :: dqv, dqc, dqr, dqi, dqs, dqg, dcv, dcp
1895 rhoq_t(:,:,:,:) = 0.0_rp
1896 rhoe_t(:,:,:) = 0.0_rp
1897 cptot_t(:,:,:) = 0.0_rp
1898 cvtot_t(:,:,:) = 0.0_rp
1904 cpa(k,i,j) = cptot0(k,i,j)
1905 cva(k,i,j) = cvtot0(k,i,j)
1906 pres(k,i,j) = pres0(k,i,j)
1907 temp(k,i,j) = temp0(k,i,j)
1932 rhoq(iq,k,i,j) = dens(k,i,j) * qtrc(k,i,j,iq)
1934 rhoq2(i_qv,k,i,j) = dens(k,i,j)*qtrc(k,i,j,i_qv)
1935 rhoq2(i_ni,k,i,j) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,i_ni) )
1936 rhoq2(i_nc,k,i,j) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,i_nc) )
1944 rrho(k,i,j) = 1.0_rp / dens(k,i,j)
1945 rhoe(k,i,j) = dens(k,i,j) * temp(k,i,j) * cva(k,i,j)
1946 wtemp(k,i,j) = max(temp(k,i,j), tem_min)
1951 if( opt_debug_tem )
call debug_tem_kij( ka, ks, ke, ia, is, ie, ja, js, je, &
1952 1, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )
1957 rho_fac = rho_0 / max(dens(k,i,j),rho_min)
1958 rho_fac_q(i_mp_qc,k,i,j) = rho_fac**gamma_v(i_mp_qc)
1959 rho_fac_q(i_mp_qr,k,i,j) = rho_fac**gamma_v(i_mp_qr)
1960 rho_fac_q(i_mp_qi,k,i,j) = (pres(k,i,j)/pre0_vt)**a_pre0_vt * (temp(k,i,j)/tem0_vt)**a_tem0_vt
1961 rho_fac_q(i_mp_qs,k,i,j) = rho_fac_q(i_mp_qi,k,i,j)
1962 rho_fac_q(i_mp_qg,k,i,j) = rho_fac_q(i_mp_qi,k,i,j)
1970 sl_plcdep(i,j) = 0.0_rp
1971 sl_plrdep(i,j) = 0.0_rp
1972 sl_pnrdep(i,j) = 0.0_rp
1980 qke_d(k,i,j) = 0.0_rp
1989 dtdt_equiv_d(k,i,j) = 0.0_rp
1999 ka, ks, ke, ia, is, ie, ja, js, je, &
2009 dtdt_equiv_d(:,:,:), &
2020 drhogqc = dt * pq(i_lcccn,k,i,j)
2021 drhognc = dt * pq(i_ncccn,k,i,j)
2022 drhogqi = dt * pq(i_liccn,k,i,j)
2023 drhogni = dt * pq(i_niccn,k,i,j)
2024 drhogqv = max( -rhoq(i_qv,k,i,j), -drhogqc-drhogqi )
2025 fac1 = drhogqv / min( -drhogqc-drhogqi, -eps )
2027 dqv = drhogqv * rrho(k,i,j)
2028 dqc = drhogqc*fac1 * rrho(k,i,j)
2029 dqi = drhogqi*fac1 * rrho(k,i,j)
2034 rhoq0_t(k,i,j,i_qv) = drhogqv / dt
2035 rhoq0_t(k,i,j,i_qc) = drhogqc*fac1 / dt
2036 rhoq0_t(k,i,j,i_qi) = drhogqi*fac1 / dt
2037 rhoq0_t(k,i,j,i_nc) = drhognc / dt
2038 rhoq0_t(k,i,j,i_ni) = drhogni / dt
2040 rhoe0_t(k,i,j) = - lhv * drhogqv / dt + lhf * drhogqi*fac1 / dt
2042 cvtot0_t(k,i,j) = dcv/dt
2043 cptot0_t(k,i,j) = dcp/dt
2054 rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe0_t(k,i,j)
2055 cvtot_t(k,i,j) = cvtot_t(k,i,j) + cvtot0_t(k,i,j)
2056 cptot_t(k,i,j) = cptot_t(k,i,j) + cptot0_t(k,i,j)
2067 rhoq(i_qv,k,i,j) = rhoq(i_qv,k,i,j) + rhoq0_t(k,i,j,i_qv)*dt
2068 rhoq(i_qc,k,i,j) = max(0.0_rp, rhoq(i_qc,k,i,j) + rhoq0_t(k,i,j,i_qc)*dt )
2069 rhoq(i_qi,k,i,j) = max(0.0_rp, rhoq(i_qi,k,i,j) + rhoq0_t(k,i,j,i_qi)*dt )
2070 rhoq(i_nc,k,i,j) = max(0.0_rp, rhoq(i_nc,k,i,j) + rhoq0_t(k,i,j,i_nc)*dt )
2071 rhoq(i_ni,k,i,j) = max(0.0_rp, rhoq(i_ni,k,i,j) + rhoq0_t(k,i,j,i_ni)*dt )
2074 rhoq(i_nc,k,i,j) = min( rhoq(i_nc,k,i,j), nc_uplim_d(1,i,j) )
2076 rhoe(k,i,j) = rhoe(k,i,j) + rhoe0_t(k,i,j)*dt
2077 cva(k,i,j) = cva(k,i,j) + cvtot0_t(k,i,j)*dt
2078 cpa(k,i,j) = cpa(k,i,j) + cptot0_t(k,i,j)*dt
2080 temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cva(k,i,j) )
2081 pres(k,i,j) = dens(k,i,j) * (cpa(k,i,j)-cva(k,i,j)) * temp(k,i,j)
2082 wtemp(k,i,j) = max( temp(k,i,j), tem_min )
2088 if( opt_debug_tem )
call debug_tem_kij( ka, ks, ke, ia, is, ie, ja, js, je, &
2089 2, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )
2103 rhoq2(i_qr,k,i,j) = rhoq(i_qr,k,i,j)
2104 rhoq2(i_nr,k,i,j) = rhoq(i_nr,k,i,j)
2105 xq(i_mp_qr,k,i,j) = max(xr_min, min(xr_max, rhoq2(i_qr,k,i,j)/(rhoq2(i_nr,k,i,j)+nr_min) ))
2107 dq_xa(i_mp_qr,k,i,j) = a_m(i_mp_qr)*xq(i_mp_qr,k,i,j)**b_m(i_mp_qr)
2108 vt_xa(i_mp_qr,1,k,i,j) = alpha_v(i_mp_qr,1)*(xq(i_mp_qr,k,i,j)**beta_v(i_mp_qr,1))*rho_fac_q(i_mp_qr,k,i,j)
2109 vt_xa(i_mp_qr,2,k,i,j) = vt_xa(i_mp_qr,1,k,i,j)
2113 rhoq2(i_qv,k,i,j) = rhoq(i_qv,k,i,j)
2114 rhoq2(i_qc,k,i,j) = rhoq(i_qc,k,i,j)
2115 rhoq2(i_qi,k,i,j) = rhoq(i_qi,k,i,j)
2116 rhoq2(i_qs,k,i,j) = rhoq(i_qs,k,i,j)
2117 rhoq2(i_qg,k,i,j) = rhoq(i_qg,k,i,j)
2119 rhoq2(i_nc,k,i,j) = rhoq(i_nc,k,i,j)
2120 rhoq2(i_ni,k,i,j) = rhoq(i_ni,k,i,j)
2121 rhoq2(i_ns,k,i,j) = rhoq(i_ns,k,i,j)
2122 rhoq2(i_ng,k,i,j) = rhoq(i_ng,k,i,j)
2127 xq(i_mp_qc,k,i,j) = min(xc_max, max(xc_min, rhoq2(i_qc,k,i,j)/(rhoq2(i_nc,k,i,j)+nc_min) ))
2128 xq(i_mp_qi,k,i,j) = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ))
2129 xq(i_mp_qs,k,i,j) = min(xs_max, max(xs_min, rhoq2(i_qs,k,i,j)/(rhoq2(i_ns,k,i,j)+ns_min) ))
2130 xq(i_mp_qg,k,i,j) = min(xg_max, max(xg_min, rhoq2(i_qg,k,i,j)/(rhoq2(i_ng,k,i,j)+ng_min) ))
2133 dq_xa(i_mp_qc,k,i,j) = a_m(i_mp_qc)*xq(i_mp_qc,k,i,j)**b_m(i_mp_qc)
2134 dq_xa(i_mp_qi,k,i,j) = a_m(i_mp_qi)*xq(i_mp_qi,k,i,j)**b_m(i_mp_qi)
2135 dq_xa(i_mp_qs,k,i,j) = a_m(i_mp_qs)*xq(i_mp_qs,k,i,j)**b_m(i_mp_qs)
2136 dq_xa(i_mp_qg,k,i,j) = a_m(i_mp_qg)*xq(i_mp_qg,k,i,j)**b_m(i_mp_qg)
2139 vt_xa(i_mp_qc,1,k,i,j) = alpha_v(i_mp_qc,1)*(xq(i_mp_qc,k,i,j)**beta_v(i_mp_qc,1))*rho_fac_q(i_mp_qc,k,i,j)
2140 vt_xa(i_mp_qi,1,k,i,j) = alpha_v(i_mp_qi,1)*(xq(i_mp_qi,k,i,j)**beta_v(i_mp_qi,1))*rho_fac_q(i_mp_qi,k,i,j)
2141 vt_xa(i_mp_qs,1,k,i,j) = alpha_v(i_mp_qs,1)*(xq(i_mp_qs,k,i,j)**beta_v(i_mp_qs,1))*rho_fac_q(i_mp_qs,k,i,j)
2142 vt_xa(i_mp_qg,1,k,i,j) = alpha_v(i_mp_qg,1)*(xq(i_mp_qg,k,i,j)**beta_v(i_mp_qg,1))*rho_fac_q(i_mp_qg,k,i,j)
2143 vt_xa(i_mp_qc,2,k,i,j) = alpha_v(i_mp_qc,2)*(xq(i_mp_qc,k,i,j)**beta_v(i_mp_qc,2))*rho_fac_q(i_mp_qc,k,i,j)
2144 vt_xa(i_mp_qi,2,k,i,j) = alpha_v(i_mp_qi,2)*(xq(i_mp_qi,k,i,j)**beta_v(i_mp_qi,2))*rho_fac_q(i_mp_qi,k,i,j)
2145 vt_xa(i_mp_qs,2,k,i,j) = alpha_v(i_mp_qs,2)*(xq(i_mp_qs,k,i,j)**beta_v(i_mp_qs,2))*rho_fac_q(i_mp_qs,k,i,j)
2146 vt_xa(i_mp_qg,2,k,i,j) = alpha_v(i_mp_qg,2)*(xq(i_mp_qg,k,i,j)**beta_v(i_mp_qg,2))*rho_fac_q(i_mp_qg,k,i,j)
2152 call moist_psat_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
2153 wtemp(:,:,:), esw(:,:,:) )
2154 call moist_psat_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
2155 wtemp(:,:,:), esi(:,:,:) )
2158 ka, ks, ke, ia, is, ie, ja, js, je, &
2165 call dep_vapor_melt_ice_kij( &
2166 ka, ks, ke, ia, is, ie, ja, js, je, &
2182 ka, ks, ke, ia, is, ie, ja, js, je, &
2184 ntmax_phase_change, &
2187 cz(:,:,:), fz(:,:,:), &
2189 dtdt_equiv_d(:,:,:), &
2214 rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe0_t(k,i,j)
2215 cvtot_t(k,i,j) = cvtot_t(k,i,j) + cvtot0_t(k,i,j)
2216 cptot_t(k,i,j) = cptot_t(k,i,j) + cptot0_t(k,i,j)
2228 rhoq(iq,k,i,j) = max(0.0_rp, rhoq(iq,k,i,j) + rhoq0_t(k,i,j,iq)*dt )
2239 rhoe(k,i,j) = rhoe(k,i,j) + rhoe0_t(k,i,j)*dt
2240 cva(k,i,j) = cva(k,i,j) + cvtot0_t(k,i,j)*dt
2241 cpa(k,i,j) = cpa(k,i,j) + cptot0_t(k,i,j)*dt
2242 temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cva(k,i,j) )
2243 pres(k,i,j) = dens(k,i,j) * (cpa(k,i,j) - cva(k,i,j)) * temp(k,i,j)
2250 if( opt_debug_tem )
call debug_tem_kij( ka, ks, ke, ia, is, ie, ja, js, je, &
2251 3, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )
2268 rhoq2(i_qv,k,i,j) = rhoq(i_qv,k,i,j)
2269 rhoq2(i_qc,k,i,j) = rhoq(i_qc,k,i,j)
2270 rhoq2(i_qr,k,i,j) = rhoq(i_qr,k,i,j)
2271 rhoq2(i_qi,k,i,j) = rhoq(i_qi,k,i,j)
2272 rhoq2(i_qs,k,i,j) = rhoq(i_qs,k,i,j)
2273 rhoq2(i_qg,k,i,j) = rhoq(i_qg,k,i,j)
2275 rhoq2(i_nc,k,i,j) = rhoq(i_nc,k,i,j)
2276 rhoq2(i_nr,k,i,j) = rhoq(i_nr,k,i,j)
2277 rhoq2(i_ni,k,i,j) = rhoq(i_ni,k,i,j)
2278 rhoq2(i_ns,k,i,j) = rhoq(i_ns,k,i,j)
2279 rhoq2(i_ng,k,i,j) = rhoq(i_ng,k,i,j)
2282 xq(i_mp_qc,k,i,j) = min(xc_max, max(xc_min, rhoq2(i_qc,k,i,j)/(rhoq2(i_nc,k,i,j)+nc_min) ) )
2283 xq(i_mp_qr,k,i,j) = min(xr_max, max(xr_min, rhoq2(i_qr,k,i,j)/(rhoq2(i_nr,k,i,j)+nr_min) ) )
2284 xq(i_mp_qi,k,i,j) = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ) )
2285 xq(i_mp_qs,k,i,j) = min(xs_max, max(xs_min, rhoq2(i_qs,k,i,j)/(rhoq2(i_ns,k,i,j)+ns_min) ) )
2286 xq(i_mp_qg,k,i,j) = min(xg_max, max(xg_min, rhoq2(i_qg,k,i,j)/(rhoq2(i_ng,k,i,j)+ng_min) ) )
2289 dq_xa(i_mp_qc,k,i,j) = 2.0_rp*a_rea(i_mp_qc)*xq(i_mp_qc,k,i,j)**b_rea(i_mp_qc)
2290 dq_xa(i_mp_qr,k,i,j) = 2.0_rp*a_rea(i_mp_qr)*xq(i_mp_qr,k,i,j)**b_rea(i_mp_qr)
2291 dq_xa(i_mp_qi,k,i,j) = 2.0_rp*a_rea(i_mp_qi)*xq(i_mp_qi,k,i,j)**b_rea(i_mp_qi)
2292 dq_xa(i_mp_qs,k,i,j) = 2.0_rp*a_rea(i_mp_qs)*xq(i_mp_qs,k,i,j)**b_rea(i_mp_qs)
2293 dq_xa(i_mp_qg,k,i,j) = 2.0_rp*a_rea(i_mp_qg)*xq(i_mp_qg,k,i,j)**b_rea(i_mp_qg)
2297 vt_xa(i_mp_qc,2,k,i,j) = alpha_v(i_mp_qc,2)*(xq(i_mp_qc,k,i,j)**beta_v(i_mp_qc,2))*rho_fac_q(i_mp_qc,k,i,j)
2298 vt_xa(i_mp_qr,2,k,i,j) = alpha_v(i_mp_qr,2)*(xq(i_mp_qr,k,i,j)**beta_v(i_mp_qr,2))*rho_fac_q(i_mp_qr,k,i,j)
2299 vt_xa(i_mp_qi,2,k,i,j) = alpha_v(i_mp_qi,2)*(xq(i_mp_qi,k,i,j)**beta_v(i_mp_qi,2))*rho_fac_q(i_mp_qi,k,i,j)
2300 vt_xa(i_mp_qs,2,k,i,j) = alpha_v(i_mp_qs,2)*(xq(i_mp_qs,k,i,j)**beta_v(i_mp_qs,2))*rho_fac_q(i_mp_qs,k,i,j)
2301 vt_xa(i_mp_qg,2,k,i,j) = alpha_v(i_mp_qg,2)*(xq(i_mp_qg,k,i,j)**beta_v(i_mp_qg,2))*rho_fac_q(i_mp_qg,k,i,j)
2308 if ( mp_doautoconversion )
then 2310 ka, ks, ke, ia, is, ie, ja, js, je, &
2322 pq(i_lcaut,k,i,j) = 0.0_rp
2323 pq(i_ncaut,k,i,j) = 0.0_rp
2324 pq(i_nraut,k,i,j) = 0.0_rp
2325 pq(i_lcacc,k,i,j) = 0.0_rp
2326 pq(i_ncacc,k,i,j) = 0.0_rp
2327 pq(i_nrslc,k,i,j) = 0.0_rp
2328 pq(i_nrbrk,k,i,j) = 0.0_rp
2336 ka, ks, ke, ia, is, ie, ja, js, je, &
2346 ka, ks, ke, ia, is, ie, ja, js, je, &
2357 profile_start(
"sn14_update_rhoq")
2363 wrm_dqc = max( dt*( pq(i_lcaut,k,i,j)+pq(i_lcacc,k,i,j) ), -rhoq2(i_qc,k,i,j) )
2364 wrm_dnc = max( dt*( pq(i_ncaut,k,i,j)+pq(i_ncacc,k,i,j) ), -rhoq2(i_nc,k,i,j) )
2365 wrm_dnr = max( dt*( pq(i_nraut,k,i,j)+pq(i_nrslc,k,i,j)+pq(i_nrbrk,k,i,j) ), -rhoq2(i_nr,k,i,j) )
2374 gc_dqc = max( dt*pac(i_lgaclc2lg,k,i,j) , min(0.0_rp, -rhoq2(i_qc,k,i,j)-wrm_dqc ))
2375 sc_dqc = max( dt*pac(i_lsaclc2ls,k,i,j) , min(0.0_rp, -rhoq2(i_qc,k,i,j)-wrm_dqc-gc_dqc ))
2376 ic_dqc = max( dt*pac(i_liaclc2li,k,i,j) , min(0.0_rp, -rhoq2(i_qc,k,i,j)-wrm_dqc-gc_dqc-sc_dqc ))
2378 gc_dnc = max( dt*pac(i_ngacnc2ng,k,i,j) , min(0.0_rp, -rhoq2(i_nc,k,i,j)-wrm_dnc ))
2379 sc_dnc = max( dt*pac(i_nsacnc2ns,k,i,j) , min(0.0_rp, -rhoq2(i_nc,k,i,j)-wrm_dnc-gc_dnc ))
2380 ic_dnc = max( dt*pac(i_niacnc2ni,k,i,j) , min(0.0_rp, -rhoq2(i_nc,k,i,j)-wrm_dnc-gc_dnc-sc_dnc ))
2383 sw = sign(0.5_rp, t00-temp(k,i,j)) + 0.5_rp
2384 rg_dqr = max( dt*pac(i_lraclg2lg, k,i,j), min(0.0_rp, -rhoq2(i_qr,k,i,j)-wrm_dqr )) * sw
2385 rg_dqg = max( dt*pac(i_lraclg2lg, k,i,j), min(0.0_rp, -rhoq2(i_qg,k,i,j) )) * ( 1.0_rp - sw )
2386 rs_dqr = max( dt*pac(i_lracls2lg_r,k,i,j), min(0.0_rp, -rhoq2(i_qr,k,i,j)-wrm_dqr-rg_dqr )) * sw
2387 ri_dqr = max( dt*pac(i_lracli2lg_r,k,i,j), min(0.0_rp, -rhoq2(i_qr,k,i,j)-wrm_dqr-rg_dqr-rs_dqr )) * sw
2389 rg_dnr = max( dt*pac(i_nracng2ng, k,i,j), min(0.0_rp, -rhoq2(i_nr,k,i,j)-wrm_dnr )) * sw
2390 rg_dng = max( dt*pac(i_nracng2ng, k,i,j), min(0.0_rp, -rhoq2(i_ng,k,i,j) )) * ( 1.0_rp - sw )
2391 rs_dnr = max( dt*pac(i_nracns2ng_r,k,i,j), min(0.0_rp, -rhoq2(i_nr,k,i,j)-wrm_dnr-rg_dnr )) * sw
2392 ri_dnr = max( dt*pac(i_nracni2ng_r,k,i,j), min(0.0_rp, -rhoq2(i_nr,k,i,j)-wrm_dnr-rg_dnr-rs_dnr )) * sw
2395 fac1 = (ri_dqr-eps)/ (dt*pac(i_lracli2lg_r,k,i,j)-eps)
2396 ri_dqi = max( dt*pac(i_lracli2lg_i,k,i,j)*fac1, min(0.0_rp, -rhoq2(i_qi,k,i,j)+ic_dqc ))
2397 ii_dqi = max( dt*pac(i_liacli2ls,k,i,j) , min(0.0_rp, -rhoq2(i_qi,k,i,j)+ic_dqc-ri_dqi ))
2398 is_dqi = max( dt*pac(i_liacls2ls,k,i,j) , min(0.0_rp, -rhoq2(i_qi,k,i,j)+ic_dqc-ri_dqi-ii_dqi ))
2400 fac4 = (ri_dnr-eps)/ (dt*pac(i_nracni2ng_r,k,i,j)-eps)
2401 ri_dni = max( dt*pac(i_nracni2ng_i,k,i,j)*fac4, min(0.0_rp, -rhoq2(i_ni,k,i,j) ))
2402 ii_dni = max( dt*pac(i_niacni2ns,k,i,j) , min(0.0_rp, -rhoq2(i_ni,k,i,j)-ri_dni ))
2403 is_dni = max( dt*pac(i_niacns2ns,k,i,j) , min(0.0_rp, -rhoq2(i_ni,k,i,j)-ri_dni-ii_dni ))
2405 fac3 = (rs_dqr-eps)/(dt*pac(i_lracls2lg_r,k,i,j)-eps)
2406 rs_dqs = max( dt*pac(i_lracls2lg_s,k,i,j)*fac3, min(0.0_rp, -rhoq2(i_qs,k,i,j)+sc_dqc+ii_dqi+is_dqi ))
2407 gs_dqs = max( dt*pac(i_lgacls2lg,k,i,j) , min(0.0_rp, -rhoq2(i_qs,k,i,j)+sc_dqc+ii_dqi+is_dqi-rs_dqs ))
2409 fac6 = (rs_dnr-eps)/(dt*pac(i_nracns2ng_r,k,i,j)-eps)
2411 rs_dns = max( dt*pac(i_nracns2ng_s,k,i,j)*fac6, min(0.0_rp, -rhoq2(i_ns,k,i,j)+0.50_rp*ii_dni+is_dni ))
2412 gs_dns = max( dt*pac(i_ngacns2ng,k,i,j) , min(0.0_rp, -rhoq2(i_ns,k,i,j)+0.50_rp*ii_dni+is_dni-rs_dns ))
2413 ss_dns = max( dt*pac(i_nsacns2ns,k,i,j) , min(0.0_rp, -rhoq2(i_ns,k,i,j)+0.50_rp*ii_dni+is_dni-rs_dns-gs_dns ))
2415 gg_dng = max( dt*pac(i_ngacng2ng,k,i,j) , min(0.0_rp, -rhoq2(i_ng,k,i,j) ))
2421 clp_dqr = (-rg_dqg-rs_dqs-ri_dqi) * (1.0_rp-sw)
2423 clp_dqs = -sc_dqc-ii_dqi-is_dqi
2424 clp_dqg = -gc_dqc -gs_dqs + (-rg_dqr-rs_dqr-rs_dqs-ri_dqr-ri_dqi) * sw
2429 clp_dns = -ii_dni*0.5_rp
2430 clp_dng = (-rs_dnr-ri_dnr) * sw
2433 clm_dqc = gc_dqc+sc_dqc+ic_dqc
2434 clm_dqr = (rg_dqr+rs_dqr+ri_dqr) * sw
2435 clm_dqi = ri_dqi+ii_dqi+is_dqi
2436 clm_dqs = rs_dqs+gs_dqs
2437 clm_dqg = rg_dqg * (1.0_rp-sw)
2439 clm_dnc = gc_dnc+sc_dnc+ic_dnc
2440 clm_dnr = (rg_dnr+rs_dnr+ri_dnr) * sw
2441 clm_dni = ri_dni+ii_dni+is_dni
2442 clm_dns = rs_dns+ss_dns+gs_dns
2443 clm_dng = gg_dng + rg_dng * (1.0_rp-sw)
2447 pco_dqi = max( dt*pq(i_licon,k,i,j), -clp_dqi )
2448 pco_dqs = max( dt*pq(i_lscon,k,i,j), -clp_dqs )
2449 pco_dqg = -pco_dqi-pco_dqs
2451 pco_dni = max( dt*pq(i_nicon,k,i,j), -clp_dni )
2452 pco_dns = max( dt*pq(i_nscon,k,i,j), -clp_dns )
2453 pco_dng = -pco_dni-pco_dns
2456 eml_dqi = max( dt*pq(i_liacm,k,i,j), min(0.0_rp, -rhoq2(i_qi,k,i,j)-(clp_dqi+clm_dqi)-pco_dqi ))
2457 eml_dqs = max( dt*pq(i_lsacm,k,i,j), min(0.0_rp, -rhoq2(i_qs,k,i,j)-(clp_dqs+clm_dqs)-pco_dqs ))
2458 eml_dqg = max( dt*(pq(i_lgacm,k,i,j)+pq(i_lgarm,k,i,j)+pq(i_lsarm,k,i,j)+pq(i_liarm,k,i,j)), &
2459 min(0.0_rp, -rhoq2(i_qg,k,i,j)-(clp_dqg+clm_dqg)-pco_dqg ))
2461 eml_dqr = -eml_dqs-eml_dqg
2463 eml_dni = max( dt*pq(i_niacm,k,i,j), min(0.0_rp, -rhoq2(i_ni,k,i,j)-(clp_dni+clm_dni)-pco_dni ))
2464 eml_dns = max( dt*pq(i_nsacm,k,i,j), min(0.0_rp, -rhoq2(i_ns,k,i,j)-(clp_dns+clm_dns)-pco_dns ))
2465 eml_dng = max( dt*(pq(i_ngacm,k,i,j)+pq(i_ngarm,k,i,j)+pq(i_nsarm,k,i,j)+pq(i_niarm,k,i,j)), &
2466 min(0.0_rp, -rhoq2(i_ng,k,i,j)-(clp_dng+clm_dng)-pco_dng ))
2468 eml_dnr = -eml_dns-eml_dng
2471 spl_dqg = max( dt*pq(i_lgspl,k,i,j), min(0.0_rp, -rhoq2(i_qg,k,i,j)-(clp_dqg+clm_dqg)-pco_dqg-eml_dqg ))
2472 spl_dqs = max( dt*pq(i_lsspl,k,i,j), min(0.0_rp, -rhoq2(i_qs,k,i,j)-(clp_dqs+clm_dqs)-pco_dqs-eml_dqs ))
2473 spl_dqi = -spl_dqg-spl_dqs
2474 fac9 = (spl_dqg-eps)/(dt*pq(i_lgspl,k,i,j)-eps)
2475 fac10 = (spl_dqs-eps)/(dt*pq(i_lsspl,k,i,j)-eps)
2476 spl_dni = dt*pq(i_nispl,k,i,j)*fac9*fac10
2479 drhogqc = (wrm_dqc + clp_dqc + clm_dqc + eml_dqc )
2480 drhognc = (wrm_dnc + clp_dnc + clm_dnc + eml_dnc )
2482 drhogqr = (wrm_dqr + clp_dqr + clm_dqr + eml_dqr )
2483 drhognr = (wrm_dnr + clp_dnr + clm_dnr + eml_dnr )
2485 drhogqi = ( clp_dqi + clm_dqi + pco_dqi + eml_dqi + spl_dqi)
2486 drhogni = ( clp_dni + clm_dni + pco_dni + eml_dni + spl_dni)
2488 drhogqs = ( clp_dqs + clm_dqs + pco_dqs + eml_dqs + spl_dqs)
2489 drhogns = ( clp_dns + clm_dns + pco_dns + eml_dns )
2491 drhogqg = ( clp_dqg + clm_dqg + pco_dqg + eml_dqg + spl_dqg)
2492 drhogng = ( clp_dng + clm_dng + pco_dng + eml_dng )
2496 rhoq0_t(k,i,j,i_qc) = drhogqc / dt
2497 rhoq0_t(k,i,j,i_nc) = drhognc / dt
2498 rhoq0_t(k,i,j,i_qr) = drhogqr / dt
2499 rhoq0_t(k,i,j,i_nr) = drhognr / dt
2500 rhoq0_t(k,i,j,i_qi) = drhogqi / dt
2501 rhoq0_t(k,i,j,i_ni) = drhogni / dt
2502 rhoq0_t(k,i,j,i_qs) = drhogqs / dt
2503 rhoq0_t(k,i,j,i_ns) = drhogns / dt
2504 rhoq0_t(k,i,j,i_qg) = drhogqg / dt
2505 rhoq0_t(k,i,j,i_ng) = drhogng / dt
2507 rhoe0_t(k,i,j) = lhf * ( drhogqi + drhogqs + drhogqg ) / dt
2509 dqc = drhogqc * rrho(k,i,j)
2510 dqr = drhogqr * rrho(k,i,j)
2511 dqi = drhogqi * rrho(k,i,j)
2512 dqs = drhogqs * rrho(k,i,j)
2513 dqg = drhogqg * rrho(k,i,j)
2518 cvtot0_t(k,i,j) = dcv/dt
2519 cptot0_t(k,i,j) = dcp/dt
2524 profile_stop(
"sn14_update_rhoq")
2535 rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe0_t(k,i,j)
2536 cvtot_t(k,i,j) = cvtot_t(k,i,j) + cvtot0_t(k,i,j)
2537 cptot_t(k,i,j) = cptot_t(k,i,j) + cptot0_t(k,i,j)
2550 rhoq(iq,k,i,j) = max(0.0_rp, rhoq(iq,k,i,j) + rhoq0_t(k,i,j,iq)*dt )
2563 rhoq_t(k,i,j,iq) = ( rhoq(iq,k,i,j) - dens(k,i,j)*qtrc(k,i,j,iq) )/dt
2571 if( opt_debug_tem )
call debug_tem_kij( ka, ks, ke, ia, is, ie, ja, js, je, &
2572 4, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )
2577 end subroutine mp_sn14
2581 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2591 integer,
intent(in) :: KA, KS, KE
2592 integer,
intent(in) :: IA, IS, IE
2593 integer,
intent(in) :: JA, JS, JE
2595 integer,
intent(in) :: point
2596 real(RP),
intent(in) :: tem(ka,ia,ja)
2597 real(RP),
intent(in) :: rho(ka,ia,ja)
2598 real(RP),
intent(in) :: pre(ka,ia,ja)
2599 real(RP),
intent(in) :: qv (ka,ia,ja)
2607 if ( tem(k,i,j) < tem_min &
2608 .OR. rho(k,i,j) < rho_min &
2609 .OR. pre(k,i,j) < 1.0_rp )
then 2611 log_info(
"ATMOS_PHY_MP_SN14_debug_tem_kij",
'(A,I3,A,4(F16.5),3(I6))') &
2612 "point: ", point,
" low tem,rho,pre:", tem(k,i,j), rho(k,i,j), pre(k,i,j), qv(k,i,j), k, i, j,
prc_myrank 2622 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2624 rho, tem, pre, qdry, &
2635 moist_psat_liq => atmos_saturation_psat_liq, &
2636 moist_psat_ice => atmos_saturation_psat_ice, &
2637 moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
2638 moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
2639 moist_dqsi_dtem_dens => atmos_saturation_dqs_dtem_dens_liq
2642 integer,
intent(in) :: KA, KS, KE
2643 integer,
intent(in) :: IA, IS, IE
2644 integer,
intent(in) :: JA, JS, JE
2646 real(RP),
intent(in) :: cz( ka,ia,ja)
2647 real(RP),
intent(in) :: fz(0:ka,ia,ja)
2648 real(RP),
intent(in) :: w (ka,ia,ja)
2649 real(RP),
intent(in) :: rho(ka,ia,ja)
2650 real(RP),
intent(in) :: tem(ka,ia,ja)
2651 real(RP),
intent(in) :: pre(ka,ia,ja)
2652 real(RP),
intent(in) :: qdry(ka,ia,ja)
2654 real(RP),
intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
2656 real(RP),
intent(in) :: cpa(ka,ia,ja)
2657 real(RP),
intent(in) :: dTdt_rad(ka,ia,ja)
2658 real(RP),
intent(in) :: qke(ka,ia,ja)
2659 real(RP),
intent(in) :: dt
2660 real(RP),
intent(in) :: CCN(ka,ia,ja)
2662 real(RP),
intent(out) :: PQ(pq_max,ka,ia,ja)
2667 real(RP),
parameter :: c_ccn_ocean= 1.00e+8_rp
2668 real(RP),
parameter :: c_ccn_land = 1.26e+9_rp
2669 real(RP),
save :: c_ccn = 1.00e+8_rp
2671 real(RP),
parameter :: kappa_ocean= 0.462_rp
2672 real(RP),
parameter :: kappa_land = 0.308_rp
2673 real(RP),
save :: kappa = 0.462_rp
2674 real(RP),
save :: c_in = 1.0_rp
2676 real(RP),
save :: nm_M92 = 1.e+3_rp
2677 real(RP),
save :: am_M92 = -0.639_rp
2678 real(RP),
save :: bm_M92 = 12.96_rp
2680 real(RP),
save :: in_max = 1000.e+3_rp
2681 real(RP),
save :: ssi_max= 0.60_rp
2682 real(RP),
save :: ssw_max= 1.1_rp
2684 logical,
save :: flag_first = .true.
2685 real(RP),
save :: qke_min = 0.03_rp
2686 real(RP),
save :: tem_ccn_low=233.150_rp
2687 real(RP),
save :: tem_in_low =173.150_rp
2688 logical,
save :: nucl_twomey = .false.
2689 logical,
save :: inucl_w = .false.
2691 namelist / param_atmos_phy_mp_sn14_nucleation / &
2694 nm_m92, am_m92, bm_m92, &
2699 nucl_twomey, inucl_w
2704 real(RP) :: esw(ka,ia,ja)
2705 real(RP) :: esi(ka,ia,ja)
2706 real(RP) :: ssw(ka,ia,ja)
2707 real(RP) :: ssi(ka,ia,ja)
2709 real(RP) :: w_dssidz(ka,ia,ja)
2711 real(RP) :: ssi_below(ka,ia,ja)
2712 real(RP) :: z_below(ka,ia,ja)
2716 real(RP) :: qsw(ka,ia,ja)
2717 real(RP) :: qsi(ka,ia,ja)
2718 real(RP) :: dqsidtem_rho(ka,ia,ja)
2719 real(RP) :: dssidt_rad(ka,ia,ja)
2721 real(RP) :: wssi, wdssi
2727 real(RP) :: sigma_w(ka,ia,ja)
2728 real(RP) :: weff(ka,ia,ja)
2729 real(RP) :: weff_max(ka,ia,ja)
2730 real(RP) :: velz(ka)
2732 real(RP) :: coef_ccn(ia,ja)
2733 real(RP) :: slope_ccn(ia,ja)
2734 real(RP) :: nc_new(ka,ia,ja)
2735 real(RP) :: nc_new_below(ka,ia,ja)
2737 real(RP) :: nc_new_max
2740 logical :: flag_nucleation(ka,ia,ja)
2742 real(RP) :: r_gravity
2743 real(RP),
parameter :: r_sqrt3=0.577350269_rp
2744 real(RP),
parameter :: eps=1.e-30_rp
2747 real(RP) :: dlcdt_max, dli_max
2748 real(RP) :: dncdt_max, dni_max
2753 if( flag_first )
then 2755 read(
io_fid_conf, nml=param_atmos_phy_mp_sn14_nucleation, end=100)
2756 100 log_nml(param_atmos_phy_mp_sn14_nucleation)
2759 if ( mp_couple_aerosol .AND. nucl_twomey )
then 2760 log_error(
"ATMOS_PHY_MP_SN14_nucleation_kij",*)
"nucl_twomey should be false when MP_couple_aerosol is true, stop" 2772 nc_uplim_d(1,i,j) = c_ccn*1.5_rp
2777 r_gravity = 1.0_rp/grav
2779 call moist_psat_liq ( ka, ks, ke, ia, is, ie, ja, js, je, &
2780 tem(:,:,:), esw(:,:,:) )
2781 call moist_psat_ice ( ka, ks, ke, ia, is, ie, ja, js, je, &
2782 tem(:,:,:), esi(:,:,:) )
2783 call moist_pres2qsat_liq ( ka, ks, ke, ia, is, ie, ja, js, je, &
2784 tem(:,:,:), pre(:,:,:), qdry(:,:,:), &
2786 call moist_pres2qsat_ice ( ka, ks, ke, ia, is, ie, ja, js, je, &
2787 tem(:,:,:), pre(:,:,:), qdry(:,:,:), &
2789 call moist_dqsi_dtem_dens( ka, ks, ke, ia, is, ie, ja, js, je, &
2790 tem(:,:,:), rho(:,:,:), &
2791 dqsidtem_rho(:,:,:) )
2794 a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
2802 pv = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)
2803 ssw(k,i,j) = min( mp_ssw_lim, ( pv/esw(k,i,j)-1.0_rp ) )*100.0_rp
2804 ssi(k,i,j) = (pv/esi(k,i,j) - 1.00_rp)
2806 ssi_below(k+1,i,j) = ssi(k,i,j)
2807 z_below(k+1,i,j) = cz(k,i,j)
2810 ssi_below(ks,i,j) = ssi(ks,i,j)
2811 z_below(ks,i,j) = cz(ks-1,i,j)
2816 coef_ccn(i,j) = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
2818 * (70.0_rp)**(kappa/(kappa + 2.0_rp))
2820 slope_ccn(i,j) = 1.5_rp*kappa/(kappa + 2.0_rp)
2823 sigma_w(k,i,j) = r_sqrt3*sqrt(max(qke(k,i,j),qke_min))
2825 sigma_w(ks-1,i,j) = sigma_w(ks,i,j)
2826 sigma_w(ke+1,i,j) = sigma_w(ke,i,j)
2829 weff(k,i,j) = w(k,i,j) - cpa(k,i,j)*r_gravity*dtdt_rad(k,i,j)
2835 if( mp_couple_aerosol )
then 2840 if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp )
then 2841 nc_new(k,i,j) = max( ccn(k,i,j), c_ccn )
2843 nc_new(k,i,j) = 0.0_rp
2851 if( nucl_twomey )
then 2858 weff_max(k,i,j) = weff(k,i,j) + sigma_w(k,i,j)
2860 if( (weff(k,i,j) > 1.e-8_rp) .AND. (ssw(k,i,j) > 1.e-10_rp) .AND. pre(k,i,j) > 300.e+2_rp )
then 2862 nc_new_max = coef_ccn(i,j)*weff_max(k,i,j)**slope_ccn(i,j)
2863 nc_new(k,i,j) = a_max*nc_new_max**b_max
2865 nc_new(k,i,j) = 0.0_rp
2875 if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp )
then 2876 nc_new(k,i,j) = c_ccn*ssw(k,i,j)**kappa
2878 nc_new(k,i,j) = 0.0_rp
2891 if( nc_new(k,i,j) > nc_uplim_d(1,i,j) )
then 2892 flag_nucleation(k,i,j) = .false.
2893 nc_new_below(k+1,i,j) = 1.e+30_rp
2894 else if( nc_new(k,i,j) > eps )
then 2895 flag_nucleation(k,i,j) = .true.
2896 nc_new_below(k+1,i,j) = nc_new(k,i,j)
2898 flag_nucleation(k,i,j) = .false.
2899 nc_new_below(k+1,i,j) = 0.0_rp
2902 nc_new_below(ks,i,j) = 0.0_rp
2915 if( mp_couple_aerosol )
then 2923 if ( flag_nucleation(k,i,j) .AND. &
2924 tem(k,i,j) > tem_ccn_low )
then 2925 dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2926 dncdt_max = dlcdt_max/xc_min
2928 dnc_new = nc_new(k,i,j)
2929 pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2930 pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2932 pq(i_ncccn,k,i,j) = 0.0_rp
2933 pq(i_lcccn,k,i,j) = 0.0_rp
2941 if( nucl_twomey )
then 2948 if ( flag_nucleation(k,i,j) .AND. &
2949 tem(k,i,j) > tem_ccn_low .AND. &
2950 nc_new(k,i,j) > rhoq(i_nc,k,i,j) )
then 2951 dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2952 dncdt_max = dlcdt_max/xc_min
2953 dnc_new = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
2954 pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2955 pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2957 pq(i_ncccn,k,i,j) = 0.0_rp
2958 pq(i_lcccn,k,i,j) = 0.0_rp
2968 if( tem(k,i,j) > tem_ccn_low .AND. &
2969 nc_new(k,i,j) > rhoq(i_nc,k,i,j) )
then 2970 dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2971 dncdt_max = dlcdt_max/xc_min
2972 dnc_new = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
2973 pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2974 pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2976 pq(i_ncccn,k,i,j) = 0.0_rp
2977 pq(i_lcccn,k,i,j) = 0.0_rp
2996 velz(k) = ( w(k,i,j) * ( cz(k+1,i,j) - fz(k,i,j) ) + w(k+1,i,j) * ( fz(k,i,j) - cz(k,i,j) ) ) / ( cz(k+1,i,j) - cz(k,i,j) )
3000 dzh = cz(k,i,j) - z_below(k,i,j)
3001 w_dssidz(k,i,j) = velz(k) * (ssi(k,i,j) - ssi_below(k,i,j))/dzh
3002 dssidt_rad(k,i,j) = -rhoq(i_qv,k,i,j)/(rho(k,i,j)*qsi(k,i,j)*qsi(k,i,j))*dqsidtem_rho(k,i,j)*dtdt_rad(k,i,j)
3003 dli_max = (rhoq(i_qv,k,i,j) - esi(k,i,j)/(rvap*tem(k,i,j)))*rdt
3004 dni_max = min( dli_max/xi_ccn, (in_max-rhoq(i_ni,k,i,j))*rdt )
3005 wdssi = min( w_dssidz(k,i,j)+dssidt_rad(k,i,j), 0.01_rp)
3006 wssi = min( ssi(k,i,j), ssi_max)
3008 if( ( wdssi > eps ) .AND. &
3009 (tem(k,i,j) < 273.15_rp ) .AND. &
3010 (rhoq(i_ni,k,i,j) < in_max ) .AND. &
3011 (wssi >= eps ) )
then 3014 pq(i_niccn,k,i,j) = min(dni_max, c_in*bm_m92*nm_m92*0.3_rp*exp(0.3_rp*bm_m92*(wssi-0.1_rp))*wdssi)
3016 pq(i_niccn,k,i,j) = min(dni_max, max(c_in*nm_m92*exp(0.3_rp*bm_m92*(wssi-0.1_rp) )-rhoq(i_ni,k,i,j),0.0_rp )*rdt )
3018 pq(i_liccn,k,i,j) = min(dli_max, pq(i_niccn,k,i,j)*xi_ccn )
3022 pq(i_niccn,k,i,j) = 0.0_rp
3023 pq(i_liccn,k,i,j) = 0.0_rp
3033 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
3035 tem, rhoq, xq, & ! in
3044 integer,
intent(in) :: KA, KS, KE
3045 integer,
intent(in) :: IA, IS, IE
3046 integer,
intent(in) :: JA, JS, JE
3048 real(RP),
intent(in) :: Pac(pac_max,ka,ia,ja)
3049 real(RP),
intent(in) :: tem(ka,ia,ja)
3050 real(RP),
intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
3051 real(RP),
intent(in) :: xq(hydro_max,ka,ia,ja)
3053 real(RP),
intent(inout):: PQ(pq_max,ka,ia,ja)
3055 logical,
save :: flag_first = .true.
3057 real(RP),
parameter :: pice = 350.0e+6_rp
3059 real(RP),
parameter :: pnc = 250.0_rp
3060 real(RP),
parameter :: rc_cr= 12.e-6_rp
3064 real(RP),
save :: xc_cr
3065 real(RP),
save :: alpha
3066 real(RP),
save :: gm, lgm
3070 real(RP) :: a0,a1,a2,a3,a4,a5
3071 real(RP) :: a6,a7,a8,a9,a10
3072 real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
3073 real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
3074 real(RP),
parameter :: eps=1.0e-30_rp
3078 real(RP) :: wn, wni, wns, wng
3081 if( flag_first )
then 3082 flag_first = .false.
3084 xc_cr = (2.0_rp*rc_cr/a_m(i_mp_qc))**(1.0_rp/b_m(i_mp_qc))
3085 alpha = (nu(i_mp_qc)+1.0_rp)/mu(i_mp_qc)
3086 gm = gammafunc(alpha)
3096 if (tem(k,i,j) > 270.16_rp)
then 3098 else if(tem(k,i,j) >= 268.16_rp)
then 3099 fp = (270.16_rp-tem(k,i,j))*0.5_rp
3100 else if(tem(k,i,j) >= 265.16_rp)
then 3101 fp = (tem(k,i,j)-265.16_rp)*0.333333333_rp
3109 x = coef_lambda(i_mp_qc)*(xc_cr/xq(i_mp_qc,k,i,j))**mu(i_mp_qc)
3111 if(x<1.e-2_rp*alpha)
then 3113 else if(x<alpha+1.0_rp)
then 3116 a1 = a0*x/(alpha+1.0_rp)
3117 a2 = a1*x/(alpha+2.0_rp)
3118 a3 = a2*x/(alpha+3.0_rp)
3119 a4 = a3*x/(alpha+4.0_rp)
3120 a5 = a4*x/(alpha+5.0_rp)
3121 a6 = a5*x/(alpha+6.0_rp)
3122 a7 = a6*x/(alpha+7.0_rp)
3123 a8 = a7*x/(alpha+8.0_rp)
3124 a9 = a8*x/(alpha+9.0_rp)
3125 a10 = a9*x/(alpha+10.0_rp)
3126 igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
3127 else if(x<alpha*1.d2)
then 3135 an1 = -(1.0_rp-alpha)
3137 d1 = 1.0_rp/(an1*d0+b1)
3142 an2 = -2.0_rp*(2.0_rp-alpha)
3144 d2 = 1.0_rp/(an2*d1+b2)
3149 igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
3154 n12 = rhoq(i_nc,k,i,j)*(1.0_rp-igm)
3156 wn = (pice + n12/((rhoq(i_qc,k,i,j)+xc_min)*pnc) )*fp
3157 wni = wn*(-pac(i_liaclc2li,k,i,j) )
3158 wns = wn*(-pac(i_lsaclc2ls,k,i,j) )
3159 wng = wn*(-pac(i_lgaclc2lg,k,i,j) )
3160 pq(i_nispl,k,i,j) = wni+wns+wng
3162 pq(i_lsspl,k,i,j) = - wns*xq(i_mp_qi,k,i,j)
3163 pq(i_lgspl,k,i,j) = - wng*xq(i_mp_qi,k,i,j)
3173 ! collection process
3174 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
3176 xq, dq_xave, vt_xave, & ! in
3177 ! rho ! [Add] 11/08/30
3181 moist_psat_ice => atmos_saturation_psat_ice
3184 integer,
intent(in) :: KA, KS, KE
3185 integer,
intent(in) :: IA, IS, IE
3186 integer,
intent(in) :: JA, JS, JE
3191 real(RP),
intent(in) :: wtem(ka,ia,ja)
3193 real(RP),
intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
3195 real(RP),
intent(in) :: xq(hydro_max,ka,ia,ja)
3197 real(RP),
intent(in) :: dq_xave(hydro_max,ka,ia,ja)
3199 real(RP),
intent(in) :: vt_xave(hydro_max,2,ka,ia,ja)
3203 real(RP),
intent(inout):: PQ(pq_max,ka,ia,ja)
3205 real(RP),
intent(out):: Pac(pac_max,ka,ia,ja)
3210 real(RP),
save :: dc0 = 15.0e-6_rp
3211 real(RP),
save :: dc1 = 40.0e-6_rp
3212 real(RP),
save :: di0 = 150.0e-6_rp
3213 real(RP),
save :: ds0 = 150.0e-6_rp
3214 real(RP),
save :: dg0 = 150.0e-6_rp
3216 real(RP),
save :: sigma_c=0.00_rp
3217 real(RP),
save :: sigma_r=0.00_rp
3218 real(RP),
save :: sigma_i=0.2_rp
3219 real(RP),
save :: sigma_s=0.2_rp
3220 real(RP),
save :: sigma_g=0.00_rp
3222 real(RP),
save :: E_im = 0.80_rp
3223 real(RP),
save :: E_sm = 0.80_rp
3224 real(RP),
save :: E_gm = 1.00_rp
3226 real(RP),
save :: E_ir=1.0_rp
3227 real(RP),
save :: E_sr=1.0_rp
3228 real(RP),
save :: E_gr=1.0_rp
3229 real(RP),
save :: E_ii=1.0_rp
3230 real(RP),
save :: E_si=1.0_rp
3231 real(RP),
save :: E_gi=1.0_rp
3232 real(RP),
save :: E_ss=1.0_rp
3233 real(RP),
save :: E_gs=1.0_rp
3234 real(RP),
save :: E_gg=1.0_rp
3237 integer,
save :: i_iconv2g=1
3238 integer,
save :: i_sconv2g=1
3240 real(RP),
save :: rho_g = 900.0_rp
3242 real(RP),
save :: cfill_i = 0.68_rp
3243 real(RP),
save :: cfill_s = 0.01_rp
3245 real(RP),
save :: di_cri = 500.e-6_rp
3247 logical,
save :: opt_stick_KS96=.false.
3248 logical,
save :: opt_stick_CO86=.false.
3249 real(RP),
parameter :: a_dec = 0.883_rp
3250 real(RP),
parameter :: b_dec = 0.093_rp
3251 real(RP),
parameter :: c_dec = 0.00348_rp
3252 real(RP),
parameter :: d_dec = 4.5185e-5_rp
3254 logical,
save :: flag_first = .true.
3255 namelist / param_atmos_phy_mp_sn14_collection / &
3256 dc0, dc1, di0, ds0, dg0, &
3257 sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
3261 e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
3262 i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
3264 real(RP) :: tem(ka,ia,ja)
3267 real(RP) :: E_c, E_r, E_i, E_s, E_g
3268 real(RP) :: E_ic, E_sc, E_gc
3270 real(RP) :: E_stick(ka,ia,ja)
3272 real(RP) :: temc, temc2, temc3
3275 real(RP) :: esi(ka,ia,ja)
3277 real(RP) :: temc_p, temc_m
3289 real(RP) :: coef_acc_LCI, coef_acc_NCI
3290 real(RP) :: coef_acc_LCS, coef_acc_NCS
3292 real(RP) :: coef_acc_LCG, coef_acc_NCG
3293 real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I
3294 real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R
3295 real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S
3296 real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R
3297 real(RP) :: coef_acc_LRG, coef_acc_NRG
3298 real(RP) :: coef_acc_LII, coef_acc_NII
3299 real(RP) :: coef_acc_LIS, coef_acc_NIS
3300 real(RP) :: coef_acc_NSS
3301 real(RP) :: coef_acc_NGG
3302 real(RP) :: coef_acc_LSG, coef_acc_NSG
3304 real(RP) :: dcdc, dcdi, dcds, dcdg
3305 real(RP) :: drdr, drdi, drds, drdg
3306 real(RP) :: didi, dids, didg
3307 real(RP) :: dsds, dsdg
3310 real(RP) :: vcvc, vcvi, vcvs, vcvg
3311 real(RP) :: vrvr, vrvi, vrvs, vrvg
3312 real(RP) :: vivi, vivs, vivg
3313 real(RP) :: vsvs, vsvg
3316 real(RP) :: wx_cri, wx_crs
3317 real(RP) :: coef_emelt
3324 if( flag_first )
then 3326 read(
io_fid_conf, nml=param_atmos_phy_mp_sn14_collection, end=100 )
3327 100 log_nml(param_atmos_phy_mp_sn14_collection)
3328 flag_first = .false.
3335 tem(k,i,j) = max( wtem(k,i,j), tem_min )
3340 call moist_psat_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
3341 tem(:,:,:), esi(:,:,:) )
3343 if( opt_stick_ks96 )
then 3348 temc = tem(k,i,j) - t00
3351 e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
3352 esi_rat = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)/esi(k,i,j)
3353 e_stick(k,i,j) = min(1.0_rp, e_dec*esi_rat)
3357 else if( opt_stick_co86 )
then 3362 temc = min(tem(k,i,j) - t00,0.0_rp)
3363 w1 = 0.035_rp*temc-0.7_rp
3364 e_stick(k,i,j) = 10._rp**w1
3373 temc_m = min(tem(k,i,j) - t00,0.0_rp)
3374 e_stick(k,i,j) = exp(0.09_rp*temc_m)
3380 profile_start(
"sn14_collection")
3387 temc_p = max(tem(k,i,j) - t00,0.0_rp)
3389 ave_dc = coef_d(i_mp_qc)*xq(i_mp_qc,k,i,j)**b_m(i_mp_qc)
3390 ave_di = coef_d(i_mp_qi)*xq(i_mp_qi,k,i,j)**b_m(i_mp_qi)
3391 ave_ds = coef_d(i_mp_qs)*xq(i_mp_qs,k,i,j)**b_m(i_mp_qs)
3392 ave_dg = coef_d(i_mp_qg)*xq(i_mp_qg,k,i,j)**b_m(i_mp_qg)
3395 e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
3396 sw = 0.5_rp - sign(0.5_rp, di0-ave_di)
3398 sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds)
3400 sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg)
3407 dcdc = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qc,k,i,j)
3408 drdr = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qr,k,i,j)
3409 didi = dq_xave(i_mp_qi,k,i,j) * dq_xave(i_mp_qi,k,i,j)
3410 dsds = dq_xave(i_mp_qs,k,i,j) * dq_xave(i_mp_qs,k,i,j)
3411 dgdg = dq_xave(i_mp_qg,k,i,j) * dq_xave(i_mp_qg,k,i,j)
3412 dcdi = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qi,k,i,j)
3413 dcds = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qs,k,i,j)
3414 dcdg = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qg,k,i,j)
3415 drdi = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qi,k,i,j)
3416 drds = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qs,k,i,j)
3417 drdg = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qg,k,i,j)
3418 dids = dq_xave(i_mp_qi,k,i,j) * dq_xave(i_mp_qs,k,i,j)
3420 dsdg = dq_xave(i_mp_qs,k,i,j) * dq_xave(i_mp_qg,k,i,j)
3422 vcvc = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qc,2,k,i,j)
3423 vrvr = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qr,2,k,i,j)
3424 vivi = vt_xave(i_mp_qi,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
3425 vsvs = vt_xave(i_mp_qs,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
3426 vgvg = vt_xave(i_mp_qg,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
3427 vcvi = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
3428 vcvs = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
3429 vcvg = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
3430 vrvi = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
3431 vrvs = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
3432 vrvg = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
3433 vivs = vt_xave(i_mp_qi,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
3435 vsvg = vt_xave(i_mp_qs,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
3444 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3445 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3446 + sigma_i + sigma_c )**0.5_rp
3448 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3449 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3450 + sigma_i + sigma_c )**0.5_rp
3451 pac(i_liaclc2li,k,i,j)= -0.25_rp*pi*e_ic*rhoq(i_ni,k,i,j)*rhoq(i_qc,k,i,j)*coef_acc_lci
3452 pac(i_niacnc2ni,k,i,j)= -0.25_rp*pi*e_ic*rhoq(i_ni,k,i,j)*rhoq(i_nc,k,i,j)*coef_acc_nci
3456 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3457 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3458 + sigma_s + sigma_c )**0.5_rp
3460 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3461 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3462 + sigma_s + sigma_c )**0.5_rp
3463 pac(i_lsaclc2ls,k,i,j)= -0.25_rp*pi*e_sc*rhoq(i_ns,k,i,j)*rhoq(i_qc,k,i,j)*coef_acc_lcs
3464 pac(i_nsacnc2ns,k,i,j)= -0.25_rp*pi*e_sc*rhoq(i_ns,k,i,j)*rhoq(i_nc,k,i,j)*coef_acc_ncs
3468 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3469 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3470 + sigma_g + sigma_c )**0.5_rp
3472 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3473 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3474 + sigma_g + sigma_c )**0.5_rp
3475 pac(i_lgaclc2lg,k,i,j)= -0.25_rp*pi*e_gc*rhoq(i_ng,k,i,j)*rhoq(i_qc,k,i,j)*coef_acc_lcg
3476 pac(i_ngacnc2ng,k,i,j)= -0.25_rp*pi*e_gc*rhoq(i_ng,k,i,j)*rhoq(i_nc,k,i,j)*coef_acc_ncg
3479 ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3480 * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3481 + sigma_g + sigma_s )**0.5_rp
3483 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3485 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3486 + sigma_g + sigma_s )**0.5_rp
3487 pac(i_lgacls2lg,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_gs*rhoq(i_ng,k,i,j)*rhoq(i_qs,k,i,j)*coef_acc_lsg
3488 pac(i_ngacns2ng,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_gs*rhoq(i_ng,k,i,j)*rhoq(i_ns,k,i,j)*coef_acc_nsg
3493 ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
3494 * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
3495 + sigma_i + sigma_s )**0.5_rp
3497 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
3498 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
3499 + sigma_i + sigma_s )**0.5_rp
3500 pac(i_liacls2ls,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_si*rhoq(i_ns,k,i,j)*rhoq(i_qi,k,i,j)*coef_acc_lis
3501 pac(i_niacns2ns,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_si*rhoq(i_ns,k,i,j)*rhoq(i_ni,k,i,j)*coef_acc_nis
3503 sw = sign(0.5_rp, t00-tem(k,i,j)) + 0.5_rp
3513 ( ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) * sw &
3514 + ( 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) ) &
3515 * sqrt( ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg ) * sw &
3516 + ( 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) &
3517 + sigma_r + sigma_g )
3518 pac(i_lraclg2lg,k,i,j) = -0.25_rp*pi*e_gr*coef_acc_lrg &
3519 * ( rhoq(i_ng,k,i,j)*rhoq(i_qr,k,i,j) * sw &
3520 + rhoq(i_nr,k,i,j)*rhoq(i_qg,k,i,j) * (1.0_rp-sw) )
3522 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) &
3523 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg &
3524 + sigma_r + sigma_g )**0.5_rp
3525 pac(i_nracng2ng,k,i,j) = -0.25_rp*pi*e_gr*rhoq(i_ng,k,i,j)*rhoq(i_nr,k,i,j)*coef_acc_nrg
3535 ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
3536 * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
3537 + sigma_r + sigma_i )**0.5_rp
3539 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
3540 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
3541 + sigma_r + sigma_i )**0.5_rp
3542 pac(i_lracli2lg_i,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_nr,k,i,j)*rhoq(i_qi,k,i,j)*coef_acc_lri_i
3543 pac(i_nracni2ng_i,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_nr,k,i,j)*rhoq(i_ni,k,i,j)*coef_acc_nri_i
3546 ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
3547 * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
3548 + sigma_r + sigma_i )**0.5_rp
3550 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
3551 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
3552 + sigma_r + sigma_i )**0.5_rp
3553 pac(i_lracli2lg_r,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_ni,k,i,j)*rhoq(i_qr,k,i,j)*coef_acc_lri_r
3554 pac(i_nracni2ng_r,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_ni,k,i,j)*rhoq(i_nr,k,i,j)*coef_acc_nri_r
3558 ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
3559 * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
3560 + sigma_r + sigma_s )**0.5_rp
3562 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
3563 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
3564 + sigma_r + sigma_s )**0.5_rp
3565 pac(i_lracls2lg_s,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_nr,k,i,j)*rhoq(i_qs,k,i,j)*coef_acc_lrs_s
3566 pac(i_nracns2ng_s,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_nr,k,i,j)*rhoq(i_ns,k,i,j)*coef_acc_nrs_s
3569 ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
3570 * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
3571 + sigma_r + sigma_s )**0.5_rp
3573 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
3574 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
3575 + sigma_r + sigma_s )**0.5_rp
3576 pac(i_lracls2lg_r,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_ns,k,i,j)*rhoq(i_qr,k,i,j)*coef_acc_lrs_r
3577 pac(i_nracns2ng_r,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_ns,k,i,j)*rhoq(i_nr,k,i,j)*coef_acc_nrs_r
3585 ( delta_b0(i_mp_qi)*didi + delta_ab1(i_mp_qi,i_mp_qi)*didi + delta_b1(i_mp_qi)*didi ) &
3586 * ( theta_b0(i_mp_qi)*vivi - theta_ab1(i_mp_qi,i_mp_qi)*vivi + theta_b1(i_mp_qi)*vivi &
3587 + sigma_i + sigma_i )**0.5_rp
3589 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qi,i_mp_qi)*didi + delta_b0(i_mp_qi)*didi ) &
3590 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qi,i_mp_qi)*vivi + theta_b0(i_mp_qi)*vivi &
3591 + sigma_i + sigma_i )**0.5_rp
3592 pac(i_liacli2ls,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_ii*rhoq(i_ni,k,i,j)*rhoq(i_qi,k,i,j)*coef_acc_lii
3593 pac(i_niacni2ns,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_ii*rhoq(i_ni,k,i,j)*rhoq(i_ni,k,i,j)*coef_acc_nii
3605 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qs,i_mp_qs)*dsds + delta_b0(i_mp_qs)*dsds ) &
3606 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qs,i_mp_qs)*vsvs + theta_b0(i_mp_qs)*vsvs &
3607 + sigma_s + sigma_s )**0.5_rp
3608 pac(i_nsacns2ns,k,i,j)= -0.125_rp*pi*e_stick(k,i,j)*e_ss*rhoq(i_ns,k,i,j)*rhoq(i_ns,k,i,j)*coef_acc_nss
3612 ( delta_b0(i_mp_qg)*dgdg + delta_ab0(i_mp_qg,i_mp_qg)*dgdg + delta_b0(i_mp_qg)*dgdg ) &
3613 * ( theta_b0(i_mp_qg)*vgvg - theta_ab0(i_mp_qg,i_mp_qg)*vgvg + theta_b0(i_mp_qg)*vgvg &
3614 + sigma_g + sigma_g )**0.5_rp
3615 pac(i_ngacng2ng,k,i,j)= -0.125_rp*pi*e_stick(k,i,j)*e_gg*rhoq(i_ng,k,i,j)*rhoq(i_ng,k,i,j)*coef_acc_ngg
3622 sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di)
3623 wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di*ave_di*ave_di/xq(i_mp_qi,k,i,j) - 1.0_rp ) * sw
3624 pq(i_licon,k,i,j) = i_iconv2g * pac(i_liaclc2li,k,i,j)/max(1.0_rp, wx_cri) * sw
3625 pq(i_nicon,k,i,j) = i_iconv2g * pq(i_licon,k,i,j)/xq(i_mp_qi,k,i,j) * sw
3628 wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds*ave_ds*ave_ds/xq(i_mp_qs,k,i,j) - 1.0_rp )
3629 pq(i_lscon,k,i,j) = i_sconv2g * (pac(i_lsaclc2ls,k,i,j))/max(1.0_rp, wx_crs)
3630 pq(i_nscon,k,i,j) = i_sconv2g * pq(i_lscon,k,i,j)/xq(i_mp_qs,k,i,j)
3639 coef_emelt = cl/lhf0*temc_p
3641 pq(i_lgacm,k,i,j) = coef_emelt*pac(i_lgaclc2lg,k,i,j)
3642 pq(i_ngacm,k,i,j) = pq(i_lgacm,k,i,j)/xq(i_mp_qg,k,i,j)
3644 pq(i_lgarm,k,i,j) = coef_emelt*pac(i_lraclg2lg,k,i,j)
3645 pq(i_ngarm,k,i,j) = pq(i_lgarm,k,i,j)/xq(i_mp_qg,k,i,j)
3647 pq(i_lsacm,k,i,j) = coef_emelt*(pac(i_lsaclc2ls,k,i,j))
3648 pq(i_nsacm,k,i,j) = pq(i_lsacm,k,i,j)/xq(i_mp_qs,k,i,j)
3650 pq(i_lsarm,k,i,j) = coef_emelt*(pac(i_lracls2lg_r,k,i,j)+pac(i_lracls2lg_s,k,i,j))
3651 pq(i_nsarm,k,i,j) = pq(i_lsarm,k,i,j)/xq(i_mp_qg,k,i,j)
3653 pq(i_liacm,k,i,j) = coef_emelt*pac(i_liaclc2li,k,i,j)
3654 pq(i_niacm,k,i,j) = pq(i_liacm,k,i,j)/xq(i_mp_qi,k,i,j)
3656 pq(i_liarm,k,i,j) = coef_emelt*(pac(i_lracli2lg_r,k,i,j)+pac(i_lracli2lg_i,k,i,j))
3657 pq(i_niarm,k,i,j) = pq(i_liarm,k,i,j)/xq(i_mp_qg,k,i,j)
3661 profile_stop(
"sn14_collection")
3669 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
3670 rhoq, xq, dq_xave, &
3675 integer,
intent(in) :: KA, KS, KE
3676 integer,
intent(in) :: IA, IS, IE
3677 integer,
intent(in) :: JA, JS, JE
3679 real(RP),
intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
3680 real(RP),
intent(in) :: xq(hydro_max,ka,ia,ja)
3681 real(RP),
intent(in) :: dq_xave(hydro_max,ka,ia,ja)
3682 real(RP),
intent(in) :: rho(ka,ia,ja)
3684 real(RP),
intent(inout) :: PQ(pq_max,ka,ia,ja)
3687 real(RP),
parameter :: kcc = 4.44e+9_rp
3688 real(RP),
parameter :: tau_min = 1.e-20_rp
3689 real(RP),
parameter :: rx_sep = 1.0_rp/x_sep
3692 real(RP),
parameter :: kcr = 5.8_rp
3693 real(RP),
parameter :: thr_acc = 5.e-5_rp
3696 real(RP),
parameter :: krr = 4.33_rp
3697 real(RP),
parameter :: kaprr = 60.7_rp
3698 real(RP),
parameter :: kbr = 1000._rp
3699 real(RP),
parameter :: kapbr = 2.3e+3_rp
3700 real(RP),
parameter :: dr_min = 0.35e-3_rp
3703 real(RP) :: coef_nuc0
3704 real(RP) :: coef_nuc1
3705 real(RP) :: coef_aut0
3706 real(RP) :: coef_aut1
3717 coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
3718 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)
3719 coef_aut0 = -kcc*coef_nuc0
3720 coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
3725 lwc = rhoq(i_qr,k,i,j)+rhoq(i_qc,k,i,j)
3726 if( lwc > xc_min )
then 3727 tau = max(tau_min, rhoq(i_qr,k,i,j)/lwc)
3731 rho_fac = sqrt(rho_0/max(rho(k,i,j),rho_min))
3734 psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3
3735 pq(i_ncaut,k,i,j) = coef_aut0*rhoq(i_qc,k,i,j)*rhoq(i_qc,k,i,j)*rho_fac*rho_fac
3737 pq(i_lcaut,k,i,j) = coef_aut1*lwc*lwc*xq(i_mp_qc,k,i,j)*xq(i_mp_qc,k,i,j) &
3738 *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac
3739 pq(i_nraut,k,i,j) = -rx_sep*pq(i_lcaut,k,i,j)
3742 psi_acc =(tau/(tau+thr_acc))**4
3743 pq(i_lcacc,k,i,j) = -kcr*rhoq(i_qc,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac*psi_acc
3744 pq(i_ncacc,k,i,j) = -kcr*rhoq(i_nc,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac*psi_acc
3747 pq(i_nrslc,k,i,j) = -krr*rhoq(i_nr,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac
3750 ddr = min(1.e-3_rp, dq_xave(i_mp_qr,k,i,j) - dr_eq )
3751 if ( dq_xave(i_mp_qr,k,i,j) < dr_min )
then 3753 else if ( dq_xave(i_mp_qr,k,i,j) <= dr_eq )
then 3756 psi_brk = exp(kapbr*ddr) - 1.0_rp
3758 pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j)
3767 subroutine dep_vapor_melt_ice_kij( &
3768 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
3769 rho, tem, pre, qd, & ! in
3772 xq, vt_xave, dq_xave, & ! in
3778 integer,
intent(in) :: KA, KS, KE
3779 integer,
intent(in) :: IA, IS, IE
3780 integer,
intent(in) :: JA, JS, JE
3783 real(RP),
intent(inout) :: PQ(pq_max,ka,ia,ja)
3785 real(RP),
intent(in) :: rho(ka,ia,ja)
3786 real(RP),
intent(in) :: tem(ka,ia,ja)
3787 real(RP),
intent(in) :: pre(ka,ia,ja)
3788 real(RP),
intent(in) :: qd(ka,ia,ja)
3789 real(RP),
intent(in) :: esw(ka,ia,ja)
3790 real(RP),
intent(in) :: esi(ka,ia,ja)
3791 real(RP),
intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
3792 real(RP),
intent(in) :: xq(hydro_max,ka,ia,ja)
3796 real(RP),
intent(in) :: vt_xave(hydro_max,2,ka,ia,ja)
3798 real(RP),
intent(in) :: dq_xave(hydro_max,ka,ia,ja)
3801 real(RP) :: temc_lim
3808 real(RP) :: nua, r_nua
3814 real(RP) :: Gwr, Gii, Gis, Gig
3819 real(RP) :: Nrers_r2, Nreis_r2
3820 real(RP) :: Nress_r2, Nregs_r2
3822 real(RP) :: Nrerl_r2, Nreil_r2
3823 real(RP) :: Nresl_r2, Nregl_r2
3824 real(RP) :: NscNrer_s, NscNrer_l
3825 real(RP) :: NscNrei_s, NscNrei_l
3826 real(RP) :: NscNres_s, NscNres_l
3827 real(RP) :: NscNreg_s, NscNreg_l
3828 real(RP) :: ventLR_s, ventLR_l
3829 real(RP) :: ventNI_s, ventNI_l, ventLI_s, ventLI_l
3830 real(RP) :: ventNS_s, ventNS_l, ventLS_s, ventLS_l
3831 real(RP) :: ventNG_s, ventNG_l, ventLG_s, ventLG_l
3833 real(RP) :: wtr, wti, wts, wtg
3834 real(RP),
parameter :: r_14=1.0_rp/1.4_rp
3835 real(RP),
parameter :: r_15=1.0_rp/1.5_rp
3838 real(RP) :: ventNI, ventLI
3839 real(RP) :: ventNS, ventLS
3840 real(RP) :: ventNG, ventLG
3842 real(RP),
parameter :: Re_max=1.e+3_rp
3843 real(RP),
parameter :: Re_min=1.e-4_rp
3855 profile_start(
"sn14_dep_vapor")
3859 temc = tem(k,i,j) - t00
3860 temc_lim= max(temc, -40._rp )
3861 rho_lim = max(rho(k,i,j),rho_min)
3862 qv = rhoq(i_qv,k,i,j)/rho_lim
3863 pre_lim = rho_lim*(qd(k,i,j)*rdry + qv*rvap)*(temc_lim+t00)
3871 dw = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(p00/pre_lim)
3872 kalfa = ka0 + temc_lim*dka_dt
3873 mua = mua0 + temc_lim*dmua_dt
3876 gw = (lhv0/kalfa/tem(k,i,j))*(lhv0/rvap/tem(k,i,j)-1.0_rp)+(rvap*tem(k,i,j)/dw/esw(k,i,j))
3877 gi = (lhs0/kalfa/tem(k,i,j))*(lhs0/rvap/tem(k,i,j)-1.0_rp)+(rvap*tem(k,i,j)/dw/esi(k,i,j))
3879 gwr = 4.0_rp*pi/cap(i_mp_qr)/gw
3880 gii = 4.0_rp*pi/cap(i_mp_qi)/gi
3881 gis = 4.0_rp*pi/cap(i_mp_qs)/gi
3882 gig = 4.0_rp*pi/cap(i_mp_qg)/gi
3885 nsc_r3 = (nua/dw)**(0.33333333_rp)
3888 nrers_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qr,1,k,i,j)*dq_xave(i_mp_qr,k,i,j)*r_nua)))
3889 nreis_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qi,1,k,i,j)*dq_xave(i_mp_qi,k,i,j)*r_nua)))
3890 nress_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qs,1,k,i,j)*dq_xave(i_mp_qs,k,i,j)*r_nua)))
3891 nregs_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qg,1,k,i,j)*dq_xave(i_mp_qg,k,i,j)*r_nua)))
3894 nrerl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qr,2,k,i,j)*dq_xave(i_mp_qr,k,i,j)*r_nua)))
3895 nreil_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qi,2,k,i,j)*dq_xave(i_mp_qi,k,i,j)*r_nua)))
3896 nresl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qs,2,k,i,j)*dq_xave(i_mp_qs,k,i,j)*r_nua)))
3897 nregl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qg,2,k,i,j)*dq_xave(i_mp_qg,k,i,j)*r_nua)))
3898 nscnrer_s=nsc_r3*nrers_r2
3899 nscnrer_l=nsc_r3*nrerl_r2
3901 nscnrei_s=nsc_r3*nreis_r2
3902 nscnrei_l=nsc_r3*nreil_r2
3904 nscnres_s=nsc_r3*nress_r2
3905 nscnres_l=nsc_r3*nresl_r2
3907 nscnreg_s=nsc_r3*nregs_r2
3908 nscnreg_l=nsc_r3*nregl_r2
3910 ventlr_s = ah_vent1(i_mp_qr,1) + bh_vent1(i_mp_qr,1)*nscnrer_s
3911 ventlr_l = ah_vent1(i_mp_qr,2) + bh_vent1(i_mp_qr,2)*nscnrer_l
3913 ventni_s = ah_vent0(i_mp_qi,1) + bh_vent0(i_mp_qi,1)*nscnrei_s
3914 ventni_l = ah_vent0(i_mp_qi,2) + bh_vent0(i_mp_qi,2)*nscnrei_l
3915 ventli_s = ah_vent1(i_mp_qi,1) + bh_vent1(i_mp_qi,1)*nscnrei_s
3916 ventli_l = ah_vent1(i_mp_qi,2) + bh_vent1(i_mp_qi,2)*nscnrei_l
3918 ventns_s = ah_vent0(i_mp_qs,1) + bh_vent0(i_mp_qs,1)*nscnres_s
3919 ventns_l = ah_vent0(i_mp_qs,2) + bh_vent0(i_mp_qs,2)*nscnres_l
3920 ventls_s = ah_vent1(i_mp_qs,1) + bh_vent1(i_mp_qs,1)*nscnres_s
3921 ventls_l = ah_vent1(i_mp_qs,2) + bh_vent1(i_mp_qs,2)*nscnres_l
3923 ventng_s = ah_vent0(i_mp_qg,1) + bh_vent0(i_mp_qg,1)*nscnreg_s
3924 ventng_l = ah_vent0(i_mp_qg,2) + bh_vent0(i_mp_qg,2)*nscnreg_l
3925 ventlg_s = ah_vent1(i_mp_qg,1) + bh_vent1(i_mp_qg,1)*nscnreg_s
3926 ventlg_l = ah_vent1(i_mp_qg,2) + bh_vent1(i_mp_qg,2)*nscnreg_l
3930 wtr = ( min(max( nscnrer_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3931 wti = ( min(max( nscnrei_s , 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3932 wts = ( min(max( nscnres_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3933 wtg = ( min(max( nscnreg_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3935 ventni = (1.0_rp-wti)*ventni_s + wti*ventni_l
3936 ventns = (1.0_rp-wts)*ventns_s + wts*ventns_l
3937 ventng = (1.0_rp-wtg)*ventng_s + wtg*ventng_l
3939 ventlr = (1.0_rp-wtr)*ventlr_s + wtr*ventlr_l
3940 ventli = (1.0_rp-wti)*ventli_s + wti*ventli_l
3941 ventls = (1.0_rp-wts)*ventls_s + wts*ventls_l
3942 ventlg = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
3961 pq(i_lcdep,k,i,j) = gwr*rhoq(i_nc,k,i,j)*dq_xave(i_mp_qc,k,i,j)*coef_deplc
3962 pq(i_lrdep,k,i,j) = gwr*rhoq(i_nr,k,i,j)*dq_xave(i_mp_qr,k,i,j)*ventlr
3963 pq(i_lidep,k,i,j) = gii*rhoq(i_ni,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventli
3964 pq(i_lsdep,k,i,j) = gis*rhoq(i_ns,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventls
3965 pq(i_lgdep,k,i,j) = gig*rhoq(i_ng,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventlg
3966 pq(i_nrdep,k,i,j) = pq(i_lrdep,k,i,j)/xq(i_mp_qr,k,i,j)
3967 pq(i_nidep,k,i,j) = 0.0_rp
3968 pq(i_nsdep,k,i,j) = pq(i_lsdep,k,i,j)/xq(i_mp_qs,k,i,j)
3969 pq(i_ngdep,k,i,j) = pq(i_lgdep,k,i,j)/xq(i_mp_qg,k,i,j)
3977 dt = kalfa/(cpvap*rho_0)
3983 gm = 2.0_rp*pi/emelt&
3984 * ( (kalfa*dt/dw)*(temc) + (dw*lhs0/rvap)*(esi(k,i,j)/tem(k,i,j)-psat0/t00) )
3989 sw = ( sign(0.5_rp,temc) + 0.5_rp ) * ( sign(0.5_rp,gm-eps) + 0.5_rp )
3993 pq(i_limlt,k,i,j) = - gm * rhoq(i_qi,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventli/xq(i_mp_qi,k,i,j) * sw
3995 pq(i_nimlt,k,i,j) = - gm * rhoq(i_ni,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventni/xq(i_mp_qi,k,i,j) * sw
3996 pq(i_lsmlt,k,i,j) = - gm * rhoq(i_qs,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventls/xq(i_mp_qs,k,i,j) * sw
3998 pq(i_nsmlt,k,i,j) = - gm * rhoq(i_ns,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventns/xq(i_mp_qs,k,i,j) * sw
3999 pq(i_lgmlt,k,i,j) = - gm * rhoq(i_qg,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventlg/xq(i_mp_qg,k,i,j) * sw
4001 pq(i_ngmlt,k,i,j) = - gm * rhoq(i_ng,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventng/xq(i_mp_qg,k,i,j) * sw
4006 profile_stop(
"sn14_dep_vapor")
4009 end subroutine dep_vapor_melt_ice_kij
4012 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
4021 integer,
intent(in) :: KA, KS, KE
4022 integer,
intent(in) :: IA, IS, IE
4023 integer,
intent(in) :: JA, JS, JE
4025 real(RP),
intent(in) :: dt
4027 real(RP),
intent(in) :: tem(ka,ia,ja)
4029 real(RP),
intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
4030 real(RP),
intent(in) :: xq(hydro_max,ka,ia,ja)
4032 real(RP),
intent(inout):: PQ(pq_max,ka,ia,ja)
4034 real(RP),
parameter :: temc_min = -65.0_rp
4035 real(RP),
parameter :: a_het = 0.2_rp
4036 real(RP),
parameter :: b_het = 0.65_rp
4038 real(RP) :: coef_m2_c
4039 real(RP) :: coef_m2_r
4041 real(RP) :: temc, temc2, temc3, temc4
4043 real(RP) :: Jhom, Jhet
4051 coef_m2_c = coef_m2(i_mp_qc)
4052 coef_m2_r = coef_m2(i_mp_qr)
4054 profile_start(
"sn14_freezing")
4058 temc = max( tem(k,i,j) - t00, temc_min )
4061 jhet = a_het*exp( -b_het*temc - 1.0_rp )
4064 if( temc < -65.0_rp )
then 4065 jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
4066 jhet = a_het*exp( 65.0_rp*b_het - 1.0_rp )
4067 else if( temc < -30.0_rp )
then 4072 - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
4073 - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
4074 else if( temc < 0.0_rp)
then 4075 jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
4087 pq(i_lchom,k,i,j) = 0.0_rp
4088 pq(i_nchom,k,i,j) = 0.0_rp
4090 #if defined(PGI) || defined(SX) 4091 tmp = min( xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp)
4092 pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
4093 pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( - tmp ) )
4095 tmp = min( xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp)
4096 pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
4097 pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( - tmp ) )
4099 pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt ) )
4100 pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( - xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt ) )
4101 pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt ) )
4102 pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( - xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt ) )
4107 profile_stop(
"sn14_freezing")
4128 integer,
intent(in) :: KA, KS, KE
4130 real(RP),
intent(in) :: RHOQ(ka,i_qc:i_ng)
4131 real(RP),
intent(in) :: DENS(ka)
4132 real(RP),
intent(in) :: TEMP(ka)
4133 real(RP),
intent(in) :: PRES(ka)
4135 real(RP),
intent(out) :: vterm(ka,
qa_mp-1)
4140 real(RP) :: rhofac_q(hydro_max)
4142 real(RP) :: rlambdar
4152 integer :: k, i, j, iq
4155 profile_start(
"sn14_terminal_vel")
4157 mud_r = 3.0_rp * nu(i_mp_qr) + 2.0_rp
4161 rhofac = rho_0 / max( dens(k), rho_min )
4164 rhofac_q(i_mp_qc) = rhofac ** gamma_v(i_mp_qc)
4165 xq = max( xqmin(i_mp_qc), min( xqmax(i_mp_qc), rhoq(k,i_qc) / ( rhoq(k,i_nc) + nqmin(i_mp_qc) ) ) )
4167 vterm(k,i_mp_qc) = -rhofac_q(i_mp_qc) * coef_vt1(i_mp_qc,1) * xq**beta_v(i_mp_qc,1)
4169 vterm(k,i_mp_nc) = -rhofac_q(i_mp_qc) * coef_vt0(i_mp_qc,1) * xq**beta_vn(i_mp_qc,1)
4172 rhofac_q(i_mp_qr) = rhofac ** gamma_v(i_mp_qr)
4173 xq = max( xqmin(i_mp_qr), min( xqmax(i_mp_qr), rhoq(k,i_qr) / ( rhoq(k,i_nr) + nqmin(i_mp_qr) ) ) )
4175 rlambdar = a_m(i_mp_qr) * xq**b_m(i_mp_qr) &
4176 * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
4177 dq = ( 4.0_rp + mud_r ) * rlambdar
4179 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
4181 velq_s = coef_vtr_ar2 * dq &
4182 * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-5.0_rp-mud_r) )
4183 velq_l = coef_vtr_ar1 - coef_vtr_br1 &
4184 * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-4.0_rp-mud_r)
4185 vterm(k,i_mp_qr) = -rhofac_q(i_mp_qr) &
4186 * ( velq_l * ( weight ) &
4187 + velq_s * ( 1.0_rp - weight ) )
4189 dq = ( 1.0_rp + mud_r ) * rlambdar
4190 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
4192 velq_s = coef_vtr_ar2 * dql &
4193 * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-2.0_rp-mud_r) )
4194 velq_l = coef_vtr_ar1 - coef_vtr_br1 &
4195 * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-1.0_rp-mud_r)
4196 vterm(k,i_mp_nr) = -rhofac_q(i_mp_qr) &
4197 * ( velq_l * ( weight ) &
4198 + velq_s * ( 1.0_rp - weight ) )
4201 rhofac_q(i_mp_qi) = ( pres(k)/pre0_vt )**a_pre0_vt * ( temp(k)/tem0_vt )**a_tem0_vt
4202 xq = max( xqmin(i_mp_qi), min( xqmax(i_mp_qi), rhoq(k,i_qi) / ( rhoq(k,i_ni) + nqmin(i_mp_qi) ) ) )
4204 tmp = a_m(i_mp_qi) * xq**b_m(i_mp_qi)
4205 dq = coef_dave_l(i_mp_qi) * tmp
4206 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_li ) ) ) )
4208 velq_s = coef_vt1(i_mp_qi,1) * xq**beta_v(i_mp_qi,1)
4209 velq_l = coef_vt1(i_mp_qi,2) * xq**beta_v(i_mp_qi,2)
4210 vterm(k,i_mp_qi) = -rhofac_q(i_mp_qi) &
4211 * ( velq_l * ( weight ) &
4212 + velq_s * ( 1.0_rp - weight ) )
4214 dq = coef_dave_n(i_mp_qi) * tmp
4215 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ni ) ) ) )
4217 velq_s = coef_vt0(i_mp_qi,1) * xq**beta_vn(i_mp_qi,1)
4218 velq_l = coef_vt0(i_mp_qi,2) * xq**beta_vn(i_mp_qi,2)
4219 vterm(k,i_mp_ni) = -rhofac_q(i_mp_qi) &
4220 * ( velq_l * ( weight ) &
4221 + velq_s * ( 1.0_rp - weight ) )
4224 rhofac_q(i_mp_qs) = rhofac_q(i_mp_qi)
4225 xq = max( xqmin(i_mp_qs), min( xqmax(i_mp_qs), rhoq(k,i_qs) / ( rhoq(k,i_ns) + nqmin(i_mp_qs) ) ) )
4227 tmp = a_m(i_mp_qs) * xq**b_m(i_mp_qs)
4228 dq = coef_dave_l(i_mp_qs) * tmp
4229 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ls ) ) ) )
4231 velq_s = coef_vt1(i_mp_qs,1) * xq**beta_v(i_mp_qs,1)
4232 velq_l = coef_vt1(i_mp_qs,2) * xq**beta_v(i_mp_qs,2)
4233 vterm(k,i_mp_qs) = -rhofac_q(i_mp_qs) &
4234 * ( velq_l * ( weight ) &
4235 + velq_s * ( 1.0_rp - weight ) )
4237 dq = coef_dave_n(i_mp_qs) * tmp
4238 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ns ) ) ) )
4240 velq_s = coef_vt0(i_mp_qs,1) * xq**beta_vn(i_mp_qs,1)
4241 velq_l = coef_vt0(i_mp_qs,2) * xq**beta_vn(i_mp_qs,2)
4242 vterm(k,i_mp_ns) = -rhofac_q(i_mp_qs) &
4243 * ( velq_l * ( weight ) &
4244 + velq_s * ( 1.0_rp - weight ) )
4247 rhofac_q(i_mp_qg) = rhofac_q(i_mp_qi)
4248 xq = max( xqmin(i_mp_qg), min( xqmax(i_mp_qg), rhoq(k,i_qg) / ( rhoq(k,i_ng) + nqmin(i_mp_qg) ) ) )
4250 tmp = a_m(i_mp_qg) * xq**b_m(i_mp_qg)
4251 dq = coef_dave_l(i_mp_qg) * tmp
4252 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_lg ) ) ) )
4254 velq_s = coef_vt1(i_mp_qg,1) * xq**beta_v(i_mp_qg,1)
4255 velq_l = coef_vt1(i_mp_qg,2) * xq**beta_v(i_mp_qg,2)
4256 vterm(k,i_mp_qg) = -rhofac_q(i_mp_qg) &
4257 * ( velq_l * ( weight ) &
4258 + velq_s * ( 1.0_rp - weight ) )
4260 dq = coef_dave_n(i_mp_qg) * tmp
4261 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ng ) ) ) )
4263 velq_s = coef_vt0(i_mp_qg,1) * xq**beta_vn(i_mp_qg,1)
4264 velq_l = coef_vt0(i_mp_qg,2) * xq**beta_vn(i_mp_qg,2)
4265 vterm(k,i_mp_ng) = -rhofac_q(i_mp_qg) &
4266 * ( velq_l * ( weight ) &
4267 + velq_s * ( 1.0_rp - weight ) )
4272 vterm(ks-1,iq) = vterm(ks,iq)
4276 profile_stop(
"sn14_terminal_vel")
4282 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
4283 ntdiv, ntmax, & ! in [Add] 10/08/03
4292 esw, esi, rhoq2, & ! in
4296 sl_PLCdep, & ! inout
4297 sl_PLRdep, sl_PNRdep, & ! inout
4312 moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
4313 moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
4314 moist_dqs_dtem_dens_liq => atmos_saturation_dqs_dtem_dens_liq, &
4315 moist_dqs_dtem_dens_ice => atmos_saturation_dqs_dtem_dens_ice, &
4316 moist_dqs_dtem_dpre_liq => atmos_saturation_dqs_dtem_dpre_liq, &
4317 moist_dqs_dtem_dpre_ice => atmos_saturation_dqs_dtem_dpre_ice
4320 integer,
intent(in) :: KA, KS, KE
4321 integer,
intent(in) :: IA, IS, IE
4322 integer,
intent(in) :: JA, JS, JE
4324 integer,
intent(in) :: ntdiv
4325 integer,
intent(in) :: ntmax
4327 real(RP),
intent(in) :: dt
4329 real(RP),
intent(in) :: cz(ka,ia,ja)
4330 real(RP),
intent(in) :: fz(ka,ia,ja)
4331 real(RP),
intent(in) :: w(ka,ia,ja)
4332 real(RP),
intent(in) :: dTdt_rad(ka,ia,ja)
4333 real(RP),
intent(in) :: rho(ka,ia,ja)
4334 real(RP),
intent(in) :: qdry(ka,ia,ja)
4335 real(RP),
intent(in) :: esw(ka,ia,ja)
4336 real(RP),
intent(in) :: esi(ka,ia,ja)
4337 real(RP),
intent(in) :: rhoq2(i_qv:i_ng,ka,ia,ja)
4339 real(RP),
intent(in) :: tem(ka,ia,ja)
4340 real(RP),
intent(in) :: pre(ka,ia,ja)
4341 real(RP),
intent(in) :: cva(ka,ia,ja)
4342 real(RP),
intent(in) :: cpa(ka,ia,ja)
4345 real(RP),
intent(inout) :: PQ(pq_max,ka,ia,ja)
4347 real(RP),
intent(inout) :: sl_PLCdep(ia,ja)
4348 real(RP),
intent(inout) :: sl_PLRdep(ia,ja), sl_PNRdep(ia,ja)
4350 real(RP),
intent(out) :: RHOQ_t(ka, ia, ja,
qa_mp)
4351 real(RP),
intent(out) :: RHOE_t(ka, ia, ja)
4352 real(RP),
intent(out) :: CPtot_t(ka,ia,ja)
4353 real(RP),
intent(out) :: CVtot_t(ka,ia,ja)
4356 real(RP),
intent(out) :: qc_evaporate(ka,ia,ja)
4360 real(RP) :: wtem(ka,ia,ja)
4365 real(RP) :: qsw(ka,ia,ja)
4366 real(RP) :: qsi(ka,ia,ja)
4367 real(RP) :: dqswdtem_rho(ka,ia,ja)
4368 real(RP) :: dqsidtem_rho(ka,ia,ja)
4369 real(RP) :: dqswdtem_pre(ka,ia,ja)
4370 real(RP) :: dqsidtem_pre(ka,ia,ja)
4371 real(RP) :: dqswdpre_tem(ka,ia,ja)
4372 real(RP) :: dqsidpre_tem(ka,ia,ja)
4377 real(RP) :: aliqliq, asolliq
4378 real(RP) :: aliqsol, asolsol
4383 real(RP) :: taucnd, r_taucnd
4384 real(RP) :: taudep, r_taudep
4385 real(RP) :: taucnd_c(ka,ia,ja), r_taucnd_c
4386 real(RP) :: taucnd_r(ka,ia,ja), r_taucnd_r
4387 real(RP) :: taudep_i(ka,ia,ja), r_taudep_i
4388 real(RP) :: taudep_s(ka,ia,ja), r_taudep_s
4389 real(RP) :: taudep_g(ka,ia,ja), r_taudep_g
4392 real(RP) :: PLR2NR, PLI2NI, PLS2NS, PLG2NG
4393 real(RP) :: coef_a_cnd, coef_b_cnd
4394 real(RP) :: coef_a_dep, coef_b_dep
4416 real(RP) :: r_xc_ccn, r_xi_ccn
4419 real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
4420 real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
4422 real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
4423 real(RP) :: r_rvaptem
4425 real(RP) :: lvsw, lvsi
4426 real(RP) :: dlvsw, dlvsi
4428 real(RP) :: dcnd, ddep
4429 real(RP) :: uplim_cnd
4430 real(RP) :: lowlim_cnd
4432 real(RP) :: uplim_dep
4433 real(RP) :: lowlim_dep
4434 real(RP) :: ssw, ssi
4435 real(RP) :: r_esw, r_esi
4436 real(RP) :: r_lvsw, r_lvsi
4438 real(RP) :: ssw_o, ssi_o
4444 real(RP),
save :: fac_cndc = 1.0_rp
4445 logical,
save :: opt_fix_taucnd_c=.false.
4446 logical,
save :: flag_first =.true.
4448 namelist / param_atmos_phy_mp_sn14_condensation / &
4449 opt_fix_taucnd_c, fac_cndc
4451 real(RP) :: fac_cndc_wrk
4453 real(RP),
parameter :: tau100day = 1.e+7_rp
4454 real(RP),
parameter :: r_tau100day = 1.e-7_rp
4455 real(RP),
parameter :: eps=1.e-30_rp
4459 integer :: i,j,k,iqw
4461 real(RP) :: dqv, dqc, dqr, dqi, dqs, dqg, dcv, dcp
4465 if( flag_first )
then 4466 flag_first = .false.
4468 read (
io_fid_conf,nml=param_atmos_phy_mp_sn14_condensation, end=100)
4469 100 log_nml(param_atmos_phy_mp_sn14_condensation)
4477 r_xc_ccn=1.0_rp/xc_ccn
4480 if( opt_fix_taucnd_c )
then 4481 fac_cndc_wrk = fac_cndc**(1.0_rp-b_m(i_mp_qc))
4485 pq(i_lcdep,k,i,j) = pq(i_lcdep,k,i,j)*fac_cndc_wrk
4489 log_info(
"ATMOS_PHY_MP_SN14_update_by_phase_change_kij",*)
"taucnd:fac_cndc_wrk=",fac_cndc_wrk
4498 wtem(k,i,j) = max( tem(k,i,j), tem_min )
4503 call moist_pres2qsat_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
4504 wtem(:,:,:), pre(:,:,:), qdry(:,:,:), &
4506 call moist_pres2qsat_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
4507 wtem(:,:,:), pre(:,:,:), qdry(:,:,:), &
4509 call moist_dqs_dtem_dens_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
4510 wtem(:,:,:), rho(:,:,:), &
4511 dqswdtem_rho(:,:,:) )
4512 call moist_dqs_dtem_dens_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
4513 wtem(:,:,:), rho(:,:,:), &
4514 dqsidtem_rho(:,:,:) )
4515 call moist_dqs_dtem_dpre_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
4516 wtem(:,:,:), pre(:,:,:), qdry(:,:,:), &
4517 dqswdtem_pre(:,:,:), dqswdpre_tem(:,:,:) )
4518 call moist_dqs_dtem_dpre_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
4519 wtem(:,:,:), pre(:,:,:), qdry(:,:,:), &
4520 dqsidtem_pre(:,:,:), dqsidpre_tem(:,:,:) )
4522 profile_start(
"sn14_update")
4526 if( cz(k,i,j) <= 25000.0_rp )
then 4531 if( pre(k,i,j) < esw(k,i,j)+1.e-10_rp )
then 4533 dqswdtem_rho(k,i,j) = 0.0_rp
4534 dqswdtem_pre(k,i,j) = 0.0_rp
4535 dqswdpre_tem(k,i,j) = 0.0_rp
4537 if( pre(k,i,j) < esi(k,i,j)+1.e-10_rp )
then 4539 dqsidtem_rho(k,i,j) = 0.0_rp
4540 dqsidtem_pre(k,i,j) = 0.0_rp
4541 dqsidpre_tem(k,i,j) = 0.0_rp
4544 r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4545 lvsw = esw(k,i,j)*r_rvaptem
4546 lvsi = esi(k,i,j)*r_rvaptem
4547 pv = rhoq2(i_qv,k,i,j)*rvap*tem(k,i,j)
4548 r_esw = 1.0_rp/esw(k,i,j)
4549 r_esi = 1.0_rp/esi(k,i,j)
4550 ssw = min( mp_ssw_lim, ( pv*r_esw-1.0_rp ) )
4551 ssi = pv*r_esi - 1.0_rp
4552 r_lvsw = 1.0_rp/lvsw
4553 r_lvsi = 1.0_rp/lvsi
4554 r_taucnd_c = pq(i_lcdep,k,i,j)*r_lvsw
4555 r_taucnd_r = pq(i_lrdep,k,i,j)*r_lvsw
4556 r_taudep_i = pq(i_lidep,k,i,j)*r_lvsi
4557 r_taudep_s = pq(i_lsdep,k,i,j)*r_lvsi
4558 r_taudep_g = pq(i_lgdep,k,i,j)*r_lvsi
4565 r_cva = 1.0_rp / cva(k,i,j)
4566 r_cpa = 1.0_rp / cpa(k,i,j)
4570 + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
4573 + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
4576 + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
4579 + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
4580 pdynliq = w2 * grav * ( r_cpa*dqswdtem_pre(k,i,j) + rho(k,i,j)*dqswdpre_tem(k,i,j) )
4581 pdynsol = w2 * grav * ( r_cpa*dqsidtem_pre(k,i,j) + rho(k,i,j)*dqsidpre_tem(k,i,j) )
4582 pradliq = -dtdt_rad(k,i,j) * dqswdtem_rho(k,i,j)
4583 pradsol = -dtdt_rad(k,i,j) * dqsidtem_rho(k,i,j)
4590 acnd = pdynliq + pradliq &
4591 - ( r_taudep_i+r_taudep_s+r_taudep_g ) * ( qsw(k,i,j) - qsi(k,i,j) )
4592 adep = pdynsol + pradsol &
4593 + ( r_taucnd_c+r_taucnd_r ) * ( qsw(k,i,j) - qsi(k,i,j) )
4595 + aliqliq*( r_taucnd_c+r_taucnd_r ) &
4596 + asolliq*( r_taudep_i+r_taudep_s+r_taudep_g )
4598 + aliqsol*( r_taucnd_c+r_taucnd_r )&
4599 + asolsol*( r_taudep_i+r_taudep_s+r_taudep_g )
4601 uplim_cnd = max( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
4602 lowlim_cnd = min( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
4603 if( r_taucnd < r_tau100day )
then 4605 pq(i_lcdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lcdep,k,i,j)*ssw_o ))
4606 pq(i_lrdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lrdep,k,i,j)*ssw_o ))
4607 pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_nrdep,k,i,j)*ssw_o )
4610 taucnd = 1.0_rp/r_taucnd
4612 coef_a_cnd = rho(k,i,j)*acnd*taucnd
4613 coef_b_cnd = rho(k,i,j)*taucnd*r_dt*(ssw_o*qsw(k,i,j)-acnd*taucnd) * ( exp(-dt*r_taucnd) - 1.0_rp )
4614 pq(i_lcdep,k,i,j) = coef_a_cnd*r_taucnd_c - coef_b_cnd*r_taucnd_c
4615 plr2nr = pq(i_nrdep,k,i,j)/(pq(i_lrdep,k,i,j)+1.e-30_rp)
4616 pq(i_lrdep,k,i,j) = coef_a_cnd*r_taucnd_r - coef_b_cnd*r_taucnd_r
4617 pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_lrdep,k,i,j)*plr2nr )
4620 uplim_dep = max( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4621 lowlim_dep = min( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4622 if( r_taudep < r_tau100day )
then 4624 pq(i_lidep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lidep,k,i,j)*ssi_o ))
4625 pq(i_lsdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lsdep,k,i,j)*ssi_o ))
4626 pq(i_lgdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lgdep,k,i,j)*ssi_o ))
4627 pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_nidep,k,i,j)*ssi_o )
4628 pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_nsdep,k,i,j)*ssi_o )
4629 pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_ngdep,k,i,j)*ssi_o )
4631 taudep = 1.0_rp/r_taudep
4633 coef_a_dep = rho(k,i,j)*adep*taudep
4634 coef_b_dep = rho(k,i,j)*taudep*r_dt*(ssi_o*qsi(k,i,j)-adep*taudep) * ( exp(-dt*r_taudep) - 1.0_rp )
4635 pli2ni = pq(i_nidep,k,i,j)/max(pq(i_lidep,k,i,j),1.e-30_rp)
4636 pls2ns = pq(i_nsdep,k,i,j)/max(pq(i_lsdep,k,i,j),1.e-30_rp)
4637 plg2ng = pq(i_ngdep,k,i,j)/max(pq(i_lgdep,k,i,j),1.e-30_rp)
4638 pq(i_lidep,k,i,j) = coef_a_dep*r_taudep_i - coef_b_dep*r_taudep_i
4639 pq(i_lsdep,k,i,j) = coef_a_dep*r_taudep_s - coef_b_dep*r_taudep_s
4640 pq(i_lgdep,k,i,j) = coef_a_dep*r_taudep_g - coef_b_dep*r_taudep_g
4641 pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_lidep,k,i,j)*pli2ni )
4642 pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_lsdep,k,i,j)*pls2ns )
4643 pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_lgdep,k,i,j)*plg2ng )
4646 sw = 0.5_rp - sign(0.5_rp, pq(i_lcdep,k,i,j)+eps)
4647 pncdep = min(0.0_rp, ((rhoq2(i_qc,k,i,j)+pq(i_lcdep,k,i,j)*dt)*r_xc_ccn - rhoq2(i_nc,k,i,j))*r_dt ) * sw
4660 r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4661 lvsw = esw(k,i,j)*r_rvaptem
4662 dlvsw = rhoq2(i_qv,k,i,j)-lvsw
4663 dcnd = dt*(pq(i_lcdep,k,i,j)+pq(i_lrdep,k,i,j))
4665 sw = ( sign(0.5_rp,dcnd) + sign(0.5_rp,dlvsw) ) &
4666 * ( 0.5_rp + sign(0.5_rp,abs(dcnd)-eps) )
4670 fac1 = min(dlvsw*sw,dcnd*sw)*sw / (abs(sw)-1.0_rp+dcnd) &
4672 dep_dqc = max( dt*pq(i_lcdep,k,i,j)*fac1, &
4673 -rhoq2(i_qc,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw)
4674 dep_dqr = max( dt*pq(i_lrdep,k,i,j)*fac1, &
4675 -rhoq2(i_qr,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw)
4694 dep_dnc = max( dt*pncdep*fac1, -rhoq2(i_nc,k,i,j) )
4695 dep_dnr = max( dt*pq(i_nrdep,k,i,j)*fac1, -rhoq2(i_nr,k,i,j) )
4697 qc_evaporate(k,i,j) = - dep_dnc
4700 lvsi = esi(k,i,j)*r_rvaptem
4701 ddep = dt*(pq(i_lidep,k,i,j)+pq(i_lsdep,k,i,j)+pq(i_lgdep,k,i,j))
4702 dlvsi = rhoq2(i_qv,k,i,j)-lvsi
4704 sw = ( sign(0.5_rp,ddep) + sign(0.5_rp,dlvsi) ) &
4705 * ( 0.5_rp + sign(0.5_rp,abs(ddep)-eps) )
4709 fac2 = min(dlvsi*sw,ddep*sw)*sw / (abs(sw)-1.0_rp+ddep) &
4711 dep_dqi = max( dt*pq(i_lidep,k,i,j) &
4712 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
4713 -rhoq2(i_qi,k,i,j) - 1e30_rp*(sw+1.0_rp) )
4714 dep_dqs = max( dt*pq(i_lsdep,k,i,j) &
4715 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
4716 -rhoq2(i_qs,k,i,j) - 1e30_rp*(sw+1.0_rp) )
4717 dep_dqg = max( dt*pq(i_lgdep,k,i,j) &
4718 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
4719 -rhoq2(i_qg,k,i,j) - 1e30_rp*(sw+1.0_rp) )
4741 dep_dni = max( dt*pq(i_nidep,k,i,j)*fac2, -rhoq2(i_ni,k,i,j) )
4742 dep_dns = max( dt*pq(i_nsdep,k,i,j)*fac2, -rhoq2(i_ns,k,i,j) )
4743 dep_dng = max( dt*pq(i_ngdep,k,i,j)*fac2, -rhoq2(i_ng,k,i,j) )
4746 frz_dqc = max( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j)), -rhoq2(i_qc,k,i,j)-dep_dqc )
4747 frz_dnc = max( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j)), -rhoq2(i_nc,k,i,j)-dep_dnc )
4748 fac3 = ( frz_dqc-eps )/( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j))-eps )
4749 fac4 = ( frz_dnc-eps )/( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j))-eps )
4750 pq(i_lchom,k,i,j) = fac3*pq(i_lchom,k,i,j)
4751 pq(i_lchet,k,i,j) = fac3*pq(i_lchet,k,i,j)
4752 pq(i_nchom,k,i,j) = fac4*pq(i_nchom,k,i,j)
4753 pq(i_nchet,k,i,j) = fac4*pq(i_nchet,k,i,j)
4757 mlt_dqi = max( dt*pq(i_limlt,k,i,j), -rhoq2(i_qi,k,i,j)-dep_dqi )
4758 mlt_dni = max( dt*pq(i_nimlt,k,i,j), -rhoq2(i_ni,k,i,j)-dep_dni )
4760 mlt_dqs = max( dt*pq(i_lsmlt,k,i,j), -rhoq2(i_qs,k,i,j)-dep_dqs )
4761 mlt_dns = max( dt*pq(i_nsmlt,k,i,j), -rhoq2(i_ns,k,i,j)-dep_dns )
4763 mlt_dqg = max( dt*pq(i_lgmlt,k,i,j), -rhoq2(i_qg,k,i,j)-dep_dqg )
4764 mlt_dng = max( dt*pq(i_ngmlt,k,i,j), -rhoq2(i_ng,k,i,j)-dep_dng )
4767 frz_dqr = max( dt*(pq(i_lrhet,k,i,j)), min(0.0_rp, -rhoq2(i_qr,k,i,j)-dep_dqr) )
4768 frz_dnr = max( dt*(pq(i_nrhet,k,i,j)), min(0.0_rp, -rhoq2(i_nr,k,i,j)-dep_dnr) )
4770 fac5 = ( frz_dqr-eps )/( dt*pq(i_lrhet,k,i,j)-eps )
4771 pq(i_lrhet,k,i,j) = fac5*pq(i_lrhet,k,i,j)
4772 fac6 = ( frz_dnr-eps )/( dt*pq(i_nrhet,k,i,j)-eps )
4773 pq(i_nrhet,k,i,j) = fac6*pq(i_nrhet,k,i,j)
4776 drhoqv = -(dep_dqc+dep_dqi+dep_dqs+dep_dqg+dep_dqr)
4778 xi = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ))
4779 sw = 0.5_rp + sign(0.5_rp,xi-x_sep)
4783 drhoqc = ( frz_dqc - mlt_dqi*(1.0_rp-sw) + dep_dqc )
4784 drhonc = ( frz_dnc - mlt_dni*(1.0_rp-sw) + dep_dnc )
4786 drhoqr = ( frz_dqr - mlt_dqg - mlt_dqs - mlt_dqi*sw + dep_dqr )
4787 drhonr = ( frz_dnr - mlt_dng - mlt_dns - mlt_dni*sw + dep_dnr )
4790 drhoqi = (-frz_dqc + mlt_dqi + dep_dqi )
4791 drhoni = (-frz_dnc + mlt_dni + dep_dni )
4794 drhoqs = ( mlt_dqs + dep_dqs )
4795 drhons = ( mlt_dns + dep_dns )
4798 drhoqg = (-frz_dqr + mlt_dqg + dep_dqg )
4799 drhong = (-frz_dnr + mlt_dng + dep_dng )
4802 rhoq_t(k,i,j,i_qv) = drhoqv / dt
4803 rhoq_t(k,i,j,i_qc) = drhoqc / dt
4804 rhoq_t(k,i,j,i_nc) = drhonc / dt
4805 rhoq_t(k,i,j,i_qr) = drhoqr / dt
4806 rhoq_t(k,i,j,i_nr) = drhonr / dt
4807 rhoq_t(k,i,j,i_qi) = drhoqi / dt
4808 rhoq_t(k,i,j,i_ni) = drhoni / dt
4809 rhoq_t(k,i,j,i_qs) = drhoqs / dt
4810 rhoq_t(k,i,j,i_ns) = drhons / dt
4811 rhoq_t(k,i,j,i_qg) = drhoqg / dt
4812 rhoq_t(k,i,j,i_ng) = drhong / dt
4814 rhoe_t(k,i,j) = ( - lhv * drhoqv + lhf * ( drhoqi+ drhoqs + drhoqg ) ) / dt
4816 rrho = 1.0_rp/rho(k,i,j)
4827 cvtot_t(k,i,j) = dcv/dt
4828 cptot_t(k,i,j) = dcp/dt
4830 dz = fz(k,i,j) - fz(k,i,j)
4831 sl_plcdep(i,j) = sl_plcdep(i,j) + dep_dqc*dz
4832 sl_plrdep(i,j) = sl_plrdep(i,j) + dep_dqr*dz
4833 sl_pnrdep(i,j) = sl_pnrdep(i,j) + dep_dnr*dz
4837 profile_stop(
"sn14_update")
character(len=h_short), dimension(qa_mp), parameter, public atmos_phy_mp_sn14_tracer_names
subroutine, public atmos_phy_mp_sn14_terminal_velocity(KA, KS, KE, DENS, TEMP, RHOQ, PRES, vterm)
ATMOS_PHY_MP_sn14_terminal_velocity Calculate terminal velocity.
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
real(rp), parameter, public const_psat0
saturate pressure of water vapor at 0C [Pa]
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
module atmosphere / saturation
module ATMOSPHERE / Physics Cloud Microphysics
real(rp), public cv_ice
CV for ice [J/kg/K].
real(rp), parameter, public const_ci
specific heat (ice) [J/kg/K]
subroutine ice_multiplication_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, Pac, tem, rhoq, xq, PQ)
real(rp), parameter, public const_dwatr
density of water [kg/m3]
real(rp), public cp_ice
CP for ice [J/kg/K].
integer, parameter, public i_hs
snow
real(rp) function, public sf_gamma(x)
Gamma function.
subroutine, public atmos_phy_mp_sn14_cloud_fraction(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC, mask_criterion, cldfrac)
ATMOS_PHY_MP_sn14_cloud_fraction Calculate Cloud Fraction.
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
subroutine mixed_phase_collection_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, wtem, rhoq, xq, dq_xave, vt_xave, PQ, Pac)
integer, parameter, public i_hr
liquid water rain
integer, parameter, public i_hi
ice water cloud
subroutine, public atmos_phy_mp_sn14_setup(KA, IA, JA)
ATMOS_PHY_MP_sn14_setup setup.
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
integer, public io_fid_conf
Config file ID.
real(rp), public const_cvvap
specific heat (water vapor, constant volume) [J/kg/K]
integer, parameter, public i_hh
hail
real(rp), public const_lhf0
latent heat of fusion at 0C [J/kg]
subroutine, public atmos_phy_mp_sn14_tendency(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, W, QTRC, PRES, TEMP, Qdry, CPtot, CVtot, CCN, dt, cz, fz, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, EVAPORATE)
ATMOS_PHY_MP_sn14_tendency calculate tendency.
subroutine, public atmos_phy_mp_sn14_qtrc2nhyd(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, Ne)
Calculate number concentration of each category.
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
real(rp), public const_undef
character(len=h_short), dimension(qa_mp), parameter, public atmos_phy_mp_sn14_tracer_units
subroutine, public atmos_phy_mp_sn14_qtrc2qhyd(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, Qe)
ATMOS_PHY_MP_sn14_qtrc2qhyd Calculate mass ratio of each category.
real(rp), parameter, public const_lhs0
latent heat of sublimation at 0C [J/kg]
real(rp), public cv_vapor
CV for vapor [J/kg/K].
module atmosphere / hydrometeor
integer, parameter, public atmos_phy_mp_sn14_ntracers
character(len=h_mid), dimension(qa_mp), parameter, public atmos_phy_mp_sn14_tracer_descriptions
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
real(rp), public const_pre00
pressure reference [Pa]
integer, parameter, public atmos_phy_mp_sn14_nwaters
integer, parameter, public qa_mp
subroutine nucleation_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, cz, fz, w, rho, tem, pre, qdry, rhoq, cpa, dTdt_rad, qke, CCN, dt, PQ)
real(rp), public const_lhf00
latent heat of fusion at 0K [J/kg]
real(rp), public const_grav
standard acceleration of gravity [m/s2]
subroutine freezing_water_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, dt, rhoq, xq, tem, PQ)
real(rp), public const_lhv00
latent heat of vaporizaion at 0K [J/kg]
integer, public prc_myrank
process num in local communicator
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
subroutine, public prc_abort
Abort Process.
integer, parameter, public i_hc
liquid water cloud
integer, parameter, public atmos_phy_mp_sn14_nices
subroutine debug_tem_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, point, tem, rho, pre, qv)
real(rp), parameter, public const_emelt
subroutine, public atmos_phy_mp_sn14_qhyd2qtrc(KA, KS, KE, IA, IS, IE, JA, JS, JE, Qe, QTRC, QNUM)
subroutine aut_acc_slc_brk_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, rhoq, xq, dq_xave, rho, PQ)
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
real(dp), parameter, public const_undef8
undefined value (REAL8)
real(rp), public const_eps
small number
subroutine, public atmos_phy_mp_sn14_effective_radius(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS0, TEMP0, QTRC0, Re)
ATMOS_PHY_MP_sn14_effective_radius Calculate Effective Radius.
real(rp), public const_pi
pi
integer, parameter, public n_hyd
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
real(rp), public cp_vapor
CP for vapor [J/kg/K].
real(rp), public cp_water
CP for water [J/kg/K].
integer, parameter, public i_hg
graupel
subroutine update_by_phase_change_kij(KA, KS, KE, IA, IS, IE, JA, JS, JE, ntdiv, ntmax, dt, cz, fz, w, dTdt_rad, rho, qdry, esw, esi, rhoq2, pre, tem, cpa, cva, PQ, sl_PLCdep, sl_PLRdep, sl_PNRdep, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, qc_evaporate)
real(rp), public cv_water
CV for water [J/kg/K].