68 private :: mp_suzuki10
76 private :: ice_nucleat
84 private :: random_setup
113 integer,
parameter :: i_qv = 1
115 character(len=3) :: namspc (8) = (/
'Qcl', &
124 character(len=27) :: lnamspc(8) = (/
'Mixing ratio of cloud bin', &
125 'Mixing ratio of colum bin', &
126 'Mixing ratio of plate bin', &
127 'Mixing ratio of dendrit bin', &
128 'Mixing ratio of snow bin', &
129 'Mixing ratio of graupel bin', &
130 'Mixing ratio of hail bin', &
131 'Mixing ratio of aerosol bin' /)
136 integer :: kphase = 0
137 integer :: iceflg = 1
139 integer,
parameter :: i_mp_qc = 1
140 integer,
parameter :: i_mp_qcl = 2
141 integer,
parameter :: i_mp_qp = 3
142 integer,
parameter :: i_mp_qd = 4
143 integer,
parameter :: i_mp_qs = 5
144 integer,
parameter :: i_mp_qg = 6
145 integer,
parameter :: i_mp_qh = 7
147 logical :: mp_doautoconversion = .true.
148 logical :: mp_couple_aerosol = .false.
151 integer :: num_start_waters
152 integer :: num_end_waters
153 integer :: num_start_ices
154 integer :: num_end_ices
157 integer,
parameter :: il = 1
158 integer,
parameter :: ic = 2
159 integer,
parameter :: ip = 3
160 integer,
parameter :: id = 4
161 integer,
parameter :: iss = 5
162 integer,
parameter :: ig = 6
163 integer,
parameter :: ih = 7
167 real(RP),
allocatable :: xctr( : )
168 real(RP),
allocatable :: xbnd( : )
169 real(RP),
allocatable :: radc( : )
170 real(RP),
allocatable :: cctr( :,: )
171 real(RP),
allocatable :: cbnd( :,: )
172 real(RP),
allocatable :: ck( :,:,:,: )
173 real(RP),
allocatable :: vt( :,: )
174 real(RP),
allocatable :: br( :,: )
175 integer,
allocatable :: ifrsl( :,:,: )
177 real(RP),
allocatable :: xactr( : )
178 real(RP),
allocatable :: xabnd( : )
179 real(RP),
allocatable :: rada( : )
181 real(RP),
allocatable :: expxctr( : )
182 real(RP),
allocatable :: expxbnd( : )
183 real(RP),
allocatable :: expxactr( : )
184 real(RP),
allocatable :: expxabnd( : )
185 real(RP),
allocatable :: rexpxctr( : )
186 real(RP),
allocatable :: rexpxbnd( : )
187 real(RP),
allocatable :: rexpxactr( : )
188 real(RP),
allocatable :: rexpxabnd( : )
194 real(RP),
allocatable,
save :: vterm(:)
197 real(RP),
parameter :: cldmin = 1.0e-10_rp
198 real(RP),
parameter :: oneovthird = 1.0_rp / 3.0_rp
199 real(RP),
parameter :: thirdovforth = 3.0_rp / 4.0_rp
200 real(RP),
parameter :: twoovthird = 2.0_rp / 3.0_rp
202 real(RP) :: rbnd = 40.e-06_rp
206 real(RP) :: rhoa = 2.25e+03_rp
207 real(RP) :: emaer = 58.0_rp
208 real(RP) :: emwtr = 18.0_rp
209 real(RP) :: rasta = 1.e-08_rp
210 real(RP) :: raend = 1.e-06_rp
211 real(RP) :: r0a = 1.e-07_rp
213 logical :: flg_regeneration = .false.
214 logical :: flg_nucl = .false.
215 logical :: flg_icenucl = .false.
216 logical :: flg_sf_aero = .false.
218 integer,
save :: rndm_flgp = 0
220 real(RP),
allocatable :: marate( : )
221 integer,
allocatable,
save :: ncld( : )
227 character(len=11),
parameter :: fname_micpara=
"micpara.dat" 229 integer(4) :: fid_micpara
232 integer,
allocatable :: blrg( :,: )
233 integer,
allocatable :: bsml( :,: )
235 integer :: mspc, mbin
236 real(RP) :: rndm(1,1,1)
239 real(RP) :: c_ccn = 100.e+6_rp
240 real(RP) :: kappa = 0.462_rp
243 real(RP) :: sigma = 7.5e-02_rp
244 real(RP) :: vhfct = 2.0_rp
246 real(RP),
parameter :: tcrit = 271.15_rp
247 integer,
allocatable :: kindx( :,: )
250 integer,
parameter :: ndat = 33, icemax = 3
251 integer,
parameter :: kdeg = 4, ldeg = 4, nspc_mk = 7
255 real(DP),
allocatable :: radc_mk( : ), xctr_mk( : ), xbnd_mk( : )
256 real(DP),
allocatable :: cctr_mk( :,: ), cbnd_mk( :,: )
257 real(DP),
allocatable :: ck_mk( :,:,:,: )
258 real(DP),
allocatable :: vt_mk( :,: )
259 real(DP),
allocatable :: br_mk( :,: )
261 real(DP) :: xmss( nspc_mk,ndat ), zcap( nspc_mk,ndat ), vtrm( nspc_mk,ndat )
262 real(DP) :: blkr( nspc_mk,ndat ), blkd( nspc_mk,ndat ), ykrn( nspc_mk,nspc_mk,ndat,ndat )
264 real(DP) :: ywll( ndat,ndat ), ywli( ndat,ndat,icemax ), ywls( ndat,ndat )
265 real(DP) :: ywlg( ndat,ndat ), ywlh( ndat,ndat )
267 real(DP) :: ywil( ndat,ndat,icemax ), ywii( ndat,ndat,icemax,icemax )
268 real(DP) :: ywis( ndat,ndat,icemax ), ywig( ndat,ndat,icemax )
269 real(DP) :: ywih( ndat,ndat,icemax )
271 real(DP) :: ywsl( ndat,ndat ), ywsi( ndat,ndat,icemax ), ywss( ndat,ndat )
272 real(DP) :: ywsg( ndat,ndat ), ywsh( ndat,ndat )
274 real(DP) :: ywgl( ndat,ndat ), ywgi( ndat,ndat,icemax ), ywgs( ndat,ndat )
275 real(DP) :: ywgg( ndat,ndat ), ywgh( ndat,ndat )
277 real(DP) :: ywhl( ndat,ndat ), ywhi( ndat,ndat,icemax ), ywhs( ndat,ndat )
278 real(DP) :: ywhg( ndat,ndat ), ywhh( ndat,ndat )
281 real(RP) :: sigma_sdf(5), r0_sdf(5), n0_sdf(5), rho_sdf(5)
291 namelist / param_atmos_phy_mp_suzuki10_bin / &
297 integer :: m, n, ierr
301 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Setup' 302 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Tracers setup for Suzuki (2010) Spectral BIN model' 305 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'READ BIN NUMBER' 308 read(
io_fid_conf,nml=param_atmos_phy_mp_suzuki10_bin,iostat=ierr)
311 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Not found namelist. Default used.' 312 elseif( ierr > 0 )
then 313 log_error(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10_bin, Check!' 317 log_nml(param_atmos_phy_mp_suzuki10_bin)
319 if( iceflg == 0 )
then 321 elseif( iceflg == 1 )
then 324 log_error(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
"ICEFLG should be 0 (warm rain) or 1 (mixed rain) check!!" 333 num_start_waters = i_qv + 1
335 num_start_ices = num_end_waters + 1
396 integer,
intent(in) :: KA
397 integer,
intent(in) :: IA
398 integer,
intent(in) :: JA
404 real(RP) :: S10_EMAER
406 logical :: S10_FLAG_REGENE
407 logical :: S10_FLAG_NUCLEAT
408 logical :: S10_FLAG_ICENUCLEAT
409 logical :: S10_FLAG_SFAERO
410 integer :: S10_RNDM_FLGP
411 integer :: S10_RNDM_MSPC
412 integer :: S10_RNDM_MBIN
414 namelist / param_atmos_phy_mp_suzuki10 / &
415 mp_doautoconversion, &
424 s10_flag_icenucleat, &
432 real(RP),
parameter :: max_term_vel = 10.0_rp
434 integer :: nnspc, nnbin
435 integer :: nn, mm, mmyu, nnyu
436 integer :: myu, nyu, i, j, k, n, ierr
440 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Setup' 441 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Suzuki (2010) Spectral BIN model' 444 allocate( xctr(
nbin ) )
445 allocate( xbnd(
nbin+1 ) )
446 allocate( radc(
nbin ) )
447 allocate( cctr(
nbin,nspc_mk ) )
448 allocate( cbnd(
nbin+1,nspc_mk ) )
449 allocate( ck( nspc_mk,nspc_mk,
nbin,
nbin ) )
450 allocate( vt( nspc_mk,
nbin ) )
451 allocate( br( nspc_mk,
nbin ) )
452 allocate( ifrsl( 2,nspc_mk,nspc_mk ) )
453 allocate( expxctr(
nbin ) )
454 allocate( expxbnd(
nbin+1 ) )
455 allocate( rexpxctr(
nbin ) )
456 allocate( rexpxbnd(
nbin+1 ) )
457 if (
nccn /= 0 )
then 458 allocate( xactr(
nccn ) )
459 allocate( xabnd(
nccn+1 ) )
460 allocate( rada(
nccn ) )
461 allocate( expxactr(
nccn ) )
462 allocate( expxabnd(
nccn+1 ) )
463 allocate( rexpxactr(
nccn ) )
464 allocate( rexpxabnd(
nccn+1 ) )
468 mspc = nspc_mk*nspc_mk
475 s10_flag_regene = flg_regeneration
476 s10_flag_nucleat = flg_nucl
477 s10_flag_icenucleat = flg_icenucl
478 s10_flag_sfaero = flg_sf_aero
479 s10_rndm_flgp = rndm_flgp
484 read(
io_fid_conf,nml=param_atmos_phy_mp_suzuki10,iostat=ierr)
487 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Not found namelist. Default used.' 488 elseif( ierr > 0 )
then 489 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10, Check!' 492 log_nml(param_atmos_phy_mp_suzuki10)
494 if ( nspc /= 1 .AND. nspc /= 7 )
then 495 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nspc should be set as 1 (warm rain) or 7 (mixed phase) check!' 504 flg_regeneration = s10_flag_regene
505 flg_nucl = s10_flag_nucleat
506 flg_icenucl = s10_flag_icenucleat
507 flg_sf_aero = s10_flag_sfaero
508 rndm_flgp = s10_rndm_flgp
517 open ( fid_micpara, file = fname_micpara, form =
'formatted', status =
'old', iostat=ierr )
520 if ( ierr == 0 )
then 522 read( fid_micpara,* ) nnspc, nnbin
524 if ( nnbin /=
nbin )
then 525 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nbin in namelist and nbin in micpara.dat is different check!' 530 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Radius of cloud *' 532 read( fid_micpara,* ) nn, xctr( n ), radc( n )
533 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
534 "Radius of ", n,
"th cloud bin (bin center)= ", radc( n ) ,
"[m]" 537 read( fid_micpara,* ) nn, xbnd( n )
539 read( fid_micpara,* ) dxmic
540 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Width of Cloud SDF= ', dxmic
546 read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
550 read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
559 read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
568 read( fid_micpara,* ) mmyu, nn, vt( myu,n )
575 read( fid_micpara,* ) mmyu, nn, br( myu,n )
579 close ( fid_micpara )
584 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'micpara.dat is created' 589 open ( fid_micpara, file = fname_micpara, form =
'formatted', status =
'old', iostat=ierr )
591 read( fid_micpara,* ) nnspc, nnbin
593 if ( nnbin /=
nbin )
then 594 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nbin in inc_tracer and nbin in micpara.dat is different check!' 599 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Radius of cloud *' 601 read( fid_micpara,* ) nn, xctr( n ), radc( n )
602 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
603 "Radius of ", n,
"th cloud bin (bin center)= ", radc( n ) ,
"[m]" 606 read( fid_micpara,* ) nn, xbnd( n )
608 read( fid_micpara,* ) dxmic
609 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Width of Cloud SDF= ', dxmic
615 read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
619 read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
628 read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
637 read( fid_micpara,* ) mmyu, nn, vt( myu,n )
644 read( fid_micpara,* ) mmyu, nn, br( myu,n )
648 close ( fid_micpara )
665 if (
nccn /= 0 )
then 667 allocate ( ncld( 1:
nccn ) )
668 xasta = log( rhoa*4.0_rp/3.0_rp*pi * ( rasta )**3 )
669 xaend = log( rhoa*4.0_rp/3.0_rp*pi * ( raend )**3 )
671 dxaer = ( xaend-xasta )/
nccn 674 xabnd( n ) = xasta + dxaer*( n-1 )
677 xactr( n ) = ( xabnd( n )+xabnd( n+1 ) )*0.50_rp
678 rada( n ) = ( exp( xactr( n ) )*thirdovforth/pi/rhoa )**( oneovthird )
679 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
680 "Radius of ", n,
"th aerosol bin (bin center)= ", rada( n ) ,
"[m]" 683 if ( flg_sf_aero )
then 684 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
"flg_sf_aero=true is not supported stop!! " 715 if( radc( n ) > rbnd )
then 720 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,ES15.7,A)')
'Radius between cloud and rain is ', radc(nbnd),
'[m]' 723 if ( rndm_flgp > 0 )
then 724 call random_setup( ia*ja*ka )
727 if ( mp_couple_aerosol .AND.
nccn /=0 )
then 728 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nccn should be 0 when MP_couple_aerosol = .true. !! stop' 732 if (
nccn /= 0 )
then 734 expxactr( n ) = exp( xactr( n ) )
735 rexpxactr( n ) = 1.0_rp / exp( xactr( n ) )
738 expxabnd( n ) = exp( xabnd( n ) )
739 rexpxabnd( n ) = 1.0_rp / exp( xabnd( n ) )
743 allocate( vterm(qa-1) )
747 vterm((myu-1)*
nbin+n) = -vt( myu,n )
751 expxctr( n ) = exp( xctr( n ) )
752 rexpxctr( n ) = 1.0_rp / exp( xctr( n ) )
755 expxbnd( n ) = exp( xbnd( n ) )
756 rexpxbnd( n ) = 1.0_rp / exp( xbnd( n ) )
760 call getrule( ifrsl,kindx )
763 sigma_sdf(1) = 0.2_rp
764 sigma_sdf(2) = 0.35_rp
765 sigma_sdf(3) = 0.35_rp
766 sigma_sdf(4) = 0.35_rp
767 sigma_sdf(5) = 0.35_rp
769 r0_sdf(2) = 2.61e-6_rp
771 r0_sdf(4) = 2.61e-6_rp
772 r0_sdf(5) = 2.61e-6_rp
773 n0_sdf(1) = 8.0e+6_rp
775 n0_sdf(3) = 3.0e+6_rp
776 n0_sdf(4) = 4.0e+6_rp
777 n0_sdf(5) = 4.0e+6_rp
780 rho_sdf(3) = 100.0_rp
781 rho_sdf(4) = 400.0_rp
782 rho_sdf(5) = 400.0_rp
785 sigma_sdf(1) = 0.2_rp
786 sigma_sdf(2) = 0.35_rp
787 sigma_sdf(3) = 0.35_rp
788 sigma_sdf(4) = 0.35_rp
789 sigma_sdf(5) = 0.35_rp
791 r0_sdf(2) = 2.61e-6_rp
793 r0_sdf(4) = 2.61e-6_rp
794 r0_sdf(5) = 2.61e-6_rp
795 n0_sdf(1) = 8.0e+6_rp
797 n0_sdf(3) = 3.0e+6_rp
798 n0_sdf(4) = 4.0e+6_rp
799 n0_sdf(5) = 4.0e+6_rp
802 rho_sdf(3) = 100.0_rp
803 rho_sdf(4) = 400.0_rp
804 rho_sdf(5) = 400.0_rp
812 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
825 atmos_saturation_pres2qsat_liq, &
826 atmos_saturation_pres2qsat_ice
829 integer,
intent(in) :: KA, KS, KE
830 integer,
intent(in) :: IA, IS, IE
831 integer,
intent(in) :: JA, JS, JE
832 integer,
intent(in) :: KIJMAX
834 real(DP),
intent(in) :: dt
835 real(RP),
intent(in) :: DENS (ka,ia,ja)
836 real(RP),
intent(in) :: PRES (ka,ia,ja)
837 real(RP),
intent(in) :: TEMP (ka,ia,ja)
838 real(RP),
intent(in) :: QTRC (ka,ia,ja,qa)
839 real(RP),
intent(in) :: QDRY (ka,ia,ja)
840 real(RP),
intent(in) :: CPtot(ka,ia,ja)
841 real(RP),
intent(in) :: CVtot(ka,ia,ja)
842 real(RP),
intent(in) :: CCN (ka,ia,ja)
844 real(RP),
intent(out) :: RHOQ_t (ka,ia,ja,qa)
845 real(RP),
intent(out) :: RHOE_t (ka,ia,ja)
846 real(RP),
intent(out) :: CPtot_t(ka,ia,ja)
847 real(RP),
intent(out) :: CVtot_t(ka,ia,ja)
848 real(RP),
intent(out) :: EVAPORATE(ka,ia,ja)
850 real(RP) :: QSAT_L(ka,ia,ja)
851 real(RP) :: QSAT_I(ka,ia,ja)
852 real(RP) :: ssliq(ka,ia,ja)
853 real(RP) :: ssice(ka,ia,ja)
855 integer :: ijk_index (kijmax,3)
856 integer :: index_cld (kijmax)
857 integer :: index_cold(kijmax)
858 integer :: index_warm(kijmax)
859 integer :: ijkcount, ijkcount_cold, ijkcount_warm
860 integer :: ijk, indirect
862 real(RP) :: DENS_ijk(kijmax)
863 real(RP) :: PRES_ijk(kijmax)
864 real(RP) :: TEMP_ijk(kijmax)
865 real(RP) :: Qdry_ijk(kijmax)
866 real(RP) :: Qvap_ijk(kijmax)
867 real(RP) :: CCN_ijk(kijmax)
868 real(RP) :: CP_ijk(kijmax)
869 real(RP) :: CV_ijk(kijmax)
870 real(RP) :: Evaporate_ijk(kijmax)
871 real(RP) :: Ghyd_ijk(
nbin,nspc,kijmax)
872 real(RP) :: Gaer_ijk(max(
nccn,1),kijmax)
878 integer :: k, i, j, m, n, iq
881 if ( nspc == 1 )
then 882 log_progress(*)
'atmosphere / physics / microphysics / SBM (Liquid water only)' 883 elseif( nspc > 1 )
then 884 log_progress(*)
'atmosphere / physics / microphysics / SBM (Mixed phase)' 887 call atmos_saturation_pres2qsat_liq( ka, ks, ke, &
890 temp(:,:,:), pres(:,:,:), qdry(:,:,:), &
893 call atmos_saturation_pres2qsat_ice( ka, ks, ke, &
896 temp(:,:,:), pres(:,:,:), qdry(:,:,:), &
902 ssliq(k,i,j) = qtrc(k,i,j,i_qv) / qsat_l(k,i,j) - 1.0_rp
903 ssice(k,i,j) = qtrc(k,i,j,i_qv) / qsat_i(k,i,j) - 1.0_rp
908 if ( nspc == 1 )
then 909 ssice(:,:,:) = 0.0_rp
964 cldsum = cldsum + qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
965 countbin = countbin + 1
969 if ( cldsum > cldmin &
970 .OR. ssliq(k,i,j) > 0.0_rp &
971 .OR. ssice(k,i,j) > 0.0_rp )
then 973 ijkcount = ijkcount + 1
975 index_cld(ijkcount) = ijk
977 dens_ijk(ijkcount) = dens(k,i,j)
978 pres_ijk(ijkcount) = pres(k,i,j)
979 temp_ijk(ijkcount) = temp(k,i,j)
980 qdry_ijk(ijkcount) = qdry(k,i,j)
981 cp_ijk(ijkcount) = cptot(k,i,j)
982 cv_ijk(ijkcount) = cvtot(k,i,j)
983 ccn_ijk(ijkcount) = ccn(k,i,j)
984 qvap_ijk(ijkcount) = qtrc(k,i,j,i_qv)
989 ghyd_ijk(n,m,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
990 countbin = countbin + 1
995 gaer_ijk(n,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxaer
996 countbin = countbin + 1
999 if ( temp(k,i,j) < tem00 .AND. nspc > 1 )
then 1000 ijkcount_cold = ijkcount_cold + 1
1001 index_cold(ijkcount_cold) = ijkcount
1003 ijkcount_warm = ijkcount_warm + 1
1004 index_warm(ijkcount_warm) = ijkcount
1011 rhoq_t(k,i,j,iq) = 0.0_rp
1013 rhoe_t(k,i,j) = 0.0_rp
1014 cptot_t(k,i,j) = 0.0_rp
1015 cvtot_t(k,i,j) = 0.0_rp
1016 evaporate(k,i,j) = 0.0_rp
1058 if ( ijkcount > 0 )
then 1062 call mp_suzuki10( ka, ia, ja, &
1066 index_cold( 1:ijkcount), &
1067 index_warm( 1:ijkcount), &
1068 dens_ijk( 1:ijkcount), &
1069 pres_ijk( 1:ijkcount), &
1070 qdry_ijk( 1:ijkcount), &
1071 ccn_ijk( 1:ijkcount), &
1072 temp_ijk( 1:ijkcount), &
1073 qvap_ijk( 1:ijkcount), &
1074 ghyd_ijk(:,:,1:ijkcount), &
1075 gaer_ijk(:, 1:ijkcount), &
1076 cp_ijk( 1:ijkcount), &
1077 cv_ijk( 1:ijkcount), &
1078 evaporate_ijk(1:ijkcount), &
1128 do ijk = 1, ijkcount
1129 indirect = index_cld(ijk)
1130 i = ijk_index(indirect,1)
1131 j = ijk_index(indirect,2)
1132 k = ijk_index(indirect,3)
1134 rhoe_t(k,i,j) = ( temp_ijk(ijk) * cv_ijk(ijk) - temp(k,i,j) * cvtot(k,i,j) ) * dens(k,i,j) / dt
1135 cptot_t(k,i,j) = ( cp_ijk(ijk) - cptot(k,i,j) ) / dt
1136 cvtot_t(k,i,j) = ( cv_ijk(ijk) - cvtot(k,i,j) ) / dt
1137 evaporate(k,i,j) = evaporate_ijk(ijk) / dt
1139 rhoq_t(k,i,j,i_qv) = ( qvap_ijk(ijk) - qtrc(k,i,j,i_qv) ) * dens(k,i,j) / dt
1144 rhoq_new = ghyd_ijk(n,m,ijk) * dxmic
1145 rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
1146 countbin = countbin + 1
1151 rhoq_new = gaer_ijk(n,ijk) * dxaer
1152 rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
1153 countbin = countbin + 1
1186 integer,
intent(in) :: KA
1188 real(RP),
intent(out) :: vterm_o(ka,qa-1)
1194 vterm_o(:,iq) = vterm(iq)
1211 integer,
intent(in) :: KA, KS, KE
1212 integer,
intent(in) :: IA, IS, IE
1213 integer,
intent(in) :: JA, JS, JE
1215 real(RP),
intent(in) :: QTRC0 (ka,ia,ja,nspc*
nbin)
1216 real(RP),
intent(in) :: mask_criterion
1217 real(RP),
intent(out) :: cldfrac(ka,ia,ja)
1220 integer :: k, i, j, iq, ihydro
1229 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1230 qhydro = qhydro + qtrc0(k,i,j,iq)
1233 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1237 elseif( nspc == 1 )
then 1242 do ihydro = 1, i_mp_qc
1243 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1244 qhydro = qhydro + qtrc0(k,i,j,iq)
1247 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1279 integer,
intent(in) :: KA, KS, KE
1280 integer,
intent(in) :: IA, IS, IE
1281 integer,
intent(in) :: JA, JS, JE
1283 real(RP),
intent(in) :: DENS0(ka,ia,ja)
1284 real(RP),
intent(in) :: TEMP0(ka,ia,ja)
1285 real(RP),
intent(in) :: QTRC0(ka,ia,ja,nspc*
nbin)
1286 real(RP),
intent(out) :: Re (ka,ia,ja,
n_hyd)
1288 real(RP),
parameter :: um2cm = 100.0_rp
1290 real(RP) :: sum0(nspc), sum2, sum3, re_tmp(nspc)
1291 integer :: i, j, k, iq, ihydro
1297 re(k,i,j,:) = 0.0_rp
1305 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1306 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1307 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1309 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1310 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1311 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1313 sum3 = max( sum3, 0.0_rp )
1314 sum2 = max( sum2, 0.0_rp )
1315 if ( sum2 /= 0.0_rp )
then 1316 re(k,i,j,
i_hc) = sum3 / sum2 * um2cm
1318 re(k,i,j,
i_hc) = 0.0_rp
1325 do iq = nbnd+1,
nbin 1327 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1328 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1329 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1331 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1332 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1333 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1335 sum3 = max( sum3, 0.0_rp )
1336 sum2 = max( sum2, 0.0_rp )
1337 if ( sum2 /= 0.0_rp )
then 1338 re(k,i,j,
i_hr) = sum3 / sum2 * um2cm
1340 re(k,i,j,
i_hr) = 0.0_rp
1348 if ( nspc > 1 )
then 1353 sum0(ihydro) = 0.0_rp
1356 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1357 sum0(ihydro) = sum0(ihydro) &
1358 + ( qtrc0(k,i,j,iq) * dens0(k,i,j) )
1360 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1361 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1362 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1364 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1365 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1366 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1368 sum3 = max( sum3, 0.0_rp )
1369 sum2 = max( sum2, 0.0_rp )
1370 if ( sum2 == 0.0_rp )
then 1371 re_tmp(ihydro) = 0.0_rp
1373 re_tmp(ihydro) = sum3 / sum2 * um2cm
1377 re(k,i,j,
i_hi) = ( re_tmp(i_mp_qcl) * sum0(i_mp_qcl) &
1378 + re_tmp(i_mp_qp ) * sum0(i_mp_qp ) &
1379 + re_tmp(i_mp_qd ) * sum0(i_mp_qd ) ) &
1380 / ( sum0(i_mp_qcl) + sum0(i_mp_qp) + sum0(i_mp_qd) + eps )
1381 re(k,i,j,
i_hs) = re_tmp(i_mp_qs)
1382 re(k,i,j,
i_hg) = re_tmp(i_mp_qg)
1383 re(k,i,j,
i_hh) = re_tmp(i_mp_qh)
1412 integer,
intent(in) :: KA, KS, KE
1413 integer,
intent(in) :: IA, IS, IE
1414 integer,
intent(in) :: JA, JS, JE
1415 real(RP),
intent(in) :: QTRC0(ka,ia,ja,
nbin*nspc)
1416 real(RP),
intent(out) :: Qe (ka,ia,ja,
n_hyd)
1418 integer :: ihydro, ibin, iq, icateg
1423 qe(:,:,:,:) = 0.0_rp
1427 iq =
nbin*(ihydro-1) + ibin
1429 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc -1) )
then 1430 if ( iq > 0 .AND. iq <= nbnd )
then 1432 elseif( iq > nbnd .AND. iq <=
nbin )
then 1435 elseif( iq >
nbin*(i_mp_qc -1) .AND. iq <=
nbin*(i_mp_qcl-1) )
then 1437 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qp -1) )
then 1439 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qs -1) )
then 1441 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qg -1) )
then 1443 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qh -1) )
then 1445 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <=
nbin*nspc )
then 1452 qe(k,i,j,icateg) = qe(k,i,j,icateg) + qtrc0(k,i,j,iq)
1480 integer,
intent(in) :: KA, KS, KE
1481 integer,
intent(in) :: IA, IS, IE
1482 integer,
intent(in) :: JA, JS, JE
1483 real(RP),
intent(in) :: DENS (ka,ia,ja)
1484 real(RP),
intent(in) :: QTRC0(ka,ia,ja,
nbin*nspc)
1485 real(RP),
intent(out) :: Ne (ka,ia,ja,
n_hyd)
1487 integer :: ihydro, ibin, iq, icateg
1492 ne(:,:,:,:) = 0.0_rp
1496 iq =
nbin*(ihydro-1) + ibin
1498 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc -1) )
then 1499 if ( iq > 0 .AND. iq <= nbnd )
then 1501 elseif( iq > nbnd .AND. iq <=
nbin )
then 1504 elseif( iq >
nbin*(i_mp_qc -1) .AND. iq <=
nbin*(i_mp_qcl-1) )
then 1506 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qp -1) )
then 1508 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qs -1) )
then 1510 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qg -1) )
then 1512 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qh -1) )
then 1514 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <=
nbin*nspc )
then 1521 ne(k,i,j,icateg) = ne(k,i,j,icateg) + dens(k,i,j) * qtrc0(k,i,j,iq) * rexpxctr(ibin)
1534 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1552 integer,
intent(in) :: KA, KS, KE
1553 integer,
intent(in) :: IA, IS, IE
1554 integer,
intent(in) :: JA, JS, JE
1556 real(RP),
intent(in) :: Qe(ka,ia,ja,
n_hyd)
1558 real(RP),
intent(out) :: QTRC(ka,ia,ja,qa-1)
1560 real(RP),
intent(in),
optional :: QNUM(ka,ia,ja,
n_hyd)
1562 real(RP) :: coef0, coef1, coef2
1563 real(RP) :: dummy(
nbin)
1564 real(RP) :: tmp_hyd, num_hyd, lambda_hyd
1566 integer :: k, i, j, iq
1568 if (
present(qnum) )
then 1569 log_warn(
"ATMOS_PHY_MP_suzuki10_qhyd2qtrc",*)
'At this moment, number concentratio is ignored' 1573 coef0 = 4.0_rp/3.0_rp*pi
1574 coef1 = 4.0_rp/3.0_rp*sqrt(pi/2.0_rp)
1576 if( nspc == 1 )
then 1584 dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
1586 - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
1587 / sigma_sdf(1) / sigma_sdf(1) &
1589 tmp_hyd = tmp_hyd + dummy(iq)
1592 coef2 = ( qe(k,i,j,
i_hc) + qe(k,i,j,
i_hr) &
1594 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1597 qtrc(k+2,i,j,(il-1)*
nbin+iq) = coef2 * dummy(iq)
1604 elseif( nspc > 1 )
then 1613 dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
1615 - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
1616 / sigma_sdf(1) / sigma_sdf(1) &
1618 tmp_hyd = tmp_hyd + dummy(iq)
1621 coef2 = ( qe(k,i,j,
i_hc) + qe(k,i,j,
i_hr) ) &
1622 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1625 qtrc(k,i,j,(il-1)*
nbin+iq) = coef2 * dummy(iq)
1631 dummy(iq) = coef1 / sigma_sdf(2) * rho_sdf(2) * radc( iq )**3 &
1633 - ( log( radc(iq) )-log( r0_sdf(2) ) )**2*0.5_rp &
1634 / sigma_sdf(2) / sigma_sdf(2) &
1636 tmp_hyd = tmp_hyd + dummy(iq)
1639 coef2 = ( qe(k,i,j,
i_hi) ) &
1640 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1643 qtrc(k,i,j,(ip-1)*
nbin+iq) = coef2 * dummy(iq)
1647 num_hyd = coef0 * n0_sdf(3) * rho_sdf(3)
1648 lambda_hyd = ( pi * rho_sdf(3) / 6.0_rp *n0_sdf(3) *
sf_gamma(4.0_rp) &
1649 / ( qe(k,i,j,
i_hs) &
1650 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hs)-eps)) &
1655 dummy(iq) = num_hyd * radc( iq )**3 &
1656 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1657 tmp_hyd = tmp_hyd + dummy(iq)
1660 coef2 = ( qe(k,i,j,
i_hs) ) &
1661 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1664 qtrc(k,i,j,(iss-1)*
nbin+iq) = coef2 * dummy(iq)
1668 num_hyd = coef0 * n0_sdf(4) * rho_sdf(4)
1669 lambda_hyd = ( pi * rho_sdf(4) / 6.0_rp *n0_sdf(4) *
sf_gamma(4.0_rp) &
1670 / ( qe(k,i,j,
i_hg) &
1671 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hg)-eps)) &
1676 dummy(iq) = num_hyd * radc( iq )**3 &
1677 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1678 tmp_hyd = tmp_hyd + dummy(iq)
1681 coef2 = ( qe(k,i,j,
i_hg) ) &
1682 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1685 qtrc(k,i,j,(ig-1)*
nbin+iq) = coef2 * dummy(iq)
1689 num_hyd = coef0 * n0_sdf(5) * rho_sdf(5)
1690 lambda_hyd = ( pi * rho_sdf(5) / 6.0_rp *n0_sdf(5) *
sf_gamma(4.0_rp) &
1691 / ( qe(k,i,j,
i_hh) &
1692 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hh)-eps)) &
1697 dummy(iq) = num_hyd * radc( iq )**3 &
1698 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1699 tmp_hyd = tmp_hyd + dummy(iq)
1702 coef2 = ( qe(k,i,j,
i_hh) ) &
1703 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1706 qtrc(k,i,j,(ih-1)*
nbin+iq) = coef2 * dummy(iq)
1715 do iq =
nbin*nspc+1, qa-1
1719 qtrc(k,i,j,iq) = 0.0_rp
1729 subroutine mp_suzuki10( &
1750 integer,
intent(in) :: KA
1751 integer,
intent(in) :: IA
1752 integer,
intent(in) :: JA
1754 integer,
intent(in) :: ijkmax
1755 integer,
intent(in) :: ijkmax_cold
1756 integer,
intent(in) :: ijkmax_warm
1757 integer ,
intent(in) :: index_cold(ijkmax)
1758 integer ,
intent(in) :: index_warm(ijkmax)
1759 real(RP),
intent(in) :: dens (ijkmax)
1760 real(RP),
intent(in) :: pres (ijkmax)
1761 real(RP),
intent(in) :: qdry (ijkmax)
1762 real(RP),
intent(in) :: ccn (ijkmax)
1763 real(RP),
intent(inout) :: temp (ijkmax)
1764 real(RP),
intent(inout) :: qvap (ijkmax)
1765 real(RP),
intent(inout) :: ghyd (
nbin,nspc,ijkmax)
1766 real(RP),
intent(inout) :: gaer (max(
nccn,1),ijkmax)
1767 real(RP),
intent(inout) :: cp (ijkmax)
1768 real(RP),
intent(inout) :: cv (ijkmax)
1769 real(RP),
intent(out) :: evaporate (ijkmax)
1770 real(DP),
intent(in) :: dt
1773 if (
nccn /= 0 )
then 1774 if ( nspc == 1 )
then 1778 call nucleata( ijkmax, &
1791 call cndevpsbla( ijkmax, &
1804 if ( mp_doautoconversion )
then 1806 call collmain( ka, ia, ja, &
1813 elseif( nspc > 1 )
then 1817 call nucleata( ijkmax, &
1830 call freezing( ijkmax, &
1853 call melting( ijkmax, &
1864 call cndevpsbla( ijkmax, &
1877 if ( mp_doautoconversion )
then 1879 call collmainf( ka, ia, ja, &
1888 elseif(
nccn == 0 )
then 1890 if ( nspc == 1 )
then 1894 call nucleat( ijkmax, &
1907 call cndevpsbl( ijkmax, &
1919 if ( mp_doautoconversion )
then 1921 call collmain( ka, ia, ja, &
1928 elseif( nspc > 1 )
then 1932 call nucleat( ijkmax, &
1945 call freezing( ijkmax, &
1955 call ice_nucleat( ijkmax, &
1968 call melting( ijkmax, &
1979 call cndevpsbl( ijkmax, &
1991 if ( mp_doautoconversion )
then 1993 call collmainf( ka, ia, ja, &
2005 end subroutine mp_suzuki10
2008 subroutine nucleat( &
2021 atmos_hydrometeor_lhv, &
2027 atmos_saturation_pres2qsat_liq
2030 integer,
intent(in) :: ijkmax
2031 real(RP),
intent(in) :: dens(ijkmax)
2032 real(RP),
intent(in) :: pres(ijkmax)
2033 real(RP),
intent(in) :: qdry(ijkmax)
2034 real(RP),
intent(in) :: ccn(ijkmax)
2035 real(RP),
intent(inout) :: temp(ijkmax)
2036 real(RP),
intent(inout) :: qvap(ijkmax)
2037 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2038 real(RP),
intent(inout) :: cp(ijkmax)
2039 real(RP),
intent(inout) :: cv(ijkmax)
2040 real(DP),
intent(in) :: dtime
2042 real(RP) :: ssliq(ijkmax)
2043 real(RP) :: qlevp(ijkmax)
2049 real(RP) :: sumnum(ijkmax)
2050 real(RP) :: gcn(
nbin,ijkmax )
2051 real(RP) :: qsat(ijkmax)
2057 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2059 if( mp_couple_aerosol )
then 2063 dmp = ccn(ijk) * expxctr( 1 )
2064 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2065 gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
2067 qvap(ijk) = qvap(ijk) - dqv
2068 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cp(ijk)
2075 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2076 temp(:), pres(:), qdry(:), &
2080 ssliq(ijk) = qvap(ijk)/qsat(ijk) - 1.0_rp
2086 if ( ssliq(ijk) > 0.0_rp )
then 2091 gcn( n,ijk ) = gc( n,il,ijk )*rexpxctr( n )
2095 sumnum(ijk) = sumnum(ijk) + gcn( n,ijk )*dxmic
2097 n_c = c_ccn * ( ssliq(ijk) * 1.e+2_rp )**( kappa )
2098 if ( n_c > sumnum(ijk) )
then 2099 dmp = ( n_c - sumnum(ijk) ) * expxctr( 1 )
2100 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2101 gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
2103 qvap(ijk) = qvap(ijk) - dqv
2104 qvap(ijk) = max( qvap(ijk),0.0_rp )
2105 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cp(ijk)
2118 end subroutine nucleat
2121 subroutine nucleata( &
2138 atmos_hydrometeor_lhv, &
2144 atmos_saturation_pres2qsat_liq, &
2145 atmos_saturation_pres2qsat_ice
2148 integer,
intent(in) :: ijkmax
2149 real(RP),
intent(in) :: dens(ijkmax)
2150 real(RP),
intent(in) :: pres(ijkmax)
2151 real(RP),
intent(in) :: qdry(ijkmax)
2152 real(RP),
intent(inout) :: temp(ijkmax)
2153 real(RP),
intent(inout) :: qvap(ijkmax)
2154 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2155 real(RP),
intent(inout) :: ga (
nccn ,ijkmax)
2156 real(RP),
intent(inout) :: cp(ijkmax)
2157 real(RP),
intent(inout) :: cv(ijkmax)
2158 real(DP),
intent(in) :: dtime
2160 real(RP) :: gan(
nccn )
2162 real(RP) :: acoef, bcoef
2165 real(RP) :: ractr, rcld, xcld, part, dmp
2166 integer :: n, nc, ncrit
2171 real(RP) :: qlevp(ijkmax)
2172 real(RP) :: qsatl(ijkmax)
2179 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2181 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2182 temp(:), pres(:), qdry(:), &
2186 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2188 if ( ssliq <= 0.0_rp ) cycle
2193 gan( n ) = ga( n,ijk )*rexpxactr( n )
2196 acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk)
2197 bcoef = vhfct* rhoa/rhow * emwtr/emaer
2201 ractr = ( expxactr( n )*thirdovforth/pi/rhoa )**( oneovthird )
2202 rcld = sqrt( 3.0_rp*bcoef*ractr*ractr*ractr / acoef )
2203 xcld = log( rhow * 4.0_rp*pi*oneovthird*rcld*rcld*rcld )
2204 if ( flg_nucl )
then 2207 ncld( n ) = int( ( xcld-xctr( 1 ) )/dxmic ) + 1
2208 ncld( n ) = min( max( ncld( n ),1 ),
nbin )
2214 if ( ssliq <= 0.0_rp )
exit 2216 acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk)
2217 rcrit = acoef*oneovthird * ( 4.0_rp/bcoef )**( oneovthird ) / ssliq**( twoovthird )
2218 xcrit = log( rhoa * 4.0_rp*pi*oneovthird * rcrit*rcrit*rcrit )
2219 ncrit = int( ( xcrit-xabnd( 1 ) )/dxaer ) + 1
2221 if ( n == ncrit )
then 2222 part = ( xabnd( ncrit+1 )-xcrit )/dxaer
2223 elseif ( n > ncrit )
then 2231 dmp = part*gan( n )*dxaer*expxctr( nc )
2232 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2233 gc( nc,il,ijk ) = gc( nc,il,ijk ) + dmp/dxmic
2234 gan( n ) = gan( n ) - dmp/dxaer*rexpxctr( nc )
2235 gan( n ) = max( gan( n ), 0.0_rp )
2237 qvap(ijk) = qvap(ijk) - dqv
2238 qvap(ijk) = max( qvap(ijk),0.0_rp )
2239 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cp(ijk)
2246 ga( n,ijk ) = gan( n )*expxactr( n )
2254 end subroutine nucleata
2257 subroutine cndevpsbl( &
2271 integer,
intent(in) :: ijkmax
2272 real(RP),
intent(in) :: dens(ijkmax)
2273 real(RP),
intent(in) :: pres(ijkmax)
2274 real(RP),
intent(in) :: qdry(ijkmax)
2275 real(RP),
intent(inout) :: temp(ijkmax)
2276 real(RP),
intent(inout) :: qvap(ijkmax)
2277 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2278 real(RP),
intent(inout) :: cp(ijkmax)
2279 real(RP),
intent(inout) :: cv(ijkmax)
2280 real(RP),
intent(out) :: evaporate(ijkmax)
2281 real(DP),
intent(in) :: dtime
2284 call liqphase( ijkmax, &
2297 call icephase( ijkmax, &
2308 call mixphase( ijkmax, &
2321 end subroutine cndevpsbl
2324 subroutine cndevpsbla( &
2339 integer,
intent(in) :: ijkmax
2340 real(RP),
intent(in) :: dens(ijkmax)
2341 real(RP),
intent(in) :: pres(ijkmax)
2342 real(RP),
intent(in) :: qdry(ijkmax)
2343 real(RP),
intent(inout) :: temp(ijkmax)
2344 real(RP),
intent(inout) :: qvap(ijkmax)
2345 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2346 real(RP),
intent(inout) :: ga (
nccn ,ijkmax)
2347 real(RP),
intent(inout) :: cp(ijkmax)
2348 real(RP),
intent(inout) :: cv(ijkmax)
2349 real(RP),
intent(out) :: evaporate(ijkmax)
2350 real(DP),
intent(in) :: dtime
2353 call liqphase( ijkmax, &
2366 if ( flg_regeneration )
then 2367 call faero( ijkmax, &
2374 call icephase( ijkmax, &
2385 call mixphase( ijkmax, &
2398 end subroutine cndevpsbla
2401 subroutine liqphase( &
2421 atmos_hydrometeor_lhv, &
2427 atmos_saturation_pres2qsat_liq
2430 integer,
intent(in) :: ijkmax
2431 real(RP),
intent(in) :: dens(ijkmax)
2432 real(RP),
intent(in) :: pres(ijkmax)
2433 real(RP),
intent(in) :: qdry(ijkmax)
2434 real(RP),
intent(inout) :: temp(ijkmax)
2435 real(RP),
intent(inout) :: qvap(ijkmax)
2436 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2437 real(RP),
intent(inout) :: cp(ijkmax)
2438 real(RP),
intent(inout) :: cv(ijkmax)
2439 real(RP),
intent(out) :: regene_gcn(ijkmax)
2440 real(DP),
intent(in) :: dtime
2442 integer :: n, myu, ncount
2443 integer :: nloop(ijkmax)
2444 real(RP) :: gclold(ijkmax)
2445 real(RP) :: gclnew(ijkmax)
2447 real(RP) :: dtcnd(ijkmax)
2448 real(RP) :: gtliq(ijkmax)
2449 real(RP) :: qlevp(ijkmax)
2450 real(RP) :: cefliq, a, sliqtnd
2451 real(RP) :: sumliq(ijkmax), umax(ijkmax)
2453 real(RP) :: gcn( -1:
nbin+2,nspc,ijkmax )
2455 real(RP),
parameter :: cflfct = 0.50_rp
2457 integer :: iflg( nspc,ijkmax )
2458 real(RP) :: csum( nspc,ijkmax )
2459 real(RP) :: f1, f2, emu, cefd, cefk, festl
2460 real(RP) :: qsatl(ijkmax)
2461 real(RP),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2462 real(RP),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2463 real(RP) :: zerosw, qvtmp
2465 integer :: ijk, nn, mm, loopflg(ijkmax)
2468 real(RP) :: uadv ( 0:
nbin+2,nspc,ijkmax )
2469 real(RP) :: flq ( 1:
nbin+1,nspc,ijkmax )
2470 real(RP) :: acoef( 0:2,0:
nbin+1,nspc,ijkmax )
2471 real(RP) :: crn ( 0:
nbin+2,nspc,ijkmax )
2472 real(RP) :: aip ( 0:
nbin+1,nspc,ijkmax )
2473 real(RP) :: aim ( 0:
nbin+1,nspc,ijkmax )
2474 real(RP) :: ai ( 0:
nbin+1,nspc,ijkmax )
2475 real(RP) :: cmins, cplus
2484 csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
2486 if( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
2490 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2492 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2493 temp(:), pres(:), qdry(:), &
2501 gclold(ijk) = gclold(ijk) + gc( n,il,ijk ) * dxmic
2506 gcn( n,il,ijk ) = gc( n,il,ijk ) * rexpxctr( n )
2508 gcn( -1,il,ijk ) = 0.0_rp
2509 gcn( 0,il,ijk ) = 0.0_rp
2510 gcn(
nbin+1,il,ijk ) = 0.0_rp
2511 gcn(
nbin+2,il,ijk ) = 0.0_rp
2514 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2516 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2517 cefd = emu/dens(ijk)
2520 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2521 f1 = rvap*temp(ijk)/festl/cefd
2522 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2523 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2525 umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq )
2526 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2527 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2528 dtcnd(ijk) = dtime / nloop(ijk)
2530 nloop(ijk) = nloop(ijk) * iflg( il,ijk )
2532 nloopmax = maxval(nloop,1)
2535 regene_gcn(:) = 0.0_rp
2538 do ncount = 1, nloopmax
2541 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
2545 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2549 do nn = 1, loopflg(ijk)
2551 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2552 cefd = emu/dens(ijk)
2554 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2555 f1 = rvap*temp(ijk)/festl/cefd
2556 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2557 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2559 sumliq(ijk) = 0.0_rp
2562 sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
2569 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2570 temp(:), pres(:), qdry(:), &
2574 do nn = 1, loopflg(ijk)
2576 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2579 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
2580 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
2581 cefliq = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)*qlevp(ijk)/cp(ijk)/rvap/temp(ijk)/temp(ijk) )
2582 a = - cefliq*sumliq(ijk)*gtliq(ijk)/dens(ijk)
2583 a = a + eps * ( 1.0_rp - zerosw )
2585 sliqtnd = zerosw * &
2586 ( ssliq * ( exp( a*dtcnd(ijk) )-1.0_rp )/( a*dtcnd(ijk) ) &
2587 * ( 0.5_rp + sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
2589 * ( 0.5_rp - sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
2591 + ssliq * ( 1.0_rp - zerosw )
2594 do mm = 1, iflg( il,ijk )
2597 uadv( n,il,ijk ) = cbnd( n,il )*rexpxbnd( n )*gtliq(ijk)*sliqtnd
2599 uadv( 0 ,il,ijk ) = 0.0_rp
2600 uadv(
nbin+2,il,ijk ) = 0.0_rp
2613 do nn = 1, loopflg(ijk)
2615 do mm = 1, iflg( myu,ijk )
2617 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
2625 do nn = 1, loopflg(ijk)
2627 do mm = 1, iflg( myu,ijk )
2629 acoef(0,n,myu,ijk) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
2630 acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
2631 acoef(2,n,myu,ijk) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
2639 do nn = 1, loopflg(ijk)
2641 do mm = 1, iflg( myu,ijk )
2643 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
2645 aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
2646 + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
2647 + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
2649 aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
2657 do nn = 1, loopflg(ijk)
2659 do mm = 1, iflg( myu,ijk )
2661 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
2663 aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
2664 - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
2665 + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
2667 aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
2675 do nn = 1, loopflg(ijk)
2677 do mm = 1, iflg( myu,ijk )
2679 ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
2680 + acoef(2,n,myu,ijk) * 2.0_rp
2682 ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
2690 do nn = 1, loopflg(ijk)
2692 do mm = 1, iflg( myu,ijk )
2694 flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
2695 - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
2703 do nn = 1, loopflg(ijk)
2705 do mm = 1, iflg( myu,ijk )
2706 regene_gcn(ijk) = regene_gcn(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
2707 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
2714 do nn = 1, loopflg(ijk)
2716 do mm = 1, iflg( myu,ijk )
2718 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
2729 do nn = 1, loopflg(ijk)
2732 gclnew(ijk) = 0.0_rp
2735 gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )
2740 gclnew(ijk) = gclnew(ijk)*dxmic
2743 cndmss = gclnew(ijk) - gclold(ijk)
2744 dqv = cndmss/dens(ijk)
2745 qvap(ijk) = qvap(ijk) - dqv
2746 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cp(ijk)
2750 gclold(ijk) = gclnew(ijk)
2762 gc( n,il,ijk ) = gcn( n,il,ijk )*expxctr( n )
2767 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2772 if ( gc( n,il,ijk ) < 0.0_rp )
then 2773 cndmss = -gc( n,il,ijk )
2774 gc( n,il,ijk ) = 0.0_rp
2775 dqv = cndmss/dens(ijk)
2776 qvap(ijk) = qvap(ijk) + dqv
2777 temp(ijk) = temp(ijk) - dqv*qlevp(ijk)/cp(ijk)
2787 end subroutine liqphase
2790 subroutine icephase( &
2809 atmos_hydrometeor_lhs, &
2815 atmos_saturation_pres2qsat_ice
2818 integer,
intent(in) :: ijkmax
2819 real(RP),
intent(in) :: dens(ijkmax)
2820 real(RP),
intent(in) :: pres(ijkmax)
2821 real(RP),
intent(in) :: qdry(ijkmax)
2822 real(RP),
intent(inout) :: temp(ijkmax)
2823 real(RP),
intent(inout) :: qvap(ijkmax)
2824 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2825 real(RP),
intent(inout) :: cp(ijkmax)
2826 real(RP),
intent(inout) :: cv(ijkmax)
2827 real(DP),
intent(in) :: dtime
2829 integer :: myu, n, ncount
2830 integer :: nloop(ijkmax)
2831 real(RP) :: gciold(ijkmax), gcinew(ijkmax)
2832 real(RP) :: dtcnd(ijkmax)
2834 real(RP) :: gtice(ijkmax), umax(ijkmax)
2835 real(RP) :: qlsbl(ijkmax)
2836 real(RP) :: cefice, d, uval, sicetnd
2837 real(RP) :: sumice(ijkmax)
2839 real(RP) :: gcn( -1:
nbin+2,nspc,ijkmax )
2840 real(RP),
parameter :: cflfct = 0.50_rp
2841 real(RP) :: dumm_regene(ijkmax)
2842 integer :: iflg( nspc,ijkmax )
2843 real(RP) :: csum( nspc,ijkmax )
2844 real(RP) :: f1, f2, emu, cefd, cefk, festi
2845 real(RP) :: qsati(ijkmax)
2846 real(RP),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2847 real(RP),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2848 integer :: ijk, mm, nn, loopflg(ijkmax)
2849 real(RP) :: zerosw, qvtmp
2854 real(RP) :: uadv ( 0:
nbin+2,nspc,ijkmax )
2855 real(RP) :: flq ( 1:
nbin+1,nspc,ijkmax )
2856 real(RP) :: acoef( 0:2,0:
nbin+1,nspc,ijkmax )
2857 real(RP) :: crn ( 0:
nbin+2,nspc,ijkmax )
2858 real(RP) :: aip ( 0:
nbin+1,nspc,ijkmax )
2859 real(RP) :: aim ( 0:
nbin+1,nspc,ijkmax )
2860 real(RP) :: ai ( 0:
nbin+1,nspc,ijkmax )
2861 real(RP) :: cmins, cplus
2868 csum( :,: ) = 0.0_rp
2873 csum( myu,ijk ) = csum( myu,ijk ) + gc( n,myu,ijk )*dxmic
2878 if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
2884 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
2886 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
2887 temp(:), pres(:), qdry(:), &
2897 gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
2904 gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
2906 gcn( -1,myu,ijk ) = 0.0_rp
2907 gcn( 0,myu,ijk ) = 0.0_rp
2908 gcn(
nbin+1,myu,ijk ) = 0.0_rp
2909 gcn(
nbin+2,myu,ijk ) = 0.0_rp
2913 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
2915 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2916 cefd = emu/dens(ijk)
2918 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2919 f1 = rvap*temp(ijk)/festi/cefd
2920 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
2921 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
2925 uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice )
2926 umax(ijk) = max( umax(ijk),uval )
2929 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2930 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2931 dtcnd(ijk) = dtime/nloop(ijk)
2933 nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) )
2936 nloopmax = maxval(nloop,1)
2940 do ncount = 1, nloopmax
2943 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
2947 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
2951 do nn = 1, loopflg(ijk)
2952 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2953 cefd = emu/dens(ijk)
2955 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2956 f1 = rvap*temp(ijk)/festi/cefd
2957 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
2958 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
2960 sumice(ijk) = 0.0_rp
2963 sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
2971 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
2972 temp(:), pres(:), qdry(:), &
2977 do nn = 1, loopflg(ijk)
2980 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
2983 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
2984 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
2986 cefice = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)*qlsbl(ijk)/cp(ijk)/rvap/temp(ijk)/temp(ijk) )
2987 d = - cefice*sumice(ijk)*gtice(ijk)/dens(ijk)
2988 d = d + eps * ( 1.0_rp - zerosw )
2990 sicetnd = zerosw * &
2991 ( ssice * ( exp( d*dtcnd(ijk) )-1.0_rp )/( d*dtcnd(ijk) ) &
2992 * ( 0.5_rp + sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
2994 * ( 0.5_rp - sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
2996 + ssice * ( 1.0_rp - zerosw )
3000 do mm = 1, iflg( myu,ijk )
3003 uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd
3005 uadv( 0, myu,ijk ) = 0.0_rp
3006 uadv(
nbin+2,myu,ijk ) = 0.0_rp
3016 do nn = 1, loopflg(ijk)
3018 do mm = 1, iflg( myu,ijk )
3020 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
3028 do nn = 1, loopflg(ijk)
3030 do mm = 1, iflg( myu,ijk )
3032 acoef(0,n,myu,ijk) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3033 acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3034 acoef(2,n,myu,ijk) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3042 do nn = 1, loopflg(ijk)
3044 do mm = 1, iflg( myu,ijk )
3046 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3048 aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3049 + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3050 + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3052 aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
3060 do nn = 1, loopflg(ijk)
3062 do mm = 1, iflg( myu,ijk )
3064 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3066 aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3067 - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3068 + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3070 aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
3078 do nn = 1, loopflg(ijk)
3080 do mm = 1, iflg( myu,ijk )
3082 ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
3083 + acoef(2,n,myu,ijk) * 2.0_rp
3085 ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
3093 do nn = 1, loopflg(ijk)
3095 do mm = 1, iflg( myu,ijk )
3097 flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
3098 - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3106 do nn = 1, loopflg(ijk)
3108 do mm = 1, iflg( myu,ijk )
3109 dumm_regene(ijk) = dumm_regene(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
3110 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
3117 do nn = 1, loopflg(ijk)
3119 do mm = 1, iflg( myu,ijk )
3121 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
3132 do nn = 1, loopflg(ijk)
3135 gcinew(ijk) = 0.0_rp
3138 gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
3143 sblmss = gcinew(ijk) - gciold(ijk)
3144 dqv = sblmss/dens(ijk)
3145 qvap(ijk) = qvap(ijk) - dqv
3146 temp(ijk) = temp(ijk) + dqv*qlsbl(ijk)/cp(ijk)
3150 gciold(ijk) = gcinew(ijk)
3163 gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
3171 end subroutine icephase
3174 subroutine mixphase( &
3193 atmos_hydrometeor_lhv, &
3194 atmos_hydrometeor_lhs, &
3202 atmos_saturation_pres2qsat_liq, &
3203 atmos_saturation_pres2qsat_ice
3206 integer,
intent(in) :: ijkmax
3207 real(RP),
intent(in) :: dens(ijkmax)
3208 real(RP),
intent(in) :: pres(ijkmax)
3209 real(RP),
intent(in) :: qdry(ijkmax)
3210 real(RP),
intent(inout) :: temp(ijkmax)
3211 real(RP),
intent(inout) :: qvap(ijkmax)
3212 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3213 real(RP),
intent(inout) :: cp(ijkmax)
3214 real(RP),
intent(inout) :: cv(ijkmax)
3215 real(DP),
intent(in) :: dtime
3217 integer :: n, myu, mm, ncount
3218 integer :: nloop(ijkmax)
3219 real(RP) :: gclold(ijkmax), gclnew(ijkmax)
3220 real(RP) :: gciold(ijkmax), gcinew(ijkmax)
3221 real(RP) :: cndmss, sblmss
3222 real(RP) :: gtliq(ijkmax), gtice(ijkmax)
3223 real(RP) :: umax(ijkmax), uval, dtcnd(ijkmax)
3224 real(RP) :: qlevp(ijkmax), qlsbl(ijkmax)
3225 real(RP) :: cef1, cef2, cef3, cef4, a, b, c, d
3226 real(RP) :: rmdplus, rmdmins, ssplus, ssmins, tplus, tmins
3227 real(RP) :: sliqtnd, sicetnd
3228 real(RP) :: ssliq, ssice
3229 real(RP) :: sumliq(ijkmax), sumice(ijkmax)
3230 real(RP) :: gcn( -1:
nbin+2,nspc,ijkmax )
3231 real(RP),
parameter :: cflfct = 0.50_rp
3232 real(RP) :: dumm_regene(ijkmax)
3233 real(RP) :: csum( nspc,ijkmax )
3234 integer :: iflg( nspc,ijkmax )
3235 real(RP) :: f1, f2, emu, cefd, cefk, festl, festi
3236 real(RP) :: qsatl(ijkmax), qsati(ijkmax)
3237 real(RP),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
3238 real(RP),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
3239 real(RP) :: qvtmp, zerosw
3240 integer :: ijk, nn, loopflg(ijkmax)
3244 real(RP) :: uadv ( 0:
nbin+2,nspc,ijkmax )
3245 real(RP) :: flq ( 1:
nbin+1,nspc,ijkmax )
3246 real(RP) :: acoef( 0:2,0:
nbin+1,nspc,ijkmax )
3247 real(RP) :: crn ( 0:
nbin+2,nspc,ijkmax )
3248 real(RP) :: aip ( 0:
nbin+1,nspc,ijkmax )
3249 real(RP) :: aim ( 0:
nbin+1,nspc,ijkmax )
3250 real(RP) :: ai ( 0:
nbin+1,nspc,ijkmax )
3251 real(RP) :: cmins, cplus
3256 dumm_regene( : ) = 0.0_rp
3258 csum( :,: ) = 0.0_rp
3263 csum( myu,ijk ) = csum( myu,ijk )+gc( n,myu,ijk )*dxmic
3268 if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
3273 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3275 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3277 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
3278 temp(:), pres(:), qdry(:), &
3280 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3281 temp(:), pres(:), qdry(:), &
3290 gclold(ijk) = gclold(ijk) + gc( n,il,ijk )*dxmic
3295 gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
3302 gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
3304 gcn( -1,myu,ijk ) = 0.0_rp
3305 gcn( 0,myu,ijk ) = 0.0_rp
3306 gcn(
nbin+1,myu,ijk ) = 0.0_rp
3307 gcn(
nbin+2,myu,ijk ) = 0.0_rp
3311 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3312 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3314 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3315 cefd = emu/dens(ijk)
3318 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3319 f1 = rvap*temp(ijk)/festl/cefd
3320 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3321 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3323 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3324 f1 = rvap*temp(ijk)/festi/cefd
3325 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3326 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3329 umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq )
3331 uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice )
3332 umax(ijk) = max( umax(ijk),uval )
3335 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
3336 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
3337 dtcnd(ijk) = dtime/nloop(ijk)
3339 nloop(ijk) = nloop(ijk) * iflg( il,ijk )
3340 nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) )
3343 nloopmax = maxval(nloop,1)
3348 do ncount = 1, nloopmax
3351 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
3356 do nn = 1, loopflg(ijk)
3359 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3361 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3364 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3365 cefd = emu/dens(ijk)
3367 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3368 f1 = rvap*temp(ijk)/festl/cefd
3369 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3370 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3371 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3372 f1 = rvap*temp(ijk)/festi/cefd
3373 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3374 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3376 sumliq(ijk) = 0.0_rp
3378 sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
3381 sumice(ijk) = 0.0_rp
3384 sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
3391 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
3392 temp(:), pres(:), qdry(:), &
3394 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3395 temp(:), pres(:), qdry(:), &
3399 do nn = 1, loopflg(ijk)
3402 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3403 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3405 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
3406 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
3407 cef1 = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cp(ijk) )
3408 cef2 = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cp(ijk) )
3409 cef3 = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cp(ijk) )
3410 cef4 = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cp(ijk) )
3412 a = - cef1*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3413 b = - cef2*sumice(ijk)*gtice(ijk)/dens(ijk)
3414 c = - cef3*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3415 d = - cef4*sumice(ijk)*gtice(ijk)/dens(ijk)
3417 b = b + eps * ( 1.0_rp - zerosw )
3419 rmdplus = ( ( a+d ) + sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
3420 rmdmins = ( ( a+d ) - sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
3422 rmdplus = rmdplus + eps * ( 1.0_rp - zerosw )
3423 rmdmins = rmdmins + eps * ( 1.0_rp - zerosw )
3426 ssplus = ( ( rmdmins-a )*ssliq - b*ssice )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
3427 ssmins = ( ( a-rmdplus )*ssliq + b*ssice )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
3429 tplus = ( exp( rmdplus*dtcnd(ijk) )-1.0_rp )/( rmdplus*dtcnd(ijk) ) &
3430 * ( 0.5_rp + sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) ) &
3432 * ( 0.5_rp - sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) )
3433 tmins = ( exp( rmdmins*dtcnd(ijk) )-1.0_rp )/( rmdmins*dtcnd(ijk) ) &
3434 * ( 0.5_rp + sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) ) &
3436 * ( 0.5_rp - sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) )
3438 sliqtnd = ssliq * ( 1.0_rp - zerosw ) &
3440 ( b*tplus*ssplus + b*tmins*ssmins )
3441 sicetnd = ssice * ( 1.0_rp - zerosw ) &
3443 ( ( rmdplus-a )*tplus*ssplus &
3444 + ( rmdmins-a )*tmins*ssmins )
3448 do mm = 1, iflg( myu,ijk )
3452 uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtliq(ijk)*sliqtnd &
3453 * ( 0.5_rp - sign( 0.5_rp,
real(myu)-1.5_RP) ) &
3454 + cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd &
3455 * ( 0.5_rp + sign( 0.5_rp,
real(myu)-1.5_rp) )
3457 uadv( 0, myu,ijk ) = 0.0_rp
3458 uadv(
nbin+2,myu,ijk ) = 0.0_rp
3467 do nn = 1, loopflg(ijk)
3469 do mm = 1, iflg( myu,ijk )
3471 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
3479 do nn = 1, loopflg(ijk)
3481 do mm = 1, iflg( myu,ijk )
3483 acoef(0,n,myu,ijk) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3484 acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3485 acoef(2,n,myu,ijk) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3493 do nn = 1, loopflg(ijk)
3495 do mm = 1, iflg( myu,ijk )
3497 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3499 aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3500 + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3501 + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3503 aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
3511 do nn = 1, loopflg(ijk)
3513 do mm = 1, iflg( myu,ijk )
3515 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3517 aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3518 - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3519 + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3521 aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
3529 do nn = 1, loopflg(ijk)
3531 do mm = 1, iflg( myu,ijk )
3533 ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
3534 + acoef(2,n,myu,ijk) * 2.0_rp
3536 ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
3544 do nn = 1, loopflg(ijk)
3546 do mm = 1, iflg( myu,ijk )
3548 flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
3549 - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3557 do nn = 1, loopflg(ijk)
3559 do mm = 1, iflg( myu,ijk )
3560 dumm_regene(ijk) = dumm_regene(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
3561 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk)+eps )
3568 do nn = 1, loopflg(ijk)
3570 do mm = 1, iflg( myu,ijk )
3572 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
3583 do nn = 1, loopflg(ijk)
3585 gclnew(ijk) = 0.0_rp
3587 gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )*dxmic
3590 gcinew(ijk) = 0.0_rp
3593 gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
3598 cndmss = gclnew(ijk) - gclold(ijk)
3599 sblmss = gcinew(ijk) - gciold(ijk)
3601 qvap(ijk) = qvap(ijk) - ( cndmss + sblmss ) / dens(ijk)
3602 temp(ijk) = temp(ijk) + ( cndmss*qlevp(ijk)+sblmss*qlsbl(ijk) ) / dens(ijk) / cp(ijk)
3606 gclold(ijk) = gclnew(ijk)
3607 gciold(ijk) = gcinew(ijk)
3618 gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
3623 if ( .NOT. flg_regeneration )
then 3624 dumm_regene(:) = 0.0_rp
3630 end subroutine mixphase
3633 subroutine ice_nucleat( &
3649 atmos_saturation_pres2qsat_ice
3657 integer,
intent(in) :: ijkmax
3658 integer,
intent(in) :: num_cold
3659 integer,
intent(in) :: index_cold(ijkmax)
3660 real(RP),
intent(in) :: dens(ijkmax)
3661 real(RP),
intent(in) :: pres(ijkmax)
3662 real(RP),
intent(in) :: qdry(ijkmax)
3663 real(RP),
intent(inout) :: temp(ijkmax)
3664 real(RP),
intent(inout) :: qvap(ijkmax)
3665 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3666 real(RP),
intent(inout) :: cp(ijkmax)
3667 real(RP),
intent(inout) :: cv(ijkmax)
3668 real(DP),
intent(in) :: dtime
3671 real(RP) :: numin, tdel, qdel
3673 real(RP),
parameter :: acoef = -0.639_rp, bcoef = 12.96_rp
3675 real(RP),
parameter :: tcolmu = 269.0_rp, tcolml = 265.0_rp
3676 real(RP),
parameter :: tdendu = 257.0_rp, tdendl = 255.0_rp
3677 real(RP),
parameter :: tplatu = 250.6_rp
3679 real(RP) :: qsati(ijkmax)
3680 integer :: ijk, indirect
3685 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3686 temp(:), pres(:), qdry(:), &
3688 do indirect = 1, num_cold
3689 ijk = index_cold(indirect)
3692 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3694 if( ssice <= 0.0_rp ) cycle
3696 numin = bcoef * exp( acoef + bcoef * ssice )
3697 numin = numin * expxctr( 1 )/dxmic
3698 numin = min( numin,qvap(ijk)*dens(ijk) )
3700 if ( temp(ijk) <= tplatu .OR. ( temp(ijk) >= tcolml .AND. temp(ijk) < tcolmu ) )
then 3701 gc( 1,ic,ijk ) = gc( 1,ic,ijk ) + numin
3703 elseif( temp(ijk) <= tdendu .AND. temp(ijk) >= tdendl )
then 3704 gc( 1,id,ijk ) = gc( 1,id,ijk ) + numin
3707 gc( 1,ip,ijk ) = gc( 1,ip,ijk ) + numin
3710 qdel = numin/dens(ijk)
3711 qvap(ijk) = qvap(ijk) - qdel
3712 tdel = numin/dens(ijk)*qlmlt/cp(ijk)
3713 temp(ijk) = temp(ijk) + tdel
3723 end subroutine ice_nucleat
3726 subroutine freezing( &
3748 integer,
intent(in) :: ijkmax
3749 integer,
intent(in) :: num_cold
3750 integer,
intent(in) :: index_cold(ijkmax)
3751 real(RP),
intent(in) :: dens(ijkmax)
3752 real(RP),
intent(inout) :: temp(ijkmax)
3753 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3754 real(RP),
intent(inout) :: cp(ijkmax)
3755 real(RP),
intent(inout) :: cv(ijkmax)
3756 real(DP),
intent(in) :: dtime
3758 integer :: nbound, n
3759 real(RP) :: xbound, tc, rate, dmp, frz, sumfrz, tdel
3760 real(RP),
parameter :: coefa = 1.0e-01_rp
3761 real(RP),
parameter :: coefb = 0.66_rp
3762 real(RP),
parameter :: rbound = 2.0e-04_rp
3766 integer :: ijk, indirect
3770 do indirect = 1, num_cold
3771 ijk = index_cold(indirect)
3774 xbound = log( rhow * 4.0_rp*pi/3.0_rp * rbound**3 )
3775 nbound = int( ( xbound-xbnd( 1 ) )/dxmic ) + 1
3778 rate = coefa*exp( -coefb*tc )
3782 dmp = rate*expxctr( n )
3783 frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
3784 frz = min( frz, gc( n,il,ijk ) )
3786 gc( n,il,ijk ) = gc( n,il,ijk ) - frz
3787 gc( n,ip,ijk ) = gc( n,ip,ijk ) + frz
3789 sumfrz = sumfrz + frz
3792 dmp = rate*expxctr( n )
3793 frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
3794 frz = min( frz, gc( n,il,ijk) )
3796 gc( n,il,ijk ) = gc( n,il,ijk ) - frz
3797 gc( n,ih,ijk ) = gc( n,ih,ijk ) + frz
3799 sumfrz = sumfrz + frz
3819 sumfrz = sumfrz*dxmic
3821 tdel = sumfrz/dens(ijk)*qlmlt/cp(ijk)
3822 temp(ijk) = temp(ijk) + tdel
3830 end subroutine freezing
3833 subroutine melting( &
3852 integer,
intent(in) :: ijkmax
3853 integer,
intent(in) :: num_warm
3854 integer,
intent(in) :: index_warm(ijkmax)
3855 real(RP),
intent(in) :: dens(ijkmax)
3856 real(RP),
intent(inout) :: temp(ijkmax)
3857 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3858 real(RP),
intent(inout) :: cp(ijkmax)
3859 real(RP),
intent(inout) :: cv(ijkmax)
3860 real(DP),
intent(in) :: dtime
3863 real(RP) :: summlt, sumice, tdel
3864 integer :: ijk, indirect
3868 do indirect = 1, num_warm
3869 ijk = index_warm(indirect)
3875 sumice = sumice + gc( n,m,ijk )
3876 gc( n,m,ijk ) = 0.0_rp
3878 gc( n,il,ijk ) = gc( n,il,ijk ) + sumice
3879 summlt = summlt + sumice
3881 summlt = summlt*dxmic
3883 tdel = - summlt/dens(ijk)*qlmlt/cp(ijk)
3884 temp(ijk) = temp(ijk) + tdel
3893 end subroutine melting
3896 subroutine collmain( &
3904 integer,
intent(in) :: KA
3905 integer,
intent(in) :: IA
3906 integer,
intent(in) :: JA
3908 integer,
intent(in) :: ijkmax
3909 real(RP),
intent(in) :: temp(ijkmax)
3910 real(RP),
intent(inout) :: ghyd(
nbin,nspc,ijkmax)
3911 real(DP),
intent(in) :: dt
3914 if ( rndm_flgp == 1 )
then 3915 call r_collcoag( ka, ia, ja, &
3922 call collcoag( ijkmax, &
3929 end subroutine collmain
3932 subroutine collmainf( &
3940 integer,
intent(in) :: KA
3941 integer,
intent(in) :: IA
3942 integer,
intent(in) :: JA
3944 integer,
intent(in) :: ijkmax
3945 real(RP),
intent(in) :: temp(ijkmax)
3946 real(RP),
intent(inout) :: ghyd(
nbin,nspc,ijkmax)
3947 real(DP),
intent(in) :: dt
3950 if ( rndm_flgp == 1 )
then 3951 call r_collcoag( ka, ia, ja, &
3958 call collcoag( ijkmax, &
3965 end subroutine collmainf
3970 subroutine collcoag( &
3977 integer,
intent(in) :: ijkmax
3978 real(RP),
intent(in) :: temp(ijkmax)
3979 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3980 real(DP),
intent(in) :: dtime
3982 integer :: i, j, k, l
3983 real(RP) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
3984 real(RP) :: gprime, gprimk, wgt, crn, sum, flux
3985 integer,
parameter :: ldeg = 2
3986 real(RP) :: acoef( 0:ldeg )
3987 real(RP),
parameter :: dmpmin = 1.e-01_rp
3988 real(RP) :: suri, surj
3990 integer :: myu, n, irsl, ilrg, isml
3991 integer :: ibnd( ijkmax )
3992 integer :: iflg( nspc,ijkmax )
3993 integer :: iexst(
nbin,nspc,ijkmax )
3994 real(RP) :: csum( nspc,ijkmax )
3995 integer :: ijk, nn, mm, pp, qq
4002 csum( :,: ) = 0.0_rp
4007 csum( myu,ijk ) = csum( myu,ijk ) + gc( n,myu,ijk )*dxmic
4011 if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
4014 if ( temp(ijk) < tcrit )
then 4022 if ( gc( n,myu,ijk ) > cldmin )
then 4023 iexst( n,myu,ijk ) = 1
4033 do nn = 1, iflg( isml,ijk )
4036 do mm = 1, iflg( ilrg,ijk )
4038 irsl = ifrsl( ibnd(ijk),isml,ilrg )
4041 do pp = 1, iexst( i,isml,ijk )
4044 do qq = 1, iexst( j,ilrg,ijk )
4051 dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
4052 dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
4054 if ( dmpi <= dmpmin )
then 4055 frci = gc( i,isml,ijk )*dmpi
4057 frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
4060 if ( dmpj <= dmpmin )
then 4061 frcj = gc( j,ilrg,ijk )*dmpj
4063 frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
4068 if ( gprime > 0.0_rp .AND. k <
nbin )
then 4070 suri = gc( i,isml,ijk )
4071 surj = gc( j,ilrg,ijk )
4072 gc( i,isml,ijk ) = gc( i,isml,ijk )-frci
4073 gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj
4074 gc( i,isml,ijk ) = max( gc( i,isml,ijk )-frci, 0.0_rp )
4075 gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk )-frcj, 0.0_rp )
4076 frci = suri - gc( i,isml,ijk )
4077 frcj = surj - gc( j,ilrg,ijk )
4080 gprimk = gc( k,irsl,ijk ) + gprime
4081 wgt = gprime / gprimk
4082 crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
4084 acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
4085 acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.50_rp
4086 acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
4090 sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
4091 *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
4095 flux = min( max( flux,0.0_rp ),gprime )
4097 gc( k,irsl,ijk ) = gprimk - flux
4098 gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
4118 end subroutine collcoag
4121 subroutine getrule( ifrsl, indx )
4123 integer,
intent(out) :: indx(
nbin,
nbin )
4124 integer,
intent(out) :: ifrsl( 2,nspc_mk,nspc_mk )
4130 xnew = log( expxctr( i )+expxctr( j ) )
4131 k = int( ( xnew-xctr( 1 ) )/dxmic ) + 1
4132 k = max( max( k,j ),i )
4139 ifrsl( 1:2,il,il ) = il
4142 ifrsl( 1:2,il,ic ) = ic
4143 ifrsl( 1,ic,il ) = ig
4144 ifrsl( 2,ic,il ) = ih
4147 ifrsl( 1:2,il,ip ) = ip
4148 ifrsl( 1,ip,il ) = ig
4149 ifrsl( 2,ip,il ) = ih
4152 ifrsl( 1:2,il,id ) = id
4153 ifrsl( 1,id,il ) = ig
4154 ifrsl( 2,id,il ) = ih
4157 ifrsl( 1:2,il,iss ) = iss
4158 ifrsl( 1,iss,il ) = ig
4159 ifrsl( 2,iss,il ) = ih
4162 ifrsl( 1:2,il,ig ) = ig
4163 ifrsl( 1,ig,il ) = ig
4164 ifrsl( 2,ig,il ) = ih
4167 ifrsl( 1:2,il,ih ) = ih
4168 ifrsl( 1,ih,il ) = ig
4169 ifrsl( 2,ih,il ) = ih
4173 ifrsl( 1:2,ic,ic ) = iss
4176 ifrsl( 1:2,ic,ip ) = iss
4177 ifrsl( 1:2,ip,ic ) = iss
4180 ifrsl( 1:2,ic,id ) = iss
4181 ifrsl( 1:2,id,ic ) = iss
4184 ifrsl( 1:2,ic,iss ) = iss
4185 ifrsl( 1:2,iss,ic ) = iss
4188 ifrsl( 1:2,ic,ig ) = ig
4189 ifrsl( 1:2,ig,ic ) = ic
4192 ifrsl( 1:2,ih,ic ) = ic
4193 ifrsl( 1,ic,ih ) = ig
4194 ifrsl( 2,ic,ih ) = ih
4198 ifrsl( 1:2,ip,ip ) = iss
4201 ifrsl( 1:2,ip,id ) = iss
4202 ifrsl( 1:2,id,ip ) = iss
4205 ifrsl( 1:2,ip,iss ) = iss
4206 ifrsl( 1:2,iss,ip ) = iss
4209 ifrsl( 1:2,ip,ig ) = ig
4210 ifrsl( 1:2,ig,ip ) = ip
4213 ifrsl( 1:2,ih,ip ) = ip
4214 ifrsl( 1,ip,ih ) = ig
4215 ifrsl( 2,ip,ih ) = ih
4219 ifrsl( 1:2,id,id ) = iss
4222 ifrsl( 1:2,id,iss ) = iss
4223 ifrsl( 1:2,iss,id ) = iss
4226 ifrsl( 1:2,id,ig ) = ig
4227 ifrsl( 1:2,ig,id ) = id
4230 ifrsl( 1:2,ih,id ) = id
4231 ifrsl( 1,id,ih ) = ig
4232 ifrsl( 2,id,ih ) = ih
4236 ifrsl( 1:2,iss,iss ) = iss
4239 ifrsl( 1:2,iss,ig ) = ig
4240 ifrsl( 1:2,ig,iss ) = iss
4243 ifrsl( 1:2,ih,iss ) = iss
4244 ifrsl( 1,iss,ih ) = ig
4245 ifrsl( 2,iss,ih ) = ih
4249 ifrsl( 1:2,ig,ig ) = ig
4252 ifrsl( 1,ig,ih ) = ig
4253 ifrsl( 1,ih,ig ) = ig
4254 ifrsl( 2,ig,ih ) = ih
4255 ifrsl( 2,ih,ig ) = ih
4258 ifrsl( 1:2,ih,ih ) = ih
4261 end subroutine getrule
4270 integer ,
intent(in) :: ijkmax
4271 real(RP),
intent(in) :: f0(ijkmax)
4272 real(RP),
intent(inout) :: ga(
nccn,ijkmax)
4283 ga(n,ijk) = ga(n,ijk) + f0(ijk) * marate(n) * expxactr(n) / dxaer
4290 end subroutine faero
4295 subroutine random_setup( mset )
4301 integer,
intent(in) :: mset
4305 real(RP) :: nbinr, tmp1
4306 real(RP) :: rans( mbin ), ranl( mbin )
4308 real(RP),
allocatable :: ranstmp( : )
4309 real(RP),
allocatable :: ranltmp( : )
4312 integer,
allocatable :: orderb( : )
4315 real(RP),
allocatable :: randnum(:,:,:)
4318 allocate( blrg( mset, mbin ) )
4319 allocate( bsml( mset, mbin ) )
4320 allocate( ranstmp( pq ) )
4321 allocate( ranltmp( pq ) )
4322 allocate( orderb( pq ) )
4323 allocate( randnum(1,1,pq) )
4325 a =
real(
nbin )*
real(
nbin-1 )*0.50_RP
4326 if ( a < mbin )
then 4327 log_error(
"ATMOS_PHY_MP_SUZUKI10_random_setup",*)
"mbin should be smaller than {nbin}_C_{2}" 4331 wgtbin = a/
real( mbin )
4332 nbinr =
real(
nbin )
4339 ranstmp( (p-1)*
nbin-(p*(p-1))/2+1 : p*
nbin-(p*(p+1))/2 ) = p
4341 ranltmp( (p-1)*
nbin-(p*(p-1))/2+q ) = p+q
4348 abq1 = randnum( 1,1,p )
4349 k = int( abq1*( pq-p-1 ) ) + p
4351 orderb( p ) = orderb( k )
4357 rans( p ) = ranstmp( orderb( p ) )
4358 ranl( p ) = ranltmp( orderb( p ) )
4360 rans( p ) = ranstmp( orderb( p-pq ) )
4361 ranl( p ) = ranltmp( orderb( p-pq ) )
4363 if ( rans( p ) >= ranl( p ) )
then 4365 rans( p ) = ranl( p )
4369 blrg( n,1:mbin ) = int( ranl( 1:mbin ) )
4370 bsml( n,1:mbin ) = int( rans( 1:mbin ) )
4373 deallocate( ranstmp )
4374 deallocate( ranltmp )
4375 deallocate( orderb )
4376 deallocate( randnum )
4378 end subroutine random_setup
4384 subroutine r_collcoag( &
4395 integer,
intent(in) :: KA
4396 integer,
intent(in) :: IA
4397 integer,
intent(in) :: JA
4399 integer,
intent(in) :: ijkmax
4400 real(RP),
intent(in) :: swgt
4401 real(RP),
intent(in) :: temp(ijkmax)
4402 real(RP),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4403 real(DP),
intent(in) :: dtime
4405 integer :: i, j, k, l
4406 real(RP) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
4407 real(RP) :: gprime, gprimk, wgt, crn, sum, flux
4408 integer,
parameter :: ldeg = 2
4409 real(RP),
parameter :: dmpmin = 1.e-01_rp, cmin = 1.e-10_rp
4410 real(RP) :: acoef( 0:ldeg )
4413 integer :: nums( mbin ), numl( mbin )
4414 real(RP),
parameter :: gt = 1.0_rp
4416 real(RP) :: nbinr, mbinr
4418 real(RP) :: tmpi, tmpj
4420 integer :: ibnd( ijkmax )
4421 integer :: iflg( nspc,ijkmax )
4422 integer :: iexst(
nbin,nspc,ijkmax )
4423 real(RP) :: csum( nspc,ijkmax )
4424 integer :: ijk, nn, mm, pp, qq, myu, n, isml, ilrg, irsl
4431 csum( :,: ) = 0.0_rp
4435 csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
4436 csum( ic,ijk ) = csum( ic,ijk ) + gc( n,ic,ijk )*dxmic
4437 csum( ip,ijk ) = csum( ip,ijk ) + gc( n,ip,ijk )*dxmic
4438 csum( id,ijk ) = csum( id,ijk ) + gc( n,id,ijk )*dxmic
4439 csum( iss,ijk ) = csum( iss,ijk ) + gc( n,iss,ijk )*dxmic
4440 csum( ig,ijk ) = csum( ig,ijk ) + gc( n,ig,ijk )*dxmic
4441 csum( ih,ijk ) = csum( ih,ijk ) + gc( n,ih,ijk )*dxmic
4443 if ( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
4444 if ( csum( ic,ijk ) > cldmin ) iflg( ic,ijk ) = 1
4445 if ( csum( ip,ijk ) > cldmin ) iflg( ip,ijk ) = 1
4446 if ( csum( id,ijk ) > cldmin ) iflg( id,ijk ) = 1
4447 if ( csum( iss,ijk ) > cldmin ) iflg( iss,ijk ) = 1
4448 if ( csum( ig,ijk ) > cldmin ) iflg( ig,ijk ) = 1
4449 if ( csum( ih,ijk ) > cldmin ) iflg( ih,ijk ) = 1
4451 if ( temp(ijk) < tcrit )
then 4459 if ( gc( n,myu,ijk ) > cldmin )
then 4460 iexst( n,myu,ijk ) = 1
4470 do nn = 1, iflg( isml,ijk )
4473 do mm = 1, iflg( ilrg,ijk )
4475 irsl = ifrsl( ibnd(ijk),isml,ilrg )
4478 det = int( rndm(1,1,1)*ia*ja*ka )
4479 nbinr =
real(
nbin )
4480 mbinr =
real( mbin )
4481 nums( 1:mbin ) = bsml( det,1:mbin )
4482 numl( 1:mbin ) = blrg( det,1:mbin )
4488 do pp = 1, iexst( i,isml,ijk )
4489 do qq = 1, iexst( j,ilrg,ijk )
4496 dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
4497 dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
4499 if ( dmpi <= dmpmin )
then 4500 frci = gc( i,isml,ijk )*dmpi
4502 frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
4505 if ( dmpj <= dmpmin )
then 4506 frcj = gc( j,ilrg,ijk )*dmpj
4508 frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
4510 tmpi = gc( i,isml,ijk )
4511 tmpj = gc( j,ilrg,ijk )
4513 gc( i,isml,ijk ) = gc( i,isml,ijk )-frci*swgt
4514 gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj*swgt
4517 gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk ), 0.0_rp )
4519 gc( i,isml,ijk ) = max( gc( i,isml,ijk ), 0.0_rp )
4521 frci = tmpi - gc( i,isml,ijk )
4522 frcj = tmpj - gc( j,ilrg,ijk )
4547 if ( gprime > 0.0_rp .AND. k <
nbin )
then 4549 gprimk = gc( k,irsl,ijk ) + gprime
4550 wgt = gprime / gprimk
4551 crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
4553 acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
4554 acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.5_rp
4555 acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
4559 sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
4560 *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
4564 flux = min( max( flux,0.0_rp ),gprime )
4566 gc( k,irsl,ijk ) = gprimk - flux
4567 gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
4586 end subroutine r_collcoag
4599 allocate( radc_mk(
nbin ) )
4600 allocate( xctr_mk(
nbin ) )
4601 allocate( xbnd_mk(
nbin+1 ) )
4602 allocate( cctr_mk( nspc_mk,
nbin ) )
4603 allocate( cbnd_mk( nspc_mk,
nbin+1 ) )
4604 allocate( ck_mk( nspc_mk,nspc_mk,
nbin,
nbin ) )
4605 allocate( vt_mk( nspc_mk,
nbin ) )
4606 allocate( br_mk( nspc_mk,
nbin ) )
4629 deallocate( radc_mk )
4630 deallocate( xctr_mk )
4631 deallocate( xbnd_mk )
4632 deallocate( cctr_mk )
4633 deallocate( cbnd_mk )
4638 end subroutine mkpara
4644 integer,
parameter :: il = 1, ic = 2, ip = 3, id = 4
4645 integer,
parameter :: is = 5, ig = 6, ih = 7
4646 integer,
parameter :: icemax = 3
4648 integer :: k, kk, i, j
4649 real(DP) :: xl( ndat ), rlec( ndat ), vrl( ndat )
4650 real(DP) :: blkradl( ndat ), blkdnsl( ndat )
4652 real(DP) :: xi( ndat,icemax ), riec( ndat,icemax ), vri( ndat,icemax )
4653 real(DP) :: blkradi( ndat,icemax ), blkdnsi( ndat,icemax )
4655 real(DP) :: xs( ndat ), rsec( ndat ), vrs( ndat )
4656 real(DP) :: blkrads( ndat ), blkdnss( ndat )
4658 real(DP) :: xg( ndat ), rgec( ndat ), vrg( ndat )
4659 real(DP) :: blkradg( ndat ), blkdnsg( ndat )
4661 real(DP) :: xh( ndat ), rhec( ndat ), vrh( ndat )
4662 real(DP) :: blkradh( ndat ), blkdnsh( ndat )
4665 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4666 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4667 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4668 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4669 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4670 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4671 data xi(1:ndat,1) / &
4672 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4673 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4674 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4675 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4676 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4677 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4678 data xi(1:ndat,2) / &
4679 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4680 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4681 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4682 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4683 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4684 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00 /
4685 data xi(1:ndat,3) / &
4686 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4687 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4688 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4689 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4690 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4691 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4693 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4694 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4695 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4696 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4697 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4698 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4700 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4701 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4702 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4703 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4704 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4705 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4707 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4708 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4709 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4710 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4711 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4712 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4713 data rlec(1:ndat) / &
4714 0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
4715 0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
4716 0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
4717 0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
4718 0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
4719 0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
4720 data riec(1:ndat,1) / &
4721 0.31936e-03_dp,0.40397e-03_dp,0.51099e-03_dp,0.64638e-03_dp,0.81764e-03_dp,0.10343e-02_dp, &
4722 0.13084e-02_dp,0.16551e-02_dp,0.20937e-02_dp,0.26486e-02_dp,0.33506e-02_dp,0.42387e-02_dp, &
4723 0.64360e-02_dp,0.81426e-02_dp,0.10302e-01_dp,0.13035e-01_dp,0.16494e-01_dp,0.20872e-01_dp, &
4724 0.26412e-01_dp,0.33426e-01_dp,0.42304e-01_dp,0.53543e-01_dp,0.67770e-01_dp,0.85783e-01_dp, &
4725 0.10859e+00_dp,0.13746e+00_dp,0.17403e+00_dp,0.22032e+00_dp,0.27895e+00_dp,0.35319e+00_dp, &
4726 0.44722e+00_dp,0.56630e+00_dp,0.71712e+00_dp /
4727 data riec(1:ndat,2) / &
4728 0.13188e-03_dp,0.16615e-03_dp,0.20953e-03_dp,0.27728e-03_dp,0.36694e-03_dp,0.48559e-03_dp, &
4729 0.64261e-03_dp,0.85040e-03_dp,0.11254e-02_dp,0.14893e-02_dp,0.19709e-02_dp,0.26082e-02_dp, &
4730 0.34515e-02_dp,0.45676e-02_dp,0.60446e-02_dp,0.79991e-02_dp,0.10586e-01_dp,0.14009e-01_dp, &
4731 0.18539e-01_dp,0.24533e-01_dp,0.32466e-01_dp,0.42964e-01_dp,0.56857e-01_dp,0.75242e-01_dp, &
4732 0.99573e-01_dp,0.13177e+00_dp,0.17438e+00_dp,0.23077e+00_dp,0.30539e+00_dp,0.40414e+00_dp, &
4733 0.53482e+00_dp,0.70775e+00_dp,0.93661e+00_dp /
4734 data riec(1:ndat,3) / 0.14998e-03_dp,0.18896e-03_dp,0.23808e-03_dp, &
4735 0.29996e-03_dp,0.37793e-03_dp,0.47616e-03_dp,0.61048e-03_dp,0.81343e-03_dp,0.10839e-02_dp, &
4736 0.14442e-02_dp,0.19243e-02_dp,0.25640e-02_dp,0.34164e-02_dp,0.45522e-02_dp,0.60656e-02_dp, &
4737 0.80820e-02_dp,0.10769e-01_dp,0.14349e-01_dp,0.19119e-01_dp,0.25475e-01_dp,0.44576e-01_dp, &
4738 0.62633e-01_dp,0.88006e-01_dp,0.12366e+00_dp,0.17375e+00_dp,0.24414e+00_dp,0.34304e+00_dp, &
4739 0.48201e+00_dp,0.67728e+00_dp,0.95164e+00_dp,0.13372e+01_dp,0.18788e+01_dp,0.26400e+01_dp /
4740 data rsec(1:ndat) / &
4741 0.92832e-03_dp,0.11696e-02_dp,0.14736e-02_dp,0.18566e-02_dp,0.23392e-02_dp,0.29472e-02_dp, &
4742 0.37133e-02_dp,0.46784e-02_dp,0.58944e-02_dp,0.74265e-02_dp,0.93569e-02_dp,0.11789e-01_dp, &
4743 0.14853e-01_dp,0.18714e-01_dp,0.23578e-01_dp,0.29706e-01_dp,0.37427e-01_dp,0.47156e-01_dp, &
4744 0.59412e-01_dp,0.74855e-01_dp,0.94311e-01_dp,0.11882e+00_dp,0.14971e+00_dp,0.18862e+00_dp, &
4745 0.23765e+00_dp,0.29942e+00_dp,0.37724e+00_dp,0.47530e+00_dp,0.59884e+00_dp,0.75449e+00_dp, &
4746 0.95060e+00_dp,0.11977e+01_dp,0.15090e+01_dp /
4747 data rgec(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
4748 0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
4749 0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
4750 0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
4751 0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
4752 0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
4753 data rhec(1:ndat) / &
4754 0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
4755 0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
4756 0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
4757 0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
4758 0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
4759 0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
4760 data vrl(1:ndat) / &
4761 0.50000e-01_dp,0.78000e-01_dp,0.12000e+00_dp,0.19000e+00_dp,0.31000e+00_dp,0.49000e+00_dp, &
4762 0.77000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp,0.48000e+01_dp,0.74000e+01_dp, &
4763 0.11000e+02_dp,0.17000e+02_dp,0.26000e+02_dp,0.37000e+02_dp,0.52000e+02_dp,0.71000e+02_dp, &
4764 0.94000e+02_dp,0.12000e+03_dp,0.16000e+03_dp,0.21000e+03_dp,0.26000e+03_dp,0.33000e+03_dp, &
4765 0.41000e+03_dp,0.48000e+03_dp,0.57000e+03_dp,0.66000e+03_dp,0.75000e+03_dp,0.82000e+03_dp, &
4766 0.88000e+03_dp,0.90000e+03_dp,0.90000e+03_dp /
4767 data vri(1:ndat,1) / 0.30000e-01_dp,0.40000e-01_dp,0.60000e-01_dp, &
4768 0.80000e-01_dp,0.11000e+00_dp,0.15000e+00_dp,0.17000e+00_dp,0.18000e+00_dp,0.20000e+00_dp, &
4769 0.25000e+00_dp,0.40000e+00_dp,0.60000e+01_dp,0.10000e+02_dp,0.15000e+02_dp,0.20000e+02_dp, &
4770 0.25000e+02_dp,0.31000e+02_dp,0.37000e+02_dp,0.41000e+02_dp,0.46000e+02_dp,0.51000e+02_dp, &
4771 0.55000e+02_dp,0.59000e+02_dp,0.62000e+02_dp,0.64000e+02_dp,0.67000e+02_dp,0.68000e+02_dp, &
4772 0.69000e+02_dp,0.70000e+02_dp,0.71000e+02_dp,0.71500e+02_dp,0.71750e+02_dp,0.72000e+02_dp /
4773 data vri(1:ndat,2) / &
4774 0.30000e-01_dp,0.40000e-01_dp,0.50000e-01_dp,0.70000e-01_dp,0.90000e-01_dp,0.12000e+00_dp, &
4775 0.50000e+00_dp,0.80000e+00_dp,0.16000e+01_dp,0.18000e+01_dp,0.20000e+01_dp,0.30000e+01_dp, &
4776 0.40000e+01_dp,0.50000e+01_dp,0.80000e+01_dp,0.13000e+02_dp,0.19000e+02_dp,0.26000e+02_dp, &
4777 0.32000e+02_dp,0.38000e+02_dp,0.47000e+02_dp,0.55000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
4778 0.77000e+02_dp,0.79000e+02_dp,0.80000e+02_dp,0.81000e+02_dp,0.81000e+02_dp,0.82000e+02_dp, &
4779 0.82000e+02_dp,0.82000e+02_dp,0.82000e+02_dp /
4780 data vri(1:ndat,3) / 0.35000e-01_dp,0.45000e-01_dp,0.55000e-01_dp, &
4781 0.75000e-01_dp,0.95000e-01_dp,0.13000e+00_dp,0.60000e+00_dp,0.90000e+00_dp,0.17000e+01_dp, &
4782 0.20000e+01_dp,0.25000e+01_dp,0.38000e+01_dp,0.50000e+01_dp,0.70000e+01_dp,0.90000e+01_dp, &
4783 0.11000e+02_dp,0.14000e+02_dp,0.17000e+02_dp,0.21000e+02_dp,0.25000e+02_dp,0.32000e+02_dp, &
4784 0.38000e+02_dp,0.44000e+02_dp,0.49000e+02_dp,0.53000e+02_dp,0.55000e+02_dp,0.58000e+02_dp, &
4785 0.59000e+02_dp,0.61000e+02_dp,0.62000e+02_dp,0.63000e+02_dp,0.64000e+02_dp,0.65000e+02_dp /
4786 data vrs(1:ndat) / &
4787 0.20000e-01_dp,0.31000e-01_dp,0.49000e-01_dp,0.77000e-01_dp,0.12000e+00_dp,0.19000e+00_dp, &
4788 0.30000e+00_dp,0.48000e+00_dp,0.76000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp, &
4789 0.48000e+01_dp,0.75000e+01_dp,0.11000e+02_dp,0.16000e+02_dp,0.21000e+02_dp,0.26000e+02_dp, &
4790 0.34000e+02_dp,0.41000e+02_dp,0.49000e+02_dp,0.57000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
4791 0.81000e+02_dp,0.87000e+02_dp,0.93000e+02_dp,0.99000e+02_dp,0.10750e+03_dp,0.11500e+03_dp, &
4792 0.12500e+03_dp,0.13500e+03_dp,0.14500e+03_dp /
4793 data vrg(1:ndat) / 0.39000e-01_dp,0.62000e-01_dp,0.97000e-01_dp, &
4794 0.15000e+00_dp,0.24000e+00_dp,0.38000e+00_dp,0.61000e+00_dp,0.96000e+00_dp,0.15000e+01_dp, &
4795 0.24000e+01_dp,0.38000e+01_dp,0.61000e+01_dp,0.96000e+01_dp,0.15000e+02_dp,0.23000e+02_dp, &
4796 0.31000e+02_dp,0.39000e+02_dp,0.49000e+02_dp,0.59000e+02_dp,0.68000e+02_dp,0.79000e+02_dp, &
4797 0.88000e+02_dp,0.10000e+03_dp,0.11000e+03_dp,0.13000e+03_dp,0.15000e+03_dp,0.17000e+03_dp, &
4798 0.20000e+03_dp,0.23000e+03_dp,0.26000e+03_dp,0.30000e+03_dp,0.35000e+03_dp,0.40000e+03_dp /
4799 data vrh(1:ndat) / &
4800 0.53000e-01_dp,0.84000e-01_dp,0.13000e+00_dp,0.21000e+00_dp,0.33000e+00_dp,0.52000e+00_dp, &
4801 0.82000e+00_dp,0.13000e+01_dp,0.21000e+01_dp,0.33000e+01_dp,0.52000e+01_dp,0.82000e+01_dp, &
4802 0.13000e+02_dp,0.20000e+02_dp,0.28000e+02_dp,0.36000e+02_dp,0.46000e+02_dp,0.56000e+02_dp, &
4803 0.67000e+02_dp,0.80000e+02_dp,0.97000e+02_dp,0.12000e+03_dp,0.14000e+03_dp,0.17000e+03_dp, &
4804 0.20000e+03_dp,0.24000e+03_dp,0.29000e+03_dp,0.35000e+03_dp,0.42000e+03_dp,0.51000e+03_dp, &
4805 0.61000e+03_dp,0.74000e+03_dp,0.89000e+03_dp /
4806 data blkradl(1:ndat) / &
4807 0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
4808 0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
4809 0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
4810 0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
4811 0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
4812 0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
4813 data blkradi(1:ndat,1) / 0.57452e-03_dp,0.72384e-03_dp,0.91199e-03_dp, &
4814 0.11490e-02_dp,0.14477e-02_dp,0.18240e-02_dp,0.22981e-02_dp,0.28954e-02_dp,0.36479e-02_dp, &
4815 0.45961e-02_dp,0.57908e-02_dp,0.72959e-02_dp,0.11572e-01_dp,0.14770e-01_dp,0.18851e-01_dp, &
4816 0.24060e-01_dp,0.30709e-01_dp,0.39194e-01_dp,0.50025e-01_dp,0.63848e-01_dp,0.81491e-01_dp, &
4817 0.10401e+00_dp,0.13275e+00_dp,0.16943e+00_dp,0.21625e+00_dp,0.27601e+00_dp,0.35228e+00_dp, &
4818 0.44962e+00_dp,0.57387e+00_dp,0.73244e+00_dp,0.93484e+00_dp,0.11932e+01_dp,0.15229e+01_dp /
4819 data blkradi(1:ndat,2) / &
4820 0.20715e-03_dp,0.26099e-03_dp,0.32912e-03_dp,0.43555e-03_dp,0.57638e-03_dp,0.76276e-03_dp, &
4821 0.10094e-02_dp,0.13358e-02_dp,0.17677e-02_dp,0.23394e-02_dp,0.30958e-02_dp,0.40969e-02_dp, &
4822 0.54216e-02_dp,0.71748e-02_dp,0.94948e-02_dp,0.12565e-01_dp,0.16628e-01_dp,0.22005e-01_dp, &
4823 0.29120e-01_dp,0.38537e-01_dp,0.50998e-01_dp,0.67488e-01_dp,0.89311e-01_dp,0.11819e+00_dp, &
4824 0.15641e+00_dp,0.20698e+00_dp,0.27391e+00_dp,0.36249e+00_dp,0.47970e+00_dp,0.63481e+00_dp, &
4825 0.84009e+00_dp,0.11117e+01_dp,0.14712e+01_dp /
4826 data blkradi(1:ndat,3) / 0.23559e-03_dp,0.29682e-03_dp,0.37397e-03_dp, &
4827 0.47118e-03_dp,0.59365e-03_dp,0.74795e-03_dp,0.95894e-03_dp,0.12777e-02_dp,0.17025e-02_dp, &
4828 0.22685e-02_dp,0.30227e-02_dp,0.40275e-02_dp,0.53665e-02_dp,0.71506e-02_dp,0.95278e-02_dp, &
4829 0.12695e-01_dp,0.16916e-01_dp,0.22539e-01_dp,0.30032e-01_dp,0.40017e-01_dp,0.70019e-01_dp, &
4830 0.98384e-01_dp,0.13824e+00_dp,0.19424e+00_dp,0.27293e+00_dp,0.38350e+00_dp,0.53885e+00_dp, &
4831 0.75714e+00_dp,0.10639e+01_dp,0.14948e+01_dp,0.21004e+01_dp,0.29513e+01_dp,0.41469e+01_dp /
4832 data blkrads(1:ndat) / &
4833 0.20715e-03_dp,0.26148e-03_dp,0.33067e-03_dp,0.41710e-03_dp,0.52691e-03_dp,0.66640e-03_dp, &
4834 0.84289e-03_dp,0.10674e-02_dp,0.13513e-02_dp,0.17129e-02_dp,0.21670e-02_dp,0.27521e-02_dp, &
4835 0.34989e-02_dp,0.44777e-02_dp,0.57347e-02_dp,0.75389e-02_dp,0.99020e-02_dp,0.13161e-01_dp, &
4836 0.17372e-01_dp,0.23337e-01_dp,0.31058e-01_dp,0.41194e-01_dp,0.55153e-01_dp,0.74854e-01_dp, &
4837 0.99806e-01_dp,0.13463e+00_dp,0.18136e+00_dp,0.24282e+00_dp,0.32955e+00_dp,0.44123e+00_dp, &
4838 0.59884e+00_dp,0.77090e+00_dp,0.99387e+00_dp /
4839 data blkradg(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
4840 0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
4841 0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
4842 0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
4843 0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
4844 0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
4845 data blkradh(1:ndat) / &
4846 0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
4847 0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
4848 0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
4849 0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
4850 0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
4851 0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
4852 data blkdnsl(1:ndat) / &
4853 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4854 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4855 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4856 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4857 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4858 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp /
4859 data blkdnsi(1:ndat,1) / 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4860 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4861 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.87368e+00_dp,0.87072e+00_dp,0.86777e+00_dp, &
4862 0.86483e+00_dp,0.86189e+00_dp,0.85897e+00_dp,0.85606e+00_dp,0.85316e+00_dp,0.85026e+00_dp, &
4863 0.84738e+00_dp,0.84451e+00_dp,0.84164e+00_dp,0.83879e+00_dp,0.83595e+00_dp,0.83311e+00_dp, &
4864 0.83029e+00_dp,0.82747e+00_dp,0.82467e+00_dp,0.82187e+00_dp,0.81908e+00_dp,0.81631e+00_dp /
4865 data blkdnsi(1:ndat,2) / &
4866 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4867 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4868 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4869 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4870 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4871 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
4872 data blkdnsi(1:ndat,3) / 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
4873 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
4874 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
4875 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.51790e+00_dp, &
4876 0.45557e+00_dp,0.40075e+00_dp,0.35252e+00_dp,0.31010e+00_dp,0.27278e+00_dp,0.23995e+00_dp, &
4877 0.21108e+00_dp,0.18567e+00_dp,0.16333e+00_dp,0.14367e+00_dp,0.12638e+00_dp,0.11118e+00_dp /
4878 data blkdnss(1:ndat) / &
4879 0.90000e+00_dp,0.89500e+00_dp,0.88500e+00_dp,0.88200e+00_dp,0.87500e+00_dp,0.86500e+00_dp, &
4880 0.85500e+00_dp,0.84200e+00_dp,0.83000e+00_dp,0.81500e+00_dp,0.80500e+00_dp,0.78600e+00_dp, &
4881 0.76500e+00_dp,0.73000e+00_dp,0.69500e+00_dp,0.61183e+00_dp,0.54000e+00_dp,0.46000e+00_dp, &
4882 0.40000e+00_dp,0.33000e+00_dp,0.28000e+00_dp,0.24000e+00_dp,0.20000e+00_dp,0.16000e+00_dp, &
4883 0.13500e+00_dp,0.11000e+00_dp,0.90000e-01_dp,0.75000e-01_dp,0.60000e-01_dp,0.50000e-01_dp, &
4884 0.40000e-01_dp,0.37500e-01_dp,0.35000e-01_dp /
4885 data blkdnsg(1:ndat) / &
4886 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4887 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4888 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4889 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4890 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4891 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp /
4892 data blkdnsh(1:ndat) / &
4893 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4894 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4895 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4896 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4897 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4898 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
4902 xmss( il,k ) = log( dble( xl( k ) )*1.e-03_dp )
4903 xmss( ic,k ) = log( dble( xi( k,1 ) )*1.e-03_dp )
4904 xmss( ip,k ) = log( dble( xi( k,2 ) )*1.e-03_dp )
4905 xmss( id,k ) = log( dble( xi( k,3 ) )*1.e-03_dp )
4906 xmss( is,k ) = log( dble( xs( k ) )*1.e-03_dp )
4907 xmss( ig,k ) = log( dble( xg( k ) )*1.e-03_dp )
4908 xmss( ih,k ) = log( dble( xh( k ) )*1.e-03_dp )
4913 zcap( il,k ) = dble(rlec( k ))*1.e-02_dp
4914 zcap( ic,k ) = dble(riec( k,1 ))*1.e-02_dp
4915 zcap( ip,k ) = dble(riec( k,2 ))*1.e-02_dp
4916 zcap( id,k ) = dble(riec( k,3 ))*1.e-02_dp
4917 zcap( is,k ) = dble(rsec( k ))*1.e-02_dp
4918 zcap( ig,k ) = dble(rgec( k ))*1.e-02_dp
4919 zcap( ih,k ) = dble(rhec( k ))*1.e-02_dp
4924 vtrm( il,k ) = dble(vrl( k ))*1.e-02_dp
4925 vtrm( ic,k ) = dble(vri( k,1 ))*1.e-02_dp
4926 vtrm( ip,k ) = dble(vri( k,2 ))*1.e-02_dp
4927 vtrm( id,k ) = dble(vri( k,3 ))*1.e-02_dp
4928 vtrm( is,k ) = dble(vrs( k ))*1.e-02_dp
4929 vtrm( ig,k ) = dble(vrg( k ))*1.e-02_dp
4930 vtrm( ih,k ) = dble(vrh( k ))*1.e-02_dp
4935 blkr( il,k ) = dble(blkradl( k ))*10.0_dp
4936 blkr( ic,k ) = dble(blkradi( k,1 ))*10.0_dp
4937 blkr( ip,k ) = dble(blkradi( k,2 ))*10.0_dp
4938 blkr( id,k ) = dble(blkradi( k,3 ))*10.0_dp
4939 blkr( is,k ) = dble(blkrads( k ))*10.0_dp
4940 blkr( ig,k ) = dble(blkradg( k ))*10.0_dp
4941 blkr( ih,k ) = dble(blkradh( k ))*10.0_dp
4946 blkd( il,k ) = dble(blkdnsl( k ))*1000._dp
4947 blkd( ic,k ) = dble(blkdnsi( k,1 ))*1000._dp
4948 blkd( ip,k ) = dble(blkdnsi( k,2 ))*1000._dp
4949 blkd( id,k ) = dble(blkdnsi( k,3 ))*1000._dp
4950 blkd( is,k ) = dble(blkdnss( k ))*1000._dp
4951 blkd( ig,k ) = dble(blkdnsg( k ))*1000._dp
4952 blkd( ih,k ) = dble(blkdnsh( k ))*1000._dp
4961 ykrn( k,kk,i,j ) = kernels( k,kk,i,j )*1.e-06_dp
4968 end subroutine rdkdat
4972 real(DP) :: xsta, xend
4974 real(DP):: pi_dp = 3.1415920_dp
4975 real(DP):: rhow_dp = 1.0e+03_dp
4977 xsta = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-06_dp )**3.0_dp )
4978 xend = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-03_dp )**3.0_dp )
4980 dxmic_mk = ( xend-xsta )/
nbin 4983 xbnd_mk( n ) = xsta + dxmic_mk*( n-1 )
4986 xctr_mk( n ) = ( xbnd_mk( n )+xbnd_mk( n+1 ) )*0.50_dp
4987 radc_mk( n ) = ( exp( xctr_mk( n ) )*3.0_dp/4.0_dp/pi_dp/rhow_dp )**( 1.0_dp/3.0_dp )
4990 end subroutine sdfgrid
4998 cctr_mk( myu,n ) = fcpc( myu,xctr_mk( n ) )
5001 cbnd_mk( myu,n ) = fcpc( myu,xbnd_mk( n ) )
5005 end subroutine getcp
5007 function fcpc( myu,x )
5009 integer,
intent(in) :: myu
5010 real(DP),
intent(in) :: x
5013 real(DP) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5016 ( ndat, kdeg, xmss( myu,: ), &
5020 ( ndat, kdeg, qknt, xmss( myu,: ), &
5024 ( ndat, kdeg, elm, zcap( myu,: ), &
5027 fcpc = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
5033 integer :: myu, nyu, i, j
5035 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Create micpara.dat' 5036 if( kphase == 0 )
then 5037 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Hydro-dynamic kernel' 5038 else if( kphase == 1 )
then 5039 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Long Kernel' 5040 else if( kphase == 2 )
then 5041 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Golovin Kernel' 5046 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
' myu, nyu :', myu, nyu
5049 ck_mk( myu,nyu,i,j ) = fckrn( myu,nyu,xctr_mk( i ),xctr_mk( j ) )
5057 end subroutine getck
5059 function fckrn( myu,nyu,x,y )
5061 integer,
intent(in) :: myu, nyu
5062 real(DP),
intent(in) :: x, y
5065 real(DP) :: qknt( ndat+kdeg ), rknt( ndat+kdeg )
5066 real(DP) :: coef( ndat,ndat )
5068 real(DP) :: xlrg, xsml, vlrg, vsml, rlrg
5069 real(DP):: pi_dp = 3.1415920_dp
5070 real(DP):: rhow_dp = 1.0e+03_dp
5072 if( kphase == 0 )
then 5074 ( ndat, kdeg, xmss( myu,: ), &
5077 rknt( : ) = qknt( : )
5080 ( ndat, ndat, kdeg, kdeg, &
5081 xmss( myu,: ), xmss( nyu,: ), &
5083 ykrn( myu,nyu,:,: ), &
5087 ( ndat, ndat, kdeg, kdeg, &
5089 xmss( myu,: ), xmss( nyu,: ), &
5091 else if( kphase == 1 )
then 5095 vlrg = (exp( xlrg ) / rhow_dp )*1.e+06_dp
5096 vsml = (exp( xsml ) / rhow_dp )*1.e+06_dp
5098 rlrg = ( exp( xlrg )/( 4.0_dp*pi_dp*rhow_dp )*3.0_dp )**(1.0_dp/3.0_dp )*1.e+06_dp
5100 if( rlrg <=50.0_dp )
then 5101 fckrn = 9.44e+03_dp*( vlrg*vlrg + vsml*vsml )
5103 fckrn = 5.78e-03_dp*( vlrg+vsml )
5105 else if( kphase == 2 )
then 5106 fckrn = 1.5_dp*( exp(x) +exp(y) )
5119 vt_mk( myu,n ) = max( fvterm( myu,xctr_mk( n ) ), 0.0_dp )
5123 end subroutine getvt
5125 function fvterm( myu,x )
5127 integer,
intent(in) :: myu
5128 real(DP),
intent(in) :: x
5131 real(DP) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5134 ( ndat, kdeg, xmss( myu,: ), &
5138 ( ndat, kdeg, qknt, xmss( myu,: ), &
5142 ( ndat, kdeg, elm, vtrm( myu,: ), &
5145 fvterm = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
5155 br_mk( myu,n ) = fbulkrad( myu, xctr_mk( n ) )
5159 end subroutine getbr
5161 function fbulkrad( myu,x )
5163 integer,
intent(in) :: myu
5164 real(DP),
intent(in) :: x
5165 real(DP) :: fbulkrad
5167 real(DP) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5170 ( ndat, kdeg, xmss( myu,: ), &
5174 ( ndat, kdeg, qknt, xmss( myu,: ), &
5178 ( ndat, kdeg, elm, blkr( myu,: ), &
5181 fbulkrad = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
5183 end function fbulkrad
5187 integer :: myu, nyu, i, j, n
5189 open ( fid_micpara, file = fname_micpara, form =
'formatted', status=
'new' )
5191 write( fid_micpara,* ) nspc_mk,
nbin 5195 xctr( n ) = xctr_mk( n )
5196 radc( n ) = radc_mk( n )
5197 write( fid_micpara,* ) n, xctr( n ), radc( n )
5200 xbnd( n ) = xbnd_mk( n )
5201 write( fid_micpara,* ) n, xbnd( n )
5203 write( fid_micpara,* ) dxmic_mk
5208 cctr( n,myu ) = cctr_mk( myu,n )
5209 write( fid_micpara,* ) myu, n, cctr( n,myu )
5212 cbnd( n,myu ) = cbnd_mk( myu,n )
5213 write( fid_micpara,* ) myu, n, cbnd( n,myu )
5222 ck( myu,nyu,i,j ) = ck_mk( myu,nyu,i,j )
5223 write( fid_micpara,* ) myu, nyu, i, j, ck( myu,nyu,i,j )
5232 vt( myu,n ) = vt_mk( myu,n )
5233 write( fid_micpara,* ) myu, n, vt( myu,n )
5240 br( myu,n ) = br_mk( myu,n )
5241 write( fid_micpara,* ) myu, n, br( myu,n )
5245 close ( fid_micpara )
5247 end subroutine paraout
5251 subroutine tinvss(n,a,dt,e,nn,iw,inder)
5255 integer,
intent(in) :: n, nn
5256 integer,
intent(inout) :: inder
5257 real(DP),
intent(inout) :: a(nn,n)
5258 real(DP),
intent(inout) :: dt
5259 real(DP),
intent(in) :: e
5260 integer,
intent(inout) :: iw( 2*n )
5261 integer :: i, j, k, kk, ij, nnn
5262 real(DP) :: work, aa, az, eps
5263 integer :: ipiv, jpiv
5270 elseif( n == 1 )
then 5272 elseif( n > 1 )
then 5279 write(6,690)
"n= ", n,
"nn= ", nn, &
5280 "n should be less than or equal to nn in TINVSS" 5281 690
format(a8,i5,5x,a4,i5,a55)
5291 if( abs(a(i,j)) <= abs(piv) )
goto 110
5299 if( abs(piv) <= eps )
goto 920
5300 if( k == 1 ) eps = abs(piv)*e
5301 if( ipiv == k )
goto 130
5309 if( jpiv == k )
goto 150
5324 if( i == k )
goto 220
5326 if( az == 0.0_dp )
goto 220
5328 a(i,j) = a(i,j)-a(k,j)*az
5338 if( ij == k )
goto 420
5346 if( ij == k )
goto 400
5359 write(*,691)
"n= ", n,
"should be positive in TINVSS" 5360 691
format(a8,i5,5x,a30)
5368 write(*,692)
'given matrix A to TINVSS is ill conditioned, or sigular withrank =', nnn, &
5369 'return with no further calculation' 5370 692
format(a,1x,i4,1x,a)
5376 if( dt == 0.0_dp )
goto 920
5377 a(1,1) = 1.0_dp/a(1,1)
5380 end subroutine tinvss
5382 subroutine getknot &
5383 ( ndat, kdeg, xdat, &
5386 integer,
intent(in) :: ndat
5387 integer,
intent(in) :: kdeg
5388 real(DP),
intent(in) :: xdat( ndat )
5390 real(DP),
intent(out) :: qknt( ndat+kdeg )
5396 qknt( i ) = xdat( 1 )
5400 qknt( i+kdeg ) = ( xdat( i )+xdat( i+kdeg ) )*0.50_dp
5404 qknt( ndat+i ) = xdat( ndat )
5409 end subroutine getknot
5411 recursive function fbspl ( ndat, inum, kdeg, qknt, xdat, x ) &
5416 integer,
intent(in) :: ndat
5417 integer,
intent(in) :: inum
5418 integer,
intent(in) :: kdeg
5419 real(DP),
intent(in) :: qknt( ndat+kdeg )
5420 real(DP),
intent(in) :: xdat( ndat )
5421 real(DP),
intent(in) :: x
5424 real(DP) :: bsp1, bsp2
5426 if ( ( inum == 1 .AND. x == xdat( 1 ) ) .OR. &
5427 ( inum == ndat .AND. x == xdat( ndat ) ) )
then 5432 if ( kdeg == 1 )
then 5433 if ( x >= qknt( inum ) .AND. x < qknt( inum+1 ) )
then 5439 if ( qknt( inum+kdeg-1 ) /= qknt( inum ) )
then 5440 bsp1 = ( x-qknt( inum ) ) &
5441 /( qknt( inum+kdeg-1 )-qknt( inum ) ) &
5442 * fbspl( ndat, inum, kdeg-1, qknt, xdat, x )
5446 if ( qknt( inum+kdeg ) /= qknt( inum+1 ) )
then 5447 bsp2 = ( qknt( inum+kdeg )-x ) &
5448 /( qknt( inum+kdeg )-qknt( inum+1 ) ) &
5449 * fbspl( ndat, inum+1, kdeg-1, qknt, xdat, x )
5458 function fpb( ndat, i, kdeg, qknt, xdat, elm, x )
5461 integer :: ndat, i, kdeg
5462 real(DP) :: qknt( ndat+kdeg ), xdat( ndat ), elm( ndat,ndat )
5470 sum = sum + elm( l,i )*fbspl( ndat, l, kdeg, qknt, xdat, x )
5477 subroutine getmatrx &
5478 ( ndat, kdeg, qknt, xdat, &
5483 integer,
intent(in) :: ndat
5484 integer,
intent(in) :: kdeg
5485 real(DP),
intent(in) :: qknt( ndat+kdeg )
5486 real(DP),
intent(in) :: xdat( ndat )
5488 real(DP),
intent(out) :: elm( ndat,ndat )
5492 integer :: iw( 2*ndat ), i, j, inder
5493 real(DP),
parameter :: eps = 0.
5497 elm( i,j ) = fbspl( ndat, j, kdeg, qknt, xdat, xdat( i ) )
5501 call tinvss( ndat, elm, dt, eps, ndat, iw, inder )
5505 end subroutine getmatrx
5507 subroutine getcoef &
5508 ( ndat, kdeg, elm, ydat, &
5511 integer,
intent(in) :: ndat
5512 integer,
intent(in) :: kdeg
5513 real(DP),
intent(in) :: elm( ndat,ndat )
5514 real(DP),
intent(in) :: ydat( ndat )
5516 real(DP),
intent(out) :: coef( ndat )
5525 sum = sum + elm( i,j )*ydat( j )
5532 end subroutine getcoef
5534 function fspline ( ndat, kdeg, coef, qknt, xdat, x )
5536 integer,
intent(in) :: ndat
5537 integer,
intent(in) :: kdeg
5538 real(DP),
intent(in) :: coef( ndat )
5539 real(DP),
intent(in) :: qknt( ndat+kdeg )
5540 real(DP),
intent(in) :: xdat( ndat )
5541 real(DP),
intent(in) :: x
5551 sum = sum + coef( i )*fbspl( ndat, i, kdeg, qknt, xdat, x )
5558 end function fspline
5560 subroutine getcoef2 &
5561 ( mdat, ndat, kdeg, ldeg, &
5562 xdat, ydat, qknt, rknt, zdat, &
5567 integer,
intent(in) :: mdat
5568 integer,
intent(in) :: ndat
5569 integer,
intent(in) :: kdeg
5570 integer,
intent(in) :: ldeg
5571 real(DP),
intent(in) :: xdat( mdat )
5572 real(DP),
intent(in) :: ydat( ndat )
5573 real(DP),
intent(in) :: qknt( mdat+kdeg )
5574 real(DP),
intent(in) :: rknt( ndat+ldeg )
5575 real(DP),
intent(in) :: zdat( mdat,ndat )
5577 real(DP),
intent(out) :: coef( mdat,ndat )
5580 real(DP) :: elmx( mdat,mdat ), elmy( ndat,ndat )
5581 integer :: iw1( 2*mdat ), iw2( 2*ndat )
5582 real(DP) :: beta( mdat,ndat ), sum, dt
5583 real(DP),
parameter :: eps = 0.0_dp
5584 integer :: i, j, k, l, inder
5588 elmx( i,j ) = fbspl( mdat, j, kdeg, qknt, xdat, xdat( i ) )
5591 call tinvss( mdat, elmx, dt, eps, mdat, iw1, inder )
5597 sum = sum + elmx( i,j )*zdat( j,l )
5605 elmy( i,j ) = fbspl( ndat, j, ldeg, rknt, ydat, ydat( i ) )
5608 call tinvss( ndat, elmy, dt, eps, ndat, iw2, inder )
5614 sum = sum + elmy( i,j )*beta( k,j )
5622 end subroutine getcoef2
5625 ( mdat, ndat, kdeg, ldeg, &
5626 coef, qknt, rknt, xdat, ydat, &
5629 integer,
intent(in) :: mdat
5630 integer,
intent(in) :: ndat
5631 integer,
intent(in) :: kdeg
5632 integer,
intent(in) :: ldeg
5633 real(DP),
intent(in) :: coef( mdat,ndat )
5634 real(DP),
intent(in) :: qknt( mdat+kdeg )
5635 real(DP),
intent(in) :: rknt( ndat+ldeg )
5636 real(DP),
intent(in) :: xdat( mdat ), ydat( ndat )
5637 real(DP),
intent(in) :: x, y
5639 real(DP) :: fspline2
5642 real(DP) :: sum, add
5648 add = coef( i,j )*fbspl( mdat, i, kdeg, qknt, xdat, x ) &
5649 *fbspl( ndat, j, ldeg, rknt, ydat, y )
5658 end function fspline2
real(rp), parameter, public const_psat0
saturate pressure of water vapor at 0C [Pa]
integer, public comm_world
communication world ID
module atmosphere / saturation
real(rp), public cv_ice
CV for ice [J/kg/K].
real(rp), parameter, public const_dwatr
density of water [kg/m3]
subroutine, public atmos_phy_mp_suzuki10_setup(KA, IA, JA)
Setup.
real(rp), public cp_ice
CP for ice [J/kg/K].
subroutine, public atmos_phy_mp_suzuki10_terminal_velocity(KA, vterm_o)
get terminal velocity
integer, parameter, public i_hs
snow
real(rp) function, public sf_gamma(x)
Gamma function.
integer, public comm_datatype
datatype of variable
integer, parameter, public i_hr
liquid water rain
integer, parameter, public i_hi
ice water cloud
subroutine, public atmos_phy_mp_suzuki10_qtrc2qhyd(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, Qe)
Calculate mass ratio of each category.
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
integer, public io_fid_conf
Config file ID.
integer, parameter, public i_hh
hail
integer, public atmos_phy_mp_suzuki10_nwaters
real(rp), parameter, public const_tmelt
real(rp), parameter, public const_dice
density of ice [kg/m3]
subroutine, public atmos_phy_mp_suzuki10_effective_radius(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS0, TEMP0, QTRC0, Re)
Calculate Effective Radius.
real(rp), public cv_vapor
CV for vapor [J/kg/K].
module atmosphere / hydrometeor
integer function, public io_get_available_fid()
search & get available file ID
subroutine, public atmos_phy_mp_suzuki10_cloud_fraction(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, mask_criterion, cldfrac)
Calculate Cloud Fraction.
integer, parameter, public prc_masterrank
master process in each communicator
integer, public atmos_phy_mp_suzuki10_ntracers
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
subroutine, public prc_abort
Abort Process.
integer, parameter, public i_hc
liquid water cloud
character(len=h_mid), dimension(:), allocatable, public atmos_phy_mp_suzuki10_tracer_descriptions
subroutine, public atmos_phy_mp_suzuki10_qtrc2nhyd(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, QTRC0, Ne)
Calculate number concentration of each category.
real(rp), parameter, public const_emelt
subroutine, public atmos_phy_mp_suzuki10_qhyd2qtrc(KA, KS, KE, IA, IS, IE, JA, JS, JE, Qe, QTRC, QNUM)
get mass ratio of each category
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
character(len=h_short), dimension(:), allocatable, public atmos_phy_mp_suzuki10_tracer_units
logical, public prc_ismaster
master process in local communicator?
real(rp), public const_eps
small number
integer, public atmos_phy_mp_suzuki10_nices
character(len=h_short), dimension(:), allocatable, public atmos_phy_mp_suzuki10_tracer_names
real(rp), public const_pi
pi
subroutine, public random_uniform(var)
Get uniform random number.
subroutine, public atmos_phy_mp_suzuki10_tracer_setup
Config.
integer, parameter, public n_hyd
module Spectran Bin Microphysics
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
integer, public atmos_phy_mp_suzuki10_nccn
real(rp), public cp_vapor
CP for vapor [J/kg/K].
real(rp), public cp_water
CP for water [J/kg/K].
subroutine, public atmos_phy_mp_suzuki10_tendency(KA, KS, KE, IA, IS, IE, JA, JS, JE, KIJMAX, dt, DENS, PRES, TEMP, QTRC, QDRY, CPtot, CVtot, CCN, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, EVAPORATE)
Cloud Microphysics.
integer, parameter, public i_hg
graupel
real(rp), public cv_water
CV for water [J/kg/K].