72 private :: mp_suzuki10
80 private :: ice_nucleat
88 private :: random_setup
118 integer,
parameter :: i_qv = 1
120 character(len=3) :: namspc (8) = (/
'Qcl', &
129 character(len=27) :: lnamspc(8) = (/
'Mixing ratio of cloud bin', &
130 'Mixing ratio of colum bin', &
131 'Mixing ratio of plate bin', &
132 'Mixing ratio of dendrit bin', &
133 'Mixing ratio of snow bin', &
134 'Mixing ratio of graupel bin', &
135 'Mixing ratio of hail bin', &
136 'Mixing ratio of aerosol bin' /)
141 integer :: kphase = 0
142 integer :: iceflg = 1
143 integer :: num_hyd = 0
145 integer,
parameter :: i_mp_qc = 1
146 integer,
parameter :: i_mp_qcl = 2
147 integer,
parameter :: i_mp_qp = 3
148 integer,
parameter :: i_mp_qd = 4
149 integer,
parameter :: i_mp_qs = 5
150 integer,
parameter :: i_mp_qg = 6
151 integer,
parameter :: i_mp_qh = 7
153 logical :: mp_doautoconversion = .true.
154 logical :: mp_couple_aerosol = .false.
157 integer :: num_start_waters
158 integer :: num_end_waters
159 integer :: num_start_ices
160 integer :: num_end_ices
163 integer,
parameter :: il = 1
164 integer,
parameter :: ic = 2
165 integer,
parameter :: ip = 3
166 integer,
parameter :: id = 4
167 integer,
parameter :: iss = 5
168 integer,
parameter :: ig = 6
169 integer,
parameter :: ih = 7
173 real(
rp),
allocatable :: xctr( : )
174 real(
rp),
allocatable :: xbnd( : )
175 real(
rp),
allocatable :: radc( : )
176 real(
rp),
allocatable :: cctr( :,: )
177 real(
rp),
allocatable :: cbnd( :,: )
178 real(
rp),
allocatable :: ck( :,:,:,: )
179 real(
rp),
allocatable :: vt( :,: )
180 real(
rp),
allocatable :: br( :,: )
181 integer,
allocatable :: ifrsl( :,:,: )
183 real(
rp),
allocatable :: xactr( : )
184 real(
rp),
allocatable :: xabnd( : )
185 real(
rp),
allocatable :: rada( : )
187 real(
rp),
allocatable :: expxctr( : )
188 real(
rp),
allocatable :: expxbnd( : )
189 real(
rp),
allocatable :: expxactr( : )
190 real(
rp),
allocatable :: expxabnd( : )
191 real(
rp),
allocatable :: rexpxctr( : )
192 real(
rp),
allocatable :: rexpxbnd( : )
193 real(
rp),
allocatable :: rexpxactr( : )
194 real(
rp),
allocatable :: rexpxabnd( : )
200 real(
rp),
allocatable,
save :: vterm(:)
204 real(
rp),
parameter :: cldmin = 1.0e-10_rp
205 real(
rp),
parameter :: oneovthird = 1.0_rp / 3.0_rp
206 real(
rp),
parameter :: thirdovforth = 3.0_rp / 4.0_rp
207 real(
rp),
parameter :: twoovthird = 2.0_rp / 3.0_rp
209 real(
rp) :: rbnd = 40.e-06_rp
213 real(
rp) :: rhoa = 2.25e+03_rp
214 real(
rp) :: emaer = 58.0_rp
215 real(
rp) :: emwtr = 18.0_rp
216 real(
rp) :: rasta = 1.e-08_rp
217 real(
rp) :: raend = 1.e-06_rp
218 real(
rp) :: r0a = 1.e-07_rp
220 logical :: flg_regeneration
222 logical :: flg_icenucl
223 logical :: flg_sf_aero
226 real(
rp),
allocatable :: marate( : )
232 character(len=11),
parameter :: fname_micpara=
"micpara.dat"
234 integer(4) :: fid_micpara
237 integer,
allocatable :: blrg( :,: )
238 integer,
allocatable :: bsml( :,: )
240 integer :: mspc, mbin
243 real(
rp) :: n0_icenucl = 1.e+3_rp
246 real(
rp) :: c_ccn = 100.e+6_rp
247 real(
rp) :: kappa = 0.462_rp
250 real(
rp) :: sigma = 7.5e-02_rp
251 real(
rp) :: vhfct = 2.0_rp
253 real(
rp),
parameter :: tcrit = 271.15_rp
254 integer,
allocatable :: kindx( :,: )
257 integer,
parameter :: ndat = 33, icemax = 3
258 integer,
parameter :: kdeg = 4, ldeg = 4, nspc_mk = 7
262 real(
dp),
allocatable :: radc_mk( : ), xctr_mk( : ), xbnd_mk( : )
263 real(
dp),
allocatable :: cctr_mk( :,: ), cbnd_mk( :,: )
264 real(
dp),
allocatable :: ck_mk( :,:,:,: )
265 real(
dp),
allocatable :: vt_mk( :,: )
266 real(
dp),
allocatable :: br_mk( :,: )
268 real(
dp) :: xmss( nspc_mk,ndat ), zcap( nspc_mk,ndat ), vtrm( nspc_mk,ndat )
269 real(
dp) :: blkr( nspc_mk,ndat ), blkd( nspc_mk,ndat ), ykrn( nspc_mk,nspc_mk,ndat,ndat )
271 real(
dp) :: ywll( ndat,ndat ), ywli( ndat,ndat,icemax ), ywls( ndat,ndat )
272 real(
dp) :: ywlg( ndat,ndat ), ywlh( ndat,ndat )
274 real(
dp) :: ywil( ndat,ndat,icemax ), ywii( ndat,ndat,icemax,icemax )
275 real(
dp) :: ywis( ndat,ndat,icemax ), ywig( ndat,ndat,icemax )
276 real(
dp) :: ywih( ndat,ndat,icemax )
278 real(
dp) :: ywsl( ndat,ndat ), ywsi( ndat,ndat,icemax ), ywss( ndat,ndat )
279 real(
dp) :: ywsg( ndat,ndat ), ywsh( ndat,ndat )
281 real(
dp) :: ywgl( ndat,ndat ), ywgi( ndat,ndat,icemax ), ywgs( ndat,ndat )
282 real(
dp) :: ywgg( ndat,ndat ), ywgh( ndat,ndat )
284 real(
dp) :: ywhl( ndat,ndat ), ywhi( ndat,ndat,icemax ), ywhs( ndat,ndat )
285 real(
dp) :: ywhg( ndat,ndat ), ywhh( ndat,ndat )
288 real(
rp) :: sigma_sdf(5), r0_sdf(5), n0_sdf(5), rho_sdf(5)
291 real(
rp),
allocatable :: flg_noninduct(:,:)
292 real(
rp),
allocatable :: ecoll(:,:,:,:), rcoll(:,:,:,:)
293 real(
rp),
allocatable,
save :: num_end(:,:,:,:)
294 real(
rp),
save :: ecoal_gsi = 0.0_rp
304 namelist / param_atmos_phy_mp_suzuki10_bin / &
310 integer :: m, n, ierr
314 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Setup'
315 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Tracers setup for Suzuki (2010) Spectral BIN model'
318 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'READ BIN NUMBER'
321 read(
io_fid_conf,nml=param_atmos_phy_mp_suzuki10_bin,iostat=ierr)
324 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Not found namelist. Default used.'
325 elseif( ierr > 0 )
then
326 log_error(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10_bin, Check!'
330 log_nml(param_atmos_phy_mp_suzuki10_bin)
332 if( iceflg == 0 )
then
334 elseif( iceflg == 1 )
then
337 log_error(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
"ICEFLG should be 0 (warm rain) or 1 (mixed rain) check!!"
346 num_hyd =
nbin * nspc
348 num_start_waters = i_qv + 1
350 num_start_ices = num_end_waters + 1
413 integer,
intent(in) :: ka
414 integer,
intent(in) :: ia
415 integer,
intent(in) :: ja
416 logical,
intent(in),
optional :: flg_lt
422 real(
rp) :: s10_emaer
424 logical :: s10_flag_regene = .false.
425 logical :: s10_flag_nucleat = .false.
426 logical :: s10_flag_icenucleat = .false.
427 logical :: s10_flag_sfaero = .false.
428 logical :: s10_flag_rndm = .false.
429 integer :: s10_rndm_mspc
430 integer :: s10_rndm_mbin
432 namelist / param_atmos_phy_mp_suzuki10 / &
433 mp_doautoconversion, &
442 s10_flag_icenucleat, &
452 real(
rp),
parameter :: max_term_vel = 10.0_rp
455 integer :: nnspc, nnbin
456 integer :: nn, mm, mmyu, nnyu
457 integer :: myu, nyu, i, j, k, n, ierr
461 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Setup'
462 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Suzuki (2010) Spectral BIN model'
465 allocate( xctr(
nbin ) )
466 allocate( xbnd(
nbin+1 ) )
467 allocate( radc(
nbin ) )
468 allocate( cctr(
nbin,nspc_mk ) )
469 allocate( cbnd(
nbin+1,nspc_mk ) )
470 allocate( ck( nspc_mk,nspc_mk,
nbin,
nbin ) )
471 allocate( vt( nspc_mk,
nbin ) )
472 allocate( br( nspc_mk,
nbin ) )
473 allocate( ifrsl( 2,nspc_mk,nspc_mk ) )
474 allocate( expxctr(
nbin ) )
475 allocate( expxbnd(
nbin+1 ) )
476 allocate( rexpxctr(
nbin ) )
477 allocate( rexpxbnd(
nbin+1 ) )
478 if (
nccn /= 0 )
then
479 allocate( xactr(
nccn ) )
480 allocate( xabnd(
nccn+1 ) )
481 allocate( rada(
nccn ) )
482 allocate( expxactr(
nccn ) )
483 allocate( expxabnd(
nccn+1 ) )
484 allocate( rexpxactr(
nccn ) )
485 allocate( rexpxabnd(
nccn+1 ) )
489 mspc = nspc_mk*nspc_mk
500 read(
io_fid_conf,nml=param_atmos_phy_mp_suzuki10,iostat=ierr)
503 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Not found namelist. Default used.'
504 elseif( ierr > 0 )
then
505 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10, Check!'
508 log_nml(param_atmos_phy_mp_suzuki10)
510 if ( nspc /= 1 .AND. nspc /= 7 )
then
511 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nspc should be set as 1 (warm rain) or 7 (mixed phase) check!'
520 flg_regeneration = s10_flag_regene
521 flg_nucl = s10_flag_nucleat
522 flg_icenucl = s10_flag_icenucleat
523 flg_sf_aero = s10_flag_sfaero
524 flg_rndm = s10_flag_rndm
533 open ( fid_micpara, file = fname_micpara, form =
'formatted', status =
'old', iostat=ierr )
536 if ( ierr == 0 )
then
538 read( fid_micpara,* ) nnspc, nnbin
540 if ( nnbin /=
nbin )
then
541 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nbin in namelist and nbin in micpara.dat is different check!'
546 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Radius of cloud *'
548 read( fid_micpara,* ) nn, xctr( n ), radc( n )
549 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
550 "Radius of ", n,
"th cloud bin (bin center)= ", radc( n ) ,
"[m]"
553 read( fid_micpara,* ) nn, xbnd( n )
555 read( fid_micpara,* ) dxmic
556 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Width of Cloud SDF= ', dxmic
562 read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
566 read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
575 read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
584 read( fid_micpara,* ) mmyu, nn, vt( myu,n )
591 read( fid_micpara,* ) mmyu, nn, br( myu,n )
595 close ( fid_micpara )
600 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'micpara.dat is created'
605 open ( fid_micpara, file = fname_micpara, form =
'formatted', status =
'old', iostat=ierr )
607 read( fid_micpara,* ) nnspc, nnbin
609 if ( nnbin /=
nbin )
then
610 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nbin in inc_tracer and nbin in micpara.dat is different check!'
615 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Radius of cloud *'
617 read( fid_micpara,* ) nn, xctr( n ), radc( n )
618 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
619 "Radius of ", n,
"th cloud bin (bin center)= ", radc( n ) ,
"[m]"
622 read( fid_micpara,* ) nn, xbnd( n )
624 read( fid_micpara,* ) dxmic
625 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Width of Cloud SDF= ', dxmic
631 read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
635 read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
644 read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
653 read( fid_micpara,* ) mmyu, nn, vt( myu,n )
660 read( fid_micpara,* ) mmyu, nn, br( myu,n )
664 close ( fid_micpara )
670 call comm_bcast(
nbin, radc(:) )
671 call comm_bcast(
nbin, xctr(:) )
672 call comm_bcast( dxmic )
673 call comm_bcast(
nbin+1, xbnd(:) )
674 call comm_bcast(
nbin, nspc_mk, cctr(:,:) )
675 call comm_bcast(
nbin+1, nspc_mk, cbnd(:,:) )
676 call comm_bcast( nspc_mk, nspc_mk,
nbin,
nbin, ck(:,:,:,:) )
677 call comm_bcast( nspc_mk,
nbin, br(:,:) )
678 call comm_bcast( nspc_mk,
nbin, vt(:,:) )
680 allocate( flg_noninduct( nspc,nspc ) )
681 allocate( ecoll( nspc,nspc,
nbin,
nbin ) )
682 allocate( rcoll( nspc,nspc,
nbin,
nbin ) )
684 flg_noninduct(:,:) = 0.0_rp
685 ecoll( :,:,:,: ) = 0.0_rp
686 rcoll( :,:,:,: ) = 0.0_rp
689 if (
present(flg_lt) )
then
699 if( ( myu >= ic .and. myu <= iss ) .and. ( nyu == ig .or. nyu == ih ) )
then
700 flg_noninduct( myu,nyu ) = 1.0_rp
710 if( vt( myu,i ) /= vt( nyu,j ) )
then
711 ecoll( myu,nyu,i,j ) = ck( myu,nyu,i,j ) &
712 / ( pi*( radc( i )+radc( j ) )**2 * abs( vt( myu,i )-vt( nyu,j ) ) )
713 ecoll( myu,nyu,i,j ) = max( min( 1.0_rp, ecoll( myu,nyu,i,j ) ),0.0_rp )
715 ecoll( myu,nyu,i,j ) = 0.0_rp
718 if( ecoal_gsi /= 0.0_rp )
then
719 ecoll( myu,nyu,i,j ) = ecoal_gsi
722 if( ecoll( myu,nyu,i,j ) /= 0.0_rp )
then
723 rcoll( myu,nyu,i,j ) = ( 1.0_rp-ecoll( myu,nyu,i,j ) ) / ecoll( myu,nyu,i,j )
724 elseif( ecoll( myu,nyu,i,j ) == 0.0_rp )
then
725 rcoll( myu,nyu,i,j ) = 1.0_rp
735 if (
nccn /= 0 )
then
737 xasta = log( rhoa*4.0_rp/3.0_rp*pi * ( rasta )**3 )
738 xaend = log( rhoa*4.0_rp/3.0_rp*pi * ( raend )**3 )
740 dxaer = ( xaend-xasta )/
nccn
744 xabnd( n ) = xasta + dxaer*( n-1 )
748 xactr( n ) = ( xabnd( n )+xabnd( n+1 ) )*0.50_rp
749 rada( n ) = ( exp( xactr( n ) )*thirdovforth/pi/rhoa )**( oneovthird )
752 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
753 "Radius of ", n,
"th aerosol bin (bin center)= ", rada( n ) ,
"[m]"
756 if ( flg_sf_aero )
then
757 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
"flg_sf_aero=true is not supported stop!! "
788 if( radc( n ) > rbnd )
then
793 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,ES15.7,A)')
'Radius between cloud and rain is ', radc(nbnd),
'[m]'
798 call random_setup( ia*ja*ka )
801 if ( mp_couple_aerosol .AND.
nccn /=0 )
then
802 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nccn should be 0 when MP_couple_aerosol = .true. !! stop'
806 if (
nccn /= 0 )
then
809 expxactr( n ) = exp( xactr( n ) )
810 rexpxactr( n ) = 1.0_rp / exp( xactr( n ) )
814 expxabnd( n ) = exp( xabnd( n ) )
815 rexpxabnd( n ) = 1.0_rp / exp( xabnd( n ) )
819 allocate( vterm(qa-1) )
826 vterm((myu-1)*
nbin+n) = -vt( myu,n )
831 expxctr( n ) = exp( xctr( n ) )
832 rexpxctr( n ) = 1.0_rp / exp( xctr( n ) )
836 expxbnd( n ) = exp( xbnd( n ) )
837 rexpxbnd( n ) = 1.0_rp / exp( xbnd( n ) )
841 call getrule( ifrsl,kindx )
844 sigma_sdf(1) = 0.2_rp
845 sigma_sdf(2) = 0.35_rp
846 sigma_sdf(3) = 0.35_rp
847 sigma_sdf(4) = 0.35_rp
848 sigma_sdf(5) = 0.35_rp
850 r0_sdf(2) = 2.61e-6_rp
852 r0_sdf(4) = 2.61e-6_rp
853 r0_sdf(5) = 2.61e-6_rp
854 n0_sdf(1) = 8.0e+6_rp
856 n0_sdf(3) = 3.0e+6_rp
857 n0_sdf(4) = 4.0e+6_rp
858 n0_sdf(5) = 4.0e+6_rp
861 rho_sdf(3) = 100.0_rp
862 rho_sdf(4) = 400.0_rp
863 rho_sdf(5) = 400.0_rp
866 sigma_sdf(1) = 0.2_rp
867 sigma_sdf(2) = 0.35_rp
868 sigma_sdf(3) = 0.35_rp
869 sigma_sdf(4) = 0.35_rp
870 sigma_sdf(5) = 0.35_rp
872 r0_sdf(2) = 2.61e-6_rp
874 r0_sdf(4) = 2.61e-6_rp
875 r0_sdf(5) = 2.61e-6_rp
876 n0_sdf(1) = 8.0e+6_rp
878 n0_sdf(3) = 3.0e+6_rp
879 n0_sdf(4) = 4.0e+6_rp
880 n0_sdf(5) = 4.0e+6_rp
883 rho_sdf(3) = 100.0_rp
884 rho_sdf(4) = 400.0_rp
885 rho_sdf(5) = 400.0_rp
909 deallocate( expxctr )
910 deallocate( expxbnd )
911 deallocate( rexpxctr )
912 deallocate( rexpxbnd )
913 if (
nccn /= 0 )
then
917 deallocate( expxactr )
918 deallocate( expxabnd )
919 deallocate( rexpxactr )
920 deallocate( rexpxabnd )
923 deallocate( flg_noninduct )
937 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
958 atmos_saturation_pres2qsat_liq, &
959 atmos_saturation_pres2qsat_ice
962 integer,
intent(in) :: ka, ks, ke
963 integer,
intent(in) :: ia, is, ie
964 integer,
intent(in) :: ja, js, je
965 integer,
intent(in) :: kijmax
967 real(
dp),
intent(in) :: dt
968 real(
rp),
intent(in) :: dens (ka,ia,ja)
969 real(
rp),
intent(in) :: pres (ka,ia,ja)
970 real(
rp),
intent(in) :: temp (ka,ia,ja)
971 real(
rp),
intent(in) :: qtrc (ka,ia,ja,qa)
972 real(
rp),
intent(in) :: qdry (ka,ia,ja)
973 real(
rp),
intent(in) :: cptot(ka,ia,ja)
974 real(
rp),
intent(in) :: cvtot(ka,ia,ja)
975 real(
rp),
intent(in) :: ccn (ka,ia,ja)
977 real(
rp),
intent(out) :: rhoq_t (ka,ia,ja,qa)
978 real(
rp),
intent(out) :: rhoe_t (ka,ia,ja)
979 real(
rp),
intent(out) :: cptot_t(ka,ia,ja)
980 real(
rp),
intent(out) :: cvtot_t(ka,ia,ja)
981 real(
rp),
intent(out) :: evaporate(ka,ia,ja)
984 logical,
intent(in),
optional :: flg_lt
985 real(
rp),
intent(in),
optional :: d0_crg, v0_crg
986 real(
rp),
intent(in),
optional :: dqcrg(ka,ia,ja)
987 real(
rp),
intent(in),
optional :: beta_crg(ka,ia,ja)
988 real(
rp),
intent(in),
optional :: qtrc_crg(ka,ia,ja,num_hyd)
989 real(
rp),
intent(out),
optional :: qsplt_in(ka,ia,ja,3)
990 real(
rp),
intent(out),
optional :: sarea(ka,ia,ja,num_hyd)
991 real(
rp),
intent(out),
optional :: rhoc_t_mp(ka,ia,ja,num_hyd)
993 real(
rp) :: qsat_l(ka,ia,ja)
994 real(
rp) :: qsat_i(ka,ia,ja)
995 real(
rp) :: ssliq(ka,ia,ja)
996 real(
rp) :: ssice(ka,ia,ja)
998 integer :: ijk_index (kijmax,3)
999 integer :: index_cld (kijmax)
1000 integer :: index_cold(kijmax)
1001 integer :: index_warm(kijmax)
1002 integer :: ijkcount, ijkcount_cold, ijkcount_warm
1003 integer :: ijk, indirect
1005 real(
rp) :: dens_ijk(kijmax)
1006 real(
rp) :: pres_ijk(kijmax)
1007 real(
rp) :: temp_ijk(kijmax)
1008 real(
rp) :: qdry_ijk(kijmax)
1009 real(
rp) :: qvap_ijk(kijmax)
1010 real(
rp) :: ccn_ijk(kijmax)
1011 real(
rp) :: cp_ijk(kijmax)
1012 real(
rp) :: cv_ijk(kijmax)
1013 real(
rp) :: evaporate_ijk(kijmax)
1014 real(
rp) :: ghyd_ijk(
nbin,nspc,kijmax)
1015 real(
rp) :: gaer_ijk(max(
nccn,1),kijmax)
1018 real(
rp) :: rhoq_new
1022 real(
rp) :: d0_crg_l, v0_crg_l
1023 real(
rp) :: gcrg_ijk(
nbin,nspc,kijmax)
1024 real(
rp) :: crg_sep_ijk(nspc,kijmax)
1025 real(
rp) :: dqcrg_ijk(kijmax)
1026 real(
rp) :: beta_crg_ijk(kijmax)
1029 integer :: k, i, j, m, n, iq
1032 if ( nspc == 1 )
then
1033 log_progress(*)
'atmosphere / physics / microphysics / SBM (Liquid water only)'
1034 elseif( nspc > 1 )
then
1035 log_progress(*)
'atmosphere / physics / microphysics / SBM (Mixed phase)'
1039 if (
present(flg_lt) )
then
1043 if ( flg_lt_l )
then
1045 crg_sep_ijk(:,:) = 0.0_rp
1046 gcrg_ijk(:,:,:) = 0.0_rp
1047 qsplt_in(:,:,:,:) = 0.0_rp
1053 call atmos_saturation_pres2qsat_liq( ka, ks, ke, &
1056 temp(:,:,:), pres(:,:,:), qdry(:,:,:), &
1059 call atmos_saturation_pres2qsat_ice( ka, ks, ke, &
1062 temp(:,:,:), pres(:,:,:), qdry(:,:,:), &
1069 ssliq(k,i,j) = qtrc(k,i,j,i_qv) / qsat_l(k,i,j) - 1.0_rp
1070 ssice(k,i,j) = qtrc(k,i,j,i_qv) / qsat_i(k,i,j) - 1.0_rp
1075 if ( nspc == 1 )
then
1077 ssice(:,:,:) = 0.0_rp
1110 ijk = 1 + (k-ks) + (ke-ks+1) * ( (i-is) + (ie-is+1) * (j-js) )
1111 ijk_index(ijk,1) = i
1112 ijk_index(ijk,2) = j
1113 ijk_index(ijk,3) = k
1127 cldsum = cldsum + qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
1128 countbin = countbin + 1
1132 if ( cldsum > cldmin &
1133 .OR. ssliq(k,i,j) > 0.0_rp &
1134 .OR. ssice(k,i,j) > 0.0_rp )
then
1135 ijkcount = ijkcount + 1
1136 ijk = 1 + (k-ks) + (ke-ks+1) * ( (i-is) + (ie-is+1) * (j-js) )
1137 index_cld(ijkcount) = ijk
1143 rhoq_t(k,i,j,iq) = 0.0_rp
1145 rhoe_t(k,i,j) = 0.0_rp
1146 cptot_t(k,i,j) = 0.0_rp
1147 cvtot_t(k,i,j) = 0.0_rp
1148 evaporate(k,i,j) = 0.0_rp
1157 do ijk = 1, ijkcount
1158 indirect = index_cld(ijk)
1159 i = ijk_index(indirect,1)
1160 j = ijk_index(indirect,2)
1161 k = ijk_index(indirect,3)
1163 dens_ijk(ijk) = dens(k,i,j)
1164 pres_ijk(ijk) = pres(k,i,j)
1165 temp_ijk(ijk) = temp(k,i,j)
1166 qdry_ijk(ijk) = qdry(k,i,j)
1167 cp_ijk(ijk) = cptot(k,i,j)
1168 cv_ijk(ijk) = cvtot(k,i,j)
1169 ccn_ijk(ijk) = ccn(k,i,j)
1170 qvap_ijk(ijk) = qtrc(k,i,j,i_qv)
1175 ghyd_ijk(n,m,ijk) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
1176 countbin = countbin + 1
1181 gaer_ijk(n,ijk) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxaer
1182 countbin = countbin + 1
1189 if ( nspc > 1 )
then
1190 do ijk = 1, ijkcount
1191 if ( temp_ijk(ijk) < tem00 )
then
1192 ijkcount_cold = ijkcount_cold + 1
1193 index_cold(ijkcount_cold) = ijk
1195 ijkcount_warm = ijkcount_warm + 1
1196 index_warm(ijkcount_warm) = ijk
1201 do ijk = 1, ijkcount
1202 index_warm(ijk) = ijk
1204 ijkcount_warm = ijkcount
1207 if ( flg_lt_l )
then
1209 do ijk = 1, ijkcount
1210 indirect = index_cld(ijk)
1211 i = ijk_index(indirect,1)
1212 j = ijk_index(indirect,2)
1213 k = ijk_index(indirect,3)
1218 gcrg_ijk(n,m,ijk) = qtrc_crg(k,i,j,countbin) * dens(k,i,j)
1219 countbin = countbin + 1
1222 beta_crg_ijk(ijk) = beta_crg(k,i,j)
1223 dqcrg_ijk(ijk) = dqcrg(k,i,j)
1261 if ( ijkcount > 0 )
then
1265 call mp_suzuki10( ka, ia, ja, &
1269 index_cold( 1:ijkcount), &
1270 index_warm( 1:ijkcount), &
1271 dens_ijk( 1:ijkcount), &
1272 pres_ijk( 1:ijkcount), &
1273 qdry_ijk( 1:ijkcount), &
1274 ccn_ijk( 1:ijkcount), &
1275 temp_ijk( 1:ijkcount), &
1276 qvap_ijk( 1:ijkcount), &
1277 ghyd_ijk(:,:,1:ijkcount), &
1278 gaer_ijk(:, 1:ijkcount), &
1279 cp_ijk( 1:ijkcount), &
1280 cv_ijk( 1:ijkcount), &
1281 evaporate_ijk(1:ijkcount), &
1284 d0_crg_l, v0_crg_l, &
1285 dqcrg_ijk( 1:ijkcount), &
1286 beta_crg_ijk( 1:ijkcount), &
1287 gcrg_ijk(:,:,1:ijkcount), &
1288 crg_sep_ijk(:,1:ijkcount) )
1338 do ijk = 1, ijkcount
1339 indirect = index_cld(ijk)
1340 i = ijk_index(indirect,1)
1341 j = ijk_index(indirect,2)
1342 k = ijk_index(indirect,3)
1344 rhoe_t(k,i,j) = ( temp_ijk(ijk) * cv_ijk(ijk) - temp(k,i,j) * cvtot(k,i,j) ) * dens(k,i,j) / dt
1345 cptot_t(k,i,j) = ( cp_ijk(ijk) - cptot(k,i,j) ) / dt
1346 cvtot_t(k,i,j) = ( cv_ijk(ijk) - cvtot(k,i,j) ) / dt
1347 evaporate(k,i,j) = evaporate_ijk(ijk) / dt
1349 rhoq_t(k,i,j,i_qv) = ( qvap_ijk(ijk) - qtrc(k,i,j,i_qv) ) * dens(k,i,j) / dt
1354 rhoq_new = ghyd_ijk(n,m,ijk) * dxmic
1355 rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
1356 countbin = countbin + 1
1361 rhoq_new = gaer_ijk(n,ijk) * dxaer
1362 rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
1363 countbin = countbin + 1
1370 rhoq_new = gcrg_ijk(n,m,ijk)
1371 rhoc_t_mp(k,i,j,countbin) = ( rhoq_new - qtrc_crg(k,i,j,countbin)*dens(k,i,j) ) / dt
1372 countbin = countbin + 1
1376 sarea(k,i,j,:) = 0.0_rp
1380 rhoq_new = ghyd_ijk(n,m,ijk) * dxmic
1381 sarea(k,i,j,countbin) = 4.0_rp*pi*radc(n)**2*rhoq_new*rexpxctr( n )
1382 countbin = countbin + 1
1386 qsplt_in(k,i,j,1) = crg_sep_ijk(ig,ijk) * dens(k,i,j) / dt
1387 qsplt_in(k,i,j,2) = ( crg_sep_ijk(ic,ijk) &
1388 + crg_sep_ijk(ip,ijk) &
1389 + crg_sep_ijk(id,ijk) ) * dens(k,i,j) / dt
1390 qsplt_in(k,i,j,3) = crg_sep_ijk(iss,ijk) * dens(k,i,j) / dt
1426 integer,
intent(in) :: ka
1428 real(
rp),
intent(out) :: vterm_o(ka,qa-1)
1434 vterm_o(:,iq) = vterm(iq)
1451 integer,
intent(in) :: ka, ks, ke
1452 integer,
intent(in) :: ia, is, ie
1453 integer,
intent(in) :: ja, js, je
1455 real(
rp),
intent(in) :: qtrc0 (ka,ia,ja,num_hyd)
1456 real(
rp),
intent(in) :: mask_criterion
1457 real(
rp),
intent(out) :: cldfrac(ka,ia,ja)
1460 integer :: k, i, j, iq, ihydro
1470 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1471 qhydro = qhydro + qtrc0(k,i,j,iq)
1474 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1478 elseif( nspc == 1 )
then
1484 do ihydro = 1, i_mp_qc
1485 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1486 qhydro = qhydro + qtrc0(k,i,j,iq)
1489 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1521 integer,
intent(in) :: ka, ks, ke
1522 integer,
intent(in) :: ia, is, ie
1523 integer,
intent(in) :: ja, js, je
1525 real(
rp),
intent(in) :: dens0(ka,ia,ja)
1526 real(
rp),
intent(in) :: temp0(ka,ia,ja)
1527 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1528 real(
rp),
intent(out) :: re (ka,ia,ja,
n_hyd)
1530 real(
rp),
parameter :: um2cm = 100.0_rp
1532 real(
rp) :: sum0(nspc), sum2, sum3, re_tmp(nspc)
1533 integer :: i, j, k, iq, ihydro
1540 re(k,i,j,:) = 0.0_rp
1548 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1549 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1550 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1552 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1553 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1554 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1556 sum3 = max( sum3, 0.0_rp )
1557 sum2 = max( sum2, 0.0_rp )
1558 if ( sum2 /= 0.0_rp )
then
1559 re(k,i,j,
i_hc) = sum3 / sum2 * um2cm
1561 re(k,i,j,
i_hc) = 0.0_rp
1568 do iq = nbnd+1,
nbin
1570 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1571 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1572 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1574 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1575 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1576 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1578 sum3 = max( sum3, 0.0_rp )
1579 sum2 = max( sum2, 0.0_rp )
1580 if ( sum2 /= 0.0_rp )
then
1581 re(k,i,j,
i_hr) = sum3 / sum2 * um2cm
1583 re(k,i,j,
i_hr) = 0.0_rp
1591 if ( nspc > 1 )
then
1597 sum0(ihydro) = 0.0_rp
1600 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1601 sum0(ihydro) = sum0(ihydro) &
1602 + ( qtrc0(k,i,j,iq) * dens0(k,i,j) )
1604 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1605 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1606 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1608 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1609 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1610 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1612 sum3 = max( sum3, 0.0_rp )
1613 sum2 = max( sum2, 0.0_rp )
1614 if ( sum2 == 0.0_rp )
then
1615 re_tmp(ihydro) = 0.0_rp
1617 re_tmp(ihydro) = sum3 / sum2 * um2cm
1621 re(k,i,j,
i_hi) = ( re_tmp(i_mp_qcl) * sum0(i_mp_qcl) &
1622 + re_tmp(i_mp_qp ) * sum0(i_mp_qp ) &
1623 + re_tmp(i_mp_qd ) * sum0(i_mp_qd ) ) &
1624 / ( sum0(i_mp_qcl) + sum0(i_mp_qp) + sum0(i_mp_qd) + eps )
1625 re(k,i,j,
i_hs) = re_tmp(i_mp_qs)
1626 re(k,i,j,
i_hg) = re_tmp(i_mp_qg)
1627 re(k,i,j,
i_hh) = re_tmp(i_mp_qh)
1656 integer,
intent(in) :: ka, ks, ke
1657 integer,
intent(in) :: ia, is, ie
1658 integer,
intent(in) :: ja, js, je
1659 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1660 real(
rp),
intent(out) :: qe (ka,ia,ja,
n_hyd)
1662 integer :: ihydro, ibin, iq, icateg
1668 qe(:,:,:,:) = 0.0_rp
1673 iq =
nbin*(ihydro-1) + ibin
1675 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc ) )
then
1676 if ( iq > 0 .AND. iq <= nbnd )
then
1678 elseif( iq > nbnd .AND. iq <=
nbin )
then
1681 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qcl) )
then
1683 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qp ) )
then
1685 elseif( iq >
nbin*(i_mp_qd -1) .AND. iq <=
nbin*(i_mp_qd ) )
then
1687 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qs ) )
then
1689 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qg ) )
then
1691 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <= num_hyd )
then
1699 qe(k,i,j,icateg) = qe(k,i,j,icateg) + qtrc0(k,i,j,iq)
1727 integer,
intent(in) :: ka, ks, ke
1728 integer,
intent(in) :: ia, is, ie
1729 integer,
intent(in) :: ja, js, je
1730 real(
rp),
intent(in) :: dens (ka,ia,ja)
1731 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1732 real(
rp),
intent(out) :: ne (ka,ia,ja,
n_hyd)
1734 integer :: ihydro, ibin, iq, icateg
1740 ne(:,:,:,:) = 0.0_rp
1745 iq =
nbin*(ihydro-1) + ibin
1747 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc ) )
then
1748 if ( iq > 0 .AND. iq <= nbnd )
then
1750 elseif( iq > nbnd .AND. iq <=
nbin )
then
1753 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qcl) )
then
1755 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qp ) )
then
1757 elseif( iq >
nbin*(i_mp_qd -1) .AND. iq <=
nbin*(i_mp_qd ) )
then
1759 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qs ) )
then
1761 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qg ) )
then
1763 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <= num_hyd )
then
1771 ne(k,i,j,icateg) = ne(k,i,j,icateg) + dens(k,i,j) * qtrc0(k,i,j,iq) * rexpxctr(ibin)
1784 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1802 integer,
intent(in) :: ka, ks, ke
1803 integer,
intent(in) :: ia, is, ie
1804 integer,
intent(in) :: ja, js, je
1806 real(
rp),
intent(in) :: qe(ka,ia,ja,
n_hyd)
1808 real(
rp),
intent(out) :: qtrc(ka,ia,ja,qa-1)
1810 real(
rp),
intent(in),
optional :: qnum(ka,ia,ja,
n_hyd)
1812 real(
rp) :: coef0, coef1, coef2
1814 real(
rp) :: tmp_hyd, num_hyd_l, lambda_hyd
1816 integer :: k, i, j, iq
1818 if (
present(qnum) )
then
1819 log_warn(
"ATMOS_PHY_MP_suzuki10_qhyd2qtrc",*)
'At this moment, number concentratio is ignored'
1823 coef0 = 4.0_rp/3.0_rp*pi
1824 coef1 = 4.0_rp/3.0_rp*sqrt(pi/2.0_rp)
1826 if( nspc == 1 )
then
1835 dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
1837 - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
1838 / sigma_sdf(1) / sigma_sdf(1) &
1840 tmp_hyd = tmp_hyd + dummy(iq)
1843 coef2 = ( qe(k,i,j,
i_hc) + qe(k,i,j,
i_hr) &
1845 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1848 qtrc(k+2,i,j,(il-1)*
nbin+iq) = coef2 * dummy(iq)
1855 elseif( nspc > 1 )
then
1865 dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
1867 - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
1868 / sigma_sdf(1) / sigma_sdf(1) &
1870 tmp_hyd = tmp_hyd + dummy(iq)
1873 coef2 = ( qe(k,i,j,
i_hc) + qe(k,i,j,
i_hr) ) &
1874 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1877 qtrc(k,i,j,(il-1)*
nbin+iq) = coef2 * dummy(iq)
1883 dummy(iq) = coef1 / sigma_sdf(2) * rho_sdf(2) * radc( iq )**3 &
1885 - ( log( radc(iq) )-log( r0_sdf(2) ) )**2*0.5_rp &
1886 / sigma_sdf(2) / sigma_sdf(2) &
1888 tmp_hyd = tmp_hyd + dummy(iq)
1891 coef2 = ( qe(k,i,j,
i_hi) ) &
1892 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1895 qtrc(k,i,j,(ic-1)*
nbin+iq) = 0.0_rp
1896 qtrc(k,i,j,(ip-1)*
nbin+iq) = coef2 * dummy(iq)
1897 qtrc(k,i,j,(id-1)*
nbin+iq) = 0.0_rp
1901 num_hyd_l = coef0 * n0_sdf(3) * rho_sdf(3)
1902 lambda_hyd = ( pi * rho_sdf(3) / 6.0_rp *n0_sdf(3) *
sf_gamma(4.0_rp) &
1903 / ( qe(k,i,j,
i_hs) &
1904 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hs)-eps)) &
1909 dummy(iq) = num_hyd_l * radc( iq )**3 &
1910 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1911 tmp_hyd = tmp_hyd + dummy(iq)
1914 coef2 = ( qe(k,i,j,
i_hs) ) &
1915 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1918 qtrc(k,i,j,(iss-1)*
nbin+iq) = coef2 * dummy(iq)
1922 num_hyd_l = coef0 * n0_sdf(4) * rho_sdf(4)
1923 lambda_hyd = ( pi * rho_sdf(4) / 6.0_rp *n0_sdf(4) *
sf_gamma(4.0_rp) &
1924 / ( qe(k,i,j,
i_hg) &
1925 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hg)-eps)) &
1930 dummy(iq) = num_hyd_l * radc( iq )**3 &
1931 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1932 tmp_hyd = tmp_hyd + dummy(iq)
1935 coef2 = ( qe(k,i,j,
i_hg) ) &
1936 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1939 qtrc(k,i,j,(ig-1)*
nbin+iq) = coef2 * dummy(iq)
1943 num_hyd_l = coef0 * n0_sdf(5) * rho_sdf(5)
1944 lambda_hyd = ( pi * rho_sdf(5) / 6.0_rp *n0_sdf(5) *
sf_gamma(4.0_rp) &
1945 / ( qe(k,i,j,
i_hh) &
1946 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hh)-eps)) &
1951 dummy(iq) = num_hyd_l * radc( iq )**3 &
1952 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1953 tmp_hyd = tmp_hyd + dummy(iq)
1956 coef2 = ( qe(k,i,j,
i_hh) ) &
1957 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1960 qtrc(k,i,j,(ih-1)*
nbin+iq) = coef2 * dummy(iq)
1970 do iq = num_hyd+1, qa-1
1974 qtrc(k,i,j,iq) = 0.0_rp
2001 integer,
intent(in) :: ka, ks, ke
2002 integer,
intent(in) :: ia, is, ie
2003 integer,
intent(in) :: ja, js, je
2004 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
2005 real(
rp),
intent(out) :: qecrg(ka,ia,ja,
n_hyd)
2007 integer :: ihydro, ibin, iq, icateg
2013 qecrg(:,:,:,:) = 0.0_rp
2018 iq =
nbin*(ihydro-1) + ibin
2020 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc ) )
then
2021 if ( iq > 0 .AND. iq <= nbnd )
then
2023 elseif( iq > nbnd .AND. iq <=
nbin )
then
2026 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qcl) )
then
2028 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qp ) )
then
2030 elseif( iq >
nbin*(i_mp_qd -1) .AND. iq <=
nbin*(i_mp_qd ) )
then
2032 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qs ) )
then
2034 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qg ) )
then
2036 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <= num_hyd )
then
2044 qecrg(k,i,j,icateg) = qecrg(k,i,j,icateg) + qtrc0(k,i,j,iq)
2055 subroutine mp_suzuki10( &
2085 integer,
intent(in) :: ka
2086 integer,
intent(in) :: ia
2087 integer,
intent(in) :: ja
2089 integer,
intent(in) :: ijkmax
2090 integer,
intent(in) :: ijkmax_cold
2091 integer,
intent(in) :: ijkmax_warm
2092 integer ,
intent(in) :: index_cold(ijkmax)
2093 integer ,
intent(in) :: index_warm(ijkmax)
2094 real(
rp),
intent(in) :: dens (ijkmax)
2095 real(
rp),
intent(in) :: pres (ijkmax)
2096 real(
rp),
intent(in) :: qdry (ijkmax)
2097 real(
rp),
intent(in) :: ccn (ijkmax)
2098 real(
rp),
intent(inout) :: temp (ijkmax)
2099 real(
rp),
intent(inout) :: qvap (ijkmax)
2100 real(
rp),
intent(inout) :: ghyd (
nbin,nspc,ijkmax)
2101 real(
rp),
intent(inout) :: gaer (max(
nccn,1),ijkmax)
2102 real(
rp),
intent(inout) :: cp (ijkmax)
2103 real(
rp),
intent(inout) :: cv (ijkmax)
2104 real(
rp),
intent(out) :: evaporate (ijkmax)
2105 real(
dp),
intent(in) :: dt
2108 logical,
intent(in) :: flg_lt
2109 real(
rp),
intent(in) :: d0_crg, v0_crg
2110 real(
rp),
intent(in) :: dqcrg(ijkmax), beta_crg(ijkmax)
2111 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2112 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
2119 crg_sep(:,:) = 0.0_rp
2123 if (
nccn /= 0 )
then
2124 if ( nspc == 1 )
then
2128 call nucleata( ijkmax, &
2141 call cndevpsbla( ijkmax, &
2156 if ( mp_doautoconversion )
then
2159 call collmain( ka, ia, ja, &
2174 elseif( nspc > 1 )
then
2178 call nucleata( ijkmax, &
2191 call freezing( ijkmax, &
2203 if( flg_icenucl )
then
2204 call ice_nucleat( ijkmax, &
2218 call melting( ijkmax, &
2231 call cndevpsbla( ijkmax, &
2246 if ( mp_doautoconversion )
then
2248 call collmainf( ka, ia, ja, &
2265 elseif(
nccn == 0 )
then
2267 if ( nspc == 1 )
then
2271 call nucleat( ijkmax, &
2284 call cndevpsbl( ijkmax, &
2298 if ( mp_doautoconversion )
then
2300 call collmain( ka, ia, ja, &
2315 elseif( nspc > 1 )
then
2319 call nucleat( ijkmax, &
2332 call freezing( ijkmax, &
2344 if( flg_icenucl )
then
2345 call ice_nucleat( ijkmax, &
2359 call melting( ijkmax, &
2372 call cndevpsbl( ijkmax, &
2386 if ( mp_doautoconversion )
then
2388 call collmainf( ka, ia, ja, &
2408 end subroutine mp_suzuki10
2411 subroutine nucleat( &
2424 atmos_hydrometeor_lhv_para, &
2430 atmos_saturation_pres2qsat_liq_para
2433 integer,
intent(in) :: ijkmax
2434 real(
rp),
intent(in) :: dens(ijkmax)
2435 real(
rp),
intent(in) :: pres(ijkmax)
2436 real(
rp),
intent(in) :: qdry(ijkmax)
2437 real(
rp),
intent(in) :: ccn(ijkmax)
2438 real(
rp),
intent(inout) :: temp(ijkmax)
2439 real(
rp),
intent(inout) :: qvap(ijkmax)
2440 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2441 real(
rp),
intent(inout) :: cp(ijkmax)
2442 real(
rp),
intent(inout) :: cv(ijkmax)
2443 real(
dp),
intent(in) :: dtime
2445 real(
rp) :: ssliq(ijkmax)
2446 real(
rp) :: qlevp(ijkmax)
2452 real(
rp) :: sumnum(ijkmax)
2453 real(
rp) :: gcn(
nbin,ijkmax )
2454 real(
rp) :: qsat(ijkmax)
2460 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2462 if( mp_couple_aerosol )
then
2467 dmp = ccn(ijk) * expxctr( 1 )
2468 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2469 gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
2471 qvap(ijk) = qvap(ijk) - dqv
2472 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
2479 call atmos_saturation_pres2qsat_liq_para( ijkmax, 1, ijkmax, &
2480 temp(:), pres(:), qdry(:), &
2485 ssliq(ijk) = qvap(ijk)/qsat(ijk) - 1.0_rp
2494 if ( ssliq(ijk) > 0.0_rp )
then
2499 gcn( n,ijk ) = gc( n,il,ijk )*rexpxctr( n )
2503 sumnum(ijk) = sumnum(ijk) + gcn( n,ijk )*dxmic
2505 n_c = c_ccn * ( ssliq(ijk) * 1.e+2_rp )**( kappa )
2506 if ( n_c > sumnum(ijk) )
then
2507 dmp = ( n_c - sumnum(ijk) ) * expxctr( 1 )
2508 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2509 gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
2511 qvap(ijk) = qvap(ijk) - dqv
2512 qvap(ijk) = max( qvap(ijk),0.0_rp )
2513 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
2526 end subroutine nucleat
2529 subroutine nucleata( &
2546 atmos_hydrometeor_lhv_para, &
2552 atmos_saturation_pres2qsat_liq_para
2555 integer,
intent(in) :: ijkmax
2556 real(
rp),
intent(in) :: dens(ijkmax)
2557 real(
rp),
intent(in) :: pres(ijkmax)
2558 real(
rp),
intent(in) :: qdry(ijkmax)
2559 real(
rp),
intent(inout) :: temp(ijkmax)
2560 real(
rp),
intent(inout) :: qvap(ijkmax)
2561 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2562 real(
rp),
intent(inout) :: ga (
nccn ,ijkmax)
2563 real(
rp),
intent(inout) :: cp(ijkmax)
2564 real(
rp),
intent(inout) :: cv(ijkmax)
2565 real(
dp),
intent(in) :: dtime
2569 real(
rp) :: acoef, bcoef
2572 real(
rp) :: ractr, rcld, xcld, part, dmp
2573 integer :: n, nc, ncrit
2574 integer :: ncld(
nccn)
2577 real(
rp) :: qlevp(ijkmax)
2578 real(
rp) :: qsatl(ijkmax)
2585 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2587 call atmos_saturation_pres2qsat_liq_para( ijkmax, 1, ijkmax, &
2588 temp(:), pres(:), qdry(:), &
2593 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2595 if ( ssliq <= 0.0_rp ) cycle
2600 gan( n ) = ga( n,ijk )*rexpxactr( n )
2603 acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk)
2604 bcoef = vhfct* rhoa/rhow * emwtr/emaer
2607 if ( flg_nucl )
then
2613 ractr = ( expxactr( n )*thirdovforth/pi/rhoa )**( oneovthird )
2614 rcld = sqrt( 3.0_rp*bcoef*ractr*ractr*ractr / acoef )
2615 xcld = log( rhow * 4.0_rp*pi*oneovthird*rcld*rcld*rcld )
2616 ncld( n ) = int( ( xcld-xctr( 1 ) )/dxmic ) + 1
2617 ncld( n ) = min( max( ncld( n ),1 ),
nbin )
2623 if ( ssliq <= 0.0_rp )
exit
2625 acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk)
2626 rcrit = acoef*oneovthird * ( 4.0_rp/bcoef )**( oneovthird ) / ssliq**( twoovthird )
2627 xcrit = log( rhoa * 4.0_rp*pi*oneovthird * rcrit*rcrit*rcrit )
2628 ncrit = int( ( xcrit-xabnd( 1 ) )/dxaer ) + 1
2630 if ( n == ncrit )
then
2631 part = ( xabnd( ncrit+1 )-xcrit )/dxaer
2632 elseif ( n > ncrit )
then
2640 dmp = part*gan( n )*dxaer*expxctr( nc )
2641 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2642 gc( nc,il,ijk ) = gc( nc,il,ijk ) + dmp/dxmic
2643 gan( n ) = gan( n ) - dmp/dxaer*rexpxctr( nc )
2644 gan( n ) = max( gan( n ), 0.0_rp )
2646 qvap(ijk) = qvap(ijk) - dqv
2647 qvap(ijk) = max( qvap(ijk),0.0_rp )
2648 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
2655 ga( n,ijk ) = gan( n )*expxactr( n )
2663 end subroutine nucleata
2666 subroutine cndevpsbl( &
2682 integer,
intent(in) :: ijkmax
2683 logical,
intent(in) :: flg_lt
2684 real(
rp),
intent(in) :: dens(ijkmax)
2685 real(
rp),
intent(in) :: pres(ijkmax)
2686 real(
rp),
intent(in) :: qdry(ijkmax)
2687 real(
rp),
intent(inout) :: temp(ijkmax)
2688 real(
rp),
intent(inout) :: qvap(ijkmax)
2689 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2690 real(
rp),
intent(inout) :: cp(ijkmax)
2691 real(
rp),
intent(inout) :: cv(ijkmax)
2692 real(
rp),
intent(out) :: evaporate(ijkmax)
2693 real(
dp),
intent(in) :: dtime
2694 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2697 call liqphase( ijkmax, &
2713 call icephase( ijkmax, &
2726 call mixphase( ijkmax, &
2742 end subroutine cndevpsbl
2745 subroutine cndevpsbla( &
2762 integer,
intent(in) :: ijkmax
2763 logical,
intent(in) :: flg_lt
2764 real(
rp),
intent(in) :: dens(ijkmax)
2765 real(
rp),
intent(in) :: pres(ijkmax)
2766 real(
rp),
intent(in) :: qdry(ijkmax)
2767 real(
rp),
intent(inout) :: temp(ijkmax)
2768 real(
rp),
intent(inout) :: qvap(ijkmax)
2769 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2770 real(
rp),
intent(inout) :: ga (
nccn ,ijkmax)
2771 real(
rp),
intent(inout) :: cp(ijkmax)
2772 real(
rp),
intent(inout) :: cv(ijkmax)
2773 real(
rp),
intent(out) :: evaporate(ijkmax)
2774 real(
dp),
intent(in) :: dtime
2775 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2778 call liqphase( ijkmax, &
2793 if ( flg_regeneration )
then
2794 call faero( ijkmax, &
2801 call icephase( ijkmax, &
2814 call mixphase( ijkmax, &
2829 end subroutine cndevpsbla
2832 subroutine liqphase( &
2854 atmos_hydrometeor_lhv_para, &
2860 atmos_saturation_pres2qsat_liq_para
2863 integer,
intent(in) :: ijkmax
2864 logical,
intent(in) :: flg_lt
2865 real(
rp),
intent(in) :: dens(ijkmax)
2866 real(
rp),
intent(in) :: pres(ijkmax)
2867 real(
rp),
intent(in) :: qdry(ijkmax)
2868 real(
rp),
intent(inout) :: temp(ijkmax)
2869 real(
rp),
intent(inout) :: qvap(ijkmax)
2870 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2871 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2872 real(
rp),
intent(inout) :: cp(ijkmax)
2873 real(
rp),
intent(inout) :: cv(ijkmax)
2874 real(
rp),
intent(out) :: regene_gcn(ijkmax)
2875 real(
dp),
intent(in) :: dtime
2877 integer :: n, myu, ncount
2878 integer :: nloop(ijkmax)
2879 real(
rp) :: gclold(ijkmax)
2880 real(
rp) :: gclnew(ijkmax)
2882 real(
rp) :: dtcnd(ijkmax)
2883 real(
rp) :: gtliq(ijkmax)
2884 real(
rp) :: qlevp(ijkmax)
2885 real(
rp) :: cefliq, a, sliqtnd
2886 real(
rp) :: sumliq(ijkmax), umax(ijkmax)
2888 real(
rp) :: gcn( -1:
nbin+2,nspc,ijkmax )
2890 real(
rp),
parameter :: cflfct = 0.50_rp
2892 integer :: iflg( nspc,ijkmax )
2894 real(
rp) :: f1, f2, emu, cefd, cefk, festl
2895 real(
rp) :: qsatl(ijkmax)
2896 real(
rp),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2897 real(
rp),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2898 real(
rp) :: zerosw, qvtmp
2900 integer :: ijk, nn, mm, loopflg(ijkmax)
2903 real(
rp) :: uadv ( 0:
nbin+2,nspc,ijkmax )
2904 real(
rp) :: crn ( 0:
nbin+2,nspc,ijkmax )
2905 real(
rp) :: flq ( 1:
nbin+1 )
2906 real(
rp) :: acoef( 0:
nbin+1,0:2 )
2907 real(
rp) :: aip ( 0:
nbin+1 )
2908 real(
rp) :: aim ( 0:
nbin+1 )
2909 real(
rp) :: ai ( 0:
nbin+1 )
2910 real(
rp) :: cmins, cplus
2914 real(
rp) :: gcrgn( -1:
nbin+2,nspc,ijkmax )
2915 real(
rp) :: flq_c( 1:
nbin+1 )
2916 real(
rp) :: acoef_c( 0:
nbin+1,0:2 )
2917 real(
rp) :: aip_c ( 0:
nbin+1 )
2918 real(
rp) :: aim_c ( 0:
nbin+1 )
2919 real(
rp) :: ai_c ( 0:
nbin+1 )
2930 csum = csum + gc( n,il,ijk )*dxmic
2932 if( csum > cldmin ) iflg( il,ijk ) = 1
2936 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2938 call atmos_saturation_pres2qsat_liq_para( ijkmax, 1, ijkmax, &
2939 temp(:), pres(:), qdry(:), &
2950 gclold(ijk) = gclold(ijk) + gc( n,il,ijk ) * dxmic
2955 gcn( n,il,ijk ) = gc( n,il,ijk ) * rexpxctr( n )
2956 gcrgn( n,il,ijk ) = gcrg( n,il,ijk )
2958 gcn( -1,il,ijk ) = 0.0_rp
2959 gcn( 0,il,ijk ) = 0.0_rp
2960 gcn(
nbin+1,il,ijk ) = 0.0_rp
2961 gcn(
nbin+2,il,ijk ) = 0.0_rp
2962 gcrgn( -1,il,ijk ) = 0.0_rp
2963 gcrgn( 0,il,ijk ) = 0.0_rp
2964 gcrgn(
nbin+1,il,ijk ) = 0.0_rp
2965 gcrgn(
nbin+2,il,ijk ) = 0.0_rp
2968 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2970 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2971 cefd = emu/dens(ijk)
2974 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2975 f1 = rvap*temp(ijk)/festl/cefd
2976 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2977 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2979 umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq )
2980 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2981 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2982 dtcnd(ijk) = dtime / nloop(ijk)
2984 nloop(ijk) = nloop(ijk) * iflg( il,ijk )
2986 nloopmax = maxval(nloop,1)
2990 regene_gcn(:) = 0.0_rp
2994 do ncount = 1, nloopmax
2998 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
3002 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3007 do nn = 1, loopflg(ijk)
3009 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3010 cefd = emu/dens(ijk)
3012 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3013 f1 = rvap*temp(ijk)/festl/cefd
3014 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3015 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3017 sumliq(ijk) = 0.0_rp
3020 sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
3027 call atmos_saturation_pres2qsat_liq_para( ijkmax, 1, ijkmax, &
3028 temp(:), pres(:), qdry(:), &
3033 do nn = 1, loopflg(ijk)
3035 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3038 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
3039 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
3040 cefliq = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)*qlevp(ijk)/cv(ijk)/rvap/temp(ijk)/temp(ijk) )
3041 a = - cefliq*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3042 a = a + eps * ( 1.0_rp - zerosw )
3044 sliqtnd = zerosw * &
3045 ( ssliq * ( exp( a*dtcnd(ijk) )-1.0_rp )/( a*dtcnd(ijk) ) &
3046 * ( 0.5_rp + sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
3048 * ( 0.5_rp - sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
3050 + ssliq * ( 1.0_rp - zerosw )
3053 do mm = 1, iflg( il,ijk )
3056 uadv( n,il,ijk ) = cbnd( n,il )*rexpxbnd( n )*gtliq(ijk)*sliqtnd
3058 uadv( 0 ,il,ijk ) = 0.0_rp
3059 uadv(
nbin+2,il,ijk ) = 0.0_rp
3073 do nn = 1, loopflg(ijk)
3075 do mm = 1, iflg( myu,ijk )
3077 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
3081 acoef(n,0) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3082 acoef(n,1) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3083 acoef(n,2) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3087 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3089 aip(n) = acoef(n,0) * ( 1.0_rp-cplus**1 ) &
3090 + acoef(n,1) * ( 1.0_rp-cplus**2 ) &
3091 + acoef(n,2) * ( 1.0_rp-cplus**3 )
3092 aip(n) = max( aip(n), 0.0_rp )
3096 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3098 aim(n) = acoef(n,0) * ( 1.0_rp-cmins**1 ) &
3099 - acoef(n,1) * ( 1.0_rp-cmins**2 ) &
3100 + acoef(n,2) * ( 1.0_rp-cmins**3 )
3101 aim(n) = max( aim(n), 0.0_rp )
3105 ai(n) = acoef(n,0) * 2.0_rp &
3106 + acoef(n,2) * 2.0_rp
3107 ai(n) = max( ai(n), aip(n)+aim(n)+cldmin )
3111 flq(n) = ( aip(n-1)/ai(n-1)*gcn( n-1,myu,ijk ) &
3112 - aim(n )/ai(n )*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3116 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1)-flq(n) )*dtcnd(ijk)/dxmic
3119 regene_gcn(ijk) = regene_gcn(ijk)+( -flq(1)*dtcnd(ijk)/dxmic ) &
3120 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
3131 do nn = 1, loopflg(ijk)
3133 do mm = 1, iflg( myu,ijk )
3135 acoef_c(n,0) = - ( gcrgn( n+1,myu,ijk )-26.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3136 acoef_c(n,1) = ( gcrgn( n+1,myu,ijk ) -gcrgn( n-1,myu,ijk ) ) / 16.0_rp
3137 acoef_c(n,2) = ( gcrgn( n+1,myu,ijk )- 2.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3141 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3143 aip_c(n) = acoef_c(n,0) * ( 1.0_rp-cplus**1 ) &
3144 + acoef_c(n,1) * ( 1.0_rp-cplus**2 ) &
3145 + acoef_c(n,2) * ( 1.0_rp-cplus**3 )
3146 aip_c(n) = max( aip_c(n), 0.0_rp )
3150 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3152 aim_c(n) = acoef_c(n,0) * ( 1.0_rp-cmins**1 ) &
3153 - acoef_c(n,1) * ( 1.0_rp-cmins**2 ) &
3154 + acoef_c(n,2) * ( 1.0_rp-cmins**3 )
3155 aim_c(n) = max( aim_c(n), 0.0_rp )
3159 ai_c(n) = acoef_c(n,0) * 2.0_rp &
3160 + acoef_c(n,2) * 2.0_rp
3161 ai_c(n) = max( ai_c(n), aip_c(n)+aim_c(n)+cldmin )
3165 flq_c(n) = ( aip_c(n-1)/ai_c(n-1)*gcrgn( n-1,myu,ijk ) &
3166 - aim_c(n )/ai_c(n )*gcrgn( n ,myu,ijk ) )/dtcnd(ijk)
3170 gcrgn( n,myu,ijk ) = gcrgn( n,myu,ijk ) - ( flq_c(n+1)-flq_c(n) )*dtcnd(ijk)
3184 do nn = 1, loopflg(ijk)
3187 gclnew(ijk) = 0.0_rp
3190 gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )
3195 gclnew(ijk) = gclnew(ijk)*dxmic
3198 cndmss = gclnew(ijk) - gclold(ijk)
3199 dqv = cndmss/dens(ijk)
3200 qvap(ijk) = qvap(ijk) - dqv
3201 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
3205 gclold(ijk) = gclnew(ijk)
3217 gc( n,il,ijk ) = gcn( n,il,ijk )*expxctr( n )
3225 gcrg( n,il,ijk ) = gcrgn( n,il,ijk )
3231 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3236 if ( gc( n,il,ijk ) < 0.0_rp )
then
3237 cndmss = -gc( n,il,ijk )
3238 gc( n,il,ijk ) = 0.0_rp
3239 gcrg( n,il,ijk ) = 0.0_rp
3240 dqv = cndmss/dens(ijk)
3241 qvap(ijk) = qvap(ijk) + dqv
3242 temp(ijk) = temp(ijk) - dqv*qlevp(ijk)/cv(ijk)
3252 end subroutine liqphase
3255 subroutine icephase( &
3276 atmos_hydrometeor_lhs_para, &
3282 atmos_saturation_pres2qsat_ice_para
3285 integer,
intent(in) :: ijkmax
3286 logical,
intent(in) :: flg_lt
3287 real(
rp),
intent(in) :: dens(ijkmax)
3288 real(
rp),
intent(in) :: pres(ijkmax)
3289 real(
rp),
intent(in) :: qdry(ijkmax)
3290 real(
rp),
intent(inout) :: temp(ijkmax)
3291 real(
rp),
intent(inout) :: qvap(ijkmax)
3292 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3293 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
3294 real(
rp),
intent(inout) :: cp(ijkmax)
3295 real(
rp),
intent(inout) :: cv(ijkmax)
3296 real(
dp),
intent(in) :: dtime
3298 integer :: myu, n, ncount
3299 integer :: nloop(ijkmax)
3300 real(
rp) :: gciold(ijkmax), gcinew(ijkmax)
3301 real(
rp) :: dtcnd(ijkmax)
3303 real(
rp) :: gtice(ijkmax), umax(ijkmax)
3304 real(
rp) :: qlsbl(ijkmax)
3305 real(
rp) :: cefice, d, uval, sicetnd
3306 real(
rp) :: sumice(ijkmax)
3308 real(
rp) :: gcn( -1:
nbin+2,nspc,ijkmax )
3309 real(
rp),
parameter :: cflfct = 0.50_rp
3311 integer :: iflg( nspc,ijkmax )
3313 real(
rp) :: f1, f2, emu, cefd, cefk, festi
3314 real(
rp) :: qsati(ijkmax)
3315 real(
rp),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
3316 real(
rp),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
3317 integer :: ijk, mm, nn, loopflg(ijkmax)
3318 real(
rp) :: zerosw, qvtmp
3323 real(
rp) :: uadv ( 0:
nbin+2,nspc,ijkmax )
3324 real(
rp) :: crn ( 0:
nbin+2,nspc,ijkmax )
3325 real(
rp) :: flq ( 1:
nbin+1 )
3326 real(
rp) :: acoef( 0:
nbin+1,0:2 )
3327 real(
rp) :: aip ( 0:
nbin+1 )
3328 real(
rp) :: aim ( 0:
nbin+1 )
3329 real(
rp) :: ai ( 0:
nbin+1 )
3330 real(
rp) :: cmins, cplus
3334 real(
rp) :: gcrgn( -1:
nbin+2,nspc,ijkmax )
3335 real(
rp) :: flq_c( 1:
nbin+1 )
3336 real(
rp) :: acoef_c( 0:
nbin+1,0:2 )
3337 real(
rp) :: aip_c ( 0:
nbin+1 )
3338 real(
rp) :: aim_c ( 0:
nbin+1 )
3339 real(
rp) :: ai_c ( 0:
nbin+1 )
3353 csum = csum + gc( n,myu,ijk )*dxmic
3355 if ( csum > cldmin ) iflg( myu,ijk ) = 1
3361 call atmos_hydrometeor_lhs_para( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3363 call atmos_saturation_pres2qsat_ice_para( ijkmax, 1, ijkmax, &
3364 temp(:), pres(:), qdry(:), &
3377 gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
3384 gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
3386 gcn( -1,myu,ijk ) = 0.0_rp
3387 gcn( 0,myu,ijk ) = 0.0_rp
3388 gcn(
nbin+1,myu,ijk ) = 0.0_rp
3389 gcn(
nbin+2,myu,ijk ) = 0.0_rp
3394 gcrgn( n,myu,ijk ) = gcrg( n,myu,ijk )
3396 gcrgn( -1,myu,ijk ) = 0.0_rp
3397 gcrgn( 0,myu,ijk ) = 0.0_rp
3398 gcrgn(
nbin+1,myu,ijk ) = 0.0_rp
3399 gcrgn(
nbin+2,myu,ijk ) = 0.0_rp
3405 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3407 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3408 cefd = emu/dens(ijk)
3410 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3411 f1 = rvap*temp(ijk)/festi/cefd
3412 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3413 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3417 uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice )
3418 umax(ijk) = max( umax(ijk),uval )
3421 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
3422 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
3423 dtcnd(ijk) = dtime/nloop(ijk)
3425 nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) )
3428 nloopmax = maxval(nloop,1)
3433 do ncount = 1, nloopmax
3437 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
3441 call atmos_hydrometeor_lhs_para( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3446 do nn = 1, loopflg(ijk)
3447 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3448 cefd = emu/dens(ijk)
3450 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3451 f1 = rvap*temp(ijk)/festi/cefd
3452 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3453 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3455 sumice(ijk) = 0.0_rp
3458 sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
3466 call atmos_saturation_pres2qsat_ice_para( ijkmax, 1, ijkmax, &
3467 temp(:), pres(:), qdry(:), &
3473 do nn = 1, loopflg(ijk)
3476 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3479 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
3480 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
3482 cefice = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)*qlsbl(ijk)/cv(ijk)/rvap/temp(ijk)/temp(ijk) )
3483 d = - cefice*sumice(ijk)*gtice(ijk)/dens(ijk)
3484 d = d + eps * ( 1.0_rp - zerosw )
3486 sicetnd = zerosw * &
3487 ( ssice * ( exp( d*dtcnd(ijk) )-1.0_rp )/( d*dtcnd(ijk) ) &
3488 * ( 0.5_rp + sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
3490 * ( 0.5_rp - sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
3492 + ssice * ( 1.0_rp - zerosw )
3496 do mm = 1, iflg( myu,ijk )
3499 uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd
3501 uadv( 0, myu,ijk ) = 0.0_rp
3502 uadv(
nbin+2,myu,ijk ) = 0.0_rp
3513 do nn = 1, loopflg(ijk)
3515 do mm = 1, iflg( myu,ijk )
3517 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
3521 acoef(n,0) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3522 acoef(n,1) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3523 acoef(n,2) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3527 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3529 aip(n) = acoef(n,0) * ( 1.0_rp-cplus**1 ) &
3530 + acoef(n,1) * ( 1.0_rp-cplus**2 ) &
3531 + acoef(n,2) * ( 1.0_rp-cplus**3 )
3532 aip(n) = max( aip(n), 0.0_rp )
3536 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3538 aim(n) = acoef(n,0) * ( 1.0_rp-cmins**1 ) &
3539 - acoef(n,1) * ( 1.0_rp-cmins**2 ) &
3540 + acoef(n,2) * ( 1.0_rp-cmins**3 )
3541 aim(n) = max( aim(n), 0.0_rp )
3545 ai(n) = acoef(n,0) * 2.0_rp &
3546 + acoef(n,2) * 2.0_rp
3547 ai(n) = max( ai(n), aip(n)+aim(n)+cldmin )
3551 flq(n) = ( aip(n-1)/ai(n-1)*gcn( n-1,myu,ijk ) &
3552 - aim(n )/ai(n )*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3556 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1)-flq(n) )*dtcnd(ijk)/dxmic
3571 do nn = 1, loopflg(ijk)
3573 do mm = 1, iflg( myu,ijk )
3575 acoef_c(n,0) = - ( gcrgn( n+1,myu,ijk )-26.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3576 acoef_c(n,1) = ( gcrgn( n+1,myu,ijk ) -gcrgn( n-1,myu,ijk ) ) / 16.0_rp
3577 acoef_c(n,2) = ( gcrgn( n+1,myu,ijk )- 2.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3581 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3583 aip_c(n) = acoef_c(n,0) * ( 1.0_rp-cplus**1 ) &
3584 + acoef_c(n,1) * ( 1.0_rp-cplus**2 ) &
3585 + acoef_c(n,2) * ( 1.0_rp-cplus**3 )
3586 aip_c(n) = max( aip_c(n), 0.0_rp )
3590 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3592 aim_c(n) = acoef_c(n,0) * ( 1.0_rp-cmins**1 ) &
3593 - acoef_c(n,1) * ( 1.0_rp-cmins**2 ) &
3594 + acoef_c(n,2) * ( 1.0_rp-cmins**3 )
3595 aim_c(n) = max( aim_c(n), 0.0_rp )
3599 ai_c(n) = acoef_c(n,0) * 2.0_rp &
3600 + acoef_c(n,2) * 2.0_rp
3601 ai_c(n) = max( ai_c(n), aip_c(n)+aim_c(n)+cldmin )
3605 flq_c(n) = ( aip_c(n-1)/ai_c(n-1)*gcrgn( n-1,myu,ijk ) &
3606 - aim_c(n )/ai_c(n )*gcrgn( n ,myu,ijk ) )/dtcnd(ijk)
3610 gcrgn( n,myu,ijk ) = gcrgn( n,myu,ijk ) - ( flq_c(n+1)-flq_c(n) )*dtcnd(ijk)
3624 do nn = 1, loopflg(ijk)
3627 gcinew(ijk) = 0.0_rp
3630 gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
3635 sblmss = gcinew(ijk) - gciold(ijk)
3636 dqv = sblmss/dens(ijk)
3637 qvap(ijk) = qvap(ijk) - dqv
3638 temp(ijk) = temp(ijk) + dqv*qlsbl(ijk)/cv(ijk)
3642 gciold(ijk) = gcinew(ijk)
3657 gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
3667 gcrg( n,myu,ijk ) = gcrgn( n,myu,ijk )
3676 end subroutine icephase
3679 subroutine mixphase( &
3700 atmos_hydrometeor_lhv_para, &
3701 atmos_hydrometeor_lhs_para, &
3709 atmos_saturation_pres2qsat_liq_para, &
3710 atmos_saturation_pres2qsat_ice_para
3713 integer,
intent(in) :: ijkmax
3714 logical,
intent(in) :: flg_lt
3715 real(
rp),
intent(in) :: dens(ijkmax)
3716 real(
rp),
intent(in) :: pres(ijkmax)
3717 real(
rp),
intent(in) :: qdry(ijkmax)
3718 real(
rp),
intent(inout) :: temp(ijkmax)
3719 real(
rp),
intent(inout) :: qvap(ijkmax)
3720 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3721 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
3722 real(
rp),
intent(inout) :: cp(ijkmax)
3723 real(
rp),
intent(inout) :: cv(ijkmax)
3724 real(
dp),
intent(in) :: dtime
3726 integer :: n, myu, mm, ncount
3727 integer :: nloop(ijkmax)
3728 real(
rp) :: gclold(ijkmax), gclnew(ijkmax)
3729 real(
rp) :: gciold(ijkmax), gcinew(ijkmax)
3730 real(
rp) :: cndmss, sblmss
3731 real(
rp) :: gtliq(ijkmax), gtice(ijkmax)
3732 real(
rp) :: umax(ijkmax), uval, dtcnd(ijkmax)
3733 real(
rp) :: qlevp(ijkmax), qlsbl(ijkmax)
3734 real(
rp) :: cef1, cef2, cef3, cef4, a, b, c, d
3735 real(
rp) :: rmdplus, rmdmins, ssplus, ssmins, tplus, tmins
3736 real(
rp) :: sliqtnd, sicetnd
3737 real(
rp) :: ssliq, ssice
3738 real(
rp) :: sumliq(ijkmax), sumice(ijkmax)
3739 real(
rp) :: gcn( -1:
nbin+2,nspc,ijkmax )
3740 real(
rp),
parameter :: cflfct = 0.50_rp
3743 integer :: iflg( nspc,ijkmax )
3744 real(
rp) :: f1, f2, emu, cefd, cefk, festl, festi
3745 real(
rp) :: qsatl(ijkmax), qsati(ijkmax)
3746 real(
rp),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
3747 real(
rp),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
3748 real(
rp) :: qvtmp, zerosw
3749 integer :: ijk, nn, loopflg(ijkmax)
3753 real(
rp) :: uadv ( 0:
nbin+2,nspc,ijkmax )
3754 real(
rp) :: crn ( 0:
nbin+2,nspc,ijkmax )
3755 real(
rp) :: flq ( 1:
nbin+1 )
3756 real(
rp) :: acoef( 0:
nbin+1,0:2 )
3757 real(
rp) :: aip ( 0:
nbin+1 )
3758 real(
rp) :: aim ( 0:
nbin+1 )
3759 real(
rp) :: ai ( 0:
nbin+1 )
3760 real(
rp) :: cmins, cplus
3764 real(
rp) :: gcrgn( -1:
nbin+2,nspc,ijkmax )
3765 real(
rp) :: flq_c( 1:
nbin+1 )
3766 real(
rp) :: acoef_c( 0:
nbin+1,0:2 )
3767 real(
rp) :: aip_c ( 0:
nbin+1 )
3768 real(
rp) :: aim_c ( 0:
nbin+1 )
3769 real(
rp) :: ai_c ( 0:
nbin+1 )
3783 csum = csum + gc( n,myu,ijk )*dxmic
3785 if ( csum > cldmin ) iflg( myu,ijk ) = 1
3791 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3793 call atmos_hydrometeor_lhs_para( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3795 call atmos_saturation_pres2qsat_liq_para( ijkmax, 1, ijkmax, &
3796 temp(:), pres(:), qdry(:), &
3798 call atmos_saturation_pres2qsat_ice_para( ijkmax, 1, ijkmax, &
3799 temp(:), pres(:), qdry(:), &
3811 gclold(ijk) = gclold(ijk) + gc( n,il,ijk )*dxmic
3816 gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
3823 gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
3824 gcrgn( n,myu,ijk ) = gcrg( n,myu,ijk )
3826 gcn( -1,myu,ijk ) = 0.0_rp
3827 gcn( 0,myu,ijk ) = 0.0_rp
3828 gcn(
nbin+1,myu,ijk ) = 0.0_rp
3829 gcn(
nbin+2,myu,ijk ) = 0.0_rp
3830 gcrgn( -1,myu,ijk ) = 0.0_rp
3831 gcrgn( 0,myu,ijk ) = 0.0_rp
3832 gcrgn(
nbin+1,myu,ijk ) = 0.0_rp
3833 gcrgn(
nbin+2,myu,ijk ) = 0.0_rp
3837 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3838 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3840 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3841 cefd = emu/dens(ijk)
3844 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3845 f1 = rvap*temp(ijk)/festl/cefd
3846 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3847 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3849 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3850 f1 = rvap*temp(ijk)/festi/cefd
3851 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3852 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3855 umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq )
3857 uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice )
3858 umax(ijk) = max( umax(ijk),uval )
3861 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
3862 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
3863 dtcnd(ijk) = dtime/nloop(ijk)
3865 nloop(ijk) = nloop(ijk) * iflg( il,ijk )
3866 nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) )
3869 nloopmax = maxval(nloop,1)
3874 do ncount = 1, nloopmax
3878 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
3882 call atmos_hydrometeor_lhv_para( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3884 call atmos_hydrometeor_lhs_para( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3890 do nn = 1, loopflg(ijk)
3893 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3894 cefd = emu/dens(ijk)
3896 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3897 f1 = rvap*temp(ijk)/festl/cefd
3898 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3899 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3900 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3901 f1 = rvap*temp(ijk)/festi/cefd
3902 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3903 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3905 sumliq(ijk) = 0.0_rp
3907 sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
3910 sumice(ijk) = 0.0_rp
3913 sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
3920 call atmos_saturation_pres2qsat_liq_para( ijkmax, 1, ijkmax, &
3921 temp(:), pres(:), qdry(:), &
3923 call atmos_saturation_pres2qsat_ice_para( ijkmax, 1, ijkmax, &
3924 temp(:), pres(:), qdry(:), &
3929 do nn = 1, loopflg(ijk)
3932 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3933 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3935 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
3936 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
3937 cef1 = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cv(ijk) )
3938 cef2 = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cv(ijk) )
3939 cef3 = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cv(ijk) )
3940 cef4 = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cv(ijk) )
3942 a = - cef1*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3943 b = - cef2*sumice(ijk)*gtice(ijk)/dens(ijk)
3944 c = - cef3*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3945 d = - cef4*sumice(ijk)*gtice(ijk)/dens(ijk)
3947 b = b + eps * ( 1.0_rp - zerosw )
3949 rmdplus = ( ( a+d ) + sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
3950 rmdmins = ( ( a+d ) - sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
3952 rmdplus = rmdplus + eps * ( 1.0_rp - zerosw )
3953 rmdmins = rmdmins + eps * ( 1.0_rp - zerosw )
3956 ssplus = ( ( rmdmins-a )*ssliq - b*ssice )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
3957 ssmins = ( ( a-rmdplus )*ssliq + b*ssice )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
3959 tplus = ( exp( rmdplus*dtcnd(ijk) )-1.0_rp )/( rmdplus*dtcnd(ijk) ) &
3960 * ( 0.5_rp + sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) ) &
3962 * ( 0.5_rp - sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) )
3963 tmins = ( exp( rmdmins*dtcnd(ijk) )-1.0_rp )/( rmdmins*dtcnd(ijk) ) &
3964 * ( 0.5_rp + sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) ) &
3966 * ( 0.5_rp - sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) )
3968 sliqtnd = ssliq * ( 1.0_rp - zerosw ) &
3970 ( b*tplus*ssplus + b*tmins*ssmins )
3971 sicetnd = ssice * ( 1.0_rp - zerosw ) &
3973 ( ( rmdplus-a )*tplus*ssplus &
3974 + ( rmdmins-a )*tmins*ssmins )
3978 do mm = 1, iflg( myu,ijk )
3982 uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtliq(ijk)*sliqtnd &
3983 * ( 0.5_rp - sign( 0.5_rp,real(myu)-1.5_rp) ) &
3984 + cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd &
3985 * ( 0.5_rp + sign( 0.5_rp,real(myu)-1.5_rp) )
3987 uadv( 0, myu,ijk ) = 0.0_rp
3988 uadv(
nbin+2,myu,ijk ) = 0.0_rp
3998 do nn = 1, loopflg(ijk)
4000 do mm = 1, iflg( myu,ijk )
4002 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
4006 acoef(n,0) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
4007 acoef(n,1) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
4008 acoef(n,2) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
4012 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
4014 aip(n) = acoef(n,0) * ( 1.0_rp-cplus**1 ) &
4015 + acoef(n,1) * ( 1.0_rp-cplus**2 ) &
4016 + acoef(n,2) * ( 1.0_rp-cplus**3 )
4017 aip(n) = max( aip(n), 0.0_rp )
4021 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
4023 aim(n) = acoef(n,0) * ( 1.0_rp-cmins**1 ) &
4024 - acoef(n,1) * ( 1.0_rp-cmins**2 ) &
4025 + acoef(n,2) * ( 1.0_rp-cmins**3 )
4026 aim(n) = max( aim(n), 0.0_rp )
4030 ai(n) = acoef(n,0) * 2.0_rp &
4031 + acoef(n,2) * 2.0_rp
4032 ai(n) = max( ai(n), aip(n)+aim(n)+cldmin )
4036 flq(n) = ( aip(n-1)/ai(n-1)*gcn( n-1,myu,ijk ) &
4037 - aim(n )/ai(n )*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
4041 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1)-flq(n) )*dtcnd(ijk)/dxmic
4063 do nn = 1, loopflg(ijk)
4065 do mm = 1, iflg( myu,ijk )
4067 acoef_c(n,0) = - ( gcrgn( n+1,myu,ijk )-26.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
4068 acoef_c(n,1) = ( gcrgn( n+1,myu,ijk ) -gcrgn( n-1,myu,ijk ) ) / 16.0_rp
4069 acoef_c(n,2) = ( gcrgn( n+1,myu,ijk )- 2.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
4073 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
4075 aip_c(n) = acoef_c(n,0) * ( 1.0_rp-cplus**1 ) &
4076 + acoef_c(n,1) * ( 1.0_rp-cplus**2 ) &
4077 + acoef_c(n,2) * ( 1.0_rp-cplus**3 )
4078 aip_c(n) = max( aip_c(n), 0.0_rp )
4082 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
4084 aim_c(n) = acoef_c(n,0) * ( 1.0_rp-cmins**1 ) &
4085 - acoef_c(n,1) * ( 1.0_rp-cmins**2 ) &
4086 + acoef_c(n,2) * ( 1.0_rp-cmins**3 )
4087 aim_c(n) = max( aim_c(n), 0.0_rp )
4091 ai_c(n) = acoef_c(n,0) * 2.0_rp &
4092 + acoef_c(n,2) * 2.0_rp
4093 ai_c(n) = max( ai_c(n), aip_c(n)+aim_c(n)+cldmin )
4097 flq_c(n) = ( aip_c(n-1)/ai_c(n-1)*gcrgn( n-1,myu,ijk ) &
4098 - aim_c(n )/ai_c(n )*gcrgn( n ,myu,ijk ) )/dtcnd(ijk)
4102 gcrgn( n,myu,ijk ) = gcrgn( n,myu,ijk ) - ( flq_c(n+1)-flq_c(n) )*dtcnd(ijk)
4115 do nn = 1, loopflg(ijk)
4117 gclnew(ijk) = 0.0_rp
4119 gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )*dxmic
4122 gcinew(ijk) = 0.0_rp
4125 gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
4130 cndmss = gclnew(ijk) - gclold(ijk)
4131 sblmss = gcinew(ijk) - gciold(ijk)
4133 qvap(ijk) = qvap(ijk) - ( cndmss + sblmss ) / dens(ijk)
4134 temp(ijk) = temp(ijk) + ( cndmss*qlevp(ijk)+sblmss*qlsbl(ijk) ) / dens(ijk) / cv(ijk)
4138 gclold(ijk) = gclnew(ijk)
4139 gciold(ijk) = gcinew(ijk)
4150 gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
4161 gcrg( n,myu,ijk ) = gcrgn( n,myu,ijk )
4174 end subroutine mixphase
4177 subroutine ice_nucleat( &
4193 atmos_saturation_pres2qsat_ice_para
4195 atmos_hydrometeor_lhs_para, &
4202 integer,
intent(in) :: ijkmax
4203 integer,
intent(in) :: num_cold
4204 integer,
intent(in) :: index_cold(ijkmax)
4205 real(
rp),
intent(in) :: dens(ijkmax)
4206 real(
rp),
intent(in) :: pres(ijkmax)
4207 real(
rp),
intent(in) :: qdry(ijkmax)
4208 real(
rp),
intent(inout) :: temp(ijkmax)
4209 real(
rp),
intent(inout) :: qvap(ijkmax)
4210 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4211 real(
rp),
intent(inout) :: cp(ijkmax)
4212 real(
rp),
intent(inout) :: cv(ijkmax)
4213 real(
dp),
intent(in) :: dtime
4216 real(
rp) :: numin, tdel, qdel
4218 real(
rp),
parameter :: acoef = -0.639_rp, bcoef = 12.96_rp
4220 real(
rp),
parameter :: tcolmu = 269.0_rp, tcolml = 265.0_rp
4221 real(
rp),
parameter :: tdendu = 257.0_rp, tdendl = 255.0_rp
4222 real(
rp),
parameter :: tplatu = 250.6_rp
4224 real(
rp) :: qsati(ijkmax), qlsbl(ijkmax)
4225 integer :: ijk, indirect
4232 call atmos_saturation_pres2qsat_ice_para( ijkmax, 1, ijkmax, &
4233 temp(:), pres(:), qdry(:), &
4236 call atmos_hydrometeor_lhs_para( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
4239 do indirect = 1, num_cold
4240 ijk = index_cold(indirect)
4243 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
4245 if( ssice <= 0.0_rp ) cycle
4250 numice = numice + gc( n,myu,ijk )*rexpxctr( n )*dxmic
4254 numin = n0_icenucl * exp( acoef + bcoef * ssice * 1.e+2_rp )
4255 if( numin > numice )
then
4257 if ( temp(ijk) <= tplatu .OR. ( temp(ijk) >= tcolml .AND. temp(ijk) < tcolmu ) )
then
4260 elseif( temp(ijk) <= tdendu .AND. temp(ijk) >= tdendl )
then
4267 numin = (numin-numice) * expxctr( 1 )
4268 numin = min( numin,qvap(ijk)*dens(ijk) )
4269 gc( 1,ispc,ijk ) = gc( 1,ispc,ijk ) + numin / dxmic
4271 tdel = numin/dens(ijk)*qlsbl(ijk)/cv(ijk)
4272 temp(ijk) = temp(ijk) + tdel
4273 qdel = numin/dens(ijk)
4274 qvap(ijk) = qvap(ijk) - qdel
4285 end subroutine ice_nucleat
4288 subroutine freezing( &
4306 atmos_hydrometeor_lhf_para, &
4313 integer,
intent(in) :: ijkmax
4314 integer,
intent(in) :: num_cold
4315 integer,
intent(in) :: index_cold(ijkmax)
4316 logical,
intent(in) :: flg_lt
4317 real(
rp),
intent(in) :: dens(ijkmax)
4318 real(
rp),
intent(inout) :: temp(ijkmax)
4319 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4320 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4321 real(
rp),
intent(inout) :: cp(ijkmax)
4322 real(
rp),
intent(inout) :: cv(ijkmax)
4323 real(
dp),
intent(in) :: dtime
4325 integer :: nbound, n
4326 real(
rp) :: xbound, tc, rate, dmp, frz, sumfrz, tdel
4327 real(
rp),
parameter :: coefa = 1.0e-01_rp
4328 real(
rp),
parameter :: coefb = 0.66_rp
4329 real(
rp),
parameter :: rbound = 2.0e-04_rp
4333 integer :: ijk, indirect
4334 real(
rp) :: qlmlt(ijkmax)
4339 call atmos_hydrometeor_lhf_para( ijkmax, 1, ijkmax, temp(:), qlmlt(:) )
4341 xbound = log( rhow * 4.0_rp*pi/3.0_rp * rbound**3 )
4342 nbound = int( ( xbound-xbnd( 1 ) )/dxmic ) + 1
4345 do indirect = 1, num_cold
4346 ijk = index_cold(indirect)
4351 rate = coefa*exp( -coefb*tc )
4356 dmp = rate*expxctr( n )
4357 frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
4358 frz = min( frz, gc( n,il,ijk ) )
4361 dcrg = frz / ( gc(n,il,ijk)+eps ) * gcrg( n,il,ijk )
4362 gcrg( n,il,ijk ) = gcrg( n,il,ijk ) - dcrg
4363 gcrg( n,ip,ijk ) = gcrg( n,ip,ijk ) + dcrg
4366 gc( n,il,ijk ) = gc( n,il,ijk ) - frz
4367 gc( n,ip,ijk ) = gc( n,ip,ijk ) + frz
4369 sumfrz = sumfrz + frz
4372 dmp = rate*expxctr( n )
4373 frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
4374 frz = min( frz, gc( n,il,ijk ) )
4377 dcrg = frz / ( gc(n,il,ijk)+eps ) * gcrg( n,il,ijk )
4378 gcrg( n,il,ijk ) = gcrg( n,il,ijk ) - dcrg
4379 gcrg( n,ih,ijk ) = gcrg( n,ih,ijk ) + dcrg
4382 gc( n,il,ijk ) = gc( n,il,ijk ) - frz
4383 gc( n,ih,ijk ) = gc( n,ih,ijk ) + frz
4385 sumfrz = sumfrz + frz
4422 sumfrz = sumfrz*dxmic
4424 tdel = sumfrz/dens(ijk)*qlmlt(ijk)/cv(ijk)
4425 temp(ijk) = temp(ijk) + tdel
4433 end subroutine freezing
4436 subroutine melting( &
4449 atmos_hydrometeor_lhf_para, &
4456 integer,
intent(in) :: ijkmax
4457 integer,
intent(in) :: num_warm
4458 integer,
intent(in) :: index_warm(ijkmax)
4459 logical,
intent(in) :: flg_lt
4460 real(
rp),
intent(in) :: dens(ijkmax)
4461 real(
rp),
intent(inout) :: temp(ijkmax)
4462 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4463 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4464 real(
rp),
intent(inout) :: cp(ijkmax)
4465 real(
rp),
intent(inout) :: cv(ijkmax)
4466 real(
dp),
intent(in) :: dtime
4469 real(
rp) :: summlt, sumice, tdel
4470 integer :: ijk, indirect
4471 real(
rp) :: qlmlt(ijkmax)
4476 call atmos_hydrometeor_lhf_para( ijkmax, 1, ijkmax, temp(:), qlmlt(:) )
4479 do indirect = 1, num_warm
4480 ijk = index_warm(indirect)
4486 sumice = sumice + gc( n,m,ijk )
4487 gc( n,m,ijk ) = 0.0_rp
4489 gc( n,il,ijk ) = gc( n,il,ijk ) + sumice
4493 dcrg = dcrg + gcrg( n,m,ijk )
4494 gcrg( n,m,ijk ) = 0.0_rp
4496 gcrg( n,il,ijk ) = gcrg( n,il,ijk ) + dcrg
4498 summlt = summlt + sumice
4500 summlt = summlt*dxmic
4502 tdel = - summlt/dens(ijk)*qlmlt(ijk)/cv(ijk)
4503 temp(ijk) = temp(ijk) + tdel
4512 end subroutine melting
4515 subroutine collmain( &
4530 integer,
intent(in) :: ka
4531 integer,
intent(in) :: ia
4532 integer,
intent(in) :: ja
4534 integer,
intent(in) :: ijkmax
4535 logical,
intent(in) :: flg_lt
4536 real(
rp),
intent(in) :: v0_crg, d0_crg
4537 real(
rp),
intent(in) :: dq(ijkmax)
4538 real(
rp),
intent(in) :: beta(ijkmax)
4539 real(
rp),
intent(in) :: temp(ijkmax)
4540 real(
rp),
intent(inout) :: ghyd(
nbin,nspc,ijkmax)
4541 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4542 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
4543 real(
dp),
intent(in) :: dt
4546 if ( flg_rndm )
then
4547 call r_collcoag( ka, ia, ja, &
4554 call collcoag( ijkmax, &
4568 end subroutine collmain
4571 subroutine collmainf( &
4586 integer,
intent(in) :: ka
4587 integer,
intent(in) :: ia
4588 integer,
intent(in) :: ja
4590 integer,
intent(in) :: ijkmax
4591 logical,
intent(in) :: flg_lt
4592 real(
rp),
intent(in) :: v0_crg, d0_crg
4593 real(
rp),
intent(in) :: dq(ijkmax)
4594 real(
rp),
intent(in) :: beta(ijkmax)
4595 real(
rp),
intent(in) :: temp(ijkmax)
4596 real(
rp),
intent(inout) :: ghyd(
nbin,nspc,ijkmax)
4597 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4598 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
4599 real(
dp),
intent(in) :: dt
4602 if ( flg_rndm )
then
4603 call r_collcoag( ka, ia, ja, &
4610 call collcoag( ijkmax, &
4624 end subroutine collmainf
4629 subroutine collcoag( &
4643 integer,
intent(in) :: ijkmax
4644 logical,
intent(in) :: flg_lt
4645 real(
rp),
intent(in) :: d0_crg, v0_crg
4646 real(
rp),
intent(in) :: dq(ijkmax)
4647 real(
rp),
intent(in) :: beta(ijkmax)
4648 real(
rp),
intent(in) :: temp(ijkmax)
4649 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4650 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4651 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
4652 real(
dp),
intent(in) :: dtime
4655 real(
rp) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
4656 real(
rp) :: gprime, gprimk, wgt, crn, sum, flux
4657 real(
rp) :: acoef0, acoef1, acoef2
4658 real(
rp),
parameter :: dmpmin = 1.e-01_rp
4659 real(
rp) :: suri, surj
4661 integer :: myu, n, irsl, ilrg, isml
4662 integer :: ibnd( ijkmax )
4663 integer :: iflg( nspc,ijkmax )
4664 integer :: iexst(
nbin,nspc,ijkmax )
4666 integer :: ijk, nn, mm, pp, qq
4669 real(
rp) :: drhoi, drhoj, drhok, alpha, dgenei, dgenej
4678 crg_sep( :,: ) = 0.0_rp
4686 csum = csum + gc( n,myu,ijk )*dxmic
4688 if ( csum > cldmin ) iflg( myu,ijk ) = 1
4691 if ( temp(ijk) < tcrit )
then
4699 if ( gc( n,myu,ijk ) > cldmin )
then
4700 iexst( n,myu,ijk ) = 1
4714 do nn = 1, iflg( isml,ijk )
4717 do mm = 1, iflg( ilrg,ijk )
4719 irsl = ifrsl( ibnd(ijk),isml,ilrg )
4721 if ( isml /= ilrg )
then
4723 else if ( isml == ilrg )
then
4728 do pp = 1, iexst( i,isml,ijk )
4733 if ( iexst(j,ilrg,ijk) == 1 )
then
4740 dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
4741 dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
4743 if ( dmpi <= dmpmin )
then
4744 frci = gc( i,isml,ijk )*dmpi
4746 frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
4749 if ( dmpj <= dmpmin )
then
4750 frcj = gc( j,ilrg,ijk )*dmpj
4752 frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
4757 if ( gprime > 0.0_rp .AND. k <
nbin )
then
4759 suri = gc( i,isml,ijk )
4760 surj = gc( j,ilrg,ijk )
4761 gc( i,isml,ijk ) = gc( i,isml,ijk )-frci
4762 gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj
4763 gc( i,isml,ijk ) = max( gc( i,isml,ijk )-frci, 0.0_rp )
4764 gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk )-frcj, 0.0_rp )
4765 frci = suri - gc( i,isml,ijk )
4766 frcj = surj - gc( j,ilrg,ijk )
4771 drhoi = frci/suri*gcrg( i,isml,ijk )
4772 drhoj = frcj/surj*gcrg( j,ilrg,ijk )
4773 drhok = drhoi + drhoj
4776 alpha = 5.0_rp * ( 2.0_rp*radc( i )/d0_crg )**2 &
4777 * abs( vt(ilrg,j)-vt(isml,i) )/v0_crg
4778 alpha = min( 10.0_rp, alpha )
4780 dgenei = frci / xi * rcoll( isml,ilrg,i,j ) * ( -dq( ijk ) ) &
4781 * alpha * beta( ijk ) * flg_noninduct( isml,ilrg )
4783 gcrg( i,isml,ijk ) = gcrg( i,isml,ijk ) + ( dgenei-drhoi )
4784 gcrg( j,ilrg,ijk ) = gcrg( j,ilrg,ijk ) + ( dgenej-drhoj )
4785 crg_sep( isml,ijk ) = crg_sep( isml,ijk ) + dgenei
4786 crg_sep( ilrg,ijk ) = crg_sep( ilrg,ijk ) + dgenej
4789 gprimk = gc( k,irsl,ijk ) + gprime
4790 wgt = gprime / gprimk
4791 crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
4793 acoef0 = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
4794 acoef1 = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.50_rp
4795 acoef2 = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
4803 sum = acoef0 / 2.0_rp * ( 1.0_rp - ( 1.0_rp-2.0_rp*crn ) ) &
4804 + acoef1 / 8.0_rp * ( 1.0_rp - ( 1.0_rp-2.0_rp*crn )**2 ) &
4805 + acoef2 / 24.0_rp * ( 1.0_rp - ( 1.0_rp-2.0_rp*crn )**3 )
4808 flux = min( max( flux,0.0_rp ),gprime )
4810 gc( k,irsl,ijk ) = gprimk - flux
4811 gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
4815 if( gprime /= 0.0_rp )
then
4816 gcrg( k,irsl,ijk ) = gcrg( k,irsl,ijk ) + drhok * ( gprime-flux )/gprime
4817 gcrg( k+1,irsl,ijk ) = gcrg( k+1,irsl,ijk ) + drhok * flux/gprime
4840 end subroutine collcoag
4843 subroutine getrule( ifrsl, indx )
4845 integer,
intent(out) :: indx(
nbin,
nbin )
4846 integer,
intent(out) :: ifrsl( 2,nspc_mk,nspc_mk )
4853 xnew = log( expxctr( i )+expxctr( j ) )
4854 k = int( ( xnew-xctr( 1 ) )/dxmic ) + 1
4855 k = max( max( k,j ),i )
4862 ifrsl( 1:2,il,il ) = il
4865 ifrsl( 1:2,il,ic ) = ic
4866 ifrsl( 1,ic,il ) = ig
4867 ifrsl( 2,ic,il ) = ih
4870 ifrsl( 1:2,il,ip ) = ip
4871 ifrsl( 1,ip,il ) = ig
4872 ifrsl( 2,ip,il ) = ih
4875 ifrsl( 1:2,il,id ) = id
4876 ifrsl( 1,id,il ) = ig
4877 ifrsl( 2,id,il ) = ih
4880 ifrsl( 1:2,il,iss ) = iss
4881 ifrsl( 1,iss,il ) = ig
4882 ifrsl( 2,iss,il ) = ih
4885 ifrsl( 1:2,il,ig ) = ig
4886 ifrsl( 1,ig,il ) = ig
4887 ifrsl( 2,ig,il ) = ih
4890 ifrsl( 1:2,il,ih ) = ih
4891 ifrsl( 1,ih,il ) = ig
4892 ifrsl( 2,ih,il ) = ih
4896 ifrsl( 1:2,ic,ic ) = iss
4899 ifrsl( 1:2,ic,ip ) = iss
4900 ifrsl( 1:2,ip,ic ) = iss
4903 ifrsl( 1:2,ic,id ) = iss
4904 ifrsl( 1:2,id,ic ) = iss
4907 ifrsl( 1:2,ic,iss ) = iss
4908 ifrsl( 1:2,iss,ic ) = iss
4911 ifrsl( 1:2,ic,ig ) = ig
4912 ifrsl( 1:2,ig,ic ) = ic
4915 ifrsl( 1:2,ih,ic ) = ic
4916 ifrsl( 1,ic,ih ) = ig
4917 ifrsl( 2,ic,ih ) = ih
4921 ifrsl( 1:2,ip,ip ) = iss
4924 ifrsl( 1:2,ip,id ) = iss
4925 ifrsl( 1:2,id,ip ) = iss
4928 ifrsl( 1:2,ip,iss ) = iss
4929 ifrsl( 1:2,iss,ip ) = iss
4932 ifrsl( 1:2,ip,ig ) = ig
4933 ifrsl( 1:2,ig,ip ) = ip
4936 ifrsl( 1:2,ih,ip ) = ip
4937 ifrsl( 1,ip,ih ) = ig
4938 ifrsl( 2,ip,ih ) = ih
4942 ifrsl( 1:2,id,id ) = iss
4945 ifrsl( 1:2,id,iss ) = iss
4946 ifrsl( 1:2,iss,id ) = iss
4949 ifrsl( 1:2,id,ig ) = ig
4950 ifrsl( 1:2,ig,id ) = id
4953 ifrsl( 1:2,ih,id ) = id
4954 ifrsl( 1,id,ih ) = ig
4955 ifrsl( 2,id,ih ) = ih
4959 ifrsl( 1:2,iss,iss ) = iss
4962 ifrsl( 1:2,iss,ig ) = ig
4963 ifrsl( 1:2,ig,iss ) = iss
4966 ifrsl( 1:2,ih,iss ) = iss
4967 ifrsl( 1,iss,ih ) = ig
4968 ifrsl( 2,iss,ih ) = ih
4972 ifrsl( 1:2,ig,ig ) = ig
4975 ifrsl( 1,ig,ih ) = ig
4976 ifrsl( 1,ih,ig ) = ig
4977 ifrsl( 2,ig,ih ) = ih
4978 ifrsl( 2,ih,ig ) = ih
4981 ifrsl( 1:2,ih,ih ) = ih
4984 end subroutine getrule
4993 integer ,
intent(in) :: ijkmax
4994 real(
rp),
intent(in) :: f0(ijkmax)
4995 real(
rp),
intent(inout) :: ga(
nccn,ijkmax)
5007 ga(n,ijk) = ga(n,ijk) + f0(ijk) * marate(n) * expxactr(n) / dxaer
5014 end subroutine faero
5019 subroutine random_setup( mset )
5025 integer,
intent(in) :: mset
5029 real(
rp) :: nbinr, tmp1
5030 real(
rp) :: rans( mbin ), ranl( mbin )
5032 real(
rp),
allocatable :: ranstmp( : )
5033 real(
rp),
allocatable :: ranltmp( : )
5036 integer,
allocatable :: orderb( : )
5039 real(
rp),
allocatable :: randnum(:,:,:)
5042 allocate( blrg( mset, mbin ) )
5043 allocate( bsml( mset, mbin ) )
5044 allocate( ranstmp( pq ) )
5045 allocate( ranltmp( pq ) )
5046 allocate( orderb( pq ) )
5047 allocate( randnum(1,1,pq) )
5049 a = real(
nbin )*real(
nbin-1 )*0.50_rp
5050 if ( a < mbin )
then
5051 log_error(
"ATMOS_PHY_MP_SUZUKI10_random_setup",*)
"mbin should be smaller than {nbin}_C_{2}"
5055 wgtbin = a/real( mbin )
5056 nbinr = real(
nbin )
5063 ranstmp( (p-1)*
nbin-(p*(p-1))/2+1 : p*
nbin-(p*(p+1))/2 ) = p
5065 ranltmp( (p-1)*
nbin-(p*(p-1))/2+q ) = p+q
5070 call random_uniform( randnum )
5072 abq1 = randnum( 1,1,p )
5073 k = int( abq1*( pq-p-1 ) ) + p
5075 orderb( p ) = orderb( k )
5081 rans( p ) = ranstmp( orderb( p ) )
5082 ranl( p ) = ranltmp( orderb( p ) )
5084 rans( p ) = ranstmp( orderb( p-pq ) )
5085 ranl( p ) = ranltmp( orderb( p-pq ) )
5087 if ( rans( p ) >= ranl( p ) )
then
5089 rans( p ) = ranl( p )
5093 blrg( n,1:mbin ) = int( ranl( 1:mbin ) )
5094 bsml( n,1:mbin ) = int( rans( 1:mbin ) )
5097 deallocate( ranstmp )
5098 deallocate( ranltmp )
5099 deallocate( orderb )
5100 deallocate( randnum )
5102 end subroutine random_setup
5108 subroutine r_collcoag( &
5119 integer,
intent(in) :: ka
5120 integer,
intent(in) :: ia
5121 integer,
intent(in) :: ja
5123 integer,
intent(in) :: ijkmax
5124 real(
rp),
intent(in) :: swgt
5125 real(
rp),
intent(in) :: temp(ijkmax)
5126 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
5127 real(
dp),
intent(in) :: dtime
5129 integer :: i, j, k, l
5130 real(
rp) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
5131 real(
rp) :: gprime, gprimk, wgt, crn, sum, flux
5132 integer,
parameter :: ldeg = 2
5133 real(
rp),
parameter :: dmpmin = 1.e-01_rp, cmin = 1.e-10_rp
5134 real(
rp) :: acoef( 0:ldeg )
5137 integer :: nums( mbin ), numl( mbin )
5138 real(
rp),
parameter :: gt = 1.0_rp
5140 real(
rp) :: nbinr, mbinr
5142 real(
rp) :: tmpi, tmpj
5144 real(
rp) :: rndm(nspc,nspc)
5146 integer :: ibnd( ijkmax )
5147 integer :: iflg( nspc,ijkmax )
5148 integer :: iexst(
nbin,nspc,ijkmax )
5150 integer :: ijk, nn, mm, pp, qq, myu, n, isml, ilrg, irsl
5165 csum = csum + gc( n,isml,ijk )*dxmic
5167 if ( csum > cldmin ) iflg( isml,ijk ) = 1
5170 if ( temp(ijk) < tcrit )
then
5178 if ( gc( n,myu,ijk ) > cldmin )
then
5179 iexst( n,myu,ijk ) = 1
5189 call random_uniform( rndm(:,:) )
5192 do nn = 1, iflg( isml,ijk )
5195 do mm = 1, iflg( ilrg,ijk )
5197 irsl = ifrsl( ibnd(ijk),isml,ilrg )
5199 det = int( rndm(ilrg,isml)*ia*ja*ka + 1 )
5200 nbinr = real(
nbin )
5201 mbinr = real( mbin )
5202 nums( 1:mbin ) = bsml( det,1:mbin )
5203 numl( 1:mbin ) = blrg( det,1:mbin )
5209 do pp = 1, iexst( i,isml,ijk )
5210 do qq = 1, iexst( j,ilrg,ijk )
5217 dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
5218 dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
5220 if ( dmpi <= dmpmin )
then
5221 frci = gc( i,isml,ijk )*dmpi
5223 frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
5226 if ( dmpj <= dmpmin )
then
5227 frcj = gc( j,ilrg,ijk )*dmpj
5229 frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
5231 tmpi = gc( i,isml,ijk )
5232 tmpj = gc( j,ilrg,ijk )
5234 gc( i,isml,ijk ) = gc( i,isml,ijk )-frci*swgt
5235 gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj*swgt
5238 gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk ), 0.0_rp )
5240 gc( i,isml,ijk ) = max( gc( i,isml,ijk ), 0.0_rp )
5242 frci = tmpi - gc( i,isml,ijk )
5243 frcj = tmpj - gc( j,ilrg,ijk )
5268 if ( gprime > 0.0_rp .AND. k <
nbin )
then
5270 gprimk = gc( k,irsl,ijk ) + gprime
5271 wgt = gprime / gprimk
5272 crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
5274 acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
5275 acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.5_rp
5276 acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
5280 sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
5281 *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
5285 flux = min( max( flux,0.0_rp ),gprime )
5287 gc( k,irsl,ijk ) = gprimk - flux
5288 gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
5307 end subroutine r_collcoag
5320 allocate( radc_mk(
nbin ) )
5321 allocate( xctr_mk(
nbin ) )
5322 allocate( xbnd_mk(
nbin+1 ) )
5323 allocate( cctr_mk( nspc_mk,
nbin ) )
5324 allocate( cbnd_mk( nspc_mk,
nbin+1 ) )
5325 allocate( ck_mk( nspc_mk,nspc_mk,
nbin,
nbin ) )
5326 allocate( vt_mk( nspc_mk,
nbin ) )
5327 allocate( br_mk( nspc_mk,
nbin ) )
5350 deallocate( radc_mk )
5351 deallocate( xctr_mk )
5352 deallocate( xbnd_mk )
5353 deallocate( cctr_mk )
5354 deallocate( cbnd_mk )
5359 end subroutine mkpara
5365 integer,
parameter :: il = 1, ic = 2, ip = 3, id = 4
5366 integer,
parameter :: is = 5, ig = 6, ih = 7
5367 integer,
parameter :: icemax = 3
5369 integer :: k, kk, i, j
5370 real(
dp) :: xl( ndat ), rlec( ndat ), vrl( ndat )
5371 real(
dp) :: blkradl( ndat ), blkdnsl( ndat )
5373 real(
dp) :: xi( ndat,icemax ), riec( ndat,icemax ), vri( ndat,icemax )
5374 real(
dp) :: blkradi( ndat,icemax ), blkdnsi( ndat,icemax )
5376 real(
dp) :: xs( ndat ), rsec( ndat ), vrs( ndat )
5377 real(
dp) :: blkrads( ndat ), blkdnss( ndat )
5379 real(
dp) :: xg( ndat ), rgec( ndat ), vrg( ndat )
5380 real(
dp) :: blkradg( ndat ), blkdnsg( ndat )
5382 real(
dp) :: xh( ndat ), rhec( ndat ), vrh( ndat )
5383 real(
dp) :: blkradh( ndat ), blkdnsh( ndat )
5386 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5387 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5388 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5389 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5390 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5391 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5392 data xi(1:ndat,1) / &
5393 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5394 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5395 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5396 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5397 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5398 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5399 data xi(1:ndat,2) / &
5400 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5401 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5402 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5403 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5404 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5405 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00 /
5406 data xi(1:ndat,3) / &
5407 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5408 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5409 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5410 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5411 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5412 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5414 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5415 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5416 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5417 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5418 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5419 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5421 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5422 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5423 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5424 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5425 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5426 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5428 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5429 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5430 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5431 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5432 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5433 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5434 data rlec(1:ndat) / &
5435 0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
5436 0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
5437 0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
5438 0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
5439 0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
5440 0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
5441 data riec(1:ndat,1) / &
5442 0.31936e-03_dp,0.40397e-03_dp,0.51099e-03_dp,0.64638e-03_dp,0.81764e-03_dp,0.10343e-02_dp, &
5443 0.13084e-02_dp,0.16551e-02_dp,0.20937e-02_dp,0.26486e-02_dp,0.33506e-02_dp,0.42387e-02_dp, &
5444 0.64360e-02_dp,0.81426e-02_dp,0.10302e-01_dp,0.13035e-01_dp,0.16494e-01_dp,0.20872e-01_dp, &
5445 0.26412e-01_dp,0.33426e-01_dp,0.42304e-01_dp,0.53543e-01_dp,0.67770e-01_dp,0.85783e-01_dp, &
5446 0.10859e+00_dp,0.13746e+00_dp,0.17403e+00_dp,0.22032e+00_dp,0.27895e+00_dp,0.35319e+00_dp, &
5447 0.44722e+00_dp,0.56630e+00_dp,0.71712e+00_dp /
5448 data riec(1:ndat,2) / &
5449 0.13188e-03_dp,0.16615e-03_dp,0.20953e-03_dp,0.27728e-03_dp,0.36694e-03_dp,0.48559e-03_dp, &
5450 0.64261e-03_dp,0.85040e-03_dp,0.11254e-02_dp,0.14893e-02_dp,0.19709e-02_dp,0.26082e-02_dp, &
5451 0.34515e-02_dp,0.45676e-02_dp,0.60446e-02_dp,0.79991e-02_dp,0.10586e-01_dp,0.14009e-01_dp, &
5452 0.18539e-01_dp,0.24533e-01_dp,0.32466e-01_dp,0.42964e-01_dp,0.56857e-01_dp,0.75242e-01_dp, &
5453 0.99573e-01_dp,0.13177e+00_dp,0.17438e+00_dp,0.23077e+00_dp,0.30539e+00_dp,0.40414e+00_dp, &
5454 0.53482e+00_dp,0.70775e+00_dp,0.93661e+00_dp /
5455 data riec(1:ndat,3) / 0.14998e-03_dp,0.18896e-03_dp,0.23808e-03_dp, &
5456 0.29996e-03_dp,0.37793e-03_dp,0.47616e-03_dp,0.61048e-03_dp,0.81343e-03_dp,0.10839e-02_dp, &
5457 0.14442e-02_dp,0.19243e-02_dp,0.25640e-02_dp,0.34164e-02_dp,0.45522e-02_dp,0.60656e-02_dp, &
5458 0.80820e-02_dp,0.10769e-01_dp,0.14349e-01_dp,0.19119e-01_dp,0.25475e-01_dp,0.44576e-01_dp, &
5459 0.62633e-01_dp,0.88006e-01_dp,0.12366e+00_dp,0.17375e+00_dp,0.24414e+00_dp,0.34304e+00_dp, &
5460 0.48201e+00_dp,0.67728e+00_dp,0.95164e+00_dp,0.13372e+01_dp,0.18788e+01_dp,0.26400e+01_dp /
5461 data rsec(1:ndat) / &
5462 0.92832e-03_dp,0.11696e-02_dp,0.14736e-02_dp,0.18566e-02_dp,0.23392e-02_dp,0.29472e-02_dp, &
5463 0.37133e-02_dp,0.46784e-02_dp,0.58944e-02_dp,0.74265e-02_dp,0.93569e-02_dp,0.11789e-01_dp, &
5464 0.14853e-01_dp,0.18714e-01_dp,0.23578e-01_dp,0.29706e-01_dp,0.37427e-01_dp,0.47156e-01_dp, &
5465 0.59412e-01_dp,0.74855e-01_dp,0.94311e-01_dp,0.11882e+00_dp,0.14971e+00_dp,0.18862e+00_dp, &
5466 0.23765e+00_dp,0.29942e+00_dp,0.37724e+00_dp,0.47530e+00_dp,0.59884e+00_dp,0.75449e+00_dp, &
5467 0.95060e+00_dp,0.11977e+01_dp,0.15090e+01_dp /
5468 data rgec(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
5469 0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
5470 0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
5471 0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
5472 0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
5473 0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
5474 data rhec(1:ndat) / &
5475 0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
5476 0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
5477 0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
5478 0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
5479 0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
5480 0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
5481 data vrl(1:ndat) / &
5482 0.50000e-01_dp,0.78000e-01_dp,0.12000e+00_dp,0.19000e+00_dp,0.31000e+00_dp,0.49000e+00_dp, &
5483 0.77000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp,0.48000e+01_dp,0.74000e+01_dp, &
5484 0.11000e+02_dp,0.17000e+02_dp,0.26000e+02_dp,0.37000e+02_dp,0.52000e+02_dp,0.71000e+02_dp, &
5485 0.94000e+02_dp,0.12000e+03_dp,0.16000e+03_dp,0.21000e+03_dp,0.26000e+03_dp,0.33000e+03_dp, &
5486 0.41000e+03_dp,0.48000e+03_dp,0.57000e+03_dp,0.66000e+03_dp,0.75000e+03_dp,0.82000e+03_dp, &
5487 0.88000e+03_dp,0.90000e+03_dp,0.90000e+03_dp /
5488 data vri(1:ndat,1) / 0.30000e-01_dp,0.40000e-01_dp,0.60000e-01_dp, &
5489 0.80000e-01_dp,0.11000e+00_dp,0.15000e+00_dp,0.17000e+00_dp,0.18000e+00_dp,0.20000e+00_dp, &
5490 0.25000e+00_dp,0.40000e+00_dp,0.60000e+01_dp,0.10000e+02_dp,0.15000e+02_dp,0.20000e+02_dp, &
5491 0.25000e+02_dp,0.31000e+02_dp,0.37000e+02_dp,0.41000e+02_dp,0.46000e+02_dp,0.51000e+02_dp, &
5492 0.55000e+02_dp,0.59000e+02_dp,0.62000e+02_dp,0.64000e+02_dp,0.67000e+02_dp,0.68000e+02_dp, &
5493 0.69000e+02_dp,0.70000e+02_dp,0.71000e+02_dp,0.71500e+02_dp,0.71750e+02_dp,0.72000e+02_dp /
5494 data vri(1:ndat,2) / &
5495 0.30000e-01_dp,0.40000e-01_dp,0.50000e-01_dp,0.70000e-01_dp,0.90000e-01_dp,0.12000e+00_dp, &
5496 0.50000e+00_dp,0.80000e+00_dp,0.16000e+01_dp,0.18000e+01_dp,0.20000e+01_dp,0.30000e+01_dp, &
5497 0.40000e+01_dp,0.50000e+01_dp,0.80000e+01_dp,0.13000e+02_dp,0.19000e+02_dp,0.26000e+02_dp, &
5498 0.32000e+02_dp,0.38000e+02_dp,0.47000e+02_dp,0.55000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
5499 0.77000e+02_dp,0.79000e+02_dp,0.80000e+02_dp,0.81000e+02_dp,0.81000e+02_dp,0.82000e+02_dp, &
5500 0.82000e+02_dp,0.82000e+02_dp,0.82000e+02_dp /
5501 data vri(1:ndat,3) / 0.35000e-01_dp,0.45000e-01_dp,0.55000e-01_dp, &
5502 0.75000e-01_dp,0.95000e-01_dp,0.13000e+00_dp,0.60000e+00_dp,0.90000e+00_dp,0.17000e+01_dp, &
5503 0.20000e+01_dp,0.25000e+01_dp,0.38000e+01_dp,0.50000e+01_dp,0.70000e+01_dp,0.90000e+01_dp, &
5504 0.11000e+02_dp,0.14000e+02_dp,0.17000e+02_dp,0.21000e+02_dp,0.25000e+02_dp,0.32000e+02_dp, &
5505 0.38000e+02_dp,0.44000e+02_dp,0.49000e+02_dp,0.53000e+02_dp,0.55000e+02_dp,0.58000e+02_dp, &
5506 0.59000e+02_dp,0.61000e+02_dp,0.62000e+02_dp,0.63000e+02_dp,0.64000e+02_dp,0.65000e+02_dp /
5507 data vrs(1:ndat) / &
5508 0.20000e-01_dp,0.31000e-01_dp,0.49000e-01_dp,0.77000e-01_dp,0.12000e+00_dp,0.19000e+00_dp, &
5509 0.30000e+00_dp,0.48000e+00_dp,0.76000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp, &
5510 0.48000e+01_dp,0.75000e+01_dp,0.11000e+02_dp,0.16000e+02_dp,0.21000e+02_dp,0.26000e+02_dp, &
5511 0.34000e+02_dp,0.41000e+02_dp,0.49000e+02_dp,0.57000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
5512 0.81000e+02_dp,0.87000e+02_dp,0.93000e+02_dp,0.99000e+02_dp,0.10750e+03_dp,0.11500e+03_dp, &
5513 0.12500e+03_dp,0.13500e+03_dp,0.14500e+03_dp /
5514 data vrg(1:ndat) / 0.39000e-01_dp,0.62000e-01_dp,0.97000e-01_dp, &
5515 0.15000e+00_dp,0.24000e+00_dp,0.38000e+00_dp,0.61000e+00_dp,0.96000e+00_dp,0.15000e+01_dp, &
5516 0.24000e+01_dp,0.38000e+01_dp,0.61000e+01_dp,0.96000e+01_dp,0.15000e+02_dp,0.23000e+02_dp, &
5517 0.31000e+02_dp,0.39000e+02_dp,0.49000e+02_dp,0.59000e+02_dp,0.68000e+02_dp,0.79000e+02_dp, &
5518 0.88000e+02_dp,0.10000e+03_dp,0.11000e+03_dp,0.13000e+03_dp,0.15000e+03_dp,0.17000e+03_dp, &
5519 0.20000e+03_dp,0.23000e+03_dp,0.26000e+03_dp,0.30000e+03_dp,0.35000e+03_dp,0.40000e+03_dp /
5520 data vrh(1:ndat) / &
5521 0.53000e-01_dp,0.84000e-01_dp,0.13000e+00_dp,0.21000e+00_dp,0.33000e+00_dp,0.52000e+00_dp, &
5522 0.82000e+00_dp,0.13000e+01_dp,0.21000e+01_dp,0.33000e+01_dp,0.52000e+01_dp,0.82000e+01_dp, &
5523 0.13000e+02_dp,0.20000e+02_dp,0.28000e+02_dp,0.36000e+02_dp,0.46000e+02_dp,0.56000e+02_dp, &
5524 0.67000e+02_dp,0.80000e+02_dp,0.97000e+02_dp,0.12000e+03_dp,0.14000e+03_dp,0.17000e+03_dp, &
5525 0.20000e+03_dp,0.24000e+03_dp,0.29000e+03_dp,0.35000e+03_dp,0.42000e+03_dp,0.51000e+03_dp, &
5526 0.61000e+03_dp,0.74000e+03_dp,0.89000e+03_dp /
5527 data blkradl(1:ndat) / &
5528 0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
5529 0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
5530 0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
5531 0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
5532 0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
5533 0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
5534 data blkradi(1:ndat,1) / 0.57452e-03_dp,0.72384e-03_dp,0.91199e-03_dp, &
5535 0.11490e-02_dp,0.14477e-02_dp,0.18240e-02_dp,0.22981e-02_dp,0.28954e-02_dp,0.36479e-02_dp, &
5536 0.45961e-02_dp,0.57908e-02_dp,0.72959e-02_dp,0.11572e-01_dp,0.14770e-01_dp,0.18851e-01_dp, &
5537 0.24060e-01_dp,0.30709e-01_dp,0.39194e-01_dp,0.50025e-01_dp,0.63848e-01_dp,0.81491e-01_dp, &
5538 0.10401e+00_dp,0.13275e+00_dp,0.16943e+00_dp,0.21625e+00_dp,0.27601e+00_dp,0.35228e+00_dp, &
5539 0.44962e+00_dp,0.57387e+00_dp,0.73244e+00_dp,0.93484e+00_dp,0.11932e+01_dp,0.15229e+01_dp /
5540 data blkradi(1:ndat,2) / &
5541 0.20715e-03_dp,0.26099e-03_dp,0.32912e-03_dp,0.43555e-03_dp,0.57638e-03_dp,0.76276e-03_dp, &
5542 0.10094e-02_dp,0.13358e-02_dp,0.17677e-02_dp,0.23394e-02_dp,0.30958e-02_dp,0.40969e-02_dp, &
5543 0.54216e-02_dp,0.71748e-02_dp,0.94948e-02_dp,0.12565e-01_dp,0.16628e-01_dp,0.22005e-01_dp, &
5544 0.29120e-01_dp,0.38537e-01_dp,0.50998e-01_dp,0.67488e-01_dp,0.89311e-01_dp,0.11819e+00_dp, &
5545 0.15641e+00_dp,0.20698e+00_dp,0.27391e+00_dp,0.36249e+00_dp,0.47970e+00_dp,0.63481e+00_dp, &
5546 0.84009e+00_dp,0.11117e+01_dp,0.14712e+01_dp /
5547 data blkradi(1:ndat,3) / 0.23559e-03_dp,0.29682e-03_dp,0.37397e-03_dp, &
5548 0.47118e-03_dp,0.59365e-03_dp,0.74795e-03_dp,0.95894e-03_dp,0.12777e-02_dp,0.17025e-02_dp, &
5549 0.22685e-02_dp,0.30227e-02_dp,0.40275e-02_dp,0.53665e-02_dp,0.71506e-02_dp,0.95278e-02_dp, &
5550 0.12695e-01_dp,0.16916e-01_dp,0.22539e-01_dp,0.30032e-01_dp,0.40017e-01_dp,0.70019e-01_dp, &
5551 0.98384e-01_dp,0.13824e+00_dp,0.19424e+00_dp,0.27293e+00_dp,0.38350e+00_dp,0.53885e+00_dp, &
5552 0.75714e+00_dp,0.10639e+01_dp,0.14948e+01_dp,0.21004e+01_dp,0.29513e+01_dp,0.41469e+01_dp /
5553 data blkrads(1:ndat) / &
5554 0.20715e-03_dp,0.26148e-03_dp,0.33067e-03_dp,0.41710e-03_dp,0.52691e-03_dp,0.66640e-03_dp, &
5555 0.84289e-03_dp,0.10674e-02_dp,0.13513e-02_dp,0.17129e-02_dp,0.21670e-02_dp,0.27521e-02_dp, &
5556 0.34989e-02_dp,0.44777e-02_dp,0.57347e-02_dp,0.75389e-02_dp,0.99020e-02_dp,0.13161e-01_dp, &
5557 0.17372e-01_dp,0.23337e-01_dp,0.31058e-01_dp,0.41194e-01_dp,0.55153e-01_dp,0.74854e-01_dp, &
5558 0.99806e-01_dp,0.13463e+00_dp,0.18136e+00_dp,0.24282e+00_dp,0.32955e+00_dp,0.44123e+00_dp, &
5559 0.59884e+00_dp,0.77090e+00_dp,0.99387e+00_dp /
5560 data blkradg(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
5561 0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
5562 0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
5563 0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
5564 0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
5565 0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
5566 data blkradh(1:ndat) / &
5567 0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
5568 0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
5569 0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
5570 0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
5571 0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
5572 0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
5573 data blkdnsl(1:ndat) / &
5574 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5575 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5576 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5577 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5578 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5579 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp /
5580 data blkdnsi(1:ndat,1) / 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5581 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5582 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.87368e+00_dp,0.87072e+00_dp,0.86777e+00_dp, &
5583 0.86483e+00_dp,0.86189e+00_dp,0.85897e+00_dp,0.85606e+00_dp,0.85316e+00_dp,0.85026e+00_dp, &
5584 0.84738e+00_dp,0.84451e+00_dp,0.84164e+00_dp,0.83879e+00_dp,0.83595e+00_dp,0.83311e+00_dp, &
5585 0.83029e+00_dp,0.82747e+00_dp,0.82467e+00_dp,0.82187e+00_dp,0.81908e+00_dp,0.81631e+00_dp /
5586 data blkdnsi(1:ndat,2) / &
5587 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5588 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5589 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5590 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5591 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5592 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
5593 data blkdnsi(1:ndat,3) / 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
5594 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
5595 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
5596 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.51790e+00_dp, &
5597 0.45557e+00_dp,0.40075e+00_dp,0.35252e+00_dp,0.31010e+00_dp,0.27278e+00_dp,0.23995e+00_dp, &
5598 0.21108e+00_dp,0.18567e+00_dp,0.16333e+00_dp,0.14367e+00_dp,0.12638e+00_dp,0.11118e+00_dp /
5599 data blkdnss(1:ndat) / &
5600 0.90000e+00_dp,0.89500e+00_dp,0.88500e+00_dp,0.88200e+00_dp,0.87500e+00_dp,0.86500e+00_dp, &
5601 0.85500e+00_dp,0.84200e+00_dp,0.83000e+00_dp,0.81500e+00_dp,0.80500e+00_dp,0.78600e+00_dp, &
5602 0.76500e+00_dp,0.73000e+00_dp,0.69500e+00_dp,0.61183e+00_dp,0.54000e+00_dp,0.46000e+00_dp, &
5603 0.40000e+00_dp,0.33000e+00_dp,0.28000e+00_dp,0.24000e+00_dp,0.20000e+00_dp,0.16000e+00_dp, &
5604 0.13500e+00_dp,0.11000e+00_dp,0.90000e-01_dp,0.75000e-01_dp,0.60000e-01_dp,0.50000e-01_dp, &
5605 0.40000e-01_dp,0.37500e-01_dp,0.35000e-01_dp /
5606 data blkdnsg(1:ndat) / &
5607 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5608 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5609 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5610 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5611 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5612 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp /
5613 data blkdnsh(1:ndat) / &
5614 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5615 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5616 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5617 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5618 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5619 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
5623 xmss( il,k ) = log( dble( xl( k ) )*1.e-03_dp )
5624 xmss( ic,k ) = log( dble( xi( k,1 ) )*1.e-03_dp )
5625 xmss( ip,k ) = log( dble( xi( k,2 ) )*1.e-03_dp )
5626 xmss( id,k ) = log( dble( xi( k,3 ) )*1.e-03_dp )
5627 xmss( is,k ) = log( dble( xs( k ) )*1.e-03_dp )
5628 xmss( ig,k ) = log( dble( xg( k ) )*1.e-03_dp )
5629 xmss( ih,k ) = log( dble( xh( k ) )*1.e-03_dp )
5634 zcap( il,k ) = dble(rlec( k ))*1.e-02_dp
5635 zcap( ic,k ) = dble(riec( k,1 ))*1.e-02_dp
5636 zcap( ip,k ) = dble(riec( k,2 ))*1.e-02_dp
5637 zcap( id,k ) = dble(riec( k,3 ))*1.e-02_dp
5638 zcap( is,k ) = dble(rsec( k ))*1.e-02_dp
5639 zcap( ig,k ) = dble(rgec( k ))*1.e-02_dp
5640 zcap( ih,k ) = dble(rhec( k ))*1.e-02_dp
5645 vtrm( il,k ) = dble(vrl( k ))*1.e-02_dp
5646 vtrm( ic,k ) = dble(vri( k,1 ))*1.e-02_dp
5647 vtrm( ip,k ) = dble(vri( k,2 ))*1.e-02_dp
5648 vtrm( id,k ) = dble(vri( k,3 ))*1.e-02_dp
5649 vtrm( is,k ) = dble(vrs( k ))*1.e-02_dp
5650 vtrm( ig,k ) = dble(vrg( k ))*1.e-02_dp
5651 vtrm( ih,k ) = dble(vrh( k ))*1.e-02_dp
5656 blkr( il,k ) = dble(blkradl( k ))*10.0_dp
5657 blkr( ic,k ) = dble(blkradi( k,1 ))*10.0_dp
5658 blkr( ip,k ) = dble(blkradi( k,2 ))*10.0_dp
5659 blkr( id,k ) = dble(blkradi( k,3 ))*10.0_dp
5660 blkr( is,k ) = dble(blkrads( k ))*10.0_dp
5661 blkr( ig,k ) = dble(blkradg( k ))*10.0_dp
5662 blkr( ih,k ) = dble(blkradh( k ))*10.0_dp
5667 blkd( il,k ) = dble(blkdnsl( k ))*1000._dp
5668 blkd( ic,k ) = dble(blkdnsi( k,1 ))*1000._dp
5669 blkd( ip,k ) = dble(blkdnsi( k,2 ))*1000._dp
5670 blkd( id,k ) = dble(blkdnsi( k,3 ))*1000._dp
5671 blkd( is,k ) = dble(blkdnss( k ))*1000._dp
5672 blkd( ig,k ) = dble(blkdnsg( k ))*1000._dp
5673 blkd( ih,k ) = dble(blkdnsh( k ))*1000._dp
5682 ykrn( k,kk,i,j ) = kernels( k,kk,i,j )*1.e-06_dp
5689 end subroutine rdkdat
5693 real(
dp) :: xsta, xend
5695 real(
dp):: pi_dp = 3.1415920_dp
5696 real(
dp):: rhow_dp = 1.0e+03_dp
5698 xsta = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-06_dp )**3.0_dp )
5699 xend = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-03_dp )**3.0_dp )
5701 dxmic_mk = ( xend-xsta )/
nbin
5704 xbnd_mk( n ) = xsta + dxmic_mk*( n-1 )
5707 xctr_mk( n ) = ( xbnd_mk( n )+xbnd_mk( n+1 ) )*0.50_dp
5708 radc_mk( n ) = ( exp( xctr_mk( n ) )*3.0_dp/4.0_dp/pi_dp/rhow_dp )**( 1.0_dp/3.0_dp )
5711 end subroutine sdfgrid
5720 cctr_mk( myu,n ) = fcpc( myu,xctr_mk( n ) )
5726 cbnd_mk( myu,n ) = fcpc( myu,xbnd_mk( n ) )
5730 end subroutine getcp
5732 function fcpc( myu,x )
5734 integer,
intent(in) :: myu
5735 real(
dp),
intent(in) :: x
5738 real(
dp) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5741 ( ndat, kdeg, myu, xmss(:,:), &
5745 ( ndat, kdeg, qknt, myu, xmss(:,:), &
5749 ( ndat, kdeg, elm, myu, zcap(:,:), &
5752 fcpc = fspline( ndat, kdeg, coef, qknt, myu, xmss(:,:), x )
5758 integer :: myu, nyu, i, j
5760 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Create micpara.dat'
5761 if( kphase == 0 )
then
5762 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Hydro-dynamic kernel'
5763 else if( kphase == 1 )
then
5764 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Long Kernel'
5765 else if( kphase == 2 )
then
5766 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Golovin Kernel'
5774 ck_mk( myu,nyu,i,j ) = fckrn( myu,nyu,xctr_mk( i ),xctr_mk( j ) )
5782 end subroutine getck
5785 function fckrn( myu,nyu,x,y )
5787 integer,
intent(in) :: myu, nyu
5788 real(
dp),
intent(in) :: x, y
5791 real(
dp) :: qknt( ndat+kdeg ), rknt( ndat+kdeg )
5792 real(
dp) :: coef( ndat,ndat )
5794 real(
dp) :: xlrg, xsml, vlrg, vsml, rlrg
5795 real(
dp):: pi_dp = 3.1415920_dp
5796 real(
dp):: rhow_dp = 1.0e+03_dp
5798 if( kphase == 0 )
then
5800 ( ndat, kdeg, myu, xmss(:,:), &
5803 rknt( : ) = qknt( : )
5806 ( ndat, ndat, kdeg, kdeg, &
5808 xmss(:,:), xmss(:,:), &
5814 ( ndat, ndat, kdeg, kdeg, &
5815 coef, qknt, rknt, myu, &
5816 xmss(:,:), xmss(:,:), &
5818 else if( kphase == 1 )
then
5822 vlrg = (exp( xlrg ) / rhow_dp )*1.e+06_dp
5823 vsml = (exp( xsml ) / rhow_dp )*1.e+06_dp
5825 rlrg = ( exp( xlrg )/( 4.0_dp*pi_dp*rhow_dp )*3.0_dp )**(1.0_dp/3.0_dp )*1.e+06_dp
5827 if( rlrg <=50.0_dp )
then
5828 fckrn = 9.44e+03_dp*( vlrg*vlrg + vsml*vsml )
5830 fckrn = 5.78e-03_dp*( vlrg+vsml )
5832 else if( kphase == 2 )
then
5833 fckrn = 1.5_dp*( exp(x) +exp(y) )
5847 vt_mk( myu,n ) = max( fvterm( myu,xctr_mk( n ) ), 0.0_dp )
5851 end subroutine getvt
5853 function fvterm( myu,x )
5855 integer,
intent(in) :: myu
5856 real(
dp),
intent(in) :: x
5859 real(
dp) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5862 ( ndat, kdeg, myu, xmss(:,:), &
5866 ( ndat, kdeg, qknt, myu, xmss(:,:), &
5870 ( ndat, kdeg, elm, myu, vtrm(:,:), &
5873 fvterm = fspline( ndat, kdeg, coef, qknt, myu, xmss(:,:), x )
5884 br_mk( myu,n ) = fbulkrad( myu, xctr_mk( n ) )
5888 end subroutine getbr
5890 function fbulkrad( myu,x )
5892 integer,
intent(in) :: myu
5893 real(
dp),
intent(in) :: x
5894 real(
dp) :: fbulkrad
5896 real(
dp) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5899 ( ndat, kdeg, myu, xmss(:,:), &
5903 ( ndat, kdeg, qknt, myu, xmss(:,:), &
5907 ( ndat, kdeg, elm, myu, blkr(:,:), &
5910 fbulkrad = fspline( ndat, kdeg, coef, qknt, myu, xmss(:,:), x )
5912 end function fbulkrad
5916 integer :: myu, nyu, i, j, n
5918 open ( fid_micpara, file = fname_micpara, form =
'formatted', status=
'new' )
5920 write( fid_micpara,* ) nspc_mk,
nbin
5924 xctr( n ) = xctr_mk( n )
5925 radc( n ) = radc_mk( n )
5926 write( fid_micpara,* ) n, xctr( n ), radc( n )
5929 xbnd( n ) = xbnd_mk( n )
5930 write( fid_micpara,* ) n, xbnd( n )
5932 write( fid_micpara,* ) dxmic_mk
5937 cctr( n,myu ) = cctr_mk( myu,n )
5938 write( fid_micpara,* ) myu, n, cctr( n,myu )
5941 cbnd( n,myu ) = cbnd_mk( myu,n )
5942 write( fid_micpara,* ) myu, n, cbnd( n,myu )
5951 ck( myu,nyu,i,j ) = ck_mk( myu,nyu,i,j )
5952 write( fid_micpara,* ) myu, nyu, i, j, ck( myu,nyu,i,j )
5961 vt( myu,n ) = vt_mk( myu,n )
5962 write( fid_micpara,* ) myu, n, vt( myu,n )
5969 br( myu,n ) = br_mk( myu,n )
5970 write( fid_micpara,* ) myu, n, br( myu,n )
5974 close ( fid_micpara )
5976 end subroutine paraout
5981 subroutine tinvss(n,a,dt,e,nn,iw,inder)
5985 integer,
intent(in) :: n, nn
5986 integer,
intent(inout) :: inder
5987 real(
dp),
intent(inout) :: a(nn,n)
5988 real(
dp),
intent(inout) :: dt
5989 real(
dp),
intent(in) :: e
5990 integer,
intent(inout) :: iw( 2*n )
5991 integer :: i, j, k, kk, ij, nnn
5992 real(
dp) :: work, aa, az, eps
5993 integer :: ipiv, jpiv
6000 elseif( n == 1 )
then
6002 elseif( n > 1 )
then
6009 write(6,690)
"n= ", n,
"nn= ", nn, &
6010 "n should be less than or equal to nn in TINVSS"
6011 690
format(a8,i5,5x,a4,i5,a55)
6021 if( abs(a(i,j)) <= abs(piv) )
goto 110
6029 if( abs(piv) <= eps )
goto 920
6030 if( k == 1 ) eps = abs(piv)*e
6031 if( ipiv == k )
goto 130
6039 if( jpiv == k )
goto 150
6054 if( i == k )
goto 220
6056 if( az == 0.0_dp )
goto 220
6058 a(i,j) = a(i,j)-a(k,j)*az
6068 if( ij == k )
goto 420
6076 if( ij == k )
goto 400
6089 write(*,691)
"n= ", n,
"should be positive in TINVSS"
6090 691
format(a8,i5,5x,a30)
6098 write(*,692)
'given matrix A to TINVSS is ill conditioned, or sigular withrank =', nnn, &
6099 'return with no further calculation'
6100 692
format(a,1x,i4,1x,a)
6106 if( dt == 0.0_dp )
goto 920
6107 a(1,1) = 1.0_dp/a(1,1)
6110 end subroutine tinvss
6113 subroutine getknot &
6114 ( ndat, kdeg, myu, xdat, &
6117 integer,
intent(in) :: ndat
6118 integer,
intent(in) :: kdeg
6119 integer,
intent(in) :: myu
6120 real(
dp),
intent(in) :: xdat(nspc_mk,ndat)
6122 real(
dp),
intent(out) :: qknt( ndat+kdeg )
6128 qknt( i ) = xdat(myu,1)
6132 qknt( i+kdeg ) = ( xdat(myu,i)+xdat(myu,i+kdeg) )*0.50_dp
6136 qknt( ndat+i ) = xdat(myu,ndat)
6141 end subroutine getknot
6143 recursive function fbspl ( ndat, inum, kdeg, qknt, myu, xdat, x ) &
6148 integer,
intent(in) :: ndat
6149 integer,
intent(in) :: inum
6150 integer,
intent(in) :: kdeg
6151 real(
dp),
intent(in) :: qknt(ndat+kdeg)
6152 integer,
intent(in) :: myu
6153 real(
dp),
intent(in) :: xdat(nspc_mk,ndat)
6154 real(
dp),
intent(in) :: x
6157 real(
dp) :: bsp1, bsp2
6159 if ( ( inum == 1 .AND. x == xdat(myu,1) ) .OR. &
6160 ( inum == ndat .AND. x == xdat(myu,ndat) ) )
then
6165 if ( kdeg == 1 )
then
6166 if ( x >= qknt( inum ) .AND. x < qknt( inum+1 ) )
then
6172 if ( qknt( inum+kdeg-1 ) /= qknt( inum ) )
then
6173 bsp1 = ( x-qknt( inum ) ) &
6174 /( qknt( inum+kdeg-1 )-qknt( inum ) ) &
6175 * fbspl( ndat, inum, kdeg-1, qknt, myu, xdat, x )
6179 if ( qknt( inum+kdeg ) /= qknt( inum+1 ) )
then
6180 bsp2 = ( qknt( inum+kdeg )-x ) &
6181 /( qknt( inum+kdeg )-qknt( inum+1 ) ) &
6182 * fbspl( ndat, inum+1, kdeg-1, qknt, myu, xdat, x )
6212 subroutine getmatrx &
6213 ( ndat, kdeg, qknt, myu, xdat, &
6218 integer,
intent(in) :: ndat
6219 integer,
intent(in) :: kdeg
6220 real(
dp),
intent(in) :: qknt(ndat+kdeg)
6221 integer,
intent(in) :: myu
6222 real(
dp),
intent(in) :: xdat(nspc_mk,ndat)
6224 real(
dp),
intent(out) :: elm(ndat,ndat)
6228 integer :: iw( 2*ndat ), i, j, inder
6229 real(
dp),
parameter :: eps = 0.
6233 elm( i,j ) = fbspl( ndat, j, kdeg, qknt, myu, xdat, xdat(myu,i) )
6237 call tinvss( ndat, elm, dt, eps, ndat, iw, inder )
6241 end subroutine getmatrx
6244 subroutine getcoef &
6245 ( ndat, kdeg, elm, myu, ydat, &
6248 integer,
intent(in) :: ndat
6249 integer,
intent(in) :: kdeg
6250 real(
dp),
intent(in) :: elm(ndat,ndat)
6251 integer,
intent(in) :: myu
6252 real(
dp),
intent(in) :: ydat(nspc_mk,ndat)
6254 real(
dp),
intent(out) :: coef(ndat)
6263 sum = sum + elm( i,j )*ydat(myu,j)
6270 end subroutine getcoef
6273 function fspline ( ndat, kdeg, coef, qknt, myu, xdat, x )
6275 integer,
intent(in) :: ndat
6276 integer,
intent(in) :: kdeg
6277 real(
dp),
intent(in) :: coef( ndat )
6278 real(
dp),
intent(in) :: qknt( ndat+kdeg )
6279 integer,
intent(in) :: myu
6280 real(
dp),
intent(in) :: xdat( ndat )
6281 real(
dp),
intent(in) :: x
6291 sum = sum + coef( i )*fbspl( ndat, i, kdeg, qknt, myu, xdat, x )
6298 end function fspline
6301 subroutine getcoef2 &
6302 ( mdat, ndat, kdeg, ldeg, myu, nyu, &
6303 xdat, ydat, qknt, rknt, zdat, &
6308 integer,
intent(in) :: mdat
6309 integer,
intent(in) :: ndat
6310 integer,
intent(in) :: kdeg
6311 integer,
intent(in) :: ldeg
6312 integer,
intent(in) :: myu
6313 integer,
intent(in) :: nyu
6314 real(
dp),
intent(in) :: xdat(nspc_mk,mdat)
6315 real(
dp),
intent(in) :: ydat(nspc_mk,ndat)
6316 real(
dp),
intent(in) :: qknt(mdat+kdeg)
6317 real(
dp),
intent(in) :: rknt(ndat+ldeg)
6318 real(
dp),
intent(in) :: zdat(nspc_mk,nspc_mk,mdat,ndat)
6320 real(
dp),
intent(out) :: coef( mdat,ndat )
6323 real(
dp) :: elmx( mdat,mdat ), elmy( ndat,ndat )
6324 integer :: iw1( 2*mdat ), iw2( 2*ndat )
6325 real(
dp) :: beta( mdat,ndat ), sum, dt
6326 real(
dp),
parameter :: eps = 0.0_dp
6327 integer :: i, j, k, l, inder
6331 elmx( i,j ) = fbspl( mdat, j, kdeg, qknt, myu, xdat, xdat(myu,i) )
6334 call tinvss( mdat, elmx, dt, eps, mdat, iw1, inder )
6340 sum = sum + elmx( i,j )*zdat(myu,nyu,j,l)
6348 elmy( i,j ) = fbspl( ndat, j, ldeg, rknt, myu, ydat, ydat(myu,i) )
6351 call tinvss( ndat, elmy, dt, eps, ndat, iw2, inder )
6357 sum = sum + elmy( i,j )*beta( k,j )
6365 end subroutine getcoef2
6369 ( mdat, ndat, kdeg, ldeg, &
6370 coef, qknt, rknt, myu, &
6374 integer,
intent(in) :: mdat
6375 integer,
intent(in) :: ndat
6376 integer,
intent(in) :: kdeg
6377 integer,
intent(in) :: ldeg
6378 real(
dp),
intent(in) :: coef( mdat,ndat )
6379 real(
dp),
intent(in) :: qknt( mdat+kdeg )
6380 real(
dp),
intent(in) :: rknt( ndat+ldeg )
6381 integer,
intent(in) :: myu
6382 real(
dp),
intent(in) :: xdat(nspc_mk,mdat), ydat(nspc_mk,ndat)
6383 real(
dp),
intent(in) :: x, y
6385 real(
dp) :: fspline2
6388 real(
dp) :: sum, add
6394 add = coef( i,j )*fbspl( mdat, i, kdeg, qknt, myu, xdat, x ) &
6395 *fbspl( ndat, j, ldeg, rknt, myu, ydat, y )
6404 end function fspline2