19 #include "inc_openmp.h" 55 integer,
public,
parameter ::
qa_mp = 3
69 'Ratio of Water Vapor mass to total mass (Specific humidity)', &
70 'Ratio of Cloud Water mass to total mass', &
71 'Ratio of Rain Water mass to total mass' /
83 private :: mp_kessler_vterm
89 integer,
private,
parameter :: i_mp_qc = 1
90 integer,
private,
parameter :: i_mp_qr = 2
91 integer,
private :: qs_mp
92 integer,
private :: qe_mp
94 logical,
private :: mp_donegative_fixer = .true.
95 logical,
private :: mp_doprecipitation = .true.
96 logical,
private :: mp_couple_aerosol = .false.
97 real(RP),
private :: mp_limit_negative = 1.0_rp
99 real(RP),
private,
allocatable :: factor_vterm(:)
101 real(RP),
private,
parameter :: re_qc = 8.e-6_rp
103 integer,
private :: mp_ntmax_sedimentation = 1
104 integer,
private :: mp_nstep_sedimentation
105 real(RP),
private :: mp_rnstep_sedimentation
106 real(DP),
private :: mp_dtsec_sedimentation
108 logical,
private :: first = .true.
123 character(len=*),
intent(in) :: mp_type
124 integer,
intent(out) :: qa
125 integer,
intent(out) :: qs
129 if(
io_l )
write(
io_fid_log,*)
'++++++ Module[Cloud Microphysics Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]' 130 if(
io_l )
write(
io_fid_log,*)
'*** Tracers for KESSLER-type 1-moment bulk 3 category' 132 if ( mp_type /=
'KESSLER' )
then 133 write(*,*)
'xxx ATMOS_PHY_MP_TYPE is not KESSLER. Check!' 145 qe_mp = qs_mp +
qa_mp - 1
168 real(RP),
parameter :: max_term_vel = 10.0_rp
171 namelist / param_atmos_phy_mp / &
172 mp_doprecipitation, &
173 mp_donegative_fixer, &
175 mp_ntmax_sedimentation, &
183 if(
io_l )
write(
io_fid_log,*)
'++++++ Module[Cloud Microphysics] / Categ[ATMOS PHYSICS] / Origin[SCALElib]' 184 if(
io_l )
write(
io_fid_log,*)
'*** KESSLER-type 1-moment bulk 3 category' 188 read(
io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
190 if(
io_l )
write(
io_fid_log,*)
'*** Not found namelist. Default used.' 191 elseif( ierr > 0 )
then 192 write(*,*)
'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!' 197 if( mp_couple_aerosol )
then 198 write(*,*)
'xxx MP_aerosol_couple should be .false. for KESSLER type MP!' 203 if(
io_l )
write(
io_fid_log,*)
'*** Enable negative fixer? : ', mp_donegative_fixer
204 if(
io_l )
write(
io_fid_log,*)
'*** Value limit of negative fixer (abs) : ', abs(mp_limit_negative)
205 if(
io_l )
write(
io_fid_log,*)
'*** Enable sedimentation (precipitation)? : ', mp_doprecipitation
208 mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
210 mp_nstep_sedimentation = mp_ntmax_sedimentation
211 mp_rnstep_sedimentation = 1.0_rp /
real(mp_ntmax_sedimentation,kind=
rp)
214 if(
io_l )
write(
io_fid_log,*)
'*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation,
'step' 215 if(
io_l )
write(
io_fid_log,*)
'*** DT of sedimentation : ', mp_dtsec_sedimentation,
'[s]' 221 allocate( factor_vterm(
ka) )
255 thermodyn_rhoe => atmos_thermodyn_rhoe, &
256 thermodyn_rhot => atmos_thermodyn_rhot, &
257 thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
264 real(RP),
intent(inout) :: dens (
ka,
ia,
ja)
265 real(RP),
intent(inout) :: momz (
ka,
ia,
ja)
266 real(RP),
intent(inout) :: momx (
ka,
ia,
ja)
267 real(RP),
intent(inout) :: momy (
ka,
ia,
ja)
268 real(RP),
intent(inout) :: rhot (
ka,
ia,
ja)
269 real(RP),
intent(inout) :: qtrc (
ka,
ia,
ja,
qa)
270 real(RP),
intent(in) :: ccn (
ka,
ia,
ja)
271 real(RP),
intent(out) :: evaporate(
ka,
ia,
ja)
272 real(RP),
intent(out) :: sflx_rain(
ia,
ja)
273 real(RP),
intent(out) :: sflx_snow(
ia,
ja)
275 real(RP) :: rhoe_t (
ka,
ia,
ja)
277 real(RP) :: rhoe (
ka,
ia,
ja)
278 real(RP) :: temp (
ka,
ia,
ja)
279 real(RP) :: pres (
ka,
ia,
ja)
284 real(RP) :: qc_t_sat (
ka,
ia,
ja)
285 real(RP) :: rho_prof (
ka)
291 if(
io_l )
write(
io_fid_log,*)
'*** Atmos physics step: Cloud microphysics(kessler)' 296 rho_prof(:) = rho_prof(:) * 1.e-3_rp
299 factor_vterm(k) = sqrt( rho_prof(
ks)/rho_prof(k) )
304 if ( mp_donegative_fixer )
then 305 call mp_negative_fixer( dens(:,:,:), &
312 call thermodyn_rhoe( rhoe(:,:,:), &
320 call mp_kessler( rhoe_t(:,:,:), &
326 vterm(:,:,:,:) = 0.0_rp
327 flx_hydro(:,:,:,:) = 0.0_rp
329 if ( mp_doprecipitation )
then 331 do step = 1, mp_nstep_sedimentation
333 call thermodyn_temp_pres_e( temp(:,:,:), &
342 call mp_kessler_vterm( vterm(:,:,:,:), &
346 call mp_precipitation(
qa_mp, &
358 mp_dtsec_sedimentation )
363 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
372 call hist_in( flx_hydro(:,:,:,i_mp_qr),
'pflux_QR',
'precipitation flux of QR',
'kg/m2/s', nohalo=.true. )
374 call hist_in( vterm(:,:,:,i_mp_qr),
'Vterm_QR',
'terminal velocity of QR',
'm/s' )
379 qc_t_sat(k,i,j) = qtrc(k,i,j,
i_qc)
384 call mp_saturation_adjustment( rhoe_t(:,:,:), &
392 flag_liquid = .true. )
397 qc_t_sat(k,i,j) = ( qtrc(k,i,j,
i_qc) - qc_t_sat(k,i,j) ) / dt_mp
402 call hist_in( qc_t_sat(:,:,:),
'Pcsat',
'QC production term by satadjust',
'kg/kg/s' )
407 evaporate(k,i,j) = max( -qc_t_sat(k,i,j), 0.0_rp )
408 evaporate(k,i,j) = evaporate(k,i,j) * dens(k,i,j) / (4.0_rp/3.0_rp*pi*dwatr*re_qc**3)
415 call thermodyn_rhot( rhot(:,:,:), &
422 if ( mp_donegative_fixer )
then 423 call mp_negative_fixer( dens(:,:,:), &
432 sflx_rain(i,j) = - flx_hydro(
ks-1,i,j,i_mp_qr)
433 sflx_snow(i,j) = 0.0_rp
442 subroutine mp_kessler( &
456 thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
460 saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq
463 real(RP),
intent(out) :: rhoe_t(
ka,
ia,
ja)
464 real(RP),
intent(out) :: qtrc_t(
ka,
ia,
ja,
qa)
465 real(RP),
intent(inout) :: rhoe0 (
ka,
ia,
ja)
466 real(RP),
intent(inout) :: qtrc0 (
ka,
ia,
ja,
qa)
467 real(RP),
intent(in) :: dens0 (
ka,
ia,
ja)
470 real(RP) :: qsat (
ka,
ia,
ja)
471 real(RP) :: temp0(
ka,
ia,
ja)
472 real(RP) :: pres0(
ka,
ia,
ja)
474 real(RP) :: dens (
ka)
475 real(RP) :: rhoe (
ka)
476 real(RP) :: temp (
ka)
477 real(RP) :: pres (
ka)
481 real(RP) :: qsatl(
ka)
482 real(RP) :: sliq (
ka)
485 real(RP) :: dq_evap(
ka)
486 real(RP) :: dq_auto(
ka)
487 real(RP) :: dq_accr(
ka)
489 real(RP) :: dqv, dqc, dqr
490 real(RP) :: vent_factor
500 call thermodyn_temp_pres_e( temp0(:,:,:), &
509 call saturation_dens2qsat_liq( qsat(:,:,:), &
530 dq_auto(k) = 1.e-3_rp * max( qc(k)-1.e-3_rp, 0.0_rp )
535 dq_accr(k) = 2.2_rp * qc(k) * qr(k)**0.875_rp
540 vent_factor = 1.6_rp + 124.9_rp * ( dens(k)*qr(k) )**0.2046_rp
542 dq_evap(k) = ( 1.0_rp-min(sliq(k),1.0_rp) ) / dens(k) * vent_factor &
543 * ( dens(k)*qr(k) )**0.525_rp / ( 5.4e5_rp + 2.55e8_rp / ( pres(k)*qsatl(k) ) )
548 dqc = (-dq_auto(k)-dq_accr(k) )
550 qtrc_t(k,i,j,i_qc) = max( dqc, -qc(k)*rdt )
554 dqr = ( dq_auto(k)+dq_accr(k)-dq_evap(k) )
556 qtrc_t(k,i,j,i_qr) = max( dqr, -qr(k)*rdt )
560 dqv = - ( qtrc_t(k,i,j,i_qc) &
561 + qtrc_t(k,i,j,i_qr) )
563 qtrc_t(k,i,j,i_qv) = max( dqv, -qv(k)*rdt )
567 rhoe_t(k,i,j) = - dens(k) * (
lhv * qtrc_t(k,i,j,i_qv) )
576 qtrc0(
ks:
ke,i,j,i_qv) = qtrc0(
ks:
ke,i,j,i_qv) + qtrc_t(
ks:
ke,i,j,i_qv) * dt
577 qtrc0(
ks:
ke,i,j,i_qc) = qtrc0(
ks:
ke,i,j,i_qc) + qtrc_t(
ks:
ke,i,j,i_qc) * dt
578 qtrc0(
ks:
ke,i,j,i_qr) = qtrc0(
ks:
ke,i,j,i_qr) + qtrc_t(
ks:
ke,i,j,i_qr) * dt
580 rhoe0(
ks:
ke,i,j) = rhoe0(
ks:
ke,i,j) + rhoe_t(
ks:
ke,i,j) * dt
587 end subroutine mp_kessler
591 subroutine mp_kessler_vterm( &
602 real(RP),
intent(in) :: dens0(
ka,
ia,
ja)
603 real(RP),
intent(in) :: qtrc0(
ka,
ia,
ja,
qa)
614 vterm(k,i,j,i_mp_qc) = 0.0_rp
622 zerosw = 0.5_rp - sign(0.5_rp, qtrc0(k,i,j,
i_qr) - 1.e-12_rp )
624 vterm(k,i,j,i_mp_qr) = - 36.34_rp * ( dens0(k,i,j) * ( qtrc0(k,i,j,
i_qr) + zerosw ) )**0.1364_rp &
625 * refstate_dens(
ks,i,j) / refstate_dens(k,i,j) * ( 1.0_rp - zerosw )
631 end subroutine mp_kessler_vterm
644 real(RP),
intent(out) :: cldfrac(
ka,
ia,
ja)
645 real(RP),
intent(in) :: qtrc (
ka,
ia,
ja,
qa)
646 real(RP),
intent(in) :: mask_criterion
649 integer :: k, i, j, iq
656 do iq = qs_mp+1, qe_mp
657 qhydro = qhydro + qtrc(k,i,j,iq)
659 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
681 real(RP),
intent(in) :: qtrc0(
ka,
ia,
ja,
qa)
682 real(RP),
intent(in) :: dens0(
ka,
ia,
ja)
683 real(RP),
intent(in) :: temp0(
ka,
ia,
ja)
685 real(RP),
parameter :: um2cm = 100.0_rp
691 if ( ih == i_hc )
then 692 re(:,:,:,ih) = 8.e-6_rp * um2cm
693 else if ( ih == i_hr )
then 694 re(:,:,:,ih) = 100.e-6_rp * um2cm
696 re(:,:,:,ih) = 0.0_rp
715 real(RP),
intent(in) :: qtrc0(
ka,
ia,
ja,
qa)
721 if ( ih == i_hc )
then 722 qe(:,:,:,ih) = qtrc0(:,:,:,i_qc)
723 else if ( ih == i_hr )
then 724 qe(:,:,:,ih) = qtrc0(:,:,:,i_qr)
726 qe(:,:,:,ih) = 0.0_rp
integer, public is
start point of inner domain: x, local
subroutine, public atmos_phy_mp_kessler(DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
Cloud Microphysics.
integer, public je
end point of inner domain: y, local
module ATMOSPHERE / Saturation adjustment
subroutine, public prc_mpistop
Abort MPI.
subroutine, public atmos_phy_mp_kessler_mixingratio(Qe, QTRC0)
Calculate mixing ratio of each category.
real(rp), parameter, public const_dwatr
density of water [kg/m3]
subroutine, public atmos_phy_mp_negative_fixer(DENS, RHOT, QTRC, I_QV, limit_negative)
Negative fixer.
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)
module ATMOSPHERE / Reference state
integer, parameter, public i_hr
liquid water rain
subroutine, public atmos_phy_mp_kessler_cloudfraction(cldfrac, QTRC, mask_criterion)
Calculate Cloud Fraction.
integer, public ke
end point of inner domain: z, local
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_kessler_name
real(rp), dimension(qa_max), public tracer_cv
subroutine, public atmos_phy_mp_kessler_effectiveradius(Re, QTRC0, DENS0, TEMP0)
Calculate Effective Radius.
real(rp), public const_undef
subroutine, public atmos_phy_mp_kessler_config(MP_TYPE, QA, QS)
Setup.
module ATMOSPHERE / Physics Cloud Microphysics - Common
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_kessler_unit
logical, public io_nml
output log or not? (for namelist, this process)
integer, public ia
of whole cells: x, local, with HALO
integer, parameter, public qa_mp
subroutine, public comm_horizontal_mean(varmean, var)
calculate horizontal mean (global total with communication)
module ATMOSPHERE / Physics Cloud Microphysics
real(rp), dimension(n_hyd), target, public atmos_phy_mp_kessler_dens
subroutine, public atmos_phy_mp_saturation_adjustment(RHOE_t, QTRC_t, RHOE0, QTRC0, DENS0, I_QV, I_QC, I_QI, flag_liquid)
Saturation adjustment.
integer, public ka
of whole cells: z, local, with HALO
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_dens
refernce density [kg/m3]
subroutine, public atmos_phy_mp_kessler_setup
Setup.
integer, public js
start point of inner domain: y, local
real(rp), public lhv
latent heat of vaporization for use [J/kg]
integer, parameter, public i_hc
liquid water cloud
integer, public ks
start point of inner domain: z, local
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
integer, public ie
end point of inner domain: x, local
module ATMOSPHERE / Thermodynamics
character(len=h_mid), dimension(qa_mp), target, public atmos_phy_mp_kessler_desc
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.
integer, public io_fid_log
Log file ID.
integer, parameter, public n_hyd
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
integer, parameter, public rp
integer, public io_fid_nml
Log file ID (only for output namelist)
real(rp), dimension(qa_max), public tracer_mass
integer, public ja
of whole cells: y, local, with HALO