55 #define PROFILE_START(name) call fapp_start(name, 1, 1) 56 #define PROFILE_STOP(name) call fapp_stop (name, 1, 1) 57 #elif defined(PROFILE_FINEPA) 58 #define PROFILE_START(name) call start_collection(name) 59 #define PROFILE_STOP(name) call stop_collection (name) 61 #define PROFILE_START(name) 62 #define PROFILE_STOP(name) 65 #include "macro_thermodyn.h" 137 integer,
public,
parameter ::
qa_mp = 11
159 'Ratio of Water Vapor mass to total mass (Specific humidity)', &
160 'Ratio of Cloud Water mass to total mass', &
161 'Ratio of Rain Water mass to total mass', &
162 'Ratio of Cloud Ice mass to total mass', &
163 'Ratio of Snow mass to total mass', &
164 'Ratio of Graupel mass to total mass', &
165 'Cloud Water Number Density', &
166 'Rain Water Number Density', &
167 'Cloud Ice Number Density', &
168 'Snow Number Density', &
169 'Graupel Number Density' /
188 private :: mp_sn14_init
190 private :: mp_terminal_velocity
197 integer,
private,
parameter :: hydro_max = 5
199 integer,
private,
parameter :: i_mp_qc = 1
200 integer,
private,
parameter :: i_mp_qr = 2
201 integer,
private,
parameter :: i_mp_qi = 3
202 integer,
private,
parameter :: i_mp_qs = 4
203 integer,
private,
parameter :: i_mp_qg = 5
204 integer,
private,
parameter :: i_mp_nc = 6
205 integer,
private,
parameter :: i_mp_nr = 7
206 integer,
private,
parameter :: i_mp_ni = 8
207 integer,
private,
parameter :: i_mp_ns = 9
208 integer,
private,
parameter :: i_mp_ng = 10
210 integer,
private :: qs_mp
211 integer,
private :: qe_mp
215 integer,
parameter :: i_lcccn = 1
216 integer,
parameter :: i_ncccn = 2
217 integer,
parameter :: i_liccn = 3
218 integer,
parameter :: i_niccn = 4
220 integer,
parameter :: i_lchom = 5
221 integer,
parameter :: i_nchom = 6
222 integer,
parameter :: i_lchet = 7
223 integer,
parameter :: i_nchet = 8
224 integer,
parameter :: i_lrhet = 9
225 integer,
parameter :: i_nrhet = 10
227 integer,
parameter :: i_limlt = 11
228 integer,
parameter :: i_nimlt = 12
229 integer,
parameter :: i_lsmlt = 13
230 integer,
parameter :: i_nsmlt = 14
231 integer,
parameter :: i_lgmlt = 15
232 integer,
parameter :: i_ngmlt = 16
234 integer,
parameter :: i_lrdep = 17
235 integer,
parameter :: i_nrdep = 18
236 integer,
parameter :: i_lidep = 19
237 integer,
parameter :: i_nidep = 20
238 integer,
parameter :: i_lsdep = 21
239 integer,
parameter :: i_nsdep = 22
240 integer,
parameter :: i_lgdep = 23
241 integer,
parameter :: i_ngdep = 24
242 integer,
parameter :: i_lcdep = 25
246 integer,
parameter :: i_lcaut = 26
247 integer,
parameter :: i_ncaut = 27
248 integer,
parameter :: i_nraut = 28
250 integer,
parameter :: i_lcacc = 29
251 integer,
parameter :: i_ncacc = 30
253 integer,
parameter :: i_nrslc = 31
254 integer,
parameter :: i_nrbrk = 32
257 integer,
parameter :: i_licon = 33
258 integer,
parameter :: i_nicon = 34
259 integer,
parameter :: i_lscon = 35
260 integer,
parameter :: i_nscon = 36
262 integer,
parameter :: i_liacm = 37
263 integer,
parameter :: i_niacm = 38
264 integer,
parameter :: i_liarm = 39
265 integer,
parameter :: i_niarm = 40
266 integer,
parameter :: i_lsacm = 41
267 integer,
parameter :: i_nsacm = 42
268 integer,
parameter :: i_lsarm = 43
269 integer,
parameter :: i_nsarm = 44
270 integer,
parameter :: i_lgacm = 45
271 integer,
parameter :: i_ngacm = 46
272 integer,
parameter :: i_lgarm = 47
273 integer,
parameter :: i_ngarm = 48
275 integer,
parameter :: i_lgspl = 49
276 integer,
parameter :: i_lsspl = 50
277 integer,
parameter :: i_nispl = 51
279 integer,
parameter :: pq_max = 51
283 integer,
parameter :: i_liaclc2li = 1
284 integer,
parameter :: i_niacnc2ni = 2
285 integer,
parameter :: i_lsaclc2ls = 3
286 integer,
parameter :: i_nsacnc2ns = 4
287 integer,
parameter :: i_lgaclc2lg = 5
288 integer,
parameter :: i_ngacnc2ng = 6
289 integer,
parameter :: i_lracli2lg_i = 7
290 integer,
parameter :: i_nracni2ng_i = 8
291 integer,
parameter :: i_lracli2lg_r = 9
292 integer,
parameter :: i_nracni2ng_r = 10
293 integer,
parameter :: i_lracls2lg_s = 11
294 integer,
parameter :: i_nracns2ng_s = 12
295 integer,
parameter :: i_lracls2lg_r = 13
296 integer,
parameter :: i_nracns2ng_r = 14
297 integer,
parameter :: i_lraclg2lg = 15
298 integer,
parameter :: i_nracng2ng = 16
299 integer,
parameter :: i_liacli2ls = 17
300 integer,
parameter :: i_niacni2ns = 18
301 integer,
parameter :: i_liacls2ls = 19
302 integer,
parameter :: i_niacns2ns = 20
303 integer,
parameter :: i_nsacns2ns = 21
304 integer,
parameter :: i_ngacng2ng = 22
305 integer,
parameter :: i_lgacls2lg = 23
306 integer,
parameter :: i_ngacns2ng = 24
308 integer,
parameter :: pac_max = 24
312 character(len=H_SHORT),
save :: wlabel(hydro_max)
315 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 /)
318 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 /)
321 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 /)
327 real(RP),
private,
parameter :: xc_min = 4.20e-15_rp
328 real(RP),
private,
parameter :: xr_min = 2.60e-10_rp
329 real(RP),
private,
parameter :: xi_min = 3.382e-13_rp
330 real(RP),
private,
parameter :: xs_min = 1.847e-12_rp
331 real(RP),
private,
parameter :: xg_min = 1.230e-10_rp
333 real(RP),
private,
parameter :: xc_max = 2.6e-10_rp
334 real(RP),
private,
parameter :: xr_max = 5.00e-6_rp
335 real(RP),
private,
parameter :: xi_max = 1.377e-6_rp
336 real(RP),
private,
parameter :: xs_max = 7.519e-6_rp
337 real(RP),
private,
parameter :: xg_max = 4.900e-5_rp
339 real(RP),
private,
parameter :: xmin_filter= xc_min
341 real(RP),
private,
parameter :: rmin_re= 1.e-6_rp
344 real(RP),
private,
parameter :: n0r_min= 2.5e+5_rp
345 real(RP),
private,
parameter :: n0r_max= 2.0e+7_rp
346 real(RP),
private,
parameter :: lambdar_min= 1.e+3_rp
347 real(RP),
private,
parameter :: lambdar_max= 1.e+4_rp
349 real(RP),
private,
parameter :: nc_min = 1.e+4_rp
350 real(RP),
private,
parameter :: nr_min = 1.0_rp
351 real(RP),
private,
parameter :: ni_min = 1.0_rp
352 real(RP),
private,
parameter :: ns_min = 1.e-4_rp
353 real(RP),
private,
parameter :: ng_min = 1.e-4_rp
355 real(RP),
private,
parameter :: lc_min = xc_min*nc_min
356 real(RP),
private,
parameter :: lr_min = xr_min*nr_min
357 real(RP),
private,
parameter :: li_min = xi_min*ni_min
358 real(RP),
private,
parameter :: ls_min = xs_min*ns_min
359 real(RP),
private,
parameter :: lg_min = xg_min*ng_min
361 real(RP),
private,
parameter :: x_sep = 2.6e-10_rp
363 real(RP),
private,
parameter :: tem_min=100.0_rp
364 real(RP),
private,
parameter :: rho_min=1.e-5_rp
365 real(RP),
private,
parameter :: rhoi = 916.70_rp
367 integer,
private,
save :: ntmax_phase_change = 1
368 integer,
private,
save :: ntmax_collection = 1
369 integer,
private,
save :: mp_ntmax_sedimentation= 1
372 real(RP),
private,
parameter :: rho_0 = 1.280_rp
374 real(RP),
allocatable,
private,
save :: nc_uplim_d(:,:,:)
377 real(RP),
private,
parameter :: ka0 = 2.428e-2_rp
379 real(RP),
private,
parameter :: dka_dt = 7.47e-5_rp
384 real(RP),
private,
parameter :: mua0 = 1.718e-5_rp
386 real(RP),
private,
parameter :: dmua_dt = 5.28e-8_rp
390 real(RP),
private,
save :: xc_ccn = 1.e-12_rp
391 real(RP),
private,
save :: xi_ccn = 1.e-12_rp
395 real(RP),
private,
save :: cap(hydro_max)
399 real(RP),
private,
save :: a_m(hydro_max)
400 real(RP),
private,
save :: b_m(hydro_max)
403 real(RP),
private,
save :: alpha_v(hydro_max,2)
404 real(RP),
private,
save :: beta_v(hydro_max,2)
405 real(RP),
private,
save :: alpha_vn(hydro_max,2)
406 real(RP),
private,
save :: beta_vn(hydro_max,2)
407 real(RP),
private,
save :: gamma_v(hydro_max)
410 real(RP),
private,
parameter :: pre0_vt = 300.e+2_rp
411 real(RP),
private,
parameter :: tem0_vt = 233.0_rp
412 real(RP),
private,
parameter :: a_pre0_vt = -0.1780_rp
413 real(RP),
private,
parameter :: a_tem0_vt = -0.3940_rp
421 real(RP),
private,
save :: nu(hydro_max)
422 real(RP),
private,
save :: mu(hydro_max)
428 real(RP),
private,
save :: a_area(hydro_max)
429 real(RP),
private,
save :: b_area(hydro_max)
430 real(RP),
private,
save :: ax_area(hydro_max)
431 real(RP),
private,
save :: bx_area(hydro_max)
434 real(RP),
private,
save :: a_rea(hydro_max)
435 real(RP),
private,
save :: b_rea(hydro_max)
436 real(RP),
private,
save :: a_rea2(hydro_max)
437 real(RP),
private,
save :: b_rea2(hydro_max)
438 real(RP),
private,
save :: a_rea3(hydro_max)
439 real(RP),
private,
save :: b_rea3(hydro_max)
441 real(RP),
private,
save :: a_d2vt(hydro_max)
442 real(RP),
private,
save :: b_d2vt(hydro_max)
446 real(RP),
private,
save :: coef_m2(hydro_max)
448 real(RP),
private,
save :: coef_d6(hydro_max)
450 real(RP),
private,
save :: coef_d3(hydro_max)
452 real(RP),
private,
save :: coef_d(hydro_max)
454 real(RP),
private,
save :: coef_d2v(hydro_max)
456 real(RP),
private,
save :: coef_md2v(hydro_max)
459 real(RP),
private,
save :: coef_r2(hydro_max)
460 real(RP),
private,
save :: coef_r3(hydro_max)
461 real(RP),
private,
save :: coef_re(hydro_max)
463 real(RP),
private,
save :: coef_rea2(hydro_max)
464 real(RP),
private,
save :: coef_rea3(hydro_max)
465 logical,
private,
save :: opt_m96_ice=.true.
466 logical,
private,
save :: opt_m96_column_ice=.false.
471 real(RP),
private,
save :: coef_vt0(hydro_max,2)
472 real(RP),
private,
save :: coef_vt1(hydro_max,2)
473 real(RP),
private,
save :: coef_deplc
474 real(RP),
private,
save :: coef_dave_n(hydro_max)
475 real(RP),
private,
save :: coef_dave_l(hydro_max)
478 real(RP),
private,
save :: d0_ni=261.76e-6_rp
479 real(RP),
private,
save :: d0_li=398.54e-6_rp
480 real(RP),
private,
parameter :: d0_ns=270.03e-6_rp
481 real(RP),
private,
parameter :: d0_ls=397.47e-6_rp
482 real(RP),
private,
parameter :: d0_ng=269.08e-6_rp
483 real(RP),
private,
parameter :: d0_lg=376.36e-6_rp
485 real(RP),
private,
parameter :: coef_vtr_ar1=9.65_rp
487 real(RP),
private,
parameter :: coef_vtr_br1=10.43_rp
488 real(RP),
private,
parameter :: coef_vtr_cr1=600.0_rp
489 real(RP),
private,
parameter :: coef_vtr_ar2=4.e+3_rp
490 real(RP),
private,
parameter :: coef_vtr_br2=12.e+3_rp
491 real(RP),
private,
parameter :: d_vtr_branch=0.745e-3_rp
493 real(RP),
private,
parameter :: dr_eq = 1.10e-3_rp
498 real(RP),
private,
save :: coef_a(hydro_max)
499 real(RP),
private,
save :: coef_lambda(hydro_max)
503 real(RP),
private,
save :: ah_vent (hydro_max,2)
504 real(RP),
private,
save :: bh_vent (hydro_max,2)
505 real(RP),
private,
save :: ah_vent0 (hydro_max,2)
506 real(RP),
private,
save :: bh_vent0 (hydro_max,2)
507 real(RP),
private,
save :: ah_vent1 (hydro_max,2)
508 real(RP),
private,
save :: bh_vent1 (hydro_max,2)
510 real(RP),
private,
save :: delta_b0 (hydro_max)
511 real(RP),
private,
save :: delta_b1 (hydro_max)
512 real(RP),
private,
save :: delta_ab0(hydro_max,hydro_max)
513 real(RP),
private,
save :: delta_ab1(hydro_max,hydro_max)
515 real(RP),
private,
save :: theta_b0 (hydro_max)
516 real(RP),
private,
save :: theta_b1 (hydro_max)
517 real(RP),
private,
save :: theta_ab0(hydro_max,hydro_max)
518 real(RP),
private,
save :: theta_ab1(hydro_max,hydro_max)
520 logical,
private,
save :: opt_debug=.false.
522 logical,
private,
save :: opt_debug_tem=.false.
523 logical,
private,
save :: opt_debug_inc=.true.
524 logical,
private,
save :: opt_debug_act=.true.
525 logical,
private,
save :: opt_debug_ree=.true.
526 logical,
private,
save :: opt_debug_bcs=.true.
528 integer,
private,
save :: mp_nstep_sedimentation
529 real(RP),
private,
save :: mp_rnstep_sedimentation
530 real(DP),
private,
save :: mp_dtsec_sedimentation
536 real(RP),
private,
allocatable,
save :: gsgam2_d (:,:,:)
537 real(RP),
private,
allocatable,
save :: gsgam2h_d(:,:,:)
538 real(RP),
private,
allocatable,
save :: gam2_d (:,:,:)
539 real(RP),
private,
allocatable,
save :: gam2h_d (:,:,:)
540 real(RP),
private,
allocatable,
save :: rgsgam2_d(:,:,:)
541 real(RP),
private,
allocatable,
save :: rgs_d (:,:,:)
542 real(RP),
private,
allocatable,
save :: rgsh_d (:,:,:)
544 logical,
private,
save :: mp_doautoconversion = .true.
545 logical,
private,
save :: mp_doprecipitation = .true.
546 logical,
private,
save :: mp_couple_aerosol = .false.
547 real(RP),
private,
save :: mp_ssw_lim = 1.e+1_rp
565 character(len=*),
intent(in) :: mp_type
566 integer,
intent(out) :: qa
567 integer,
intent(out) :: qs
573 if(
io_l )
write(
io_fid_log,*)
'++++++ Module[Cloud Microphysics Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]' 574 if(
io_l )
write(
io_fid_log,*)
'*** Tracers for Seiki and Nakajima (2014) 2-moment bulk 6 category' 576 if ( mp_type /=
'SN14' )
then 577 write(*,*)
'xxx ATMOS_PHY_MP_TYPE is not SN14. Check!' 595 qe_mp = qs_mp +
qa_mp - 1
627 namelist / param_atmos_phy_mp / &
628 mp_doautoconversion, &
629 mp_doprecipitation, &
632 mp_ntmax_sedimentation
634 real(RP),
parameter :: max_term_vel = 10.0_rp
640 if(
io_l )
write(
io_fid_log,*)
'++++++ Module[Cloud Microphysics] / Categ[ATMOS PHYSICS] / Origin[SCALElib]' 641 if(
io_l )
write(
io_fid_log,*)
'*** Seiki and Nakajima (2014) 2-moment bulk 6 category' 645 read(
io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
647 if(
io_l )
write(
io_fid_log,*)
'*** Not found namelist. Default used.' 648 elseif( ierr > 0 )
then 649 write(*,*)
'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!' 665 wlabel(5) =
"GRAUPEL" 669 allocate(nc_uplim_d(1,
ia,
ja))
670 nc_uplim_d(:,:,:) = 150.e6_rp
673 mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
675 mp_nstep_sedimentation = mp_ntmax_sedimentation
676 mp_rnstep_sedimentation = 1.0_rp /
real(mp_ntmax_sedimentation,kind=
rp)
680 if(
io_l )
write(
io_fid_log,*)
'*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation,
' step(s)' 681 if(
io_l )
write(
io_fid_log,*)
'*** DT of sedimentation [s] : ', mp_dtsec_sedimentation
684 allocate( gsgam2_d(
ka,
ia,
ja) )
685 allocate( gsgam2h_d(
ka,
ia,
ja) )
686 allocate( gam2_d(
ka,
ia,
ja) )
687 allocate( gam2h_d(
ka,
ia,
ja) )
688 allocate( rgsgam2_d(
ka,
ia,
ja) )
689 allocate( rgs_d(
ka,
ia,
ja) )
690 allocate( rgsh_d(
ka,
ia,
ja) )
691 gsgam2_d(:,:,:) = 1.0_rp
692 gsgam2h_d(:,:,:) = 1.0_rp
693 gam2_d(:,:,:) = 1.0_rp
694 gam2h_d(:,:,:) = 1.0_rp
695 rgsgam2_d(:,:,:) = 1.0_rp
696 rgs_d(:,:,:) = 1.0_rp
697 rgsh_d(:,:,:) = 1.0_rp
721 real(RP),
intent(inout) :: dens (
ka,
ia,
ja)
722 real(RP),
intent(inout) :: momz (
ka,
ia,
ja)
723 real(RP),
intent(inout) :: momx (
ka,
ia,
ja)
724 real(RP),
intent(inout) :: momy (
ka,
ia,
ja)
725 real(RP),
intent(inout) :: rhot (
ka,
ia,
ja)
726 real(RP),
intent(inout) :: qtrc (
ka,
ia,
ja,
qa)
727 real(RP),
intent(in) :: ccn (
ka,
ia,
ja)
728 real(RP),
intent(out) :: evaporate(
ka,
ia,
ja)
729 real(RP),
intent(out) :: sflx_rain(
ia,
ja)
730 real(RP),
intent(out) :: sflx_snow(
ia,
ja)
733 if(
io_l )
write(
io_fid_log,*)
'*** Atmos physics step: Cloud microphysics(SN14)' 741 call mp_sn14( dens, &
762 subroutine mp_sn14_init
769 real(RP) :: w1(hydro_max)
770 real(RP) :: w2(hydro_max)
771 real(RP) :: w3(hydro_max)
772 real(RP) :: w4(hydro_max)
773 real(RP) :: w5(hydro_max)
774 real(RP) :: w6(hydro_max)
775 real(RP) :: w7(hydro_max)
776 real(RP) :: w8(hydro_max)
779 real(RP) :: ar_ice_fix = 0.7_rp
780 real(RP) :: wcap1, wcap2
782 logical :: flag_vent0(hydro_max), flag_vent1(hydro_max)
784 integer :: iw,
ia, ib
787 namelist /nm_mp_sn14_init/ &
794 ntmax_phase_change, &
797 namelist /nm_mp_sn14_particles/ &
798 a_m, b_m, alpha_v, beta_v, gamma_v, &
800 a_area, b_area, cap, &
802 opt_m96_column_ice, &
805 real(RP),
parameter :: eps_gamma=1.e-30_rp
809 alpha_v(:,:) = undef8
811 alpha_vn(:,:) = undef8
812 beta_vn(:,:) = undef8
830 coef_dave_n(:) = undef8
831 coef_dave_l(:) = undef8
836 coef_md2v(:) = undef8
840 coef_rea2(:) = undef8
841 coef_rea3(:) = undef8
844 coef_lambda(:) = undef8
845 coef_vt0(:,:) = undef8
846 coef_vt1(:,:) = undef8
849 delta_ab0(:,:) = undef8
850 delta_ab1(:,:) = undef8
853 theta_ab0(:,:) = undef8
854 theta_ab1(:,:) = undef8
856 ah_vent(:,:) = undef8
857 ah_vent0(:,:) = undef8
858 ah_vent1(:,:) = undef8
859 bh_vent(:,:) = undef8
860 bh_vent0(:,:) = undef8
861 bh_vent1(:,:) = undef8
868 if(
io_l )
write(
io_fid_log,*)
'*** Not found namelist. Default used.' 869 elseif( ierr > 0 )
then 870 write(*,*)
'xxx Not appropriate names in namelist nm_mp_sn14_init. Check!' 879 a_area(i_mp_qc) = pi/4.0_rp
880 a_area(i_mp_qr) = pi/4.0_rp
881 a_area(i_mp_qi) = 0.65_rp*1.e-4_rp*100.0_rp**(2.00_rp)
882 a_area(i_mp_qs) = 0.2285_rp*1.e-4_rp*100.0_rp**(1.88_rp)
883 a_area(i_mp_qg) = 0.50_rp*1.e-4_rp*100.0_rp**(2.0_rp)
884 b_area(i_mp_qc) = 2.0_rp
885 b_area(i_mp_qr) = 2.0_rp
886 b_area(i_mp_qi) = 2.0_rp
887 b_area(i_mp_qs) = 1.88_rp
888 b_area(i_mp_qg) = 2.0_rp
894 a_m(i_mp_qc) = 0.124_rp
895 a_m(i_mp_qr) = 0.124_rp
896 a_m(i_mp_qi) = 0.217_rp
897 a_m(i_mp_qs) = 8.156_rp
898 a_m(i_mp_qg) = 0.190_rp
899 b_m(i_mp_qc) = 1.0_rp/3.0_rp
900 b_m(i_mp_qr) = 1.0_rp/3.0_rp
901 b_m(i_mp_qi) = 0.302_rp
902 b_m(i_mp_qs) = 0.526_rp
903 b_m(i_mp_qg) = 0.323_rp
907 alpha_v(i_mp_qc,:)= 3.75e+5_rp
908 alpha_v(i_mp_qr,:)= 159.0_rp
909 alpha_v(i_mp_qi,:)= 317.0_rp
910 alpha_v(i_mp_qs,:)= 27.70_rp
911 alpha_v(i_mp_qg,:)= 40.0_rp
912 beta_v(i_mp_qc,:) = 2.0_rp/3.0_rp
913 beta_v(i_mp_qr,:) = 0.266_rp
914 beta_v(i_mp_qi,:) = 0.363_rp
915 beta_v(i_mp_qs,:) = 0.216_rp
916 beta_v(i_mp_qg,:) = 0.230_rp
917 gamma_v(i_mp_qc) = 1.0_rp
919 gamma_v(i_mp_qr) = 1.0_rp/2.0_rp
920 gamma_v(i_mp_qi) = 1.0_rp/2.0_rp
921 gamma_v(i_mp_qs) = 1.0_rp/2.0_rp
922 gamma_v(i_mp_qg) = 1.0_rp/2.0_rp
931 nu(i_mp_qr) = -1.0_rp/3.0_rp
937 mu(i_mp_qr) = 1.0_rp/3.0_rp
938 mu(i_mp_qi) = 1.0_rp/3.0_rp
939 mu(i_mp_qs) = 1.0_rp/3.0_rp
940 mu(i_mp_qg) = 1.0_rp/3.0_rp
949 cap(i_mp_qc) = 2.0_rp
950 cap(i_mp_qr) = 2.0_rp
952 cap(i_mp_qs) = 2.0_rp
953 cap(i_mp_qg) = 2.0_rp
955 alpha_vn(:,:) = alpha_v(:,:)
956 beta_vn(:,:) = beta_v(:,:)
964 read(
io_fid_conf,nml=nm_mp_sn14_particles,iostat=ierr)
967 if(
io_l )
write(
io_fid_log,*)
'*** Not found namelist. Default used.' 968 elseif( ierr > 0 )
then 969 write(*,*)
'xxx Not appropriate names in namelist nm_mp_sn14_particles. Check!' 976 if( opt_m96_ice )
then 980 a_area(i_mp_qi) = 0.120284936_rp
981 a_area(i_mp_qs) = 0.131488_rp
982 a_area(i_mp_qg) = 0.5_rp
983 b_area(i_mp_qi) = 1.850000_rp
984 b_area(i_mp_qs) = 1.880000_rp
985 b_area(i_mp_qg) = 2.0_rp
986 a_m(i_mp_qi) = 1.23655360084766_rp
987 a_m(i_mp_qs) = a_m(i_mp_qi)
988 a_m(i_mp_qg) = 0.346111225718402_rp
989 b_m(i_mp_qi) = 0.408329930583912_rp
990 b_m(i_mp_qs) = b_m(i_mp_qi)
991 b_m(i_mp_qg) = 0.357142857142857_rp
993 if( opt_m96_column_ice )
then 996 a_area(i_mp_qi)= (0.684_rp*1.e-4_rp)*10.0_rp**(2.0_rp*2.00_rp)
997 b_area(i_mp_qi)= 2.0_rp
998 a_m(i_mp_qi) = 0.19834046116844_rp
999 b_m(i_mp_qi) = 0.343642611683849_rp
1002 wcap1 = sqrt(1.0_rp-ar_ice_fix**2)
1003 wcap2 = log( (1.0_rp+wcap1)/ar_ice_fix )
1004 cap(i_mp_qi) = 2.0_rp*wcap2/wcap1
1013 alpha_v(i_mp_qi,:) =(/ 5798.60107421875_rp, 167.347076416016_rp/)
1014 alpha_vn(i_mp_qi,:) =(/ 12408.177734375_rp, 421.799865722656_rp/)
1015 if( opt_m96_column_ice )
then 1016 alpha_v(i_mp_qi,:) = (/2901.0_rp, 32.20_rp/)
1017 alpha_vn(i_mp_qi,:) = (/9675.2_rp, 64.16_rp/)
1019 alpha_v(i_mp_qs,:) =(/ 15173.3916015625_rp, 305.678619384766_rp/)
1020 alpha_vn(i_mp_qs,:) =(/ 29257.1601562500_rp, 817.985717773438_rp/)
1021 alpha_v(i_mp_qg,:) =(/ 15481.6904296875_rp, 311.642242431641_rp/)
1022 alpha_vn(i_mp_qg,:) =(/ 27574.6562500000_rp, 697.536132812500_rp/)
1024 beta_v(i_mp_qi,:) =(/ 0.504873454570770_rp, 0.324817866086960_rp/)
1025 beta_vn(i_mp_qi,:) =(/ 0.548495233058929_rp, 0.385287821292877_rp/)
1026 if( opt_m96_column_ice )
then 1027 beta_v(i_mp_qi,:) =(/ 0.465552181005478_rp, 0.223826110363007_rp/)
1028 beta_vn(i_mp_qi,:) =(/ 0.530453503131866_rp, 0.273761242628098_rp/)
1030 beta_v(i_mp_qs,:) =(/ 0.528109610080719_rp, 0.329863965511322_rp/)
1031 beta_vn(i_mp_qs,:) =(/ 0.567154467105865_rp, 0.393876969814301_rp/)
1032 beta_v(i_mp_qg,:) =(/ 0.534656763076782_rp, 0.330253750085831_rp/)
1033 beta_vn(i_mp_qg,:) =(/ 0.570551633834839_rp, 0.387124240398407_rp/)
1037 ax_area(:) = a_area(:)*a_m(:)**b_area(:)
1038 bx_area(:) = b_area(:)*b_m(:)
1042 a_rea(:) = sqrt(ax_area(:)/pi)
1043 b_rea(:) = bx_area(:)/2.0_rp
1044 a_rea2(:) = a_rea(:)**2
1045 b_rea2(:) = b_rea(:)*2.0_rp
1046 a_rea3(:) = a_rea(:)**3
1047 b_rea3(:) = b_rea(:)*3.0_rp
1049 a_d2vt(:)=alpha_v(:,2)*(1.0_rp/alpha_v(:,2))**(beta_v(:,2)/b_m(:))
1050 b_d2vt(:)=(beta_v(:,2)/b_m(:))
1070 w1(iw) = gammafunc( (n+nu(iw)+1.0_rp)/mu(iw) )
1071 w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1072 w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1073 coef_m2(iw) = w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**n
1075 w4(iw) = gammafunc( (b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1076 coef_d(iw) = a_m(iw) * w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**b_m(iw)
1077 w5(iw) = gammafunc( (2.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
1078 w6(iw) = gammafunc( (3.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
1079 coef_d2v(iw) = a_m(iw) * w6(iw)/w5(iw)* ( w2(iw)/w3(iw) )**b_m(iw)
1080 coef_md2v(iw)= w5(iw)/w2(iw)* ( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw)+beta_v(iw,2))
1082 w7(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1083 coef_d3(iw) = a_m(iw)**3 * w7(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
1084 w8(iw) = gammafunc( (6.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1085 coef_d6(iw) = a_m(iw)**6 * w8(iw)/w2(iw)*( w2(iw)/w3(iw) )**(6.0_rp*b_m(iw))
1088 coef_deplc = coef_d(i_mp_qc)/a_m(i_mp_qc)
1094 w1(iw) = gammafunc( (2.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1095 w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1096 w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1098 w4(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
1100 coef_r2(iw)=w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw))
1101 coef_r3(iw)=w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
1102 coef_re(iw)=coef_r3(iw)/coef_r2(iw)
1109 w1(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
1110 w2(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
1111 w3(iw) = gammafunc( (b_rea2(iw)+nu(iw)+1.0_rp)/mu(iw) )
1112 w4(iw) = gammafunc( (b_rea3(iw)+nu(iw)+1.0_rp)/mu(iw) )
1114 coef_rea2(iw) = w3(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea2(iw)
1115 coef_rea3(iw) = w4(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea3(iw)
1121 w1(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1122 w2(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1123 coef_a(iw) = mu(iw)/w1(iw)
1125 coef_lambda(iw) = (w1(iw)/w2(iw))**(-mu(iw))
1133 w1(iw) = gammafunc( (beta_vn(iw,
ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
1134 w2(iw) = gammafunc( ( nu(iw) + 1.0_rp + n)/mu(iw) )
1135 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1136 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1138 coef_vt0(iw,
ia)=alpha_vn(iw,
ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_vn(iw,
ia)
1140 w1(iw) = gammafunc( (beta_v(iw,
ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
1141 w2(iw) = gammafunc( ( nu(iw) + 1.0_rp + n)/mu(iw) )
1143 coef_vt1(iw,
ia)=alpha_v(iw,
ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_v(iw,
ia)
1148 w1(iw) = gammafunc( ( b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
1149 w2(iw) = gammafunc( (1.0_rp + b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
1150 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1151 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1152 coef_dave_n(iw) = (w1(iw)/w3(iw))*(w3(iw)/w4(iw))** b_m(iw)
1153 coef_dave_l(iw) = (w2(iw)/w3(iw))*(w3(iw)/w4(iw))**(1.0_rp+b_m(iw))
1157 ah_vent(i_mp_qc,1:2) = (/1.0000_rp,1.0000_rp/)
1158 ah_vent(i_mp_qr,1:2) = (/1.0000_rp,0.780_rp/)
1159 ah_vent(i_mp_qi,1:2) = (/1.0000_rp,0.860_rp/)
1160 ah_vent(i_mp_qs,1:2) = (/1.0000_rp,0.780_rp/)
1161 ah_vent(i_mp_qg,1:2) = (/1.0000_rp,0.780_rp/)
1162 bh_vent(i_mp_qc,1:2) = (/0.0000_rp,0.0000_rp/)
1163 bh_vent(i_mp_qr,1:2) = (/0.108_rp,0.308_rp/)
1164 bh_vent(i_mp_qi,1:2) = (/0.140_rp,0.280_rp/)
1165 bh_vent(i_mp_qs,1:2) = (/0.108_rp,0.308_rp/)
1166 bh_vent(i_mp_qg,1:2) = (/0.108_rp,0.308_rp/)
1170 if( (nu(iw) + b_m(iw) + n) > eps_gamma )
then 1171 w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
1172 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1173 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1174 ah_vent0(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1175 ah_vent0(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1176 flag_vent0(iw)=.true.
1178 ah_vent0(iw,1)= 1.0_rp
1179 ah_vent0(iw,2)= 1.0_rp
1180 flag_vent0(iw)=.false.
1183 if( (nu(iw) + b_m(iw) + n) > eps_gamma )
then 1184 w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
1185 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1186 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1187 ah_vent1(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1188 ah_vent1(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
1189 flag_vent1(iw)=.true.
1191 ah_vent1(iw,1)= 1.0_rp
1192 ah_vent1(iw,2)= 1.0_rp
1193 flag_vent1(iw)=.true.
1198 if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )
then 1199 flag_vent0(iw)=.false.
1201 if(flag_vent0(iw))
then 1202 w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
1203 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1204 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1206 w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
1207 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)
1208 w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
1209 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)
1211 bh_vent0(iw,1) = 0.0_rp
1212 bh_vent0(iw,2) = 0.0_rp
1216 if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )
then 1217 flag_vent1(iw)=.false.
1219 if(flag_vent1(iw))
then 1220 w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
1221 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1222 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1224 w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
1225 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)
1227 w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
1228 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)
1230 bh_vent1(iw,1) = 0.0_rp
1231 bh_vent1(iw,2) = 0.0_rp
1240 w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1241 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1242 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1243 delta_b0(iw) = w1(iw)/w2(iw) &
1244 *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
1246 w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1247 delta_b1(iw) = w1(iw)/w2(iw) &
1248 *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
1254 w1(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1255 w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1256 w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1257 w4(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1259 w5(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1267 delta_ab0(
ia,ib) = 2.0_rp*(w1(ib)/w2(ib))*(w4(
ia)/w2(
ia)) &
1268 * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
1269 * ( w2(
ia)/w3(
ia) )**(b_rea(
ia) )
1271 delta_ab1(
ia,ib) = 2.0_rp*(w5(ib)/w2(ib))*(w4(
ia)/w2(
ia)) &
1272 * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
1273 * ( w2(
ia)/w3(
ia) )**(b_rea(
ia) )
1281 w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1282 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1283 w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1284 w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1285 theta_b0(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
1287 w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1288 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1289 theta_b1(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
1296 w1(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1297 w2(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1298 w3(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1299 w4(iw) = gammafunc( ( 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp )/mu(iw) )
1301 w5(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
1302 w6(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
1304 w7(iw) = gammafunc( (beta_v(iw,2) + b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1305 w8(iw) = gammafunc( ( b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
1310 theta_ab0(
ia,ib) = 2.0_rp * (w1(ib)/w2(ib))*(w3(
ia)/w4(
ia)) &
1311 * (w5(
ia)/w6(
ia))**beta_v(
ia,2) &
1312 * (w5(ib)/w6(ib))**beta_v(ib,2)
1313 theta_ab1(
ia,ib) = 2.0_rp * (w7(ib)/w8(ib))*(w3(
ia)/w4(
ia)) &
1314 * (w5(
ia)/w6(
ia))**beta_v(
ia,2) &
1315 * (w5(ib)/w6(ib))**beta_v(ib,2)
1321 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_m2 ",coef_m2(:)
1324 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_d3 ",coef_d3(:)
1325 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_d6 ",coef_d6(:)
1326 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_d2v ",coef_d2v(:)
1327 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_md2v ",coef_md2v(:)
1331 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_r2 ",coef_r2(:)
1332 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_r3 ",coef_r3(:)
1333 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_re ",coef_re(:)
1337 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"ax_area ",ax_area(:)
1338 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"bx_area ",bx_area(:)
1344 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_rea2 ",coef_rea2(:)
1345 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_rea3 ",coef_rea3(:)
1346 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_vt0 ",coef_vt0(:,1)
1347 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_vt1 ",coef_vt1(:,1)
1349 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"coef_lambda ",coef_lambda(:)
1351 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"ah_vent0 sml",ah_vent0(:,1)
1352 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"ah_vent0 lrg",ah_vent0(:,2)
1353 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"ah_vent1 sml",ah_vent1(:,1)
1354 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"ah_vent1 lrg",ah_vent1(:,2)
1355 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"bh_vent0 sml",bh_vent0(:,1)
1356 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"bh_vent0 lrg",bh_vent0(:,2)
1357 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"bh_vent1 sml",bh_vent1(:,1)
1358 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"bh_vent1 lrg",bh_vent1(:,2)
1360 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"delta_b0 ",delta_b0(:)
1361 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"delta_b1 ",delta_b1(:)
1362 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"theta_b0 ",theta_b0(:)
1363 if(
io_l )
write(
io_fid_log,
'(a,100ES16.6)')
"theta_b1 ",theta_b1(:)
1366 if(
io_l )
write(
io_fid_log,
'(a,a10,a,100ES16.6)')
"delta0(a,b)=(",trim(wlabel(
ia)),
",b)=",(delta_ab0(
ia,ib),ib=1,hydro_max)
1369 if(
io_l )
write(
io_fid_log,
'(a,a10,a,100ES16.6)')
"delta1(a,b)=(",trim(wlabel(
ia)),
",b)=",(delta_ab1(
ia,ib),ib=1,hydro_max)
1372 if(
io_l )
write(
io_fid_log,
'(a,a10,a,100ES16.6)')
"theta0(a,b)=(",trim(wlabel(
ia)),
",b)=",(theta_ab0(
ia,ib),ib=1,hydro_max)
1375 if(
io_l )
write(
io_fid_log,
'(a,a10,a,100ES16.6)')
"theta1(a,b)=(",trim(wlabel(
ia)),
",b)=",(theta_ab1(
ia,ib),ib=1,hydro_max)
1379 end subroutine mp_sn14_init
1381 subroutine mp_sn14 ( &
1403 moist_psat_liq => atmos_saturation_psat_liq, &
1404 moist_psat_ice => atmos_saturation_psat_ice
1411 real(RP),
intent(inout) :: dens (
ka,
ia,
ja)
1412 real(RP),
intent(inout) :: momz (
ka,
ia,
ja)
1413 real(RP),
intent(inout) :: momx (
ka,
ia,
ja)
1414 real(RP),
intent(inout) :: momy (
ka,
ia,
ja)
1415 real(RP),
intent(inout) :: rhot (
ka,
ia,
ja)
1416 real(RP),
intent(inout) :: qtrc (
ka,
ia,
ja,
qa)
1417 real(RP),
intent(in) :: ccn (
ka,
ia,
ja)
1418 real(RP),
intent(out) :: evaporate(
ka,
ia,
ja)
1419 real(RP),
intent(out) :: sflx_rain(
ia,
ja)
1420 real(RP),
intent(out) :: sflx_snow(
ia,
ja)
1426 real(RP) :: pott (
ka,
ia,
ja)
1430 real(RP) :: qdry(
ka,
ia,
ja)
1431 real(RP) :: temp(
ka,
ia,
ja)
1432 real(RP) :: pres(
ka,
ia,
ja)
1433 real(RP) :: rhoe(
ka,
ia,
ja)
1434 real(RP) :: velz(
ka,
ia,
ja)
1438 real(RP) :: xq(hydro_max,
ka,
ia,
ja)
1440 real(RP) :: dq_xa(hydro_max,
ka,
ia,
ja)
1441 real(RP) :: vt_xa(hydro_max,2,
ka,
ia,
ja)
1443 real(RP) :: wtemp(
ka,
ia,
ja)
1444 real(RP) :: esw(
ka,
ia,
ja)
1445 real(RP) :: esi(
ka,
ia,
ja)
1448 real(RP) :: rho_fac_q(hydro_max,
ka,
ia,
ja)
1449 real(RP) :: cva(
ka,
ia,
ja)
1450 real(RP) :: cpa(
ka,
ia,
ja)
1453 real(RP) :: drhogqc, drhognc
1454 real(RP) :: drhogqr, drhognr
1455 real(RP) :: drhogqi, drhogni
1456 real(RP) :: drhogqs, drhogns
1457 real(RP) :: drhogqg, drhogng
1460 real(RP) :: pq(pq_max,
ka,
ia,
ja)
1462 real(RP) :: wrm_dqc, wrm_dnc
1463 real(RP) :: wrm_dqr, wrm_dnr
1466 real(RP) :: pac(pac_max,
ka,
ia,
ja)
1468 real(RP) :: gc_dqc, gc_dnc
1469 real(RP) :: sc_dqc, sc_dnc
1470 real(RP) :: ic_dqc, ic_dnc
1471 real(RP) :: rg_dqg, rg_dng
1472 real(RP) :: rg_dqr, rg_dnr
1473 real(RP) :: rs_dqr, rs_dnr, rs_dqs, rs_dns
1474 real(RP) :: ri_dqr, ri_dnr
1475 real(RP) :: ri_dqi, ri_dni
1476 real(RP) :: ii_dqi, ii_dni
1477 real(RP) :: is_dqi, is_dni, ss_dns
1478 real(RP) :: gs_dqs, gs_dns, gg_dng
1480 real(RP) :: clp_dqc, clp_dnc, clm_dqc, clm_dnc
1481 real(RP) :: clp_dqr, clp_dnr, clm_dqr, clm_dnr
1482 real(RP) :: clp_dqi, clp_dni, clm_dqi, clm_dni
1483 real(RP) :: clp_dqs, clp_dns, clm_dqs, clm_dns
1484 real(RP) :: clp_dqg, clp_dng, clm_dqg, clm_dng
1485 real(RP) :: fac1, fac3, fac4, fac6, fac7, fac9, fac10
1487 real(RP) :: pco_dqi, pco_dni
1488 real(RP) :: pco_dqs, pco_dns
1489 real(RP) :: pco_dqg, pco_dng
1491 real(RP) :: eml_dqc, eml_dnc
1492 real(RP) :: eml_dqr, eml_dnr
1493 real(RP) :: eml_dqi, eml_dni
1494 real(RP) :: eml_dqs, eml_dns
1495 real(RP) :: eml_dqg, eml_dng
1497 real(RP) :: spl_dqi, spl_dni
1498 real(RP) :: spl_dqg, spl_dqs
1500 real(RP) :: rrho(
ka,
ia,
ja)
1505 real(RP) :: dtdt_equiv_d(
ka,
ia,
ja)
1512 real(RP) :: sl_plcdep(
ia,
ja)
1513 real(RP) :: sl_plrdep(
ia,
ja), sl_pnrdep(
ia,
ja)
1515 real(RP) :: qke_d(
ka,
ia,
ja)
1517 real(RP),
parameter :: eps = 1.e-19_rp
1518 real(RP),
parameter :: eps_qv = 1.e-19_rp
1519 real(RP),
parameter :: eps_rhoge = 1.e-19_rp
1520 real(RP),
parameter :: eps_rhog = 1.e-19_rp
1533 integer :: k, i, j, iq
1549 dt =
real(dt_dp,kind=
rp)
1562 rhoq(iq,k,i,j) = dens(k,i,j) * qtrc(k,i,j,iq)
1564 rhoq2(
i_qv,k,i,j) = dens(k,i,j)*qtrc(k,i,j,
i_qv)
1565 rhoq2(
i_ni,k,i,j) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,
i_ni) )
1566 rhoq2(
i_nc,k,i,j) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,
i_nc) )
1573 velz(
ks-1,i,j) = 0.0_rp
1575 velz(k,i,j) = momz(k,i,j) / ( dens(k,i,j) + dens(k+1,i,j) ) * 2.0_rp
1577 velz(
ke,i,j) = 0.0_rp
1584 rrho(k,i,j) = 1.0_rp / dens(k,i,j)
1585 pott(k,i,j) = rhot(k,i,j) * rrho(k,i,j)
1586 calc_qdry( qdry(k,i,j), qtrc,
tracer_mass, k, i, j, iq )
1587 calc_cv( cva(k,i,j), qdry(k,i,j), qtrc, k, i, j, iq, cvdry,
tracer_cv )
1588 calc_r( rmoist, qdry(k,i,j), qtrc, k, i, j, iq, rdry,
tracer_r )
1589 cpa(k,i,j) = cva(k,i,j) + rmoist
1590 calc_pre( pres(k,i,j), dens(k,i,j), pott(k,i,j), rmoist, cpa(k,i,j), p00 )
1591 temp(k,i,j) = pres(k,i,j) / ( dens(k,i,j) * rmoist )
1592 rhoe(k,i,j) = dens(k,i,j) * temp(k,i,j) * cva(k,i,j)
1593 wtemp(k,i,j) = max(temp(k,i,j), tem_min)
1598 if( opt_debug_tem )
call debug_tem_kij( 1, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,
i_qv) )
1603 rho_fac = rho_0 / max(dens(k,i,j),rho_min)
1604 rho_fac_q(i_mp_qc,k,i,j) = rho_fac**gamma_v(i_mp_qc)
1605 rho_fac_q(i_mp_qr,k,i,j) = rho_fac**gamma_v(i_mp_qr)
1606 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
1607 rho_fac_q(i_mp_qs,k,i,j) = rho_fac_q(i_mp_qi,k,i,j)
1608 rho_fac_q(i_mp_qg,k,i,j) = rho_fac_q(i_mp_qi,k,i,j)
1616 sl_plcdep(i,j) = 0.0_rp
1617 sl_plrdep(i,j) = 0.0_rp
1618 sl_pnrdep(i,j) = 0.0_rp
1626 qke_d(k,i,j) = 0.0_rp
1635 dtdt_equiv_d(k,i,j) = 0.0_rp
1646 dens, wtemp, pres, &
1659 drhogqc = dt * pq(i_lcccn,k,i,j)
1660 drhognc = dt * pq(i_ncccn,k,i,j)
1661 drhogqi = dt * pq(i_liccn,k,i,j)
1662 drhogni = dt * pq(i_niccn,k,i,j)
1663 drhogqv = max( -rhoq(
i_qv,k,i,j), -drhogqc-drhogqi )
1664 fac1 = drhogqv / min( -drhogqc-drhogqi, -eps )
1666 rhoq(
i_qv,k,i,j) = rhoq(
i_qv,k,i,j) + drhogqv
1667 rhoq(
i_qc,k,i,j) = max(0.0_rp, rhoq(
i_qc,k,i,j) + drhogqc*fac1)
1668 rhoq(
i_qi,k,i,j) = max(0.0_rp, rhoq(
i_qi,k,i,j) + drhogqi*fac1)
1669 rhoq(
i_nc,k,i,j) = max(0.0_rp, rhoq(
i_nc,k,i,j) + drhognc)
1670 rhoq(
i_ni,k,i,j) = max(0.0_rp, rhoq(
i_ni,k,i,j) + drhogni)
1673 rhoq(
i_nc,k,i,j) = min( rhoq(
i_nc,k,i,j), nc_uplim_d(1,i,j) )
1675 rhoe(k,i,j) = rhoe(k,i,j) - lhv * drhogqv + lhf * drhogqi*fac1
1677 qtrc(k,i,j,
i_qv) = rhoq(
i_qv,k,i,j) * rrho(k,i,j)
1678 qtrc(k,i,j,
i_qc) = rhoq(
i_qc,k,i,j) * rrho(k,i,j)
1679 qtrc(k,i,j,
i_qi) = rhoq(
i_qi,k,i,j) * rrho(k,i,j)
1680 qtrc(k,i,j,
i_nc) = rhoq(
i_nc,k,i,j) * rrho(k,i,j)
1681 qtrc(k,i,j,
i_ni) = rhoq(
i_ni,k,i,j) * rrho(k,i,j)
1683 calc_qdry( qdry(k,i,j), qtrc,
tracer_mass, k, i, j, iq )
1684 calc_cv( cva(k,i,j), qdry(k,i,j), qtrc, k, i, j, iq, cvdry,
tracer_cv )
1685 calc_r( rmoist, qdry(k,i,j), qtrc, k, i, j, iq, rdry,
tracer_r )
1686 temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cva(k,i,j) )
1687 pres(k,i,j) = dens(k,i,j) * rmoist * temp(k,i,j)
1688 wtemp(k,i,j) = max( temp(k,i,j), tem_min )
1694 if( opt_debug_tem )
call debug_tem_kij( 2, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,
i_qv) )
1708 rhoq2(
i_qr,k,i,j) = rhoq(
i_qr,k,i,j)
1709 rhoq2(
i_nr,k,i,j) = rhoq(
i_nr,k,i,j)
1710 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) ))
1712 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)
1713 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)
1714 vt_xa(i_mp_qr,2,k,i,j) = vt_xa(i_mp_qr,1,k,i,j)
1718 rhoq2(
i_qv,k,i,j) = rhoq(
i_qv,k,i,j)
1719 rhoq2(
i_qc,k,i,j) = rhoq(
i_qc,k,i,j)
1720 rhoq2(
i_qi,k,i,j) = rhoq(
i_qi,k,i,j)
1721 rhoq2(
i_qs,k,i,j) = rhoq(
i_qs,k,i,j)
1722 rhoq2(
i_qg,k,i,j) = rhoq(
i_qg,k,i,j)
1724 rhoq2(
i_nc,k,i,j) = rhoq(
i_nc,k,i,j)
1725 rhoq2(
i_ni,k,i,j) = rhoq(
i_ni,k,i,j)
1726 rhoq2(
i_ns,k,i,j) = rhoq(
i_ns,k,i,j)
1727 rhoq2(
i_ng,k,i,j) = rhoq(
i_ng,k,i,j)
1732 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) ))
1733 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) ))
1734 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) ))
1735 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) ))
1738 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)
1739 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)
1740 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)
1741 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)
1744 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)
1745 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)
1746 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)
1747 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)
1748 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)
1749 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)
1750 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)
1751 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)
1757 call moist_psat_liq( esw, wtemp )
1758 call moist_psat_ice( esi, wtemp )
1767 dens, wtemp, pres, qdry, &
1778 ntdiv, ntmax_phase_change, &
1794 sl_plrdep, sl_pnrdep )
1797 if( opt_debug_tem )
call debug_tem_kij( 3, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,
i_qv) )
1813 rhoq2(
i_qv,k,i,j) = rhoq(
i_qv,k,i,j)
1814 rhoq2(
i_qc,k,i,j) = rhoq(
i_qc,k,i,j)
1815 rhoq2(
i_qr,k,i,j) = rhoq(
i_qr,k,i,j)
1816 rhoq2(
i_qi,k,i,j) = rhoq(
i_qi,k,i,j)
1817 rhoq2(
i_qs,k,i,j) = rhoq(
i_qs,k,i,j)
1818 rhoq2(
i_qg,k,i,j) = rhoq(
i_qg,k,i,j)
1820 rhoq2(
i_nc,k,i,j) = rhoq(
i_nc,k,i,j)
1821 rhoq2(
i_nr,k,i,j) = rhoq(
i_nr,k,i,j)
1822 rhoq2(
i_ni,k,i,j) = rhoq(
i_ni,k,i,j)
1823 rhoq2(
i_ns,k,i,j) = rhoq(
i_ns,k,i,j)
1824 rhoq2(
i_ng,k,i,j) = rhoq(
i_ng,k,i,j)
1827 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) ) )
1828 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) ) )
1829 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) ) )
1830 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) ) )
1831 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) ) )
1834 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)
1835 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)
1836 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)
1837 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)
1838 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)
1842 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)
1843 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)
1844 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)
1845 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)
1846 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)
1853 if ( mp_doautoconversion )
then 1863 pq(i_lcaut,k,i,j) = 0.0_rp
1864 pq(i_ncaut,k,i,j) = 0.0_rp
1865 pq(i_nraut,k,i,j) = 0.0_rp
1866 pq(i_lcacc,k,i,j) = 0.0_rp
1867 pq(i_ncacc,k,i,j) = 0.0_rp
1868 pq(i_nrslc,k,i,j) = 0.0_rp
1869 pq(i_nrbrk,k,i,j) = 0.0_rp
1891 profile_start(
"sn14_update_rhoq")
1897 wrm_dqc = max( dt*( pq(i_lcaut,k,i,j)+pq(i_lcacc,k,i,j) ), -rhoq2(
i_qc,k,i,j) )
1898 wrm_dnc = max( dt*( pq(i_ncaut,k,i,j)+pq(i_ncacc,k,i,j) ), -rhoq2(
i_nc,k,i,j) )
1899 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) )
1908 gc_dqc = max( dt*pac(i_lgaclc2lg,k,i,j) , min(0.0_rp, -rhoq2(
i_qc,k,i,j)-wrm_dqc ))
1909 sc_dqc = max( dt*pac(i_lsaclc2ls,k,i,j) , min(0.0_rp, -rhoq2(
i_qc,k,i,j)-wrm_dqc-gc_dqc ))
1910 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 ))
1912 gc_dnc = max( dt*pac(i_ngacnc2ng,k,i,j) , min(0.0_rp, -rhoq2(
i_nc,k,i,j)-wrm_dnc ))
1913 sc_dnc = max( dt*pac(i_nsacnc2ns,k,i,j) , min(0.0_rp, -rhoq2(
i_nc,k,i,j)-wrm_dnc-gc_dnc ))
1914 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 ))
1917 sw = sign(0.5_rp, t00-temp(k,i,j)) + 0.5_rp
1918 rg_dqr = max( dt*pac(i_lraclg2lg, k,i,j), min(0.0_rp, -rhoq2(
i_qr,k,i,j)-wrm_dqr )) * sw
1919 rg_dqg = max( dt*pac(i_lraclg2lg, k,i,j), min(0.0_rp, -rhoq2(
i_qg,k,i,j) )) * ( 1.0_rp - sw )
1920 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
1921 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
1923 rg_dnr = max( dt*pac(i_nracng2ng, k,i,j), min(0.0_rp, -rhoq2(
i_nr,k,i,j)-wrm_dnr )) * sw
1924 rg_dng = max( dt*pac(i_nracng2ng, k,i,j), min(0.0_rp, -rhoq2(
i_ng,k,i,j) )) * ( 1.0_rp - sw )
1925 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
1926 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
1929 fac1 = (ri_dqr-eps)/ (dt*pac(i_lracli2lg_r,k,i,j)-eps)
1930 ri_dqi = max( dt*pac(i_lracli2lg_i,k,i,j)*fac1, min(0.0_rp, -rhoq2(
i_qi,k,i,j)+ic_dqc ))
1931 ii_dqi = max( dt*pac(i_liacli2ls,k,i,j) , min(0.0_rp, -rhoq2(
i_qi,k,i,j)+ic_dqc-ri_dqi ))
1932 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 ))
1934 fac4 = (ri_dnr-eps)/ (dt*pac(i_nracni2ng_r,k,i,j)-eps)
1935 ri_dni = max( dt*pac(i_nracni2ng_i,k,i,j)*fac4, min(0.0_rp, -rhoq2(
i_ni,k,i,j) ))
1936 ii_dni = max( dt*pac(i_niacni2ns,k,i,j) , min(0.0_rp, -rhoq2(
i_ni,k,i,j)-ri_dni ))
1937 is_dni = max( dt*pac(i_niacns2ns,k,i,j) , min(0.0_rp, -rhoq2(
i_ni,k,i,j)-ri_dni-ii_dni ))
1939 fac3 = (rs_dqr-eps)/(dt*pac(i_lracls2lg_r,k,i,j)-eps)
1940 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 ))
1941 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 ))
1943 fac6 = (rs_dnr-eps)/(dt*pac(i_nracns2ng_r,k,i,j)-eps)
1945 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 ))
1946 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 ))
1947 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 ))
1949 gg_dng = max( dt*pac(i_ngacng2ng,k,i,j) , min(0.0_rp, -rhoq2(
i_ng,k,i,j) ))
1955 clp_dqr = (-rg_dqg-rs_dqs-ri_dqi) * (1.0_rp-sw)
1957 clp_dqs = -sc_dqc-ii_dqi-is_dqi
1958 clp_dqg = -gc_dqc -gs_dqs + (-rg_dqr-rs_dqr-rs_dqs-ri_dqr-ri_dqi) * sw
1963 clp_dns = -ii_dni*0.5_rp
1964 clp_dng = (-rs_dnr-ri_dnr) * sw
1967 clm_dqc = gc_dqc+sc_dqc+ic_dqc
1968 clm_dqr = (rg_dqr+rs_dqr+ri_dqr) * sw
1969 clm_dqi = ri_dqi+ii_dqi+is_dqi
1970 clm_dqs = rs_dqs+gs_dqs
1971 clm_dqg = rg_dqg * (1.0_rp-sw)
1973 clm_dnc = gc_dnc+sc_dnc+ic_dnc
1974 clm_dnr = (rg_dnr+rs_dnr+ri_dnr) * sw
1975 clm_dni = ri_dni+ii_dni+is_dni
1976 clm_dns = rs_dns+ss_dns+gs_dns
1977 clm_dng = gg_dng + rg_dng * (1.0_rp-sw)
1981 pco_dqi = max( dt*pq(i_licon,k,i,j), -clp_dqi )
1982 pco_dqs = max( dt*pq(i_lscon,k,i,j), -clp_dqs )
1983 pco_dqg = -pco_dqi-pco_dqs
1985 pco_dni = max( dt*pq(i_nicon,k,i,j), -clp_dni )
1986 pco_dns = max( dt*pq(i_nscon,k,i,j), -clp_dns )
1987 pco_dng = -pco_dni-pco_dns
1990 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 ))
1991 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 ))
1992 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)), &
1993 min(0.0_rp, -rhoq2(
i_qg,k,i,j)-(clp_dqg+clm_dqg)-pco_dqg ))
1995 eml_dqr = -eml_dqs-eml_dqg
1997 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 ))
1998 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 ))
1999 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)), &
2000 min(0.0_rp, -rhoq2(
i_ng,k,i,j)-(clp_dng+clm_dng)-pco_dng ))
2002 eml_dnr = -eml_dns-eml_dng
2005 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 ))
2006 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 ))
2007 spl_dqi = -spl_dqg-spl_dqs
2008 fac9 = (spl_dqg-eps)/(dt*pq(i_lgspl,k,i,j)-eps)
2009 fac10 = (spl_dqs-eps)/(dt*pq(i_lsspl,k,i,j)-eps)
2010 spl_dni = dt*pq(i_nispl,k,i,j)*fac9*fac10
2013 drhogqc = (wrm_dqc + clp_dqc + clm_dqc + eml_dqc )
2014 drhognc = (wrm_dnc + clp_dnc + clm_dnc + eml_dnc )
2016 drhogqr = (wrm_dqr + clp_dqr + clm_dqr + eml_dqr )
2017 drhognr = (wrm_dnr + clp_dnr + clm_dnr + eml_dnr )
2019 drhogqi = ( clp_dqi + clm_dqi + pco_dqi + eml_dqi + spl_dqi)
2020 drhogni = ( clp_dni + clm_dni + pco_dni + eml_dni + spl_dni)
2022 drhogqs = ( clp_dqs + clm_dqs + pco_dqs + eml_dqs + spl_dqs)
2023 drhogns = ( clp_dns + clm_dns + pco_dns + eml_dns )
2025 drhogqg = ( clp_dqg + clm_dqg + pco_dqg + eml_dqg + spl_dqg)
2026 drhogng = ( clp_dng + clm_dng + pco_dng + eml_dng )
2030 rhoq(
i_qc,k,i,j) = max(0.0_rp, rhoq(
i_qc,k,i,j) + drhogqc )
2031 rhoq(
i_nc,k,i,j) = max(0.0_rp, rhoq(
i_nc,k,i,j) + drhognc )
2032 rhoq(
i_qr,k,i,j) = max(0.0_rp, rhoq(
i_qr,k,i,j) + drhogqr )
2033 rhoq(
i_nr,k,i,j) = max(0.0_rp, rhoq(
i_nr,k,i,j) + drhognr )
2034 rhoq(
i_qi,k,i,j) = max(0.0_rp, rhoq(
i_qi,k,i,j) + drhogqi )
2035 rhoq(
i_ni,k,i,j) = max(0.0_rp, rhoq(
i_ni,k,i,j) + drhogni )
2036 rhoq(
i_qs,k,i,j) = max(0.0_rp, rhoq(
i_qs,k,i,j) + drhogqs )
2037 rhoq(
i_ns,k,i,j) = max(0.0_rp, rhoq(
i_ns,k,i,j) + drhogns )
2038 rhoq(
i_qg,k,i,j) = max(0.0_rp, rhoq(
i_qg,k,i,j) + drhogqg )
2039 rhoq(
i_ng,k,i,j) = max(0.0_rp, rhoq(
i_ng,k,i,j) + drhogng )
2043 rhoe(k,i,j) = rhoe(k,i,j) + lhf * ( drhogqi + drhogqs + drhogqg )
2047 profile_stop(
"sn14_update_rhoq")
2058 qtrc(k,i,j,iq) = rhoq(iq,k,i,j) * rrho(k,i,j)
2061 calc_qdry( qdry(k,i,j), qtrc,
tracer_mass, k, i, j, iq )
2062 calc_cv( cva(k,i,j), qdry(k,i,j), qtrc, k, i, j, iq, cvdry,
tracer_cv )
2063 calc_r( rmoist, qdry(k,i,j), qtrc, k, i, j, iq, rdry,
tracer_r )
2064 cpa(k,i,j) = cva(k,i,j) + rmoist
2065 temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cva(k,i,j) )
2066 pres(k,i,j) = dens(k,i,j) * rmoist * temp(k,i,j)
2067 rhot(k,i,j) = temp(k,i,j) * ( p00 / pres(k,i,j) )**(rmoist/cpa(k,i,j)) &
2074 if( opt_debug_tem )
call debug_tem_kij( 4, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,
i_qv) )
2093 vterm(:,:,:,:) = 0.0_rp
2094 flx_hydro(:,:,:,:) = 0.0_rp
2096 if ( mp_doprecipitation )
then 2098 do step = 1, mp_nstep_sedimentation
2100 call mp_terminal_velocity( vterm(:,:,:,:), &
2106 call mp_precipitation(
qa_mp, &
2118 mp_dtsec_sedimentation )
2123 flx_hydro(k,i,j,i_mp_qc) = flx_hydro(k,i,j,i_mp_qc) + pflux(k,i,j,i_mp_qc) * mp_rnstep_sedimentation
2124 flx_hydro(k,i,j,i_mp_qr) = flx_hydro(k,i,j,i_mp_qr) + pflux(k,i,j,i_mp_qr) * mp_rnstep_sedimentation
2125 flx_hydro(k,i,j,i_mp_qi) = flx_hydro(k,i,j,i_mp_qi) + pflux(k,i,j,i_mp_qi) * mp_rnstep_sedimentation
2126 flx_hydro(k,i,j,i_mp_qs) = flx_hydro(k,i,j,i_mp_qs) + pflux(k,i,j,i_mp_qs) * mp_rnstep_sedimentation
2127 flx_hydro(k,i,j,i_mp_qg) = flx_hydro(k,i,j,i_mp_qg) + pflux(k,i,j,i_mp_qg) * mp_rnstep_sedimentation
2128 flx_hydro(k,i,j,i_mp_nc) = flx_hydro(k,i,j,i_mp_nc) + pflux(k,i,j,i_mp_nc) * mp_rnstep_sedimentation
2129 flx_hydro(k,i,j,i_mp_nr) = flx_hydro(k,i,j,i_mp_nr) + pflux(k,i,j,i_mp_nr) * mp_rnstep_sedimentation
2130 flx_hydro(k,i,j,i_mp_ni) = flx_hydro(k,i,j,i_mp_ni) + pflux(k,i,j,i_mp_ni) * mp_rnstep_sedimentation
2131 flx_hydro(k,i,j,i_mp_ns) = flx_hydro(k,i,j,i_mp_ns) + pflux(k,i,j,i_mp_ns) * mp_rnstep_sedimentation
2132 flx_hydro(k,i,j,i_mp_ng) = flx_hydro(k,i,j,i_mp_ng) + pflux(k,i,j,i_mp_ng) * mp_rnstep_sedimentation
2141 call hist_in( flx_hydro(:,:,:,i_mp_qc),
'pflux_QC',
'precipitation flux of QC',
'kg/m2/s', nohalo=.true. )
2142 call hist_in( flx_hydro(:,:,:,i_mp_qr),
'pflux_QR',
'precipitation flux of QR',
'kg/m2/s', nohalo=.true. )
2143 call hist_in( flx_hydro(:,:,:,i_mp_qi),
'pflux_QI',
'precipitation flux of QI',
'kg/m2/s', nohalo=.true. )
2144 call hist_in( flx_hydro(:,:,:,i_mp_qs),
'pflux_QS',
'precipitation flux of QS',
'kg/m2/s', nohalo=.true. )
2145 call hist_in( flx_hydro(:,:,:,i_mp_qg),
'pflux_QG',
'precipitation flux of QG',
'kg/m2/s', nohalo=.true. )
2146 call hist_in( flx_hydro(:,:,:,i_mp_nc),
'pflux_QC',
'precipitation flux of NC',
'kg/m2/s', nohalo=.true. )
2147 call hist_in( flx_hydro(:,:,:,i_mp_nr),
'pflux_QR',
'precipitation flux of NR',
'kg/m2/s', nohalo=.true. )
2148 call hist_in( flx_hydro(:,:,:,i_mp_ni),
'pflux_QI',
'precipitation flux of NI',
'kg/m2/s', nohalo=.true. )
2149 call hist_in( flx_hydro(:,:,:,i_mp_ns),
'pflux_QS',
'precipitation flux of NS',
'kg/m2/s', nohalo=.true. )
2150 call hist_in( flx_hydro(:,:,:,i_mp_ng),
'pflux_QG',
'precipitation flux of NG',
'kg/m2/s', nohalo=.true. )
2154 sflx_rain(i,j) = - ( flx_hydro(
ks-1,i,j,i_mp_qc) &
2155 + flx_hydro(
ks-1,i,j,i_mp_qr) )
2156 sflx_snow(i,j) = - ( flx_hydro(
ks-1,i,j,i_mp_qi) &
2157 + flx_hydro(
ks-1,i,j,i_mp_qs) &
2158 + flx_hydro(
ks-1,i,j,i_mp_qg) )
2165 end subroutine mp_sn14
2178 integer,
intent(in) :: point
2179 real(RP),
intent(in) :: tem(KA,IA,JA)
2180 real(RP),
intent(in) :: rho(KA,IA,JA)
2181 real(RP),
intent(in) :: pre(KA,IA,JA)
2182 real(RP),
intent(in) :: qv (KA,IA,JA)
2190 if ( tem(k,i,j) < tem_min &
2191 .OR. rho(k,i,j) < rho_min &
2192 .OR. pre(k,i,j) < 1.0_rp )
then 2195 "*** 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 2219 moist_psat_liq => atmos_saturation_psat_liq, &
2220 moist_psat_ice => atmos_saturation_psat_ice, &
2221 moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
2222 moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
2226 real(RP),
intent(in) :: z(KA)
2227 real(RP),
intent(in) :: velz(KA,IA,JA)
2228 real(RP),
intent(in) :: rho(KA,IA,JA)
2229 real(RP),
intent(in) :: tem(KA,IA,JA)
2230 real(RP),
intent(in) :: pre(KA,IA,JA)
2232 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
2233 real(RP),
intent(out) :: PQ(PQ_MAX,KA,IA,JA)
2235 real(RP),
intent(in) :: cpa(KA,IA,JA)
2236 real(RP),
intent(in) :: dTdt_rad(KA,IA,JA)
2237 real(RP),
intent(in) :: qke(KA,IA,JA)
2238 real(RP),
intent(in) :: dt
2239 real(RP),
intent(in) :: CCN(KA,IA,JA)
2244 real(RP),
parameter :: c_ccn_ocean= 1.00e+8_rp
2245 real(RP),
parameter :: c_ccn_land = 1.26e+9_rp
2246 real(RP),
save :: c_ccn = 1.00e+8_rp
2248 real(RP),
parameter :: kappa_ocean= 0.462_rp
2249 real(RP),
parameter :: kappa_land = 0.308_rp
2250 real(RP),
save :: kappa = 0.462_rp
2251 real(RP),
save :: c_in = 1.0_rp
2253 real(RP),
save :: nm_M92 = 1.e+3_rp
2254 real(RP),
save :: am_M92 = -0.639_rp
2255 real(RP),
save :: bm_M92 = 12.96_rp
2257 real(RP),
save :: in_max = 1000.e+3_rp
2258 real(RP),
save :: ssi_max= 0.60_rp
2259 real(RP),
save :: ssw_max= 1.1_rp
2261 logical,
save :: flag_first = .true.
2262 real(RP),
save :: qke_min = 0.03_rp
2263 real(RP),
save :: tem_ccn_low=233.150_rp
2264 real(RP),
save :: tem_in_low =173.150_rp
2265 logical,
save :: nucl_twomey = .false.
2266 logical,
save :: inucl_w = .false.
2268 namelist /nm_mp_sn14_nucleation/ &
2271 nm_m92, am_m92, bm_m92, &
2276 nucl_twomey, inucl_w
2281 real(RP) :: esw(KA,IA,JA)
2282 real(RP) :: esi(KA,IA,JA)
2283 real(RP) :: ssw(KA,IA,JA)
2284 real(RP) :: ssi(KA,IA,JA)
2286 real(RP) :: w_dssidz(KA,IA,JA)
2288 real(RP) :: ssi_below(KA,IA,JA)
2289 real(RP) :: z_below(KA,IA,JA)
2293 real(RP) :: qsw(KA,IA,JA)
2294 real(RP) :: qsi(KA,IA,JA)
2295 real(RP) :: dqsidtem_rho(KA,IA,JA)
2296 real(RP) :: dssidt_rad(KA,IA,JA)
2298 real(RP) :: wssi, wdssi
2304 real(RP) :: sigma_w(KA,IA,JA)
2305 real(RP) :: weff(KA,IA,JA)
2306 real(RP) :: weff_max(KA,IA,JA)
2308 real(RP) :: coef_ccn(IA,JA)
2309 real(RP) :: slope_ccn(IA,JA)
2310 real(RP) :: nc_new(KA,IA,JA)
2311 real(RP) :: nc_new_below(KA,IA,JA)
2313 real(RP) :: nc_new_max
2316 logical :: flag_nucleation(KA,IA,JA)
2318 real(RP) :: r_gravity
2319 real(RP),
parameter :: r_sqrt3=0.577350269_rp
2320 real(RP),
parameter :: eps=1.e-30_rp
2323 real(RP) :: dlcdt_max, dli_max
2324 real(RP) :: dncdt_max, dni_max
2329 if( flag_first )
then 2331 read(
io_fid_conf, nml=nm_mp_sn14_nucleation, end=100)
2335 if ( mp_couple_aerosol .AND. nucl_twomey )
then 2336 write(*,*)
"xxx [mp_sn14/nucleation] nucl_twomey should be false when MP_couple_aerosol is true, stop" 2348 nc_uplim_d(1,i,j) = c_ccn*1.5_rp
2353 r_gravity = 1.0_rp/grav
2355 call moist_psat_liq ( esw, tem )
2356 call moist_psat_ice ( esi, tem )
2357 call moist_pres2qsat_liq( qsw, tem, pre )
2358 call moist_pres2qsat_ice( qsi, tem, pre )
2359 call moist_dqsi_dtem_rho( dqsidtem_rho, tem, rho )
2362 a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
2370 pv = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)
2371 ssw(k,i,j) = min( mp_ssw_lim, ( pv/esw(k,i,j)-1.0_rp ) )*100.0_rp
2372 ssi(k,i,j) = (pv/esi(k,i,j) - 1.00_rp)
2374 ssi_below(k+1,i,j) = ssi(k,i,j)
2375 z_below(k+1,i,j) = z(k)
2378 ssi_below(
ks,i,j) = ssi(
ks,i,j)
2379 z_below(
ks,i,j) = z(
ks-1)
2384 coef_ccn(i,j) = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
2386 * (70.0_rp)**(kappa/(kappa + 2.0_rp))
2388 slope_ccn(i,j) = 1.5_rp*kappa/(kappa + 2.0_rp)
2391 sigma_w(k,i,j) = r_sqrt3*sqrt(max(qke(k,i,j),qke_min))
2393 sigma_w(
ks-1,i,j) = sigma_w(
ks,i,j)
2394 sigma_w(
ke+1,i,j) = sigma_w(
ke,i,j)
2397 weff(k,i,j) = 0.5_rp*(velz(k-1,i,j) + velz(k,i,j)) - cpa(k,i,j)*r_gravity*dtdt_rad(k,i,j)
2403 if( mp_couple_aerosol )
then 2408 if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp )
then 2409 nc_new(k,i,j) = max( ccn(k,i,j), c_ccn )
2411 nc_new(k,i,j) = 0.0_rp
2419 if( nucl_twomey )
then 2426 weff_max(k,i,j) = weff(k,i,j) + sigma_w(k,i,j)
2428 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 2430 nc_new_max = coef_ccn(i,j)*weff_max(k,i,j)**slope_ccn(i,j)
2431 nc_new(k,i,j) = a_max*nc_new_max**b_max
2433 nc_new(k,i,j) = 0.0_rp
2443 if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp )
then 2444 nc_new(k,i,j) = c_ccn*ssw(k,i,j)**kappa
2446 nc_new(k,i,j) = 0.0_rp
2459 if( nc_new(k,i,j) > nc_uplim_d(1,i,j) )
then 2460 flag_nucleation(k,i,j) = .false.
2461 nc_new_below(k+1,i,j) = 1.e+30_rp
2462 else if( nc_new(k,i,j) > eps )
then 2463 flag_nucleation(k,i,j) = .true.
2464 nc_new_below(k+1,i,j) = nc_new(k,i,j)
2466 flag_nucleation(k,i,j) = .false.
2467 nc_new_below(k+1,i,j) = 0.0_rp
2470 nc_new_below(
ks,i,j) = 0.0_rp
2483 if( mp_couple_aerosol )
then 2491 if ( flag_nucleation(k,i,j) .AND. &
2492 tem(k,i,j) > tem_ccn_low )
then 2493 dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2494 dncdt_max = dlcdt_max/xc_min
2496 dnc_new = nc_new(k,i,j)
2497 pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2498 pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2500 pq(i_ncccn,k,i,j) = 0.0_rp
2501 pq(i_lcccn,k,i,j) = 0.0_rp
2509 if( nucl_twomey )
then 2516 if ( flag_nucleation(k,i,j) .AND. &
2517 tem(k,i,j) > tem_ccn_low .AND. &
2518 nc_new(k,i,j) > rhoq(
i_nc,k,i,j) )
then 2519 dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2520 dncdt_max = dlcdt_max/xc_min
2521 dnc_new = nc_new(k,i,j)-rhoq(
i_nc,k,i,j)
2522 pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2523 pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2525 pq(i_ncccn,k,i,j) = 0.0_rp
2526 pq(i_lcccn,k,i,j) = 0.0_rp
2536 if( tem(k,i,j) > tem_ccn_low .AND. &
2537 nc_new(k,i,j) > rhoq(
i_nc,k,i,j) )
then 2538 dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2539 dncdt_max = dlcdt_max/xc_min
2540 dnc_new = nc_new(k,i,j)-rhoq(
i_nc,k,i,j)
2541 pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2542 pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2544 pq(i_ncccn,k,i,j) = 0.0_rp
2545 pq(i_lcccn,k,i,j) = 0.0_rp
2564 dz = z(k) - z_below(k,i,j)
2565 w_dssidz(k,i,j) = velz(k,i,j)*(ssi(k,i,j) - ssi_below(k,i,j))/dz
2566 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)
2567 dli_max = (rhoq(i_qv,k,i,j) - esi(k,i,j)/(rvap*tem(k,i,j)))*rdt
2568 dni_max = min( dli_max/xi_ccn, (in_max-rhoq(
i_ni,k,i,j))*rdt )
2569 wdssi = min( w_dssidz(k,i,j)+dssidt_rad(k,i,j), 0.01_rp)
2570 wssi = min( ssi(k,i,j), ssi_max)
2572 if( ( wdssi > eps ) .AND. &
2573 (tem(k,i,j) < 273.15_rp ) .AND. &
2574 (rhoq(
i_ni,k,i,j) < in_max ) .AND. &
2575 (wssi >= eps ) )
then 2578 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)
2580 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 )
2582 pq(i_liccn,k,i,j) = min(dli_max, pq(i_niccn,k,i,j)*xi_ccn )
2586 pq(i_niccn,k,i,j) = 0.0_rp
2587 pq(i_liccn,k,i,j) = 0.0_rp
2608 real(RP),
intent(inout):: PQ(PQ_MAX,KA,IA,JA)
2610 real(RP),
intent(in) :: Pac(Pac_MAX,KA,IA,JA)
2611 real(RP),
intent(in) :: tem(KA,IA,JA)
2612 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
2613 real(RP),
intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
2615 logical,
save :: flag_first = .true.
2617 real(RP),
parameter :: pice = 350.0e+6_rp
2619 real(RP),
parameter :: pnc = 250.0_rp
2620 real(RP),
parameter :: rc_cr= 12.e-6_rp
2624 real(RP),
save :: xc_cr
2625 real(RP),
save :: alpha
2626 real(RP),
save :: gm, lgm
2630 real(RP) :: a0,a1,a2,a3,a4,a5
2631 real(RP) :: a6,a7,a8,a9,a10
2632 real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
2633 real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
2634 real(RP),
parameter :: eps=1.0e-30_rp
2638 real(RP) :: wn, wni, wns, wng
2641 if( flag_first )
then 2642 flag_first = .false.
2644 xc_cr = (2.0_rp*rc_cr/a_m(i_mp_qc))**(1.0_rp/b_m(i_mp_qc))
2645 alpha = (nu(i_mp_qc)+1.0_rp)/mu(i_mp_qc)
2646 gm = gammafunc(alpha)
2656 if (tem(k,i,j) > 270.16_rp)
then 2658 else if(tem(k,i,j) >= 268.16_rp)
then 2659 fp = (270.16_rp-tem(k,i,j))*0.5_rp
2660 else if(tem(k,i,j) >= 265.16_rp)
then 2661 fp = (tem(k,i,j)-265.16_rp)*0.333333333_rp
2669 x = coef_lambda(i_mp_qc)*(xc_cr/xq(i_mp_qc,k,i,j))**mu(i_mp_qc)
2671 if(x<1.e-2_rp*alpha)
then 2673 else if(x<alpha+1.0_rp)
then 2676 a1 = a0*x/(alpha+1.0_rp)
2677 a2 = a1*x/(alpha+2.0_rp)
2678 a3 = a2*x/(alpha+3.0_rp)
2679 a4 = a3*x/(alpha+4.0_rp)
2680 a5 = a4*x/(alpha+5.0_rp)
2681 a6 = a5*x/(alpha+6.0_rp)
2682 a7 = a6*x/(alpha+7.0_rp)
2683 a8 = a7*x/(alpha+8.0_rp)
2684 a9 = a8*x/(alpha+9.0_rp)
2685 a10 = a9*x/(alpha+10.0_rp)
2686 igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
2687 else if(x<alpha*1.d2)
then 2695 an1 = -(1.0_rp-alpha)
2697 d1 = 1.0_rp/(an1*d0+b1)
2702 an2 = -2.0_rp*(2.0_rp-alpha)
2704 d2 = 1.0_rp/(an2*d1+b2)
2709 igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
2714 n12 = rhoq(
i_nc,k,i,j)*(1.0_rp-igm)
2716 wn = (pice + n12/((rhoq(
i_qc,k,i,j)+xc_min)*pnc) )*fp
2717 wni = wn*(-pac(i_liaclc2li,k,i,j) )
2718 wns = wn*(-pac(i_lsaclc2ls,k,i,j) )
2719 wng = wn*(-pac(i_lgaclc2lg,k,i,j) )
2720 pq(i_nispl,k,i,j) = wni+wns+wng
2722 pq(i_lsspl,k,i,j) = - wns*xq(i_mp_qi,k,i,j)
2723 pq(i_lgspl,k,i,j) = - wng*xq(i_mp_qi,k,i,j)
2733 ! collection process
2736 xq, dq_xave, vt_xave )
2741 moist_psat_ice => atmos_saturation_psat_ice
2747 real(RP),
intent(out):: Pac(Pac_MAX,KA,IA,JA)
2749 real(RP),
intent(inout):: PQ(PQ_MAX,KA,IA,JA)
2751 real(RP),
intent(in) :: wtem(KA,IA,JA)
2753 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
2755 real(RP),
intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
2757 real(RP),
intent(in) :: dq_xave(HYDRO_MAX,KA,IA,JA)
2759 real(RP),
intent(in) :: vt_xave(HYDRO_MAX,2,KA,IA,JA)
2766 real(RP),
save :: dc0 = 15.0e-6_rp
2767 real(RP),
save :: dc1 = 40.0e-6_rp
2768 real(RP),
save :: di0 = 150.0e-6_rp
2769 real(RP),
save :: ds0 = 150.0e-6_rp
2770 real(RP),
save :: dg0 = 150.0e-6_rp
2772 real(RP),
save :: sigma_c=0.00_rp
2773 real(RP),
save :: sigma_r=0.00_rp
2774 real(RP),
save :: sigma_i=0.2_rp
2775 real(RP),
save :: sigma_s=0.2_rp
2776 real(RP),
save :: sigma_g=0.00_rp
2778 real(RP),
save :: E_im = 0.80_rp
2779 real(RP),
save :: E_sm = 0.80_rp
2780 real(RP),
save :: E_gm = 1.00_rp
2782 real(RP),
save :: E_ir=1.0_rp
2783 real(RP),
save :: E_sr=1.0_rp
2784 real(RP),
save :: E_gr=1.0_rp
2785 real(RP),
save :: E_ii=1.0_rp
2786 real(RP),
save :: E_si=1.0_rp
2787 real(RP),
save :: E_gi=1.0_rp
2788 real(RP),
save :: E_ss=1.0_rp
2789 real(RP),
save :: E_gs=1.0_rp
2790 real(RP),
save :: E_gg=1.0_rp
2793 integer,
save :: i_iconv2g=1
2794 integer,
save :: i_sconv2g=1
2796 real(RP),
save :: rho_g = 900.0_rp
2798 real(RP),
save :: cfill_i = 0.68_rp
2799 real(RP),
save :: cfill_s = 0.01_rp
2801 real(RP),
save :: di_cri = 500.e-6_rp
2803 logical,
save :: opt_stick_KS96=.false.
2804 logical,
save :: opt_stick_CO86=.false.
2805 real(RP),
parameter :: a_dec = 0.883_rp
2806 real(RP),
parameter :: b_dec = 0.093_rp
2807 real(RP),
parameter :: c_dec = 0.00348_rp
2808 real(RP),
parameter :: d_dec = 4.5185e-5_rp
2810 logical,
save :: flag_first = .true.
2811 namelist /nm_mp_sn14_collection/ &
2812 dc0, dc1, di0, ds0, dg0, &
2813 sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
2817 e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
2818 i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
2820 real(RP) :: tem(KA,IA,JA)
2823 real(RP) :: E_c, E_r, E_i, E_s, E_g
2824 real(RP) :: E_ic, E_sc, E_gc
2826 real(RP) :: E_stick(KA,IA,JA)
2828 real(RP) :: temc, temc2, temc3
2831 real(RP) :: esi(KA,IA,JA)
2833 real(RP) :: temc_p, temc_m
2845 real(RP) :: coef_acc_LCI, coef_acc_NCI
2846 real(RP) :: coef_acc_LCS, coef_acc_NCS
2848 real(RP) :: coef_acc_LCG, coef_acc_NCG
2849 real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I
2850 real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R
2851 real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S
2852 real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R
2853 real(RP) :: coef_acc_LRG, coef_acc_NRG
2854 real(RP) :: coef_acc_LII, coef_acc_NII
2855 real(RP) :: coef_acc_LIS, coef_acc_NIS
2856 real(RP) :: coef_acc_NSS
2857 real(RP) :: coef_acc_NGG
2858 real(RP) :: coef_acc_LSG, coef_acc_NSG
2860 real(RP) :: dcdc, dcdi, dcds, dcdg
2861 real(RP) :: drdr, drdi, drds, drdg
2862 real(RP) :: didi, dids, didg
2863 real(RP) :: dsds, dsdg
2866 real(RP) :: vcvc, vcvi, vcvs, vcvg
2867 real(RP) :: vrvr, vrvi, vrvs, vrvg
2868 real(RP) :: vivi, vivs, vivg
2869 real(RP) :: vsvs, vsvg
2872 real(RP) :: wx_cri, wx_crs
2873 real(RP) :: coef_emelt
2880 if( flag_first )
then 2882 read(
io_fid_conf, nml=nm_mp_sn14_collection, end=100 )
2884 flag_first = .false.
2891 tem(k,i,j) = max( wtem(k,i,j), tem_min )
2896 call moist_psat_ice( esi, tem )
2898 if( opt_stick_ks96 )
then 2903 temc = tem(k,i,j) - t00
2906 e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
2907 esi_rat = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)/esi(k,i,j)
2908 e_stick(k,i,j) = min(1.0_rp, e_dec*esi_rat)
2912 else if( opt_stick_co86 )
then 2917 temc = min(tem(k,i,j) - t00,0.0_rp)
2918 w1 = 0.035_rp*temc-0.7_rp
2919 e_stick(k,i,j) = 10._rp**w1
2928 temc_m = min(tem(k,i,j) - t00,0.0_rp)
2929 e_stick(k,i,j) = exp(0.09_rp*temc_m)
2935 profile_start(
"sn14_collection")
2942 temc_p = max(tem(k,i,j) - t00,0.0_rp)
2944 ave_dc = coef_d(i_mp_qc)*xq(i_mp_qc,k,i,j)**b_m(i_mp_qc)
2945 ave_di = coef_d(i_mp_qi)*xq(i_mp_qi,k,i,j)**b_m(i_mp_qi)
2946 ave_ds = coef_d(i_mp_qs)*xq(i_mp_qs,k,i,j)**b_m(i_mp_qs)
2947 ave_dg = coef_d(i_mp_qg)*xq(i_mp_qg,k,i,j)**b_m(i_mp_qg)
2950 e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
2951 sw = 0.5_rp - sign(0.5_rp, di0-ave_di)
2953 sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds)
2955 sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg)
2962 dcdc = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qc,k,i,j)
2963 drdr = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qr,k,i,j)
2964 didi = dq_xave(i_mp_qi,k,i,j) * dq_xave(i_mp_qi,k,i,j)
2965 dsds = dq_xave(i_mp_qs,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2966 dgdg = dq_xave(i_mp_qg,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2967 dcdi = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qi,k,i,j)
2968 dcds = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2969 dcdg = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2970 drdi = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qi,k,i,j)
2971 drds = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2972 drdg = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2973 dids = dq_xave(i_mp_qi,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2975 dsdg = dq_xave(i_mp_qs,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2977 vcvc = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qc,2,k,i,j)
2978 vrvr = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qr,2,k,i,j)
2979 vivi = vt_xave(i_mp_qi,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
2980 vsvs = vt_xave(i_mp_qs,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2981 vgvg = vt_xave(i_mp_qg,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2982 vcvi = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
2983 vcvs = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2984 vcvg = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2985 vrvi = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
2986 vrvs = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2987 vrvg = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2988 vivs = vt_xave(i_mp_qi,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2990 vsvg = vt_xave(i_mp_qs,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2999 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3000 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3001 + sigma_i + sigma_c )**0.5_rp
3003 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3004 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3005 + sigma_i + sigma_c )**0.5_rp
3006 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
3007 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
3011 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3012 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3013 + sigma_s + sigma_c )**0.5_rp
3015 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3016 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3017 + sigma_s + sigma_c )**0.5_rp
3018 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
3019 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
3023 ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3024 * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3025 + sigma_g + sigma_c )**0.5_rp
3027 ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3028 * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3029 + sigma_g + sigma_c )**0.5_rp
3030 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
3031 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
3034 ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3035 * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3036 + sigma_g + sigma_s )**0.5_rp
3038 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3040 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3041 + sigma_g + sigma_s )**0.5_rp
3042 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
3043 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
3048 ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
3049 * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
3050 + sigma_i + sigma_s )**0.5_rp
3052 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
3053 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
3054 + sigma_i + sigma_s )**0.5_rp
3055 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
3056 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
3058 sw = sign(0.5_rp, t00-tem(k,i,j)) + 0.5_rp
3068 ( ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) * sw &
3069 + ( 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) ) &
3070 * sqrt( ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg ) * sw &
3071 + ( 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) &
3072 + sigma_r + sigma_g )
3073 pac(i_lraclg2lg,k,i,j) = -0.25_rp*pi*e_gr*coef_acc_lrg &
3074 * ( rhoq(i_ng,k,i,j)*rhoq(
i_qr,k,i,j) * sw &
3075 + rhoq(
i_nr,k,i,j)*rhoq(
i_qg,k,i,j) * (1.0_rp-sw) )
3077 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) &
3078 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg &
3079 + sigma_r + sigma_g )**0.5_rp
3080 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
3090 ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
3091 * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
3092 + sigma_r + sigma_i )**0.5_rp
3094 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
3095 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
3096 + sigma_r + sigma_i )**0.5_rp
3097 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
3098 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
3101 ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
3102 * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
3103 + sigma_r + sigma_i )**0.5_rp
3105 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
3106 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
3107 + sigma_r + sigma_i )**0.5_rp
3108 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
3109 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
3113 ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
3114 * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
3115 + sigma_r + sigma_s )**0.5_rp
3117 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
3118 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
3119 + sigma_r + sigma_s )**0.5_rp
3120 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
3121 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
3124 ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
3125 * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
3126 + sigma_r + sigma_s )**0.5_rp
3128 ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
3129 * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
3130 + sigma_r + sigma_s )**0.5_rp
3131 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
3132 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
3140 ( delta_b0(i_mp_qi)*didi + delta_ab1(i_mp_qi,i_mp_qi)*didi + delta_b1(i_mp_qi)*didi ) &
3141 * ( theta_b0(i_mp_qi)*vivi - theta_ab1(i_mp_qi,i_mp_qi)*vivi + theta_b1(i_mp_qi)*vivi &
3142 + sigma_i + sigma_i )**0.5_rp
3144 ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qi,i_mp_qi)*didi + delta_b0(i_mp_qi)*didi ) &
3145 * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qi,i_mp_qi)*vivi + theta_b0(i_mp_qi)*vivi &
3146 + sigma_i + sigma_i )**0.5_rp
3147 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
3148 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
3160 ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qs,i_mp_qs)*dsds + delta_b0(i_mp_qs)*dsds ) &
3161 * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qs,i_mp_qs)*vsvs + theta_b0(i_mp_qs)*vsvs &
3162 + sigma_s + sigma_s )**0.5_rp
3163 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
3167 ( delta_b0(i_mp_qg)*dgdg + delta_ab0(i_mp_qg,i_mp_qg)*dgdg + delta_b0(i_mp_qg)*dgdg ) &
3168 * ( theta_b0(i_mp_qg)*vgvg - theta_ab0(i_mp_qg,i_mp_qg)*vgvg + theta_b0(i_mp_qg)*vgvg &
3169 + sigma_g + sigma_g )**0.5_rp
3170 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
3177 sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di)
3178 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
3179 pq(i_licon,k,i,j) = i_iconv2g * pac(i_liaclc2li,k,i,j)/max(1.0_rp, wx_cri) * sw
3180 pq(i_nicon,k,i,j) = i_iconv2g * pq(i_licon,k,i,j)/xq(i_mp_qi,k,i,j) * sw
3183 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 )
3184 pq(i_lscon,k,i,j) = i_sconv2g * (pac(i_lsaclc2ls,k,i,j))/max(1.0_rp, wx_crs)
3185 pq(i_nscon,k,i,j) = i_sconv2g * pq(i_lscon,k,i,j)/xq(i_mp_qs,k,i,j)
3194 coef_emelt = cl/lhf0*temc_p
3196 pq(i_lgacm,k,i,j) = coef_emelt*pac(i_lgaclc2lg,k,i,j)
3197 pq(i_ngacm,k,i,j) = pq(i_lgacm,k,i,j)/xq(i_mp_qg,k,i,j)
3199 pq(i_lgarm,k,i,j) = coef_emelt*pac(i_lraclg2lg,k,i,j)
3200 pq(i_ngarm,k,i,j) = pq(i_lgarm,k,i,j)/xq(i_mp_qg,k,i,j)
3202 pq(i_lsacm,k,i,j) = coef_emelt*(pac(i_lsaclc2ls,k,i,j))
3203 pq(i_nsacm,k,i,j) = pq(i_lsacm,k,i,j)/xq(i_mp_qs,k,i,j)
3205 pq(i_lsarm,k,i,j) = coef_emelt*(pac(i_lracls2lg_r,k,i,j)+pac(i_lracls2lg_s,k,i,j))
3206 pq(i_nsarm,k,i,j) = pq(i_lsarm,k,i,j)/xq(i_mp_qg,k,i,j)
3208 pq(i_liacm,k,i,j) = coef_emelt*pac(i_liaclc2li,k,i,j)
3209 pq(i_niacm,k,i,j) = pq(i_liacm,k,i,j)/xq(i_mp_qi,k,i,j)
3211 pq(i_liarm,k,i,j) = coef_emelt*(pac(i_lracli2lg_r,k,i,j)+pac(i_lracli2lg_i,k,i,j))
3212 pq(i_niarm,k,i,j) = pq(i_liarm,k,i,j)/xq(i_mp_qg,k,i,j)
3216 profile_stop(
"sn14_collection")
3225 rhoq, xq, dq_xave, &
3231 real(RP),
intent(inout) :: PQ(PQ_MAX,KA,IA,JA)
3233 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3234 real(RP),
intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
3235 real(RP),
intent(in) :: dq_xave(HYDRO_MAX,KA,IA,JA)
3236 real(RP),
intent(in) :: rho(KA,IA,JA)
3239 real(RP),
parameter :: kcc = 4.44e+9_rp
3240 real(RP),
parameter :: tau_min = 1.e-20_rp
3241 real(RP),
parameter :: rx_sep = 1.0_rp/x_sep
3244 real(RP),
parameter :: kcr = 5.8_rp
3245 real(RP),
parameter :: thr_acc = 5.e-5_rp
3248 real(RP),
parameter :: krr = 4.33_rp
3249 real(RP),
parameter :: kaprr = 60.7_rp
3250 real(RP),
parameter :: kbr = 1000._rp
3251 real(RP),
parameter :: kapbr = 2.3e+3_rp
3252 real(RP),
parameter :: dr_min = 0.35e-3_rp
3255 real(RP) :: coef_nuc0
3256 real(RP) :: coef_nuc1
3257 real(RP) :: coef_aut0
3258 real(RP) :: coef_aut1
3269 coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
3270 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)
3271 coef_aut0 = -kcc*coef_nuc0
3272 coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
3277 lwc = rhoq(
i_qr,k,i,j)+rhoq(
i_qc,k,i,j)
3278 if( lwc > xc_min )
then 3279 tau = max(tau_min, rhoq(
i_qr,k,i,j)/lwc)
3283 rho_fac = sqrt(rho_0/max(rho(k,i,j),rho_min))
3286 psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3
3287 pq(i_ncaut,k,i,j) = coef_aut0*rhoq(
i_qc,k,i,j)*rhoq(
i_qc,k,i,j)*rho_fac*rho_fac
3289 pq(i_lcaut,k,i,j) = coef_aut1*lwc*lwc*xq(i_mp_qc,k,i,j)*xq(i_mp_qc,k,i,j) &
3290 *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac
3291 pq(i_nraut,k,i,j) = -rx_sep*pq(i_lcaut,k,i,j)
3294 psi_acc =(tau/(tau+thr_acc))**4
3295 pq(i_lcacc,k,i,j) = -kcr*rhoq(
i_qc,k,i,j)*rhoq(
i_qr,k,i,j)*rho_fac*psi_acc
3296 pq(i_ncacc,k,i,j) = -kcr*rhoq(
i_nc,k,i,j)*rhoq(
i_qr,k,i,j)*rho_fac*psi_acc
3299 pq(i_nrslc,k,i,j) = -krr*rhoq(
i_nr,k,i,j)*rhoq(
i_qr,k,i,j)*rho_fac
3302 ddr = min(1.e-3_rp, dq_xave(i_mp_qr,k,i,j) - dr_eq )
3303 if (dq_xave(i_mp_qr,k,i,j) < dr_min )
then 3305 pq(i_nrbrk,k,i,j) = 0.0_rp
3306 else if (dq_xave(i_mp_qr,k,i,j) <= dr_eq )
then 3307 psi_brk = kbr*ddr + 1.0_rp
3308 pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j)
3310 psi_brk = 2.0_rp*exp(kapbr*ddr) - 1.0_rp
3311 pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j)
3323 rho, tem, pre, qd, & ! in
3326 xq, vt_xave, dq_xave )
3334 real(RP),
intent(inout) :: PQ(PQ_MAX,KA,IA,JA)
3336 real(RP),
intent(in) :: rho(KA,IA,JA)
3337 real(RP),
intent(in) :: tem(KA,IA,JA)
3338 real(RP),
intent(in) :: pre(KA,IA,JA)
3339 real(RP),
intent(in) :: qd(KA,IA,JA)
3340 real(RP),
intent(in) :: esw(KA,IA,JA)
3341 real(RP),
intent(in) :: esi(KA,IA,JA)
3342 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3343 real(RP),
intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
3347 real(RP),
intent(in) :: vt_xave(HYDRO_MAX,2,KA,IA,JA)
3349 real(RP),
intent(in) :: dq_xave(HYDRO_MAX,KA,IA,JA)
3352 real(RP) :: temc_lim
3359 real(RP) :: nua, r_nua
3365 real(RP) :: Gwr, Gii, Gis, Gig
3370 real(RP) :: Nrers_r2, Nreis_r2
3371 real(RP) :: Nress_r2, Nregs_r2
3373 real(RP) :: Nrerl_r2, Nreil_r2
3374 real(RP) :: Nresl_r2, Nregl_r2
3375 real(RP) :: NscNrer_s, NscNrer_l
3376 real(RP) :: NscNrei_s, NscNrei_l
3377 real(RP) :: NscNres_s, NscNres_l
3378 real(RP) :: NscNreg_s, NscNreg_l
3379 real(RP) :: ventLR_s, ventLR_l
3380 real(RP) :: ventNI_s, ventNI_l, ventLI_s, ventLI_l
3381 real(RP) :: ventNS_s, ventNS_l, ventLS_s, ventLS_l
3382 real(RP) :: ventNG_s, ventNG_l, ventLG_s, ventLG_l
3384 real(RP) :: wtr, wti, wts, wtg
3385 real(RP),
parameter :: r_14=1.0_rp/1.4_rp
3386 real(RP),
parameter :: r_15=1.0_rp/1.5_rp
3389 real(RP) :: ventNI, ventLI
3390 real(RP) :: ventNS, ventLS
3391 real(RP) :: ventNG, ventLG
3393 real(RP),
parameter :: Re_max=1.e+3_rp
3394 real(RP),
parameter :: Re_min=1.e-4_rp
3406 profile_start(
"sn14_dep_vapor")
3410 temc = tem(k,i,j) - t00
3411 temc_lim= max(temc, -40._rp )
3412 rho_lim = max(rho(k,i,j),rho_min)
3413 qv = rhoq(i_qv,k,i,j)/rho_lim
3414 pre_lim = rho_lim*(qd(k,i,j)*rdry + qv*rvap)*(temc_lim+t00)
3422 dw = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(p00/pre_lim)
3423 kalfa = ka0 + temc_lim*dka_dt
3424 mua = mua0 + temc_lim*dmua_dt
3427 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))
3428 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))
3430 gwr = 4.0_rp*pi/cap(i_mp_qr)/gw
3431 gii = 4.0_rp*pi/cap(i_mp_qi)/gi
3432 gis = 4.0_rp*pi/cap(i_mp_qs)/gi
3433 gig = 4.0_rp*pi/cap(i_mp_qg)/gi
3436 nsc_r3 = (nua/dw)**(0.33333333_rp)
3439 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)))
3440 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)))
3441 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)))
3442 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)))
3445 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)))
3446 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)))
3447 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)))
3448 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)))
3449 nscnrer_s=nsc_r3*nrers_r2
3450 nscnrer_l=nsc_r3*nrerl_r2
3452 nscnrei_s=nsc_r3*nreis_r2
3453 nscnrei_l=nsc_r3*nreil_r2
3455 nscnres_s=nsc_r3*nress_r2
3456 nscnres_l=nsc_r3*nresl_r2
3458 nscnreg_s=nsc_r3*nregs_r2
3459 nscnreg_l=nsc_r3*nregl_r2
3461 ventlr_s = ah_vent1(i_mp_qr,1) + bh_vent1(i_mp_qr,1)*nscnrer_s
3462 ventlr_l = ah_vent1(i_mp_qr,2) + bh_vent1(i_mp_qr,2)*nscnrer_l
3464 ventni_s = ah_vent0(i_mp_qi,1) + bh_vent0(i_mp_qi,1)*nscnrei_s
3465 ventni_l = ah_vent0(i_mp_qi,2) + bh_vent0(i_mp_qi,2)*nscnrei_l
3466 ventli_s = ah_vent1(i_mp_qi,1) + bh_vent1(i_mp_qi,1)*nscnrei_s
3467 ventli_l = ah_vent1(i_mp_qi,2) + bh_vent1(i_mp_qi,2)*nscnrei_l
3469 ventns_s = ah_vent0(i_mp_qs,1) + bh_vent0(i_mp_qs,1)*nscnres_s
3470 ventns_l = ah_vent0(i_mp_qs,2) + bh_vent0(i_mp_qs,2)*nscnres_l
3471 ventls_s = ah_vent1(i_mp_qs,1) + bh_vent1(i_mp_qs,1)*nscnres_s
3472 ventls_l = ah_vent1(i_mp_qs,2) + bh_vent1(i_mp_qs,2)*nscnres_l
3474 ventng_s = ah_vent0(i_mp_qg,1) + bh_vent0(i_mp_qg,1)*nscnreg_s
3475 ventng_l = ah_vent0(i_mp_qg,2) + bh_vent0(i_mp_qg,2)*nscnreg_l
3476 ventlg_s = ah_vent1(i_mp_qg,1) + bh_vent1(i_mp_qg,1)*nscnreg_s
3477 ventlg_l = ah_vent1(i_mp_qg,2) + bh_vent1(i_mp_qg,2)*nscnreg_l
3481 wtr = ( min(max( nscnrer_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3482 wti = ( min(max( nscnrei_s , 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3483 wts = ( min(max( nscnres_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3484 wtg = ( min(max( nscnreg_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15
3486 ventni = (1.0_rp-wti)*ventni_s + wti*ventni_l
3487 ventns = (1.0_rp-wts)*ventns_s + wts*ventns_l
3488 ventng = (1.0_rp-wtg)*ventng_s + wtg*ventng_l
3490 ventlr = (1.0_rp-wtr)*ventlr_s + wtr*ventlr_l
3491 ventli = (1.0_rp-wti)*ventli_s + wti*ventli_l
3492 ventls = (1.0_rp-wts)*ventls_s + wts*ventls_l
3493 ventlg = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
3512 pq(i_lcdep,k,i,j) = gwr*rhoq(
i_nc,k,i,j)*dq_xave(i_mp_qc,k,i,j)*coef_deplc
3513 pq(i_lrdep,k,i,j) = gwr*rhoq(
i_nr,k,i,j)*dq_xave(i_mp_qr,k,i,j)*ventlr
3514 pq(i_lidep,k,i,j) = gii*rhoq(
i_ni,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventli
3515 pq(i_lsdep,k,i,j) = gis*rhoq(
i_ns,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventls
3516 pq(i_lgdep,k,i,j) = gig*rhoq(i_ng,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventlg
3517 pq(i_nrdep,k,i,j) = pq(i_lrdep,k,i,j)/xq(i_mp_qr,k,i,j)
3518 pq(i_nidep,k,i,j) = 0.0_rp
3519 pq(i_nsdep,k,i,j) = pq(i_lsdep,k,i,j)/xq(i_mp_qs,k,i,j)
3520 pq(i_ngdep,k,i,j) = pq(i_lgdep,k,i,j)/xq(i_mp_qg,k,i,j)
3528 dt = kalfa/(cpvap*rho_0)
3534 gm = 2.0_rp*pi/emelt&
3535 * ( (kalfa*dt/dw)*(temc) + (dw*lhs0/rvap)*(esi(k,i,j)/tem(k,i,j)-psat0/t00) )
3540 sw = ( sign(0.5_rp,temc) + 0.5_rp ) * ( sign(0.5_rp,gm-eps) + 0.5_rp )
3544 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
3546 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
3547 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
3549 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
3550 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
3552 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
3557 profile_stop(
"sn14_dep_vapor")
3573 real(RP),
intent(in) :: dt
3574 real(RP),
intent(inout):: PQ(PQ_MAX,KA,IA,JA)
3576 real(RP),
intent(in) :: tem(KA,IA,JA)
3578 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3579 real(RP),
intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
3581 real(RP),
parameter :: temc_min = -65.0_rp
3582 real(RP),
parameter :: a_het = 0.2_rp
3583 real(RP),
parameter :: b_het = 0.65_rp
3585 real(RP) :: coef_m2_c
3586 real(RP) :: coef_m2_r
3588 real(RP) :: temc, temc2, temc3, temc4
3590 real(RP) :: Jhom, Jhet
3598 coef_m2_c = coef_m2(i_mp_qc)
3599 coef_m2_r = coef_m2(i_mp_qr)
3601 profile_start(
"sn14_freezing")
3605 temc = max( tem(k,i,j) - t00, temc_min )
3608 jhet = a_het*exp( -b_het*temc - 1.0_rp )
3611 if( temc < -65.0_rp )
then 3612 jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
3613 jhet = a_het*exp( 65.0_rp*b_het - 1.0_rp )
3614 else if( temc < -30.0_rp )
then 3619 - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
3620 - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
3621 else if( temc < 0.0_rp)
then 3622 jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
3634 pq(i_lchom,k,i,j) = 0.0_rp
3635 pq(i_nchom,k,i,j) = 0.0_rp
3637 #if defined(PGI) || defined(SX) 3638 tmp = min( xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp)
3639 pq(i_lchet,k,i,j) = -rdt*rhoq(
i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
3640 pq(i_nchet,k,i,j) = -rdt*rhoq(
i_nc,k,i,j)*( 1.0_rp - exp( - tmp ) )
3642 tmp = min( xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp)
3643 pq(i_lrhet,k,i,j) = -rdt*rhoq(
i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
3644 pq(i_nrhet,k,i,j) = -rdt*rhoq(
i_nr,k,i,j)*( 1.0_rp - exp( - tmp ) )
3646 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 ) )
3647 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 ) )
3648 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 ) )
3649 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 ) )
3654 profile_stop(
"sn14_freezing")
3659 subroutine mp_terminal_velocity( &
3671 real(RP),
intent(out) :: velw(KA,IA,JA,QA_MP-1)
3672 real(RP),
intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3673 real(RP),
intent(in) :: DENS(KA,IA,JA)
3674 real(RP),
intent(in) :: temp(KA,IA,JA)
3675 real(RP),
intent(in) :: pres(KA,IA,JA)
3680 real(RP) :: rhofac_q(HYDRO_MAX)
3682 real(RP) :: rlambdar
3692 integer :: k, i, j, iq
3695 profile_start(
"sn14_terminal_vel")
3697 mud_r = 3.0_rp * nu(i_mp_qr) + 2.0_rp
3703 rhofac = rho_0 / max( dens(k,i,j), rho_min )
3706 rhofac_q(i_mp_qc) = rhofac ** gamma_v(i_mp_qc)
3707 xq = max( xqmin(i_mp_qc), min( xqmax(i_mp_qc), rhoq(
i_qc,k,i,j) / ( rhoq(
i_nc,k,i,j) + nqmin(i_mp_qc) ) ) )
3709 velw(k,i,j,i_mp_qc) = -rhofac_q(i_mp_qc) * coef_vt1(i_mp_qc,1) * xq**beta_v(i_mp_qc,1)
3711 velw(k,i,j,i_mp_nc) = -rhofac_q(i_mp_qc) * coef_vt0(i_mp_qc,1) * xq**beta_vn(i_mp_qc,1)
3714 rhofac_q(i_mp_qr) = rhofac ** gamma_v(i_mp_qr)
3715 xq = max( xqmin(i_mp_qr), min( xqmax(i_mp_qr), rhoq(
i_qr,k,i,j) / ( rhoq(
i_nr,k,i,j) + nqmin(i_mp_qr) ) ) )
3717 rlambdar = a_m(i_mp_qr) * xq**b_m(i_mp_qr) &
3718 * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
3719 dq = ( 4.0_rp + mud_r ) * rlambdar
3721 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
3723 velq_s = coef_vtr_ar2 * dq &
3724 * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-5.0_rp-mud_r) )
3725 velq_l = coef_vtr_ar1 - coef_vtr_br1 &
3726 * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-4.0_rp-mud_r)
3727 velw(k,i,j,i_mp_qr) = -rhofac_q(i_mp_qr) &
3728 * ( velq_l * ( weight ) &
3729 + velq_s * ( 1.0_rp - weight ) )
3731 dq = ( 1.0_rp + mud_r ) * rlambdar
3732 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
3734 velq_s = coef_vtr_ar2 * dql &
3735 * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-2.0_rp-mud_r) )
3736 velq_l = coef_vtr_ar1 - coef_vtr_br1 &
3737 * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-1.0_rp-mud_r)
3738 velw(k,i,j,i_mp_nr) = -rhofac_q(i_mp_qr) &
3739 * ( velq_l * ( weight ) &
3740 + velq_s * ( 1.0_rp - weight ) )
3743 rhofac_q(i_mp_qi) = ( pres(k,i,j)/pre0_vt )**a_pre0_vt * ( temp(k,i,j)/tem0_vt )**a_tem0_vt
3744 xq = max( xqmin(i_mp_qi), min( xqmax(i_mp_qi), rhoq(
i_qi,k,i,j) / ( rhoq(
i_ni,k,i,j) + nqmin(i_mp_qi) ) ) )
3746 tmp = a_m(i_mp_qi) * xq**b_m(i_mp_qi)
3747 dq = coef_dave_l(i_mp_qi) * tmp
3748 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_li ) ) ) )
3750 velq_s = coef_vt1(i_mp_qi,1) * xq**beta_v(i_mp_qi,1)
3751 velq_l = coef_vt1(i_mp_qi,2) * xq**beta_v(i_mp_qi,2)
3752 velw(k,i,j,i_mp_qi) = -rhofac_q(i_mp_qi) &
3753 * ( velq_l * ( weight ) &
3754 + velq_s * ( 1.0_rp - weight ) )
3756 dq = coef_dave_n(i_mp_qi) * tmp
3757 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ni ) ) ) )
3759 velq_s = coef_vt0(i_mp_qi,1) * xq**beta_vn(i_mp_qi,1)
3760 velq_l = coef_vt0(i_mp_qi,2) * xq**beta_vn(i_mp_qi,2)
3761 velw(k,i,j,i_mp_ni) = -rhofac_q(i_mp_qi) &
3762 * ( velq_l * ( weight ) &
3763 + velq_s * ( 1.0_rp - weight ) )
3766 rhofac_q(i_mp_qs) = rhofac_q(i_mp_qi)
3767 xq = max( xqmin(i_mp_qs), min( xqmax(i_mp_qs), rhoq(
i_qs,k,i,j) / ( rhoq(
i_ns,k,i,j) + nqmin(i_mp_qs) ) ) )
3769 tmp = a_m(i_mp_qs) * xq**b_m(i_mp_qs)
3770 dq = coef_dave_l(i_mp_qs) * tmp
3771 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ls ) ) ) )
3773 velq_s = coef_vt1(i_mp_qs,1) * xq**beta_v(i_mp_qs,1)
3774 velq_l = coef_vt1(i_mp_qs,2) * xq**beta_v(i_mp_qs,2)
3775 velw(k,i,j,i_mp_qs) = -rhofac_q(i_mp_qs) &
3776 * ( velq_l * ( weight ) &
3777 + velq_s * ( 1.0_rp - weight ) )
3779 dq = coef_dave_n(i_mp_qs) * tmp
3780 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ns ) ) ) )
3782 velq_s = coef_vt0(i_mp_qs,1) * xq**beta_vn(i_mp_qs,1)
3783 velq_l = coef_vt0(i_mp_qs,2) * xq**beta_vn(i_mp_qs,2)
3784 velw(k,i,j,i_mp_ns) = -rhofac_q(i_mp_qs) &
3785 * ( velq_l * ( weight ) &
3786 + velq_s * ( 1.0_rp - weight ) )
3789 rhofac_q(i_mp_qg) = rhofac_q(i_mp_qi)
3790 xq = max( xqmin(i_mp_qg), min( xqmax(i_mp_qg), rhoq(
i_qg,k,i,j) / ( rhoq(i_ng,k,i,j) + nqmin(i_mp_qg) ) ) )
3792 tmp = a_m(i_mp_qg) * xq**b_m(i_mp_qg)
3793 dq = coef_dave_l(i_mp_qg) * tmp
3794 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_lg ) ) ) )
3796 velq_s = coef_vt1(i_mp_qg,1) * xq**beta_v(i_mp_qg,1)
3797 velq_l = coef_vt1(i_mp_qg,2) * xq**beta_v(i_mp_qg,2)
3798 velw(k,i,j,i_mp_qg) = -rhofac_q(i_mp_qg) &
3799 * ( velq_l * ( weight ) &
3800 + velq_s * ( 1.0_rp - weight ) )
3802 dq = coef_dave_n(i_mp_qg) * tmp
3803 weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ng ) ) ) )
3805 velq_s = coef_vt0(i_mp_qg,1) * xq**beta_vn(i_mp_qg,1)
3806 velq_l = coef_vt0(i_mp_qg,2) * xq**beta_vn(i_mp_qg,2)
3807 velw(k,i,j,i_mp_ng) = -rhofac_q(i_mp_qg) &
3808 * ( velq_l * ( weight ) &
3809 + velq_s * ( 1.0_rp - weight ) )
3818 velw(
ks-1,i,j,iq) = velw(
ks,i,j,iq)
3824 profile_stop(
"sn14_terminal_vel")
3827 end subroutine mp_terminal_velocity
3830 ntdiv, ntmax, & ! in [Add] 10/08/03
3842 esw, esi, rhoq2, & ! in
3844 qc_evaporate, & ! in
3846 sl_PLRdep, sl_PNRdep )
3854 moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
3855 moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
3862 integer,
intent(in) :: ntdiv
3863 integer,
intent(in) :: ntmax
3865 real(RP),
intent(in) :: dt
3866 real(RP),
intent(in) :: gsgam2(KA,IA,JA)
3867 real(RP),
intent(in) :: z(KA)
3868 real(RP),
intent(in) :: dz(KA)
3869 real(RP),
intent(in) :: velz(KA,IA,JA)
3870 real(RP),
intent(in) :: dTdt_rad(KA,IA,JA)
3871 real(RP),
intent(in) :: rho(KA,IA,JA)
3872 real(RP),
intent(inout) :: rhoe(KA,IA,JA)
3873 real(RP),
intent(inout) :: rhoq(I_QV:I_NG,KA,IA,JA)
3874 real(RP),
intent(inout) :: q(KA,IA,JA,QA)
3875 real(RP),
intent(inout) :: tem(KA,IA,JA)
3876 real(RP),
intent(inout) :: pre(KA,IA,JA)
3877 real(RP),
intent(out) :: cva(KA,IA,JA)
3878 real(RP),
intent(in) :: esw(KA,IA,JA)
3879 real(RP),
intent(in) :: esi(KA,IA,JA)
3880 real(RP),
intent(in) :: rhoq2(I_QV:I_NG,KA,IA,JA)
3882 real(RP),
intent(inout) :: PQ(PQ_MAX,KA,IA,JA)
3883 real(RP),
intent(out) :: qc_evaporate(KA,IA,JA)
3885 real(RP),
intent(inout) :: sl_PLCdep(IA,JA)
3886 real(RP),
intent(inout) :: sl_PLRdep(IA,JA), sl_PNRdep(IA,JA)
3892 real(RP) :: wtem(KA,IA,JA)
3898 real(RP) :: qsw(KA,IA,JA)
3899 real(RP) :: qsi(KA,IA,JA)
3900 real(RP) :: dqswdtem_rho(KA,IA,JA)
3901 real(RP) :: dqsidtem_rho(KA,IA,JA)
3902 real(RP) :: dqswdtem_pre(KA,IA,JA)
3903 real(RP) :: dqsidtem_pre(KA,IA,JA)
3904 real(RP) :: dqswdpre_tem(KA,IA,JA)
3905 real(RP) :: dqsidpre_tem(KA,IA,JA)
3910 real(RP) :: aliqliq, asolliq
3911 real(RP) :: aliqsol, asolsol
3916 real(RP) :: taucnd, r_taucnd
3917 real(RP) :: taudep, r_taudep
3918 real(RP) :: taucnd_c(KA,IA,JA), r_taucnd_c
3919 real(RP) :: taucnd_r(KA,IA,JA), r_taucnd_r
3920 real(RP) :: taudep_i(KA,IA,JA), r_taudep_i
3921 real(RP) :: taudep_s(KA,IA,JA), r_taudep_s
3922 real(RP) :: taudep_g(KA,IA,JA), r_taudep_g
3925 real(RP) :: PLR2NR, PLI2NI, PLS2NS, PLG2NG
3926 real(RP) :: coef_a_cnd, coef_b_cnd
3927 real(RP) :: coef_a_dep, coef_b_dep
3949 real(RP) :: r_xc_ccn, r_xi_ccn
3952 real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
3953 real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
3955 real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
3956 real(RP) :: r_rvaptem
3958 real(RP) :: lvsw, lvsi
3959 real(RP) :: dlvsw, dlvsi
3961 real(RP) :: dcnd, ddep
3962 real(RP) :: uplim_cnd
3963 real(RP) :: lowlim_cnd
3965 real(RP) :: uplim_dep
3966 real(RP) :: lowlim_dep
3967 real(RP) :: ssw, ssi
3968 real(RP) :: r_esw, r_esi
3969 real(RP) :: r_lvsw, r_lvsi
3971 real(RP) :: ssw_o, ssi_o
3977 real(RP),
save :: fac_cndc = 1.0_rp
3978 logical,
save :: opt_fix_taucnd_c=.false.
3979 logical,
save :: flag_first =.true.
3981 namelist /nm_mp_sn14_condensation/ &
3982 opt_fix_taucnd_c, fac_cndc
3984 real(RP) :: fac_cndc_wrk
3986 real(RP),
parameter :: tau100day = 1.e+7_rp
3987 real(RP),
parameter :: r_tau100day = 1.e-7_rp
3988 real(RP),
parameter :: eps=1.e-30_rp
3990 integer :: i,j,k,iqw
3995 if( flag_first )
then 3996 flag_first = .false.
3998 read (
io_fid_conf,nml=nm_mp_sn14_condensation, end=100)
4007 r_xc_ccn=1.0_rp/xc_ccn
4010 if( opt_fix_taucnd_c )
then 4011 fac_cndc_wrk = fac_cndc**(1.0_rp-b_m(i_mp_qc))
4015 pq(i_lcdep,k,i,j) = pq(i_lcdep,k,i,j)*fac_cndc_wrk
4019 if(
io_l )
write(
io_fid_log,*)
"taucnd:fac_cndc_wrk=",fac_cndc_wrk
4028 wtem(k,i,j) = max( tem(k,i,j), tem_min )
4033 call moist_pres2qsat_liq ( qsw, wtem, pre )
4034 call moist_pres2qsat_ice ( qsi, wtem, pre )
4035 call moist_dqsw_dtem_rho ( dqswdtem_rho, wtem, rho )
4036 call moist_dqsi_dtem_rho ( dqsidtem_rho, wtem, rho )
4037 call moist_dqsw_dtem_dpre( dqswdtem_pre, dqswdpre_tem, wtem, pre )
4038 call moist_dqsi_dtem_dpre( dqsidtem_pre, dqsidpre_tem, wtem, pre )
4040 profile_start(
"sn14_update")
4044 if( z(k) <= 25000.0_rp )
then 4045 w = 0.5_rp*(velz(k-1,i,j) + velz(k,i,j))
4049 if( pre(k,i,j) < esw(k,i,j)+1.e-10_rp )
then 4051 dqswdtem_rho(k,i,j) = 0.0_rp
4052 dqswdtem_pre(k,i,j) = 0.0_rp
4053 dqswdpre_tem(k,i,j) = 0.0_rp
4055 if( pre(k,i,j) < esi(k,i,j)+1.e-10_rp )
then 4057 dqsidtem_rho(k,i,j) = 0.0_rp
4058 dqsidtem_pre(k,i,j) = 0.0_rp
4059 dqsidpre_tem(k,i,j) = 0.0_rp
4062 r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4063 lvsw = esw(k,i,j)*r_rvaptem
4064 lvsi = esi(k,i,j)*r_rvaptem
4065 pv = rhoq2(i_qv,k,i,j)*rvap*tem(k,i,j)
4066 r_esw = 1.0_rp/esw(k,i,j)
4067 r_esi = 1.0_rp/esi(k,i,j)
4068 ssw = min( mp_ssw_lim, ( pv*r_esw-1.0_rp ) )
4069 ssi = pv*r_esi - 1.0_rp
4070 r_lvsw = 1.0_rp/lvsw
4071 r_lvsi = 1.0_rp/lvsi
4072 r_taucnd_c = pq(i_lcdep,k,i,j)*r_lvsw
4073 r_taucnd_r = pq(i_lrdep,k,i,j)*r_lvsw
4074 r_taudep_i = pq(i_lidep,k,i,j)*r_lvsi
4075 r_taudep_s = pq(i_lsdep,k,i,j)*r_lvsi
4076 r_taudep_g = pq(i_lgdep,k,i,j)*r_lvsi
4084 calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry,
tracer_cv )
4085 calc_cp( cpa, qdry, q, k, i, j, iqw, cpdry,
tracer_cp )
4086 calc_r( rmoist, qdry, q, k, i, j, iqw, rdry,
tracer_r )
4087 r_cva = 1.0_rp / cva(k,i,j)
4088 r_cpa = 1.0_rp / cpa
4092 + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
4095 + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
4098 + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
4101 + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
4102 pdynliq = w * grav * ( r_cpa*dqswdtem_pre(k,i,j) + rho(k,i,j)*dqswdpre_tem(k,i,j) )
4103 pdynsol = w * grav * ( r_cpa*dqsidtem_pre(k,i,j) + rho(k,i,j)*dqsidpre_tem(k,i,j) )
4104 pradliq = -dtdt_rad(k,i,j) * dqswdtem_rho(k,i,j)
4105 pradsol = -dtdt_rad(k,i,j) * dqsidtem_rho(k,i,j)
4112 acnd = pdynliq + pradliq &
4113 - ( r_taudep_i+r_taudep_s+r_taudep_g ) * ( qsw(k,i,j) - qsi(k,i,j) )
4114 adep = pdynsol + pradsol &
4115 + ( r_taucnd_c+r_taucnd_r ) * ( qsw(k,i,j) - qsi(k,i,j) )
4117 + aliqliq*( r_taucnd_c+r_taucnd_r ) &
4118 + asolliq*( r_taudep_i+r_taudep_s+r_taudep_g )
4120 + aliqsol*( r_taucnd_c+r_taucnd_r )&
4121 + asolsol*( r_taudep_i+r_taudep_s+r_taudep_g )
4123 uplim_cnd = max( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
4124 lowlim_cnd = min( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
4125 if( r_taucnd < r_tau100day )
then 4127 pq(i_lcdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lcdep,k,i,j)*ssw_o ))
4128 pq(i_lrdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lrdep,k,i,j)*ssw_o ))
4129 pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_nrdep,k,i,j)*ssw_o )
4132 taucnd = 1.0_rp/r_taucnd
4134 coef_a_cnd = rho(k,i,j)*acnd*taucnd
4135 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 )
4136 pq(i_lcdep,k,i,j) = coef_a_cnd*r_taucnd_c - coef_b_cnd*r_taucnd_c
4137 plr2nr = pq(i_nrdep,k,i,j)/(pq(i_lrdep,k,i,j)+1.e-30_rp)
4138 pq(i_lrdep,k,i,j) = coef_a_cnd*r_taucnd_r - coef_b_cnd*r_taucnd_r
4139 pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_lrdep,k,i,j)*plr2nr )
4142 uplim_dep = max( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4143 lowlim_dep = min( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4144 if( r_taudep < r_tau100day )
then 4146 pq(i_lidep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lidep,k,i,j)*ssi_o ))
4147 pq(i_lsdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lsdep,k,i,j)*ssi_o ))
4148 pq(i_lgdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lgdep,k,i,j)*ssi_o ))
4149 pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_nidep,k,i,j)*ssi_o )
4150 pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_nsdep,k,i,j)*ssi_o )
4151 pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_ngdep,k,i,j)*ssi_o )
4153 taudep = 1.0_rp/r_taudep
4155 coef_a_dep = rho(k,i,j)*adep*taudep
4156 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 )
4157 pli2ni = pq(i_nidep,k,i,j)/max(pq(i_lidep,k,i,j),1.e-30_rp)
4158 pls2ns = pq(i_nsdep,k,i,j)/max(pq(i_lsdep,k,i,j),1.e-30_rp)
4159 plg2ng = pq(i_ngdep,k,i,j)/max(pq(i_lgdep,k,i,j),1.e-30_rp)
4160 pq(i_lidep,k,i,j) = coef_a_dep*r_taudep_i - coef_b_dep*r_taudep_i
4161 pq(i_lsdep,k,i,j) = coef_a_dep*r_taudep_s - coef_b_dep*r_taudep_s
4162 pq(i_lgdep,k,i,j) = coef_a_dep*r_taudep_g - coef_b_dep*r_taudep_g
4163 pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_lidep,k,i,j)*pli2ni )
4164 pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_lsdep,k,i,j)*pls2ns )
4165 pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_lgdep,k,i,j)*plg2ng )
4168 sw = 0.5_rp - sign(0.5_rp, pq(i_lcdep,k,i,j)+eps)
4169 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
4182 r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4183 lvsw = esw(k,i,j)*r_rvaptem
4184 dlvsw = rhoq2(i_qv,k,i,j)-lvsw
4185 dcnd = dt*(pq(i_lcdep,k,i,j)+pq(i_lrdep,k,i,j))
4187 sw = ( sign(0.5_rp,dcnd) + sign(0.5_rp,dlvsw) ) &
4188 * ( 0.5_rp + sign(0.5_rp,abs(dcnd)-eps) )
4192 fac1 = min(dlvsw*sw,dcnd*sw)*sw / (abs(sw)-1.0_rp+dcnd) &
4194 dep_dqc = max( dt*pq(i_lcdep,k,i,j)*fac1, &
4195 -rhoq2(
i_qc,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw)
4196 dep_dqr = max( dt*pq(i_lrdep,k,i,j)*fac1, &
4197 -rhoq2(
i_qr,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw)
4216 dep_dnc = max( dt*pncdep*fac1, -rhoq2(
i_nc,k,i,j) )
4217 dep_dnr = max( dt*pq(i_nrdep,k,i,j)*fac1, -rhoq2(
i_nr,k,i,j) )
4219 qc_evaporate(k,i,j) = - dep_dnc
4222 lvsi = esi(k,i,j)*r_rvaptem
4223 ddep = dt*(pq(i_lidep,k,i,j)+pq(i_lsdep,k,i,j)+pq(i_lgdep,k,i,j))
4224 dlvsi = rhoq2(i_qv,k,i,j)-lvsi
4226 sw = ( sign(0.5_rp,ddep) + sign(0.5_rp,dlvsi) ) &
4227 * ( 0.5_rp + sign(0.5_rp,abs(ddep)-eps) )
4231 fac2 = min(dlvsi*sw,ddep*sw)*sw / (abs(sw)-1.0_rp+ddep) &
4233 dep_dqi = max( dt*pq(i_lidep,k,i,j) &
4234 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
4235 -rhoq2(
i_qi,k,i,j) - 1e30_rp*(sw+1.0_rp) )
4236 dep_dqs = max( dt*pq(i_lsdep,k,i,j) &
4237 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
4238 -rhoq2(
i_qs,k,i,j) - 1e30_rp*(sw+1.0_rp) )
4239 dep_dqg = max( dt*pq(i_lgdep,k,i,j) &
4240 * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), &
4241 -rhoq2(
i_qg,k,i,j) - 1e30_rp*(sw+1.0_rp) )
4263 dep_dni = max( dt*pq(i_nidep,k,i,j)*fac2, -rhoq2(
i_ni,k,i,j) )
4264 dep_dns = max( dt*pq(i_nsdep,k,i,j)*fac2, -rhoq2(
i_ns,k,i,j) )
4265 dep_dng = max( dt*pq(i_ngdep,k,i,j)*fac2, -rhoq2(i_ng,k,i,j) )
4268 frz_dqc = max( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j)), -rhoq2(
i_qc,k,i,j)-dep_dqc )
4269 frz_dnc = max( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j)), -rhoq2(
i_nc,k,i,j)-dep_dnc )
4270 fac3 = ( frz_dqc-eps )/( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j))-eps )
4271 fac4 = ( frz_dnc-eps )/( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j))-eps )
4272 pq(i_lchom,k,i,j) = fac3*pq(i_lchom,k,i,j)
4273 pq(i_lchet,k,i,j) = fac3*pq(i_lchet,k,i,j)
4274 pq(i_nchom,k,i,j) = fac4*pq(i_nchom,k,i,j)
4275 pq(i_nchet,k,i,j) = fac4*pq(i_nchet,k,i,j)
4279 mlt_dqi = max( dt*pq(i_limlt,k,i,j), -rhoq2(
i_qi,k,i,j)-dep_dqi )
4280 mlt_dni = max( dt*pq(i_nimlt,k,i,j), -rhoq2(
i_ni,k,i,j)-dep_dni )
4282 mlt_dqs = max( dt*pq(i_lsmlt,k,i,j), -rhoq2(
i_qs,k,i,j)-dep_dqs )
4283 mlt_dns = max( dt*pq(i_nsmlt,k,i,j), -rhoq2(
i_ns,k,i,j)-dep_dns )
4285 mlt_dqg = max( dt*pq(i_lgmlt,k,i,j), -rhoq2(
i_qg,k,i,j)-dep_dqg )
4286 mlt_dng = max( dt*pq(i_ngmlt,k,i,j), -rhoq2(i_ng,k,i,j)-dep_dng )
4289 frz_dqr = max( dt*(pq(i_lrhet,k,i,j)), min(0.0_rp, -rhoq2(
i_qr,k,i,j)-dep_dqr) )
4290 frz_dnr = max( dt*(pq(i_nrhet,k,i,j)), min(0.0_rp, -rhoq2(
i_nr,k,i,j)-dep_dnr) )
4292 fac5 = ( frz_dqr-eps )/( dt*pq(i_lrhet,k,i,j)-eps )
4293 pq(i_lrhet,k,i,j) = fac5*pq(i_lrhet,k,i,j)
4294 fac6 = ( frz_dnr-eps )/( dt*pq(i_nrhet,k,i,j)-eps )
4295 pq(i_nrhet,k,i,j) = fac6*pq(i_nrhet,k,i,j)
4298 drhoqv = -(dep_dqc+dep_dqi+dep_dqs+dep_dqg+dep_dqr)
4300 rhoq(i_qv,k,i,j) = max(0.0_rp, rhoq(i_qv,k,i,j) + drhoqv )
4302 rhoe(k,i,j) = rhoe(k,i,j) - lhv * drhoqv
4304 xi = min(xi_max, max(xi_min, rhoq2(
i_qi,k,i,j)/(rhoq2(
i_ni,k,i,j)+ni_min) ))
4305 sw = 0.5_rp + sign(0.5_rp,xi-x_sep)
4309 drhoqc = ( frz_dqc - mlt_dqi*(1.0_rp-sw) + dep_dqc )
4310 drhonc = ( frz_dnc - mlt_dni*(1.0_rp-sw) + dep_dnc )
4312 drhoqr = ( frz_dqr - mlt_dqg - mlt_dqs - mlt_dqi*sw + dep_dqr )
4313 drhonr = ( frz_dnr - mlt_dng - mlt_dns - mlt_dni*sw + dep_dnr )
4315 rhoq(
i_qc,k,i,j) = max(0.0_rp, rhoq(
i_qc,k,i,j) + drhoqc )
4316 rhoq(
i_nc,k,i,j) = max(0.0_rp, rhoq(
i_nc,k,i,j) + drhonc )
4317 rhoq(
i_qr,k,i,j) = max(0.0_rp, rhoq(
i_qr,k,i,j) + drhoqr )
4318 rhoq(
i_nr,k,i,j) = max(0.0_rp, rhoq(
i_nr,k,i,j) + drhonr )
4321 drhoqi = (-frz_dqc + mlt_dqi + dep_dqi )
4322 drhoni = (-frz_dnc + mlt_dni + dep_dni )
4324 rhoq(
i_qi,k,i,j) = max(0.0_rp, rhoq(
i_qi,k,i,j) + drhoqi )
4325 rhoq(
i_ni,k,i,j) = max(0.0_rp, rhoq(
i_ni,k,i,j) + drhoni )
4327 rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqi
4330 drhoqs = ( mlt_dqs + dep_dqs )
4331 drhons = ( mlt_dns + dep_dns )
4333 rhoq(
i_qs,k,i,j) = max(0.0_rp, rhoq(
i_qs,k,i,j) + drhoqs )
4334 rhoq(
i_ns,k,i,j) = max(0.0_rp, rhoq(
i_ns,k,i,j) + drhons )
4336 rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqs
4339 drhoqg = (-frz_dqr + mlt_dqg + dep_dqg )
4340 drhong = (-frz_dnr + mlt_dng + dep_dng )
4342 rhoq(
i_qg,k,i,j) = max(0.0_rp, rhoq(
i_qg,k,i,j) + drhoqg )
4343 rhoq(i_ng,k,i,j) = max(0.0_rp, rhoq(i_ng,k,i,j) + drhong )
4345 rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqg
4348 rrho = 1.0_rp/rho(k,i,j)
4350 q(k,i,j,i_qv) = rhoq(i_qv,k,i,j) * rrho
4351 q(k,i,j,
i_qc) = rhoq(
i_qc,k,i,j) * rrho
4352 q(k,i,j,
i_qr) = rhoq(
i_qr,k,i,j) * rrho
4353 q(k,i,j,
i_qi) = rhoq(
i_qi,k,i,j) * rrho
4354 q(k,i,j,
i_qs) = rhoq(
i_qs,k,i,j) * rrho
4355 q(k,i,j,
i_qg) = rhoq(
i_qg,k,i,j) * rrho
4356 q(k,i,j,
i_nc) = rhoq(
i_nc,k,i,j) * rrho
4357 q(k,i,j,
i_nr) = rhoq(
i_nr,k,i,j) * rrho
4358 q(k,i,j,
i_ni) = rhoq(
i_ni,k,i,j) * rrho
4359 q(k,i,j,
i_ns) = rhoq(
i_ns,k,i,j) * rrho
4360 q(k,i,j,i_ng) = rhoq(i_ng,k,i,j) * rrho
4363 calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry,
tracer_cv )
4364 calc_r( rmoist, qdry, q, k, i, j, iqw, rdry,
tracer_r )
4365 tem(k,i,j) = rhoe(k,i,j) / ( rho(k,i,j) * cva(k,i,j) )
4366 pre(k,i,j) = rho(k,i,j) * rmoist * tem(k,i,j)
4368 sl_plcdep(i,j) = sl_plcdep(i,j) + dep_dqc*dz(k)*gsgam2(k,i,j)
4369 sl_plrdep(i,j) = sl_plrdep(i,j) + dep_dqr*dz(k)*gsgam2(k,i,j)
4370 sl_pnrdep(i,j) = sl_pnrdep(i,j) + dep_dnr*dz(k)*gsgam2(k,i,j)
4374 profile_stop(
"sn14_update")
4385 real(RP),
intent(inout) :: DENS(KA,IA,JA)
4386 real(RP),
intent(inout) :: QTRC(KA,IA,JA,QA)
4388 real(RP) :: diffq(KA,IA,JA)
4391 integer :: k, i, j, iq
4396 r_xmin = 1.0_rp / xmin_filter
4403 diffq(k,i,j) = qtrc(k,i,j,
i_qv) &
4404 + qtrc(k,i,j,
i_qc) &
4405 + qtrc(k,i,j,
i_qr) &
4406 + qtrc(k,i,j,
i_qi) &
4407 + qtrc(k,i,j,
i_qs) &
4414 qtrc(k,i,j,iq) = max(0.0_rp, qtrc(k,i,j,iq))
4421 dens(k,i,j) = dens(k,i,j) &
4423 + qtrc(k,i,j,
i_qv) &
4424 + qtrc(k,i,j,
i_qc) &
4425 + qtrc(k,i,j,
i_qr) &
4426 + qtrc(k,i,j,
i_qi) &
4427 + qtrc(k,i,j,
i_qs) &
4428 + qtrc(k,i,j,
i_qg) &
4436 if ( qtrc(k,i,j,
i_nc) > qtrc(k,i,j,
i_qc)*r_xmin )
then 4437 qtrc(k,i,j,
i_nc) = qtrc(k,i,j,
i_qc)*r_xmin
4441 if ( qtrc(k,i,j,
i_nr) > qtrc(k,i,j,
i_qr)*r_xmin )
then 4442 qtrc(k,i,j,
i_nr) = qtrc(k,i,j,
i_qr)*r_xmin
4446 if ( qtrc(k,i,j,
i_ni) > qtrc(k,i,j,
i_qi)*r_xmin )
then 4447 qtrc(k,i,j,
i_ni) = qtrc(k,i,j,
i_qi)*r_xmin
4451 if ( qtrc(k,i,j,
i_ns) > qtrc(k,i,j,
i_qs)*r_xmin )
then 4452 qtrc(k,i,j,
i_ns) = qtrc(k,i,j,
i_qs)*r_xmin
4456 if ( qtrc(k,i,j,
i_ng) > qtrc(k,i,j,
i_qg)*r_xmin )
then 4457 qtrc(k,i,j,
i_ng) = qtrc(k,i,j,
i_qg)*r_xmin
4481 real(RP),
intent(out) :: cldfrac(
ka,
ia,
ja)
4482 real(RP),
intent(in) :: qtrc (
ka,
ia,
ja,
qa)
4483 real(RP),
intent(in) :: mask_criterion
4486 integer :: k, i, j, iq
4494 qhydro = qhydro + qtrc(k,i,j,iq)
4496 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
4519 real(RP),
intent(in) :: qtrc0(
ka,
ia,
ja,
qa)
4520 real(RP),
intent(in) :: dens0(
ka,
ia,
ja)
4521 real(RP),
intent(in) :: temp0(
ka,
ia,
ja)
4530 real(RP) :: dc_ave(
ka,
ia,
ja)
4531 real(RP) :: dr_ave(
ka,
ia,
ja)
4539 real(RP) :: coef_fuetal1998
4541 real(RP),
parameter :: r2m_min=1.e-12_rp
4542 real(RP),
parameter :: um2cm = 100.0_rp
4544 real(RP) :: limitsw, zerosw
4552 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) ) )
4553 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) ) )
4554 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) ) )
4555 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) ) )
4556 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) ) )
4565 dc_ave(k,i,j) = a_m(i_mp_qc) * xc(k,i,j)**b_m(i_mp_qc)
4566 dr_ave(k,i,j) = a_m(i_mp_qr) * xr(k,i,j)**b_m(i_mp_qr)
4575 rc = 0.5_rp * dc_ave(k,i,j)
4576 limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
4577 re(k,i,j,i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
4586 rr = 0.5_rp * dr_ave(k,i,j)
4587 limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
4588 re(k,i,j,i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
4596 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)
4597 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)
4598 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)
4604 coef_fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
4608 ri3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k,i,j)
4609 rs3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k,i,j)
4610 rg3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k,i,j)
4619 zerosw = 0.5_rp - sign(0.5_rp, ri2m(k,i,j) - r2m_min )
4620 re(k,i,j,i_hi) = ri3m(k,i,j) / ( ri2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4629 zerosw = 0.5_rp - sign(0.5_rp, rs2m(k,i,j) - r2m_min )
4630 re(k,i,j,i_hs) = rs3m(k,i,j) / ( rs2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4639 zerosw = 0.5_rp - sign(0.5_rp, rg2m(k,i,j) - r2m_min )
4640 re(k,i,j,i_hg) = rg3m(k,i,j) / ( rg2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4645 re(:,:,:,i_hg+1:) = 0.0_rp
4662 real(RP),
intent(in) :: qtrc0(
ka,
ia,
ja,
qa)
4664 integer :: ihydro, iqa
4668 qe(:,:,:,i_hc) = qtrc0(:,:,:,i_qc)
4669 qe(:,:,:,i_hr) = qtrc0(:,:,:,i_qr)
4670 qe(:,:,:,i_hi) = qtrc0(:,:,:,i_qi)
4671 qe(:,:,:,i_hs) = qtrc0(:,:,:,i_qs)
4672 qe(:,:,:,i_hg) = qtrc0(:,:,:,i_qg)
4673 qe(:,:,:,i_hg+1:) = 0.0_rp
character(len=h_mid), dimension(qa_mp), target, public atmos_phy_mp_sn14_desc
integer, public is
start point of inner domain: x, local
subroutine nucleation_kij(z, velz, rho, tem, pre, rhoq, PQ, cpa, dTdt_rad, qke, CCN, dt)
subroutine, public atmos_phy_mp_sn14_effectiveradius(Re, QTRC0, DENS0, TEMP0)
Calculate Effective Radius.
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]
integer, public je
end point of inner domain: y, local
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
module ATMOSPHERE / Saturation adjustment
subroutine, public prc_mpistop
Abort MPI.
module ATMOSPHERE / Physics Cloud Microphysics
subroutine, public atmos_phy_mp_sn14_setup
Setup Cloud Microphysics.
real(rp), parameter, public const_ci
specific heat (ice) [J/kg/K]
real(rp), parameter, public const_dwatr
density of water [kg/m3]
real(rp), dimension(qa_max), public tracer_r
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
logical, public io_l
output log or not? (this process)
integer, parameter, public i_hs
snow
real(rp) function, public sf_gamma(x)
Gamma function.
real(rp), dimension(:), allocatable, public grid_cz
center coordinate [m]: z, local=global
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
subroutine freezing_water_kij(dt, PQ, rhoq, xq, tem)
integer, parameter, public i_hr
liquid water rain
integer, parameter, public i_hi
ice water cloud
integer, public ke
end point of inner domain: z, local
subroutine ice_multiplication_kij(PQ, Pac, tem, rhoq, xq)
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
subroutine dep_vapor_melt_ice_kij(PQ, rho, tem, pre, qd, rhoq, esw, esi, xq, vt_xave, dq_xave)
real(rp), public const_cvvap
specific heat (water vapor, constant volume) [J/kg/K]
real(rp), dimension(n_hyd), target, public atmos_phy_mp_sn14_dens
real(rp), public const_lhf0
latent heat of fusion at 0C [J/kg]
real(rp), parameter, public const_dice
density of ice [kg/m3]
real(rp), dimension(qa_max), public tracer_cv
subroutine, public atmos_saturation_dqsw_dtem_rho(dqsdtem, temp, dens)
subroutine, public atmos_saturation_dqsw_dtem_dpre(dqsdtem, dqsdpre, temp, pres)
real(rp), public dz
length in the main region [m]: z
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
real(rp), public const_undef
subroutine mixed_phase_collection_kij(Pac, PQ, wtem, rhoq, xq, dq_xave, vt_xave)
real(rp), parameter, public const_lhs0
latent heat of sublimation at 0C [J/kg]
subroutine, public atmos_phy_mp_sn14_cloudfraction(cldfrac, QTRC, mask_criterion)
Calculate Cloud Fraction.
module ATMOSPHERE / Physics Cloud Microphysics - Common
subroutine, public atmos_saturation_dqsi_dtem_rho(dqsdtem, temp, dens)
real(rp), dimension(qa_max), public tracer_cp
logical, public io_nml
output log or not? (for namelist, this process)
subroutine debug_tem_kij(point, tem, rho, pre, qv)
integer, public ia
of whole cells: x, local, with HALO
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_sn14_name
integer, public ka
of whole cells: z, local, with HALO
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 qa_mp
real(rp), public const_lhf00
latent heat of fusion at 0K [J/kg]
real(rp), public const_grav
standard acceleration of gravity [m/s2]
integer, public js
start point of inner domain: y, local
subroutine update_by_phase_change_kij(ntdiv, ntmax, dt, gsgam2, z, dz, velz, dTdt_rad, rho, rhoe, rhoq, q, tem, pre, cva, esw, esi, rhoq2, PQ, qc_evaporate, sl_PLCdep, sl_PLRdep, sl_PNRdep)
real(rp), public const_lhv00
latent heat of vaporizaion at 0K [J/kg]
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
subroutine, public atmos_saturation_dqsi_dtem_dpre(dqsdtem, dqsdpre, temp, pres)
integer, parameter, public i_hc
liquid water cloud
integer, public ks
start point of inner domain: z, local
real(rp), parameter, public const_emelt
integer, public prc_myrank
process num in local communicator
subroutine, public atmos_phy_mp_sn14_config(MP_TYPE, QA, QS)
Configure.
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
real(dp), parameter, public const_undef8
undefined value (REAL8)
subroutine, public atmos_phy_mp_sn14_mixingratio(Qe, QTRC0)
Calculate mixing ratio of each category.
subroutine, public atmos_phy_mp_sn14(DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
Cloud Microphysics.
subroutine aut_acc_slc_brk_kij(PQ, rhoq, xq, dq_xave, rho)
integer, public ie
end point of inner domain: x, local
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_sn14_unit
real(rp), public const_eps
small number
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
subroutine, public atmos_hydrometeor_regist(Q0, NV, NL, NI, NAME, DESC, UNIT, ADVC)
Regist tracer.
real(rp), public const_pi
pi
subroutine, public atmos_phy_mp_precipitation(QA_MP, QS_MP, qflx, vterm, DENS, MOMZ, MOMX, MOMY, RHOE, QTRC, temp, CVq, dt, vt_fixed)
precipitation transport
integer, public io_fid_conf
Config file ID.
subroutine, public tracer_regist(QS, NQ, NAME, DESC, UNIT, CV, CP, R, ADVC, MASS)
Regist tracer.
integer, public io_fid_log
Log file ID.
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.
integer, parameter, public rp
integer, parameter, public i_hg
graupel
integer, public io_fid_nml
Log file ID (only for output namelist)
real(rp), dimension(qa_max), public tracer_mass
subroutine mp_negativefilter(DENS, QTRC)
integer, public ja
of whole cells: y, local, with HALO