71 private :: mp_suzuki10
79 private :: ice_nucleat
87 private :: random_setup
116 integer,
parameter :: i_qv = 1
118 character(len=3) :: namspc (8) = (/
'Qcl', &
127 character(len=27) :: lnamspc(8) = (/
'Mixing ratio of cloud bin', &
128 'Mixing ratio of colum bin', &
129 'Mixing ratio of plate bin', &
130 'Mixing ratio of dendrit bin', &
131 'Mixing ratio of snow bin', &
132 'Mixing ratio of graupel bin', &
133 'Mixing ratio of hail bin', &
134 'Mixing ratio of aerosol bin' /)
139 integer :: kphase = 0
140 integer :: iceflg = 1
141 integer :: num_hyd = 0
143 integer,
parameter :: i_mp_qc = 1
144 integer,
parameter :: i_mp_qcl = 2
145 integer,
parameter :: i_mp_qp = 3
146 integer,
parameter :: i_mp_qd = 4
147 integer,
parameter :: i_mp_qs = 5
148 integer,
parameter :: i_mp_qg = 6
149 integer,
parameter :: i_mp_qh = 7
151 logical :: mp_doautoconversion = .true.
152 logical :: mp_couple_aerosol = .false.
155 integer :: num_start_waters
156 integer :: num_end_waters
157 integer :: num_start_ices
158 integer :: num_end_ices
161 integer,
parameter :: il = 1
162 integer,
parameter :: ic = 2
163 integer,
parameter :: ip = 3
164 integer,
parameter :: id = 4
165 integer,
parameter :: iss = 5
166 integer,
parameter :: ig = 6
167 integer,
parameter :: ih = 7
171 real(
rp),
allocatable :: xctr( : )
172 real(
rp),
allocatable :: xbnd( : )
173 real(
rp),
allocatable :: radc( : )
174 real(
rp),
allocatable :: cctr( :,: )
175 real(
rp),
allocatable :: cbnd( :,: )
176 real(
rp),
allocatable :: ck( :,:,:,: )
177 real(
rp),
allocatable :: vt( :,: )
178 real(
rp),
allocatable :: br( :,: )
179 integer,
allocatable :: ifrsl( :,:,: )
181 real(
rp),
allocatable :: xactr( : )
182 real(
rp),
allocatable :: xabnd( : )
183 real(
rp),
allocatable :: rada( : )
185 real(
rp),
allocatable :: expxctr( : )
186 real(
rp),
allocatable :: expxbnd( : )
187 real(
rp),
allocatable :: expxactr( : )
188 real(
rp),
allocatable :: expxabnd( : )
189 real(
rp),
allocatable :: rexpxctr( : )
190 real(
rp),
allocatable :: rexpxbnd( : )
191 real(
rp),
allocatable :: rexpxactr( : )
192 real(
rp),
allocatable :: rexpxabnd( : )
198 real(
rp),
allocatable,
save :: vterm(:)
201 real(
rp),
parameter :: cldmin = 1.0e-10_rp
202 real(
rp),
parameter :: oneovthird = 1.0_rp / 3.0_rp
203 real(
rp),
parameter :: thirdovforth = 3.0_rp / 4.0_rp
204 real(
rp),
parameter :: twoovthird = 2.0_rp / 3.0_rp
206 real(
rp) :: rbnd = 40.e-06_rp
210 real(
rp) :: rhoa = 2.25e+03_rp
211 real(
rp) :: emaer = 58.0_rp
212 real(
rp) :: emwtr = 18.0_rp
213 real(
rp) :: rasta = 1.e-08_rp
214 real(
rp) :: raend = 1.e-06_rp
215 real(
rp) :: r0a = 1.e-07_rp
217 logical :: flg_regeneration
219 logical :: flg_icenucl
220 logical :: flg_sf_aero
223 real(
rp),
allocatable :: marate( : )
224 integer,
allocatable,
save :: ncld( : )
230 character(len=11),
parameter :: fname_micpara=
"micpara.dat"
232 integer(4) :: fid_micpara
235 integer,
allocatable :: blrg( :,: )
236 integer,
allocatable :: bsml( :,: )
238 integer :: mspc, mbin
239 real(
rp) :: rndm(1,1,1)
242 real(
rp) :: n0_icenucl = 1.e+3_rp
245 real(
rp) :: c_ccn = 100.e+6_rp
246 real(
rp) :: kappa = 0.462_rp
249 real(
rp) :: sigma = 7.5e-02_rp
250 real(
rp) :: vhfct = 2.0_rp
252 real(
rp),
parameter :: tcrit = 271.15_rp
253 integer,
allocatable :: kindx( :,: )
256 integer,
parameter :: ndat = 33, icemax = 3
257 integer,
parameter :: kdeg = 4, ldeg = 4, nspc_mk = 7
261 real(
dp),
allocatable :: radc_mk( : ), xctr_mk( : ), xbnd_mk( : )
262 real(
dp),
allocatable :: cctr_mk( :,: ), cbnd_mk( :,: )
263 real(
dp),
allocatable :: ck_mk( :,:,:,: )
264 real(
dp),
allocatable :: vt_mk( :,: )
265 real(
dp),
allocatable :: br_mk( :,: )
267 real(
dp) :: xmss( nspc_mk,ndat ), zcap( nspc_mk,ndat ), vtrm( nspc_mk,ndat )
268 real(
dp) :: blkr( nspc_mk,ndat ), blkd( nspc_mk,ndat ), ykrn( nspc_mk,nspc_mk,ndat,ndat )
270 real(
dp) :: ywll( ndat,ndat ), ywli( ndat,ndat,icemax ), ywls( ndat,ndat )
271 real(
dp) :: ywlg( ndat,ndat ), ywlh( ndat,ndat )
273 real(
dp) :: ywil( ndat,ndat,icemax ), ywii( ndat,ndat,icemax,icemax )
274 real(
dp) :: ywis( ndat,ndat,icemax ), ywig( ndat,ndat,icemax )
275 real(
dp) :: ywih( ndat,ndat,icemax )
277 real(
dp) :: ywsl( ndat,ndat ), ywsi( ndat,ndat,icemax ), ywss( ndat,ndat )
278 real(
dp) :: ywsg( ndat,ndat ), ywsh( ndat,ndat )
280 real(
dp) :: ywgl( ndat,ndat ), ywgi( ndat,ndat,icemax ), ywgs( ndat,ndat )
281 real(
dp) :: ywgg( ndat,ndat ), ywgh( ndat,ndat )
283 real(
dp) :: ywhl( ndat,ndat ), ywhi( ndat,ndat,icemax ), ywhs( ndat,ndat )
284 real(
dp) :: ywhg( ndat,ndat ), ywhh( ndat,ndat )
287 real(
rp) :: sigma_sdf(5), r0_sdf(5), n0_sdf(5), rho_sdf(5)
290 real(
rp),
allocatable :: flg_noninduct(:,:)
291 real(
rp),
allocatable :: ecoll(:,:,:,:), rcoll(:,:,:,:)
292 real(
rp),
allocatable,
save :: num_end(:,:,:,:)
293 real(
rp),
save :: ecoal_gsi = 0.0_rp
303 namelist / param_atmos_phy_mp_suzuki10_bin / &
309 integer :: m, n, ierr
313 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Setup'
314 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Tracers setup for Suzuki (2010) Spectral BIN model'
317 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'READ BIN NUMBER'
320 read(
io_fid_conf,nml=param_atmos_phy_mp_suzuki10_bin,iostat=ierr)
323 log_info(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Not found namelist. Default used.'
324 elseif( ierr > 0 )
then
325 log_error(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10_bin, Check!'
329 log_nml(param_atmos_phy_mp_suzuki10_bin)
331 if( iceflg == 0 )
then
333 elseif( iceflg == 1 )
then
336 log_error(
"ATMOS_PHY_MP_suzuki10_tracer_setup",*)
"ICEFLG should be 0 (warm rain) or 1 (mixed rain) check!!"
345 num_hyd =
nbin * nspc
347 num_start_waters = i_qv + 1
349 num_start_ices = num_end_waters + 1
410 integer,
intent(in) :: ka
411 integer,
intent(in) :: ia
412 integer,
intent(in) :: ja
413 logical,
intent(in),
optional :: flg_lt
419 real(
rp) :: s10_emaer
421 logical :: s10_flag_regene = .false.
422 logical :: s10_flag_nucleat = .false.
423 logical :: s10_flag_icenucleat = .false.
424 logical :: s10_flag_sfaero = .false.
425 logical :: s10_flag_rndm = .false.
426 integer :: s10_rndm_mspc
427 integer :: s10_rndm_mbin
429 namelist / param_atmos_phy_mp_suzuki10 / &
430 mp_doautoconversion, &
439 s10_flag_icenucleat, &
449 real(
rp),
parameter :: max_term_vel = 10.0_rp
452 integer :: nnspc, nnbin
453 integer :: nn, mm, mmyu, nnyu
454 integer :: myu, nyu, i, j, k, n, ierr
458 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Setup'
459 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Suzuki (2010) Spectral BIN model'
462 allocate( xctr(
nbin ) )
463 allocate( xbnd(
nbin+1 ) )
464 allocate( radc(
nbin ) )
465 allocate( cctr(
nbin,nspc_mk ) )
466 allocate( cbnd(
nbin+1,nspc_mk ) )
467 allocate( ck( nspc_mk,nspc_mk,
nbin,
nbin ) )
468 allocate( vt( nspc_mk,
nbin ) )
469 allocate( br( nspc_mk,
nbin ) )
470 allocate( ifrsl( 2,nspc_mk,nspc_mk ) )
471 allocate( expxctr(
nbin ) )
472 allocate( expxbnd(
nbin+1 ) )
473 allocate( rexpxctr(
nbin ) )
474 allocate( rexpxbnd(
nbin+1 ) )
475 if (
nccn /= 0 )
then
476 allocate( xactr(
nccn ) )
477 allocate( xabnd(
nccn+1 ) )
478 allocate( rada(
nccn ) )
479 allocate( expxactr(
nccn ) )
480 allocate( expxabnd(
nccn+1 ) )
481 allocate( rexpxactr(
nccn ) )
482 allocate( rexpxabnd(
nccn+1 ) )
486 mspc = nspc_mk*nspc_mk
497 read(
io_fid_conf,nml=param_atmos_phy_mp_suzuki10,iostat=ierr)
500 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Not found namelist. Default used.'
501 elseif( ierr > 0 )
then
502 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10, Check!'
505 log_nml(param_atmos_phy_mp_suzuki10)
507 if ( nspc /= 1 .AND. nspc /= 7 )
then
508 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nspc should be set as 1 (warm rain) or 7 (mixed phase) check!'
517 flg_regeneration = s10_flag_regene
518 flg_nucl = s10_flag_nucleat
519 flg_icenucl = s10_flag_icenucleat
520 flg_sf_aero = s10_flag_sfaero
521 flg_rndm = s10_flag_rndm
530 open ( fid_micpara, file = fname_micpara, form =
'formatted', status =
'old', iostat=ierr )
533 if ( ierr == 0 )
then
535 read( fid_micpara,* ) nnspc, nnbin
537 if ( nnbin /=
nbin )
then
538 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nbin in namelist and nbin in micpara.dat is different check!'
543 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Radius of cloud *'
545 read( fid_micpara,* ) nn, xctr( n ), radc( n )
546 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
547 "Radius of ", n,
"th cloud bin (bin center)= ", radc( n ) ,
"[m]"
550 read( fid_micpara,* ) nn, xbnd( n )
552 read( fid_micpara,* ) dxmic
553 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Width of Cloud SDF= ', dxmic
559 read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
563 read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
572 read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
581 read( fid_micpara,* ) mmyu, nn, vt( myu,n )
588 read( fid_micpara,* ) mmyu, nn, br( myu,n )
592 close ( fid_micpara )
597 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'micpara.dat is created'
602 open ( fid_micpara, file = fname_micpara, form =
'formatted', status =
'old', iostat=ierr )
604 read( fid_micpara,* ) nnspc, nnbin
606 if ( nnbin /=
nbin )
then
607 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nbin in inc_tracer and nbin in micpara.dat is different check!'
612 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Radius of cloud *'
614 read( fid_micpara,* ) nn, xctr( n ), radc( n )
615 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
616 "Radius of ", n,
"th cloud bin (bin center)= ", radc( n ) ,
"[m]"
619 read( fid_micpara,* ) nn, xbnd( n )
621 read( fid_micpara,* ) dxmic
622 log_info(
"ATMOS_PHY_MP_suzuki10_setup",*)
'Width of Cloud SDF= ', dxmic
628 read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
632 read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
641 read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
650 read( fid_micpara,* ) mmyu, nn, vt( myu,n )
657 read( fid_micpara,* ) mmyu, nn, br( myu,n )
661 close ( fid_micpara )
667 call comm_bcast( radc(:),
nbin )
668 call comm_bcast( xctr(:),
nbin )
669 call comm_bcast( dxmic )
670 call comm_bcast( xbnd(:),
nbin+1 )
671 call comm_bcast( cctr(:,:),
nbin, nspc_mk )
672 call comm_bcast( cbnd(:,:),
nbin+1, nspc_mk )
673 call comm_bcast( ck(:,:,:,:), nspc_mk, nspc_mk,
nbin,
nbin )
674 call comm_bcast( br(:,:), nspc_mk,
nbin )
675 call comm_bcast( vt(:,:), nspc_mk,
nbin )
677 allocate( flg_noninduct( nspc,nspc ) )
678 allocate( ecoll( nspc,nspc,
nbin,
nbin ) )
679 allocate( rcoll( nspc,nspc,
nbin,
nbin ) )
680 flg_noninduct(:,:) = 0.0_rp
681 ecoll( :,:,:,: ) = 0.0_rp
682 rcoll( :,:,:,: ) = 0.0_rp
684 if (
present(flg_lt) )
then
693 if( ( myu >= ic .and. myu <= iss ) .and. ( nyu == ig .or. nyu == ih ) )
then
694 flg_noninduct( myu,nyu ) = 1.0_rp
703 if( vt( myu,i ) /= vt( nyu,j ) )
then
704 ecoll( myu,nyu,i,j ) = ck( myu,nyu,i,j ) &
705 / ( pi*( radc( i )+radc( j ) )**2 * abs( vt( myu,i )-vt( nyu,j ) ) )
706 ecoll( myu,nyu,i,j ) = max( min( 1.0_rp, ecoll( myu,nyu,i,j ) ),0.0_rp )
708 ecoll( myu,nyu,i,j ) = 0.0_rp
711 if( ecoal_gsi /= 0.0_rp )
then
712 ecoll( myu,nyu,i,j ) = ecoal_gsi
715 if( ecoll( myu,nyu,i,j ) /= 0.0_rp )
then
716 rcoll( myu,nyu,i,j ) = ( 1.0_rp-ecoll( myu,nyu,i,j ) ) / ecoll( myu,nyu,i,j )
717 elseif( ecoll( myu,nyu,i,j ) == 0.0_rp )
then
718 rcoll( myu,nyu,i,j ) = 1.0_rp
728 if (
nccn /= 0 )
then
730 allocate ( ncld( 1:
nccn ) )
731 xasta = log( rhoa*4.0_rp/3.0_rp*pi * ( rasta )**3 )
732 xaend = log( rhoa*4.0_rp/3.0_rp*pi * ( raend )**3 )
734 dxaer = ( xaend-xasta )/
nccn
737 xabnd( n ) = xasta + dxaer*( n-1 )
740 xactr( n ) = ( xabnd( n )+xabnd( n+1 ) )*0.50_rp
741 rada( n ) = ( exp( xactr( n ) )*thirdovforth/pi/rhoa )**( oneovthird )
742 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
743 "Radius of ", n,
"th aerosol bin (bin center)= ", rada( n ) ,
"[m]"
746 if ( flg_sf_aero )
then
747 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
"flg_sf_aero=true is not supported stop!! "
778 if( radc( n ) > rbnd )
then
783 log_info(
"ATMOS_PHY_MP_suzuki10_setup",
'(A,ES15.7,A)')
'Radius between cloud and rain is ', radc(nbnd),
'[m]'
788 call random_setup( ia*ja*ka )
791 if ( mp_couple_aerosol .AND.
nccn /=0 )
then
792 log_error(
"ATMOS_PHY_MP_suzuki10_setup",*)
'nccn should be 0 when MP_couple_aerosol = .true. !! stop'
796 if (
nccn /= 0 )
then
798 expxactr( n ) = exp( xactr( n ) )
799 rexpxactr( n ) = 1.0_rp / exp( xactr( n ) )
802 expxabnd( n ) = exp( xabnd( n ) )
803 rexpxabnd( n ) = 1.0_rp / exp( xabnd( n ) )
807 allocate( vterm(qa-1) )
811 vterm((myu-1)*
nbin+n) = -vt( myu,n )
815 expxctr( n ) = exp( xctr( n ) )
816 rexpxctr( n ) = 1.0_rp / exp( xctr( n ) )
819 expxbnd( n ) = exp( xbnd( n ) )
820 rexpxbnd( n ) = 1.0_rp / exp( xbnd( n ) )
824 call getrule( ifrsl,kindx )
827 sigma_sdf(1) = 0.2_rp
828 sigma_sdf(2) = 0.35_rp
829 sigma_sdf(3) = 0.35_rp
830 sigma_sdf(4) = 0.35_rp
831 sigma_sdf(5) = 0.35_rp
833 r0_sdf(2) = 2.61e-6_rp
835 r0_sdf(4) = 2.61e-6_rp
836 r0_sdf(5) = 2.61e-6_rp
837 n0_sdf(1) = 8.0e+6_rp
839 n0_sdf(3) = 3.0e+6_rp
840 n0_sdf(4) = 4.0e+6_rp
841 n0_sdf(5) = 4.0e+6_rp
844 rho_sdf(3) = 100.0_rp
845 rho_sdf(4) = 400.0_rp
846 rho_sdf(5) = 400.0_rp
849 sigma_sdf(1) = 0.2_rp
850 sigma_sdf(2) = 0.35_rp
851 sigma_sdf(3) = 0.35_rp
852 sigma_sdf(4) = 0.35_rp
853 sigma_sdf(5) = 0.35_rp
855 r0_sdf(2) = 2.61e-6_rp
857 r0_sdf(4) = 2.61e-6_rp
858 r0_sdf(5) = 2.61e-6_rp
859 n0_sdf(1) = 8.0e+6_rp
861 n0_sdf(3) = 3.0e+6_rp
862 n0_sdf(4) = 4.0e+6_rp
863 n0_sdf(5) = 4.0e+6_rp
866 rho_sdf(3) = 100.0_rp
867 rho_sdf(4) = 400.0_rp
868 rho_sdf(5) = 400.0_rp
876 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
897 atmos_saturation_pres2qsat_liq, &
898 atmos_saturation_pres2qsat_ice
901 integer,
intent(in) :: ka, ks, ke
902 integer,
intent(in) :: ia, is, ie
903 integer,
intent(in) :: ja, js, je
904 integer,
intent(in) :: kijmax
906 real(
dp),
intent(in) :: dt
907 real(
rp),
intent(in) :: dens (ka,ia,ja)
908 real(
rp),
intent(in) :: pres (ka,ia,ja)
909 real(
rp),
intent(in) :: temp (ka,ia,ja)
910 real(
rp),
intent(in) :: qtrc (ka,ia,ja,qa)
911 real(
rp),
intent(in) :: qdry (ka,ia,ja)
912 real(
rp),
intent(in) :: cptot(ka,ia,ja)
913 real(
rp),
intent(in) :: cvtot(ka,ia,ja)
914 real(
rp),
intent(in) :: ccn (ka,ia,ja)
916 real(
rp),
intent(out) :: rhoq_t (ka,ia,ja,qa)
917 real(
rp),
intent(out) :: rhoe_t (ka,ia,ja)
918 real(
rp),
intent(out) :: cptot_t(ka,ia,ja)
919 real(
rp),
intent(out) :: cvtot_t(ka,ia,ja)
920 real(
rp),
intent(out) :: evaporate(ka,ia,ja)
923 logical,
intent(in),
optional :: flg_lt
924 real(
rp),
intent(in),
optional :: d0_crg, v0_crg
925 real(
rp),
intent(in),
optional :: dqcrg(ka,ia,ja)
926 real(
rp),
intent(in),
optional :: beta_crg(ka,ia,ja)
927 real(
rp),
intent(in),
optional :: qtrc_crg(ka,ia,ja,num_hyd)
928 real(
rp),
intent(out),
optional :: qsplt_in(ka,ia,ja,3)
929 real(
rp),
intent(out),
optional :: sarea(ka,ia,ja,num_hyd)
930 real(
rp),
intent(out),
optional :: rhoc_t_mp(ka,ia,ja,num_hyd)
932 real(
rp) :: qsat_l(ka,ia,ja)
933 real(
rp) :: qsat_i(ka,ia,ja)
934 real(
rp) :: ssliq(ka,ia,ja)
935 real(
rp) :: ssice(ka,ia,ja)
937 integer :: ijk_index (kijmax,3)
938 integer :: index_cld (kijmax)
939 integer :: index_cold(kijmax)
940 integer :: index_warm(kijmax)
941 integer :: ijkcount, ijkcount_cold, ijkcount_warm
942 integer :: ijk, indirect
944 real(
rp) :: dens_ijk(kijmax)
945 real(
rp) :: pres_ijk(kijmax)
946 real(
rp) :: temp_ijk(kijmax)
947 real(
rp) :: qdry_ijk(kijmax)
948 real(
rp) :: qvap_ijk(kijmax)
949 real(
rp) :: ccn_ijk(kijmax)
950 real(
rp) :: cp_ijk(kijmax)
951 real(
rp) :: cv_ijk(kijmax)
952 real(
rp) :: evaporate_ijk(kijmax)
953 real(
rp) :: ghyd_ijk(
nbin,nspc,kijmax)
954 real(
rp) :: gaer_ijk(max(
nccn,1),kijmax)
961 real(
rp) :: gcrg_ijk(
nbin,nspc,kijmax)
962 real(
rp) :: crg_sep_ijk(nspc,kijmax)
963 real(
rp) :: dqcrg_ijk(kijmax)
964 real(
rp) :: beta_crg_ijk(kijmax)
967 integer :: k, i, j, m, n, iq
970 if ( nspc == 1 )
then
971 log_progress(*)
'atmosphere / physics / microphysics / SBM (Liquid water only)'
972 elseif( nspc > 1 )
then
973 log_progress(*)
'atmosphere / physics / microphysics / SBM (Mixed phase)'
976 if (
present(flg_lt) )
then
978 crg_sep_ijk(:,:) = 0.0_rp
979 qsplt_in(:,:,:,:) = 0.0_rp
984 call atmos_saturation_pres2qsat_liq( ka, ks, ke, &
987 temp(:,:,:), pres(:,:,:), qdry(:,:,:), &
990 call atmos_saturation_pres2qsat_ice( ka, ks, ke, &
993 temp(:,:,:), pres(:,:,:), qdry(:,:,:), &
999 ssliq(k,i,j) = qtrc(k,i,j,i_qv) / qsat_l(k,i,j) - 1.0_rp
1000 ssice(k,i,j) = qtrc(k,i,j,i_qv) / qsat_i(k,i,j) - 1.0_rp
1005 if ( nspc == 1 )
then
1006 ssice(:,:,:) = 0.0_rp
1039 ijk_index(ijk,1) = i
1040 ijk_index(ijk,2) = j
1041 ijk_index(ijk,3) = k
1045 gcrg_ijk(:,:,:) = 0.0_rp
1061 cldsum = cldsum + qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
1062 countbin = countbin + 1
1066 if ( cldsum > cldmin &
1067 .OR. ssliq(k,i,j) > 0.0_rp &
1068 .OR. ssice(k,i,j) > 0.0_rp )
then
1070 ijkcount = ijkcount + 1
1072 index_cld(ijkcount) = ijk
1074 dens_ijk(ijkcount) = dens(k,i,j)
1075 pres_ijk(ijkcount) = pres(k,i,j)
1076 temp_ijk(ijkcount) = temp(k,i,j)
1077 qdry_ijk(ijkcount) = qdry(k,i,j)
1078 cp_ijk(ijkcount) = cptot(k,i,j)
1079 cv_ijk(ijkcount) = cvtot(k,i,j)
1080 ccn_ijk(ijkcount) = ccn(k,i,j)
1081 qvap_ijk(ijkcount) = qtrc(k,i,j,i_qv)
1086 ghyd_ijk(n,m,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
1087 countbin = countbin + 1
1092 gaer_ijk(n,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxaer
1093 countbin = countbin + 1
1096 if ( temp(k,i,j) < tem00 .AND. nspc > 1 )
then
1097 ijkcount_cold = ijkcount_cold + 1
1098 index_cold(ijkcount_cold) = ijkcount
1100 ijkcount_warm = ijkcount_warm + 1
1101 index_warm(ijkcount_warm) = ijkcount
1104 if ( flg_lt_l )
then
1108 gcrg_ijk(n,m,ijkcount) = qtrc_crg(k,i,j,countbin) * dens(k,i,j)
1109 countbin = countbin + 1
1112 beta_crg_ijk(ijkcount) = beta_crg(k,i,j)
1113 dqcrg_ijk(ijkcount) = dqcrg(k,i,j)
1120 rhoq_t(k,i,j,iq) = 0.0_rp
1122 rhoe_t(k,i,j) = 0.0_rp
1123 cptot_t(k,i,j) = 0.0_rp
1124 cvtot_t(k,i,j) = 0.0_rp
1125 evaporate(k,i,j) = 0.0_rp
1167 if ( ijkcount > 0 )
then
1171 if ( flg_lt_l )
then
1173 call mp_suzuki10( ka, ia, ja, &
1177 index_cold( 1:ijkcount), &
1178 index_warm( 1:ijkcount), &
1179 dens_ijk( 1:ijkcount), &
1180 pres_ijk( 1:ijkcount), &
1181 qdry_ijk( 1:ijkcount), &
1182 ccn_ijk( 1:ijkcount), &
1183 temp_ijk( 1:ijkcount), &
1184 qvap_ijk( 1:ijkcount), &
1185 ghyd_ijk(:,:,1:ijkcount), &
1186 gaer_ijk(:, 1:ijkcount), &
1187 cp_ijk( 1:ijkcount), &
1188 cv_ijk( 1:ijkcount), &
1189 evaporate_ijk(1:ijkcount), &
1191 flg_lt_l, d0_crg, v0_crg, &
1192 dqcrg_ijk( 1:ijkcount), &
1193 beta_crg_ijk( 1:ijkcount), &
1194 gcrg_ijk(:,:,1:ijkcount), &
1195 crg_sep_ijk(:,1:ijkcount) )
1197 call mp_suzuki10( ka, ia, ja, &
1201 index_cold( 1:ijkcount), &
1202 index_warm( 1:ijkcount), &
1203 dens_ijk( 1:ijkcount), &
1204 pres_ijk( 1:ijkcount), &
1205 qdry_ijk( 1:ijkcount), &
1206 ccn_ijk( 1:ijkcount), &
1207 temp_ijk( 1:ijkcount), &
1208 qvap_ijk( 1:ijkcount), &
1209 ghyd_ijk(:,:,1:ijkcount), &
1210 gaer_ijk(:, 1:ijkcount), &
1211 cp_ijk( 1:ijkcount), &
1212 cv_ijk( 1:ijkcount), &
1213 evaporate_ijk(1:ijkcount), &
1264 do ijk = 1, ijkcount
1265 indirect = index_cld(ijk)
1266 i = ijk_index(indirect,1)
1267 j = ijk_index(indirect,2)
1268 k = ijk_index(indirect,3)
1270 rhoe_t(k,i,j) = ( temp_ijk(ijk) * cv_ijk(ijk) - temp(k,i,j) * cvtot(k,i,j) ) * dens(k,i,j) / dt
1271 cptot_t(k,i,j) = ( cp_ijk(ijk) - cptot(k,i,j) ) / dt
1272 cvtot_t(k,i,j) = ( cv_ijk(ijk) - cvtot(k,i,j) ) / dt
1273 evaporate(k,i,j) = evaporate_ijk(ijk) / dt
1275 rhoq_t(k,i,j,i_qv) = ( qvap_ijk(ijk) - qtrc(k,i,j,i_qv) ) * dens(k,i,j) / dt
1280 rhoq_new = ghyd_ijk(n,m,ijk) * dxmic
1281 rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
1282 countbin = countbin + 1
1287 rhoq_new = gaer_ijk(n,ijk) * dxaer
1288 rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
1289 countbin = countbin + 1
1296 rhoq_new = gcrg_ijk(n,m,ijk)
1297 rhoc_t_mp(k,i,j,countbin) = ( rhoq_new - qtrc_crg(k,i,j,countbin)*dens(k,i,j) ) / dt
1298 countbin = countbin + 1
1302 sarea(k,i,j,:) = 0.0_rp
1306 rhoq_new = ghyd_ijk(n,m,ijk) * dxmic
1307 sarea(k,i,j,countbin) = 4.0_rp*pi*radc(n)**2*rhoq_new*rexpxctr( n )
1308 countbin = countbin + 1
1312 qsplt_in(k,i,j,1) = crg_sep_ijk(ig,ijk) * dens(k,i,j) / dt
1313 qsplt_in(k,i,j,2) = ( crg_sep_ijk(ic,ijk) &
1314 + crg_sep_ijk(ip,ijk) &
1315 + crg_sep_ijk(id,ijk) ) * dens(k,i,j) / dt
1316 qsplt_in(k,i,j,3) = crg_sep_ijk(iss,ijk) * dens(k,i,j) / dt
1351 integer,
intent(in) :: ka
1353 real(
rp),
intent(out) :: vterm_o(ka,qa-1)
1359 vterm_o(:,iq) = vterm(iq)
1376 integer,
intent(in) :: ka, ks, ke
1377 integer,
intent(in) :: ia, is, ie
1378 integer,
intent(in) :: ja, js, je
1380 real(
rp),
intent(in) :: qtrc0 (ka,ia,ja,num_hyd)
1381 real(
rp),
intent(in) :: mask_criterion
1382 real(
rp),
intent(out) :: cldfrac(ka,ia,ja)
1385 integer :: k, i, j, iq, ihydro
1394 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1395 qhydro = qhydro + qtrc0(k,i,j,iq)
1398 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1402 elseif( nspc == 1 )
then
1407 do ihydro = 1, i_mp_qc
1408 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1409 qhydro = qhydro + qtrc0(k,i,j,iq)
1412 cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1444 integer,
intent(in) :: ka, ks, ke
1445 integer,
intent(in) :: ia, is, ie
1446 integer,
intent(in) :: ja, js, je
1448 real(
rp),
intent(in) :: dens0(ka,ia,ja)
1449 real(
rp),
intent(in) :: temp0(ka,ia,ja)
1450 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1451 real(
rp),
intent(out) :: re (ka,ia,ja,
n_hyd)
1453 real(
rp),
parameter :: um2cm = 100.0_rp
1455 real(
rp) :: sum0(nspc), sum2, sum3, re_tmp(nspc)
1456 integer :: i, j, k, iq, ihydro
1462 re(k,i,j,:) = 0.0_rp
1470 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1471 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1472 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1474 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1475 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1476 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1478 sum3 = max( sum3, 0.0_rp )
1479 sum2 = max( sum2, 0.0_rp )
1480 if ( sum2 /= 0.0_rp )
then
1481 re(k,i,j,
i_hc) = sum3 / sum2 * um2cm
1483 re(k,i,j,
i_hc) = 0.0_rp
1490 do iq = nbnd+1,
nbin
1492 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1493 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1494 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1496 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1497 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1498 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1500 sum3 = max( sum3, 0.0_rp )
1501 sum2 = max( sum2, 0.0_rp )
1502 if ( sum2 /= 0.0_rp )
then
1503 re(k,i,j,
i_hr) = sum3 / sum2 * um2cm
1505 re(k,i,j,
i_hr) = 0.0_rp
1513 if ( nspc > 1 )
then
1518 sum0(ihydro) = 0.0_rp
1521 do iq =
nbin*(ihydro-1)+1,
nbin*ihydro
1522 sum0(ihydro) = sum0(ihydro) &
1523 + ( qtrc0(k,i,j,iq) * dens0(k,i,j) )
1525 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1526 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1527 * radc( iq-
nbin*(ihydro-1) )**3.0_rp )
1529 + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) &
1530 * rexpxctr( iq-
nbin*(ihydro-1) ) &
1531 * radc( iq-
nbin*(ihydro-1) )**2.0_rp )
1533 sum3 = max( sum3, 0.0_rp )
1534 sum2 = max( sum2, 0.0_rp )
1535 if ( sum2 == 0.0_rp )
then
1536 re_tmp(ihydro) = 0.0_rp
1538 re_tmp(ihydro) = sum3 / sum2 * um2cm
1542 re(k,i,j,
i_hi) = ( re_tmp(i_mp_qcl) * sum0(i_mp_qcl) &
1543 + re_tmp(i_mp_qp ) * sum0(i_mp_qp ) &
1544 + re_tmp(i_mp_qd ) * sum0(i_mp_qd ) ) &
1545 / ( sum0(i_mp_qcl) + sum0(i_mp_qp) + sum0(i_mp_qd) + eps )
1546 re(k,i,j,
i_hs) = re_tmp(i_mp_qs)
1547 re(k,i,j,
i_hg) = re_tmp(i_mp_qg)
1548 re(k,i,j,
i_hh) = re_tmp(i_mp_qh)
1577 integer,
intent(in) :: ka, ks, ke
1578 integer,
intent(in) :: ia, is, ie
1579 integer,
intent(in) :: ja, js, je
1580 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1581 real(
rp),
intent(out) :: qe (ka,ia,ja,
n_hyd)
1583 integer :: ihydro, ibin, iq, icateg
1588 qe(:,:,:,:) = 0.0_rp
1592 iq =
nbin*(ihydro-1) + ibin
1594 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc ) )
then
1595 if ( iq > 0 .AND. iq <= nbnd )
then
1597 elseif( iq > nbnd .AND. iq <=
nbin )
then
1600 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qcl) )
then
1602 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qp ) )
then
1604 elseif( iq >
nbin*(i_mp_qd -1) .AND. iq <=
nbin*(i_mp_qd ) )
then
1606 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qs ) )
then
1608 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qg ) )
then
1610 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <= num_hyd )
then
1617 qe(k,i,j,icateg) = qe(k,i,j,icateg) + qtrc0(k,i,j,iq)
1645 integer,
intent(in) :: ka, ks, ke
1646 integer,
intent(in) :: ia, is, ie
1647 integer,
intent(in) :: ja, js, je
1648 real(
rp),
intent(in) :: dens (ka,ia,ja)
1649 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1650 real(
rp),
intent(out) :: ne (ka,ia,ja,
n_hyd)
1652 integer :: ihydro, ibin, iq, icateg
1657 ne(:,:,:,:) = 0.0_rp
1661 iq =
nbin*(ihydro-1) + ibin
1663 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc ) )
then
1664 if ( iq > 0 .AND. iq <= nbnd )
then
1666 elseif( iq > nbnd .AND. iq <=
nbin )
then
1669 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qcl) )
then
1671 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qp ) )
then
1673 elseif( iq >
nbin*(i_mp_qd -1) .AND. iq <=
nbin*(i_mp_qd ) )
then
1675 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qs ) )
then
1677 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qg ) )
then
1679 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <= num_hyd )
then
1686 ne(k,i,j,icateg) = ne(k,i,j,icateg) + dens(k,i,j) * qtrc0(k,i,j,iq) * rexpxctr(ibin)
1699 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1717 integer,
intent(in) :: ka, ks, ke
1718 integer,
intent(in) :: ia, is, ie
1719 integer,
intent(in) :: ja, js, je
1721 real(
rp),
intent(in) :: qe(ka,ia,ja,
n_hyd)
1723 real(
rp),
intent(out) :: qtrc(ka,ia,ja,qa-1)
1725 real(
rp),
intent(in),
optional :: qnum(ka,ia,ja,
n_hyd)
1727 real(
rp) :: coef0, coef1, coef2
1729 real(
rp) :: tmp_hyd, num_hyd_l, lambda_hyd
1731 integer :: k, i, j, iq
1733 if (
present(qnum) )
then
1734 log_warn(
"ATMOS_PHY_MP_suzuki10_qhyd2qtrc",*)
'At this moment, number concentratio is ignored'
1738 coef0 = 4.0_rp/3.0_rp*pi
1739 coef1 = 4.0_rp/3.0_rp*sqrt(pi/2.0_rp)
1741 if( nspc == 1 )
then
1749 dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
1751 - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
1752 / sigma_sdf(1) / sigma_sdf(1) &
1754 tmp_hyd = tmp_hyd + dummy(iq)
1757 coef2 = ( qe(k,i,j,
i_hc) + qe(k,i,j,
i_hr) &
1759 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1762 qtrc(k+2,i,j,(il-1)*
nbin+iq) = coef2 * dummy(iq)
1769 elseif( nspc > 1 )
then
1778 dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
1780 - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
1781 / sigma_sdf(1) / sigma_sdf(1) &
1783 tmp_hyd = tmp_hyd + dummy(iq)
1786 coef2 = ( qe(k,i,j,
i_hc) + qe(k,i,j,
i_hr) ) &
1787 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1790 qtrc(k,i,j,(il-1)*
nbin+iq) = coef2 * dummy(iq)
1796 dummy(iq) = coef1 / sigma_sdf(2) * rho_sdf(2) * radc( iq )**3 &
1798 - ( log( radc(iq) )-log( r0_sdf(2) ) )**2*0.5_rp &
1799 / sigma_sdf(2) / sigma_sdf(2) &
1801 tmp_hyd = tmp_hyd + dummy(iq)
1804 coef2 = ( qe(k,i,j,
i_hi) ) &
1805 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1808 qtrc(k,i,j,(ip-1)*
nbin+iq) = coef2 * dummy(iq)
1812 num_hyd_l = coef0 * n0_sdf(3) * rho_sdf(3)
1813 lambda_hyd = ( pi * rho_sdf(3) / 6.0_rp *n0_sdf(3) *
sf_gamma(4.0_rp) &
1814 / ( qe(k,i,j,
i_hs) &
1815 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hs)-eps)) &
1820 dummy(iq) = num_hyd_l * radc( iq )**3 &
1821 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1822 tmp_hyd = tmp_hyd + dummy(iq)
1825 coef2 = ( qe(k,i,j,
i_hs) ) &
1826 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1829 qtrc(k,i,j,(iss-1)*
nbin+iq) = coef2 * dummy(iq)
1833 num_hyd_l = coef0 * n0_sdf(4) * rho_sdf(4)
1834 lambda_hyd = ( pi * rho_sdf(4) / 6.0_rp *n0_sdf(4) *
sf_gamma(4.0_rp) &
1835 / ( qe(k,i,j,
i_hg) &
1836 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hg)-eps)) &
1841 dummy(iq) = num_hyd_l * radc( iq )**3 &
1842 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1843 tmp_hyd = tmp_hyd + dummy(iq)
1846 coef2 = ( qe(k,i,j,
i_hg) ) &
1847 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1850 qtrc(k,i,j,(ig-1)*
nbin+iq) = coef2 * dummy(iq)
1854 num_hyd_l = coef0 * n0_sdf(5) * rho_sdf(5)
1855 lambda_hyd = ( pi * rho_sdf(5) / 6.0_rp *n0_sdf(5) *
sf_gamma(4.0_rp) &
1856 / ( qe(k,i,j,
i_hh) &
1857 + (0.50_rp-sign(0.50_rp,qe(k,i,j,
i_hh)-eps)) &
1862 dummy(iq) = num_hyd_l * radc( iq )**3 &
1863 * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
1864 tmp_hyd = tmp_hyd + dummy(iq)
1867 coef2 = ( qe(k,i,j,
i_hh) ) &
1868 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
1871 qtrc(k,i,j,(ih-1)*
nbin+iq) = coef2 * dummy(iq)
1880 do iq = num_hyd+1, qa-1
1884 qtrc(k,i,j,iq) = 0.0_rp
1911 integer,
intent(in) :: ka, ks, ke
1912 integer,
intent(in) :: ia, is, ie
1913 integer,
intent(in) :: ja, js, je
1914 real(
rp),
intent(in) :: qtrc0(ka,ia,ja,num_hyd)
1915 real(
rp),
intent(out) :: qecrg(ka,ia,ja,
n_hyd)
1917 integer :: ihydro, ibin, iq, icateg
1922 qecrg(:,:,:,:) = 0.0_rp
1926 iq =
nbin*(ihydro-1) + ibin
1928 if ( iq > 0 .AND. iq <=
nbin*(i_mp_qc ) )
then
1929 if ( iq > 0 .AND. iq <= nbnd )
then
1931 elseif( iq > nbnd .AND. iq <=
nbin )
then
1934 elseif( iq >
nbin*(i_mp_qcl-1) .AND. iq <=
nbin*(i_mp_qcl) )
then
1936 elseif( iq >
nbin*(i_mp_qp -1) .AND. iq <=
nbin*(i_mp_qp ) )
then
1938 elseif( iq >
nbin*(i_mp_qd -1) .AND. iq <=
nbin*(i_mp_qd ) )
then
1940 elseif( iq >
nbin*(i_mp_qs -1) .AND. iq <=
nbin*(i_mp_qs ) )
then
1942 elseif( iq >
nbin*(i_mp_qg -1) .AND. iq <=
nbin*(i_mp_qg ) )
then
1944 elseif( iq >
nbin*(i_mp_qh -1) .AND. iq <= num_hyd )
then
1951 qecrg(k,i,j,icateg) = qecrg(k,i,j,icateg) + qtrc0(k,i,j,iq)
1962 subroutine mp_suzuki10( &
1992 integer,
intent(in) :: ka
1993 integer,
intent(in) :: ia
1994 integer,
intent(in) :: ja
1996 integer,
intent(in) :: ijkmax
1997 integer,
intent(in) :: ijkmax_cold
1998 integer,
intent(in) :: ijkmax_warm
1999 integer ,
intent(in) :: index_cold(ijkmax)
2000 integer ,
intent(in) :: index_warm(ijkmax)
2001 real(
rp),
intent(in) :: dens (ijkmax)
2002 real(
rp),
intent(in) :: pres (ijkmax)
2003 real(
rp),
intent(in) :: qdry (ijkmax)
2004 real(
rp),
intent(in) :: ccn (ijkmax)
2005 real(
rp),
intent(inout) :: temp (ijkmax)
2006 real(
rp),
intent(inout) :: qvap (ijkmax)
2007 real(
rp),
intent(inout) :: ghyd (
nbin,nspc,ijkmax)
2008 real(
rp),
intent(inout) :: gaer (max(
nccn,1),ijkmax)
2009 real(
rp),
intent(inout) :: cp (ijkmax)
2010 real(
rp),
intent(inout) :: cv (ijkmax)
2011 real(
rp),
intent(out) :: evaporate (ijkmax)
2012 real(
dp),
intent(in) :: dt
2015 logical,
intent(in),
optional :: flg_lt
2016 real(
rp),
intent(in),
optional :: d0_crg, v0_crg
2017 real(
rp),
intent(in),
optional :: dqcrg(ijkmax), beta_crg(ijkmax)
2018 real(
rp),
intent(inout),
optional :: gcrg(
nbin,nspc,ijkmax)
2019 real(
rp),
intent(out),
optional :: crg_sep(nspc,ijkmax)
2022 real(
rp) :: gcrg_l(
nbin,nspc,ijkmax), crg_sep_l(nspc,ijkmax)
2023 real(
rp) :: csum(il,ijkmax)
2025 real(
rp) :: v0_crg_l, d0_crg_l, tcrglimit_l
2028 if (
present(flg_lt) )
then
2035 gcrg_l(:,:,:) = gcrg(:,:,:)
2038 crg_sep_l(:,:) = 0.0_rp
2040 gcrg_l(:,:,:) = 0.0_rp
2041 d0_crg_l = 100.e-6_rp
2043 crg_sep_l(:,:) = 0.0_rp
2046 if (
nccn /= 0 )
then
2047 if ( nspc == 1 )
then
2051 call nucleata( ijkmax, &
2064 call cndevpsbla( ijkmax, &
2078 if ( mp_doautoconversion )
then
2081 call collmain( ka, ia, ja, &
2095 crg_sep(:,:) = crg_sep_l(:,:)
2100 elseif( nspc > 1 )
then
2104 call nucleata( ijkmax, &
2117 call freezing( ijkmax, &
2129 if( flg_icenucl )
then
2130 call ice_nucleat( ijkmax, &
2144 call melting( ijkmax, &
2157 call cndevpsbla( ijkmax, &
2171 if ( mp_doautoconversion )
then
2173 call collmainf( ka, ia, ja, &
2187 crg_sep(:,:) = crg_sep_l(:,:)
2194 elseif(
nccn == 0 )
then
2196 if ( nspc == 1 )
then
2200 call nucleat( ijkmax, &
2213 call cndevpsbl( ijkmax, &
2226 if ( mp_doautoconversion )
then
2228 call collmain( ka, ia, ja, &
2242 crg_sep(:,:) = crg_sep_l(:,:)
2247 elseif( nspc > 1 )
then
2251 call nucleat( ijkmax, &
2264 call freezing( ijkmax, &
2276 if( flg_icenucl )
then
2277 call ice_nucleat( ijkmax, &
2291 call melting( ijkmax, &
2304 call cndevpsbl( ijkmax, &
2317 if ( mp_doautoconversion )
then
2319 call collmainf( ka, ia, ja, &
2333 crg_sep(:,:) = crg_sep_l(:,:)
2343 gcrg(:,:,:) = gcrg_l(:,:,:)
2347 end subroutine mp_suzuki10
2350 subroutine nucleat( &
2363 atmos_hydrometeor_lhv, &
2369 atmos_saturation_pres2qsat_liq
2372 integer,
intent(in) :: ijkmax
2373 real(
rp),
intent(in) :: dens(ijkmax)
2374 real(
rp),
intent(in) :: pres(ijkmax)
2375 real(
rp),
intent(in) :: qdry(ijkmax)
2376 real(
rp),
intent(in) :: ccn(ijkmax)
2377 real(
rp),
intent(inout) :: temp(ijkmax)
2378 real(
rp),
intent(inout) :: qvap(ijkmax)
2379 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2380 real(
rp),
intent(inout) :: cp(ijkmax)
2381 real(
rp),
intent(inout) :: cv(ijkmax)
2382 real(
dp),
intent(in) :: dtime
2384 real(
rp) :: ssliq(ijkmax)
2385 real(
rp) :: qlevp(ijkmax)
2391 real(
rp) :: sumnum(ijkmax)
2392 real(
rp) :: gcn(
nbin,ijkmax )
2393 real(
rp) :: qsat(ijkmax)
2399 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2401 if( mp_couple_aerosol )
then
2405 dmp = ccn(ijk) * expxctr( 1 )
2406 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2407 gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
2409 qvap(ijk) = qvap(ijk) - dqv
2410 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
2417 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2418 temp(:), pres(:), qdry(:), &
2422 ssliq(ijk) = qvap(ijk)/qsat(ijk) - 1.0_rp
2428 if ( ssliq(ijk) > 0.0_rp )
then
2433 gcn( n,ijk ) = gc( n,il,ijk )*rexpxctr( n )
2437 sumnum(ijk) = sumnum(ijk) + gcn( n,ijk )*dxmic
2439 n_c = c_ccn * ( ssliq(ijk) * 1.e+2_rp )**( kappa )
2440 if ( n_c > sumnum(ijk) )
then
2441 dmp = ( n_c - sumnum(ijk) ) * expxctr( 1 )
2442 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2443 gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
2445 qvap(ijk) = qvap(ijk) - dqv
2446 qvap(ijk) = max( qvap(ijk),0.0_rp )
2447 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
2460 end subroutine nucleat
2463 subroutine nucleata( &
2480 atmos_hydrometeor_lhv, &
2486 atmos_saturation_pres2qsat_liq, &
2487 atmos_saturation_pres2qsat_ice
2490 integer,
intent(in) :: ijkmax
2491 real(
rp),
intent(in) :: dens(ijkmax)
2492 real(
rp),
intent(in) :: pres(ijkmax)
2493 real(
rp),
intent(in) :: qdry(ijkmax)
2494 real(
rp),
intent(inout) :: temp(ijkmax)
2495 real(
rp),
intent(inout) :: qvap(ijkmax)
2496 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2497 real(
rp),
intent(inout) :: ga (
nccn ,ijkmax)
2498 real(
rp),
intent(inout) :: cp(ijkmax)
2499 real(
rp),
intent(inout) :: cv(ijkmax)
2500 real(
dp),
intent(in) :: dtime
2504 real(
rp) :: acoef, bcoef
2507 real(
rp) :: ractr, rcld, xcld, part, dmp
2508 integer :: n, nc, ncrit
2513 real(
rp) :: qlevp(ijkmax)
2514 real(
rp) :: qsatl(ijkmax)
2521 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2523 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2524 temp(:), pres(:), qdry(:), &
2528 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2530 if ( ssliq <= 0.0_rp ) cycle
2535 gan( n ) = ga( n,ijk )*rexpxactr( n )
2538 acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk)
2539 bcoef = vhfct* rhoa/rhow * emwtr/emaer
2543 ractr = ( expxactr( n )*thirdovforth/pi/rhoa )**( oneovthird )
2544 rcld = sqrt( 3.0_rp*bcoef*ractr*ractr*ractr / acoef )
2545 xcld = log( rhow * 4.0_rp*pi*oneovthird*rcld*rcld*rcld )
2546 if ( flg_nucl )
then
2549 ncld( n ) = int( ( xcld-xctr( 1 ) )/dxmic ) + 1
2550 ncld( n ) = min( max( ncld( n ),1 ),
nbin )
2556 if ( ssliq <= 0.0_rp )
exit
2558 acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk)
2559 rcrit = acoef*oneovthird * ( 4.0_rp/bcoef )**( oneovthird ) / ssliq**( twoovthird )
2560 xcrit = log( rhoa * 4.0_rp*pi*oneovthird * rcrit*rcrit*rcrit )
2561 ncrit = int( ( xcrit-xabnd( 1 ) )/dxaer ) + 1
2563 if ( n == ncrit )
then
2564 part = ( xabnd( ncrit+1 )-xcrit )/dxaer
2565 elseif ( n > ncrit )
then
2573 dmp = part*gan( n )*dxaer*expxctr( nc )
2574 dmp = min( dmp,qvap(ijk)*dens(ijk) )
2575 gc( nc,il,ijk ) = gc( nc,il,ijk ) + dmp/dxmic
2576 gan( n ) = gan( n ) - dmp/dxaer*rexpxctr( nc )
2577 gan( n ) = max( gan( n ), 0.0_rp )
2579 qvap(ijk) = qvap(ijk) - dqv
2580 qvap(ijk) = max( qvap(ijk),0.0_rp )
2581 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
2588 ga( n,ijk ) = gan( n )*expxactr( n )
2596 end subroutine nucleata
2599 subroutine cndevpsbl( &
2614 integer,
intent(in) :: ijkmax
2615 real(
rp),
intent(in) :: dens(ijkmax)
2616 real(
rp),
intent(in) :: pres(ijkmax)
2617 real(
rp),
intent(in) :: qdry(ijkmax)
2618 real(
rp),
intent(inout) :: temp(ijkmax)
2619 real(
rp),
intent(inout) :: qvap(ijkmax)
2620 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2621 real(
rp),
intent(inout) :: cp(ijkmax)
2622 real(
rp),
intent(inout) :: cv(ijkmax)
2623 real(
rp),
intent(out) :: evaporate(ijkmax)
2624 real(
dp),
intent(in) :: dtime
2625 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2628 call liqphase( ijkmax, &
2642 call icephase( ijkmax, &
2654 call mixphase( ijkmax, &
2668 end subroutine cndevpsbl
2671 subroutine cndevpsbla( &
2687 integer,
intent(in) :: ijkmax
2688 real(
rp),
intent(in) :: dens(ijkmax)
2689 real(
rp),
intent(in) :: pres(ijkmax)
2690 real(
rp),
intent(in) :: qdry(ijkmax)
2691 real(
rp),
intent(inout) :: temp(ijkmax)
2692 real(
rp),
intent(inout) :: qvap(ijkmax)
2693 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2694 real(
rp),
intent(inout) :: ga (
nccn ,ijkmax)
2695 real(
rp),
intent(inout) :: cp(ijkmax)
2696 real(
rp),
intent(inout) :: cv(ijkmax)
2697 real(
rp),
intent(out) :: evaporate(ijkmax)
2698 real(
dp),
intent(in) :: dtime
2699 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2702 call liqphase( ijkmax, &
2716 if ( flg_regeneration )
then
2717 call faero( ijkmax, &
2724 call icephase( ijkmax, &
2736 call mixphase( ijkmax, &
2750 end subroutine cndevpsbla
2753 subroutine liqphase( &
2774 atmos_hydrometeor_lhv, &
2780 atmos_saturation_pres2qsat_liq
2783 integer,
intent(in) :: ijkmax
2784 real(
rp),
intent(in) :: dens(ijkmax)
2785 real(
rp),
intent(in) :: pres(ijkmax)
2786 real(
rp),
intent(in) :: qdry(ijkmax)
2787 real(
rp),
intent(inout) :: temp(ijkmax)
2788 real(
rp),
intent(inout) :: qvap(ijkmax)
2789 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
2790 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
2791 real(
rp),
intent(inout) :: cp(ijkmax)
2792 real(
rp),
intent(inout) :: cv(ijkmax)
2793 real(
rp),
intent(out) :: regene_gcn(ijkmax)
2794 real(
dp),
intent(in) :: dtime
2796 integer :: n, myu, ncount
2797 integer :: nloop(ijkmax)
2798 real(
rp) :: gclold(ijkmax)
2799 real(
rp) :: gclnew(ijkmax)
2801 real(
rp) :: dtcnd(ijkmax)
2802 real(
rp) :: gtliq(ijkmax)
2803 real(
rp) :: qlevp(ijkmax)
2804 real(
rp) :: cefliq, a, sliqtnd
2805 real(
rp) :: sumliq(ijkmax), umax(ijkmax)
2807 real(
rp) :: gcn( -1:
nbin+2,nspc,ijkmax )
2809 real(
rp),
parameter :: cflfct = 0.50_rp
2811 integer :: iflg( nspc,ijkmax )
2812 real(
rp) :: csum( nspc,ijkmax )
2813 real(
rp) :: f1, f2, emu, cefd, cefk, festl
2814 real(
rp) :: qsatl(ijkmax)
2815 real(
rp),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2816 real(
rp),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2817 real(
rp) :: zerosw, qvtmp
2819 integer :: ijk, nn, mm, loopflg(ijkmax)
2822 real(
rp) :: uadv ( 0:
nbin+2,nspc,ijkmax )
2823 real(
rp) :: flq ( 1:
nbin+1,nspc,ijkmax )
2824 real(
rp) :: acoef( 0:2,0:
nbin+1,nspc,ijkmax )
2825 real(
rp) :: crn ( 0:
nbin+2,nspc,ijkmax )
2826 real(
rp) :: aip ( 0:
nbin+1,nspc,ijkmax )
2827 real(
rp) :: aim ( 0:
nbin+1,nspc,ijkmax )
2828 real(
rp) :: ai ( 0:
nbin+1,nspc,ijkmax )
2829 real(
rp) :: cmins, cplus
2833 real(
rp) :: gcrgn( -1:
nbin+2,nspc,ijkmax )
2834 real(
rp) :: flq_c( 1:
nbin+1,nspc,ijkmax )
2835 real(
rp) :: acoef_c( 0:2,0:
nbin+1,nspc,ijkmax )
2836 real(
rp) :: aip_c ( 0:
nbin+1,nspc,ijkmax )
2837 real(
rp) :: aim_c ( 0:
nbin+1,nspc,ijkmax )
2838 real(
rp) :: ai_c ( 0:
nbin+1,nspc,ijkmax )
2846 csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
2848 if( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
2853 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2855 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2856 temp(:), pres(:), qdry(:), &
2864 gclold(ijk) = gclold(ijk) + gc( n,il,ijk ) * dxmic
2869 gcn( n,il,ijk ) = gc( n,il,ijk ) * rexpxctr( n )
2870 gcrgn( n,il,ijk ) = gcrg( n,il,ijk )
2872 gcn( -1,il,ijk ) = 0.0_rp
2873 gcn( 0,il,ijk ) = 0.0_rp
2874 gcn(
nbin+1,il,ijk ) = 0.0_rp
2875 gcn(
nbin+2,il,ijk ) = 0.0_rp
2876 gcrgn( -1,il,ijk ) = 0.0_rp
2877 gcrgn( 0,il,ijk ) = 0.0_rp
2878 gcrgn(
nbin+1,il,ijk ) = 0.0_rp
2879 gcrgn(
nbin+2,il,ijk ) = 0.0_rp
2882 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2884 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2885 cefd = emu/dens(ijk)
2888 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2889 f1 = rvap*temp(ijk)/festl/cefd
2890 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2891 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2893 umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq )
2894 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2895 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2896 dtcnd(ijk) = dtime / nloop(ijk)
2898 nloop(ijk) = nloop(ijk) * iflg( il,ijk )
2900 nloopmax = maxval(nloop,1)
2903 regene_gcn(:) = 0.0_rp
2906 do ncount = 1, nloopmax
2909 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
2913 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
2917 do nn = 1, loopflg(ijk)
2919 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2920 cefd = emu/dens(ijk)
2922 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
2923 f1 = rvap*temp(ijk)/festl/cefd
2924 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2925 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2927 sumliq(ijk) = 0.0_rp
2930 sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
2937 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
2938 temp(:), pres(:), qdry(:), &
2942 do nn = 1, loopflg(ijk)
2944 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
2947 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
2948 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
2949 cefliq = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)*qlevp(ijk)/cv(ijk)/rvap/temp(ijk)/temp(ijk) )
2950 a = - cefliq*sumliq(ijk)*gtliq(ijk)/dens(ijk)
2951 a = a + eps * ( 1.0_rp - zerosw )
2953 sliqtnd = zerosw * &
2954 ( ssliq * ( exp( a*dtcnd(ijk) )-1.0_rp )/( a*dtcnd(ijk) ) &
2955 * ( 0.5_rp + sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
2957 * ( 0.5_rp - sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
2959 + ssliq * ( 1.0_rp - zerosw )
2962 do mm = 1, iflg( il,ijk )
2965 uadv( n,il,ijk ) = cbnd( n,il )*rexpxbnd( n )*gtliq(ijk)*sliqtnd
2967 uadv( 0 ,il,ijk ) = 0.0_rp
2968 uadv(
nbin+2,il,ijk ) = 0.0_rp
2981 do nn = 1, loopflg(ijk)
2983 do mm = 1, iflg( myu,ijk )
2985 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
2993 do nn = 1, loopflg(ijk)
2995 do mm = 1, iflg( myu,ijk )
2997 acoef(0,n,myu,ijk) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
2998 acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
2999 acoef(2,n,myu,ijk) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3000 acoef_c(0,n,myu,ijk) = - ( gcrgn( n+1,myu,ijk )-26.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3001 acoef_c(1,n,myu,ijk) = ( gcrgn( n+1,myu,ijk ) -gcrgn( n-1,myu,ijk ) ) / 16.0_rp
3002 acoef_c(2,n,myu,ijk) = ( gcrgn( n+1,myu,ijk )- 2.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3010 do nn = 1, loopflg(ijk)
3012 do mm = 1, iflg( myu,ijk )
3014 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3016 aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3017 + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3018 + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3020 aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
3022 aip_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3023 + acoef_c(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3024 + acoef_c(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3026 aip_c(n,myu,ijk) = max( aip_c(n,myu,ijk), 0.0_rp )
3034 do nn = 1, loopflg(ijk)
3036 do mm = 1, iflg( myu,ijk )
3038 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3040 aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3041 - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3042 + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3044 aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
3046 aim_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3047 - acoef_c(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3048 + acoef_c(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3050 aim_c(n,myu,ijk) = max( aim_c(n,myu,ijk), 0.0_rp )
3058 do nn = 1, loopflg(ijk)
3060 do mm = 1, iflg( myu,ijk )
3062 ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
3063 + acoef(2,n,myu,ijk) * 2.0_rp
3065 ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
3067 ai_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * 2.0_rp &
3068 + acoef_c(2,n,myu,ijk) * 2.0_rp
3070 ai_c(n,myu,ijk) = max( ai_c(n,myu,ijk), aip_c(n,myu,ijk)+aim_c(n,myu,ijk)+cldmin )
3078 do nn = 1, loopflg(ijk)
3080 do mm = 1, iflg( myu,ijk )
3082 flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
3083 - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3084 flq_c(n,myu,ijk) = ( aip_c(n-1,myu,ijk)/ai_c(n-1,myu,ijk)*gcrgn( n-1,myu,ijk ) &
3085 - aim_c(n ,myu,ijk)/ai_c(n ,myu,ijk)*gcrgn( n ,myu,ijk ) )/dtcnd(ijk)
3093 do nn = 1, loopflg(ijk)
3095 do mm = 1, iflg( myu,ijk )
3096 regene_gcn(ijk) = regene_gcn(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
3097 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
3104 do nn = 1, loopflg(ijk)
3106 do mm = 1, iflg( myu,ijk )
3108 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
3109 gcrgn( n,myu,ijk ) = gcrgn( n,myu,ijk ) - ( flq_c(n+1,myu,ijk)-flq_c(n,myu,ijk) )*dtcnd(ijk)
3120 do nn = 1, loopflg(ijk)
3123 gclnew(ijk) = 0.0_rp
3126 gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )
3131 gclnew(ijk) = gclnew(ijk)*dxmic
3134 cndmss = gclnew(ijk) - gclold(ijk)
3135 dqv = cndmss/dens(ijk)
3136 qvap(ijk) = qvap(ijk) - dqv
3137 temp(ijk) = temp(ijk) + dqv*qlevp(ijk)/cv(ijk)
3141 gclold(ijk) = gclnew(ijk)
3153 gc( n,il,ijk ) = gcn( n,il,ijk )*expxctr( n )
3154 gcrg( n,il,ijk ) = gcrgn( n,il,ijk )
3159 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3164 if ( gc( n,il,ijk ) < 0.0_rp )
then
3165 cndmss = -gc( n,il,ijk )
3166 gc( n,il,ijk ) = 0.0_rp
3167 gcrg( n,il,ijk ) = 0.0_rp
3168 dqv = cndmss/dens(ijk)
3169 qvap(ijk) = qvap(ijk) + dqv
3170 temp(ijk) = temp(ijk) - dqv*qlevp(ijk)/cv(ijk)
3180 end subroutine liqphase
3183 subroutine icephase( &
3203 atmos_hydrometeor_lhs, &
3209 atmos_saturation_pres2qsat_ice
3212 integer,
intent(in) :: ijkmax
3213 real(
rp),
intent(in) :: dens(ijkmax)
3214 real(
rp),
intent(in) :: pres(ijkmax)
3215 real(
rp),
intent(in) :: qdry(ijkmax)
3216 real(
rp),
intent(inout) :: temp(ijkmax)
3217 real(
rp),
intent(inout) :: qvap(ijkmax)
3218 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3219 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
3220 real(
rp),
intent(inout) :: cp(ijkmax)
3221 real(
rp),
intent(inout) :: cv(ijkmax)
3222 real(
dp),
intent(in) :: dtime
3224 integer :: myu, n, ncount
3225 integer :: nloop(ijkmax)
3226 real(
rp) :: gciold(ijkmax), gcinew(ijkmax)
3227 real(
rp) :: dtcnd(ijkmax)
3229 real(
rp) :: gtice(ijkmax), umax(ijkmax)
3230 real(
rp) :: qlsbl(ijkmax)
3231 real(
rp) :: cefice, d, uval, sicetnd
3232 real(
rp) :: sumice(ijkmax)
3234 real(
rp) :: gcn( -1:
nbin+2,nspc,ijkmax )
3235 real(
rp),
parameter :: cflfct = 0.50_rp
3236 real(
rp) :: dumm_regene(ijkmax)
3237 integer :: iflg( nspc,ijkmax )
3238 real(
rp) :: csum( nspc,ijkmax )
3239 real(
rp) :: f1, f2, emu, cefd, cefk, festi
3240 real(
rp) :: qsati(ijkmax)
3241 real(
rp),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
3242 real(
rp),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
3243 integer :: ijk, mm, nn, loopflg(ijkmax)
3244 real(
rp) :: zerosw, qvtmp
3249 real(
rp) :: uadv ( 0:
nbin+2,nspc,ijkmax )
3250 real(
rp) :: flq ( 1:
nbin+1,nspc,ijkmax )
3251 real(
rp) :: acoef( 0:2,0:
nbin+1,nspc,ijkmax )
3252 real(
rp) :: crn ( 0:
nbin+2,nspc,ijkmax )
3253 real(
rp) :: aip ( 0:
nbin+1,nspc,ijkmax )
3254 real(
rp) :: aim ( 0:
nbin+1,nspc,ijkmax )
3255 real(
rp) :: ai ( 0:
nbin+1,nspc,ijkmax )
3256 real(
rp) :: cmins, cplus
3260 real(
rp) :: gcrgn( -1:
nbin+2,nspc,ijkmax )
3261 real(
rp) :: flq_c( 1:
nbin+1,nspc,ijkmax )
3262 real(
rp) :: acoef_c( 0:2,0:
nbin+1,nspc,ijkmax )
3263 real(
rp) :: aip_c ( 0:
nbin+1,nspc,ijkmax )
3264 real(
rp) :: aim_c ( 0:
nbin+1,nspc,ijkmax )
3265 real(
rp) :: ai_c ( 0:
nbin+1,nspc,ijkmax )
3271 csum( :,: ) = 0.0_rp
3276 csum( myu,ijk ) = csum( myu,ijk ) + gc( n,myu,ijk )*dxmic
3281 if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
3287 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3289 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3290 temp(:), pres(:), qdry(:), &
3300 gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
3307 gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
3308 gcrgn( n,myu,ijk ) = gcrg( n,myu,ijk )
3310 gcn( -1,myu,ijk ) = 0.0_rp
3311 gcn( 0,myu,ijk ) = 0.0_rp
3312 gcn(
nbin+1,myu,ijk ) = 0.0_rp
3313 gcn(
nbin+2,myu,ijk ) = 0.0_rp
3314 gcrgn( -1,myu,ijk ) = 0.0_rp
3315 gcrgn( 0,myu,ijk ) = 0.0_rp
3316 gcrgn(
nbin+1,myu,ijk ) = 0.0_rp
3317 gcrgn(
nbin+2,myu,ijk ) = 0.0_rp
3321 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3323 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3324 cefd = emu/dens(ijk)
3326 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3327 f1 = rvap*temp(ijk)/festi/cefd
3328 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3329 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3333 uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice )
3334 umax(ijk) = max( umax(ijk),uval )
3337 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
3338 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
3339 dtcnd(ijk) = dtime/nloop(ijk)
3341 nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) )
3344 nloopmax = maxval(nloop,1)
3348 do ncount = 1, nloopmax
3351 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
3355 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3359 do nn = 1, loopflg(ijk)
3360 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3361 cefd = emu/dens(ijk)
3363 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3364 f1 = rvap*temp(ijk)/festi/cefd
3365 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3366 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3368 sumice(ijk) = 0.0_rp
3371 sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
3379 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3380 temp(:), pres(:), qdry(:), &
3385 do nn = 1, loopflg(ijk)
3388 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3391 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
3392 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
3394 cefice = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)*qlsbl(ijk)/cv(ijk)/rvap/temp(ijk)/temp(ijk) )
3395 d = - cefice*sumice(ijk)*gtice(ijk)/dens(ijk)
3396 d = d + eps * ( 1.0_rp - zerosw )
3398 sicetnd = zerosw * &
3399 ( ssice * ( exp( d*dtcnd(ijk) )-1.0_rp )/( d*dtcnd(ijk) ) &
3400 * ( 0.5_rp + sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
3402 * ( 0.5_rp - sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
3404 + ssice * ( 1.0_rp - zerosw )
3408 do mm = 1, iflg( myu,ijk )
3411 uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd
3413 uadv( 0, myu,ijk ) = 0.0_rp
3414 uadv(
nbin+2,myu,ijk ) = 0.0_rp
3424 do nn = 1, loopflg(ijk)
3426 do mm = 1, iflg( myu,ijk )
3428 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
3436 do nn = 1, loopflg(ijk)
3438 do mm = 1, iflg( myu,ijk )
3440 acoef(0,n,myu,ijk) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3441 acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3442 acoef(2,n,myu,ijk) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3443 acoef_c(0,n,myu,ijk) = - ( gcrgn( n+1,myu,ijk )-26.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3444 acoef_c(1,n,myu,ijk) = ( gcrgn( n+1,myu,ijk ) -gcrgn( n-1,myu,ijk ) ) / 16.0_rp
3445 acoef_c(2,n,myu,ijk) = ( gcrgn( n+1,myu,ijk )- 2.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3453 do nn = 1, loopflg(ijk)
3455 do mm = 1, iflg( myu,ijk )
3457 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3459 aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3460 + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3461 + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3463 aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
3465 aip_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3466 + acoef_c(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3467 + acoef_c(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3469 aip_c(n,myu,ijk) = max( aip_c(n,myu,ijk), 0.0_rp )
3477 do nn = 1, loopflg(ijk)
3479 do mm = 1, iflg( myu,ijk )
3481 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3483 aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3484 - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3485 + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3487 aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
3489 aim_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3490 - acoef_c(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3491 + acoef_c(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3493 aim_c(n,myu,ijk) = max( aim_c(n,myu,ijk), 0.0_rp )
3501 do nn = 1, loopflg(ijk)
3503 do mm = 1, iflg( myu,ijk )
3505 ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
3506 + acoef(2,n,myu,ijk) * 2.0_rp
3508 ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
3510 ai_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * 2.0_rp &
3511 + acoef_c(2,n,myu,ijk) * 2.0_rp
3513 ai_c(n,myu,ijk) = max( ai_c(n,myu,ijk), aip_c(n,myu,ijk)+aim_c(n,myu,ijk)+cldmin )
3521 do nn = 1, loopflg(ijk)
3523 do mm = 1, iflg( myu,ijk )
3525 flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
3526 - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3528 flq_c(n,myu,ijk) = ( aip_c(n-1,myu,ijk)/ai_c(n-1,myu,ijk)*gcrgn( n-1,myu,ijk ) &
3529 - aim_c(n ,myu,ijk)/ai_c(n ,myu,ijk)*gcrgn( n ,myu,ijk ) )/dtcnd(ijk)
3537 do nn = 1, loopflg(ijk)
3539 do mm = 1, iflg( myu,ijk )
3540 dumm_regene(ijk) = dumm_regene(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
3541 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
3548 do nn = 1, loopflg(ijk)
3550 do mm = 1, iflg( myu,ijk )
3552 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
3553 gcrgn( n,myu,ijk ) = gcrgn( n,myu,ijk ) - ( flq_c(n+1,myu,ijk)-flq_c(n,myu,ijk) )*dtcnd(ijk)
3564 do nn = 1, loopflg(ijk)
3567 gcinew(ijk) = 0.0_rp
3570 gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
3575 sblmss = gcinew(ijk) - gciold(ijk)
3576 dqv = sblmss/dens(ijk)
3577 qvap(ijk) = qvap(ijk) - dqv
3578 temp(ijk) = temp(ijk) + dqv*qlsbl(ijk)/cv(ijk)
3582 gciold(ijk) = gcinew(ijk)
3595 gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
3596 gcrg( n,myu,ijk ) = gcrgn( n,myu,ijk )
3604 end subroutine icephase
3607 subroutine mixphase( &
3627 atmos_hydrometeor_lhv, &
3628 atmos_hydrometeor_lhs, &
3636 atmos_saturation_pres2qsat_liq, &
3637 atmos_saturation_pres2qsat_ice
3640 integer,
intent(in) :: ijkmax
3641 real(
rp),
intent(in) :: dens(ijkmax)
3642 real(
rp),
intent(in) :: pres(ijkmax)
3643 real(
rp),
intent(in) :: qdry(ijkmax)
3644 real(
rp),
intent(inout) :: temp(ijkmax)
3645 real(
rp),
intent(inout) :: qvap(ijkmax)
3646 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
3647 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
3648 real(
rp),
intent(inout) :: cp(ijkmax)
3649 real(
rp),
intent(inout) :: cv(ijkmax)
3650 real(
dp),
intent(in) :: dtime
3652 integer :: n, myu, mm, ncount
3653 integer :: nloop(ijkmax)
3654 real(
rp) :: gclold(ijkmax), gclnew(ijkmax)
3655 real(
rp) :: gciold(ijkmax), gcinew(ijkmax)
3656 real(
rp) :: cndmss, sblmss
3657 real(
rp) :: gtliq(ijkmax), gtice(ijkmax)
3658 real(
rp) :: umax(ijkmax), uval, dtcnd(ijkmax)
3659 real(
rp) :: qlevp(ijkmax), qlsbl(ijkmax)
3660 real(
rp) :: cef1, cef2, cef3, cef4, a, b, c, d
3661 real(
rp) :: rmdplus, rmdmins, ssplus, ssmins, tplus, tmins
3662 real(
rp) :: sliqtnd, sicetnd
3663 real(
rp) :: ssliq, ssice
3664 real(
rp) :: sumliq(ijkmax), sumice(ijkmax)
3665 real(
rp) :: gcn( -1:
nbin+2,nspc,ijkmax )
3666 real(
rp),
parameter :: cflfct = 0.50_rp
3667 real(
rp) :: dumm_regene(ijkmax)
3668 real(
rp) :: csum( nspc,ijkmax )
3669 integer :: iflg( nspc,ijkmax )
3670 real(
rp) :: f1, f2, emu, cefd, cefk, festl, festi
3671 real(
rp) :: qsatl(ijkmax), qsati(ijkmax)
3672 real(
rp),
parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
3673 real(
rp),
parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
3674 real(
rp) :: qvtmp, zerosw
3675 integer :: ijk, nn, loopflg(ijkmax)
3679 real(
rp) :: uadv ( 0:
nbin+2,nspc,ijkmax )
3680 real(
rp) :: flq ( 1:
nbin+1,nspc,ijkmax )
3681 real(
rp) :: acoef( 0:2,0:
nbin+1,nspc,ijkmax )
3682 real(
rp) :: crn ( 0:
nbin+2,nspc,ijkmax )
3683 real(
rp) :: aip ( 0:
nbin+1,nspc,ijkmax )
3684 real(
rp) :: aim ( 0:
nbin+1,nspc,ijkmax )
3685 real(
rp) :: ai ( 0:
nbin+1,nspc,ijkmax )
3686 real(
rp) :: cmins, cplus
3690 real(
rp) :: gcrgn( -1:
nbin+2,nspc,ijkmax )
3691 real(
rp) :: flq_c( 1:
nbin+1,nspc,ijkmax )
3692 real(
rp) :: acoef_c( 0:2,0:
nbin+1,nspc,ijkmax )
3693 real(
rp) :: aip_c ( 0:
nbin+1,nspc,ijkmax )
3694 real(
rp) :: aim_c ( 0:
nbin+1,nspc,ijkmax )
3695 real(
rp) :: ai_c ( 0:
nbin+1,nspc,ijkmax )
3699 dumm_regene( : ) = 0.0_rp
3701 csum( :,: ) = 0.0_rp
3706 csum( myu,ijk ) = csum( myu,ijk )+gc( n,myu,ijk )*dxmic
3711 if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
3717 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3719 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3721 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
3722 temp(:), pres(:), qdry(:), &
3724 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3725 temp(:), pres(:), qdry(:), &
3734 gclold(ijk) = gclold(ijk) + gc( n,il,ijk )*dxmic
3739 gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
3746 gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
3747 gcrgn( n,myu,ijk ) = gcrg( n,myu,ijk )
3749 gcn( -1,myu,ijk ) = 0.0_rp
3750 gcn( 0,myu,ijk ) = 0.0_rp
3751 gcn(
nbin+1,myu,ijk ) = 0.0_rp
3752 gcn(
nbin+2,myu,ijk ) = 0.0_rp
3753 gcrgn( -1,myu,ijk ) = 0.0_rp
3754 gcrgn( 0,myu,ijk ) = 0.0_rp
3755 gcrgn(
nbin+1,myu,ijk ) = 0.0_rp
3756 gcrgn(
nbin+2,myu,ijk ) = 0.0_rp
3760 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3761 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3763 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3764 cefd = emu/dens(ijk)
3767 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3768 f1 = rvap*temp(ijk)/festl/cefd
3769 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3770 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3772 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3773 f1 = rvap*temp(ijk)/festi/cefd
3774 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3775 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3778 umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq )
3780 uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice )
3781 umax(ijk) = max( umax(ijk),uval )
3784 dtcnd(ijk) = cflfct*dxmic/umax(ijk)
3785 nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
3786 dtcnd(ijk) = dtime/nloop(ijk)
3788 nloop(ijk) = nloop(ijk) * iflg( il,ijk )
3789 nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) )
3792 nloopmax = maxval(nloop,1)
3797 do ncount = 1, nloopmax
3800 loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) )
3805 do nn = 1, loopflg(ijk)
3808 call atmos_hydrometeor_lhv( ijkmax, 1, ijkmax, temp(:), qlevp(:) )
3810 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
3813 emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
3814 cefd = emu/dens(ijk)
3816 festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3817 f1 = rvap*temp(ijk)/festl/cefd
3818 f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
3819 gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
3820 festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/tem00 - 1.0_rp/temp(ijk) ) )
3821 f1 = rvap*temp(ijk)/festi/cefd
3822 f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
3823 gtice(ijk) = 4.0_rp*pi/( f1+f2 )
3825 sumliq(ijk) = 0.0_rp
3827 sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
3830 sumice(ijk) = 0.0_rp
3833 sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
3840 call atmos_saturation_pres2qsat_liq( ijkmax, 1, ijkmax, &
3841 temp(:), pres(:), qdry(:), &
3843 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
3844 temp(:), pres(:), qdry(:), &
3848 do nn = 1, loopflg(ijk)
3851 ssliq = qvap(ijk)/qsatl(ijk) - 1.0_rp
3852 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
3854 zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps )
3855 qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
3856 cef1 = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cv(ijk) )
3857 cef2 = ( ssliq+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cv(ijk) )
3858 cef3 = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cv(ijk) )
3859 cef4 = ( ssice+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cv(ijk) )
3861 a = - cef1*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3862 b = - cef2*sumice(ijk)*gtice(ijk)/dens(ijk)
3863 c = - cef3*sumliq(ijk)*gtliq(ijk)/dens(ijk)
3864 d = - cef4*sumice(ijk)*gtice(ijk)/dens(ijk)
3866 b = b + eps * ( 1.0_rp - zerosw )
3868 rmdplus = ( ( a+d ) + sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
3869 rmdmins = ( ( a+d ) - sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
3871 rmdplus = rmdplus + eps * ( 1.0_rp - zerosw )
3872 rmdmins = rmdmins + eps * ( 1.0_rp - zerosw )
3875 ssplus = ( ( rmdmins-a )*ssliq - b*ssice )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
3876 ssmins = ( ( a-rmdplus )*ssliq + b*ssice )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
3878 tplus = ( exp( rmdplus*dtcnd(ijk) )-1.0_rp )/( rmdplus*dtcnd(ijk) ) &
3879 * ( 0.5_rp + sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) ) &
3881 * ( 0.5_rp - sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) )
3882 tmins = ( exp( rmdmins*dtcnd(ijk) )-1.0_rp )/( rmdmins*dtcnd(ijk) ) &
3883 * ( 0.5_rp + sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) ) &
3885 * ( 0.5_rp - sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) )
3887 sliqtnd = ssliq * ( 1.0_rp - zerosw ) &
3889 ( b*tplus*ssplus + b*tmins*ssmins )
3890 sicetnd = ssice * ( 1.0_rp - zerosw ) &
3892 ( ( rmdplus-a )*tplus*ssplus &
3893 + ( rmdmins-a )*tmins*ssmins )
3897 do mm = 1, iflg( myu,ijk )
3901 uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtliq(ijk)*sliqtnd &
3902 * ( 0.5_rp - sign( 0.5_rp,real(myu)-1.5_rp) ) &
3903 + cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd &
3904 * ( 0.5_rp + sign( 0.5_rp,real(myu)-1.5_rp) )
3906 uadv( 0, myu,ijk ) = 0.0_rp
3907 uadv(
nbin+2,myu,ijk ) = 0.0_rp
3916 do nn = 1, loopflg(ijk)
3918 do mm = 1, iflg( myu,ijk )
3920 crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
3928 do nn = 1, loopflg(ijk)
3930 do mm = 1, iflg( myu,ijk )
3932 acoef(0,n,myu,ijk) = - ( gcn( n+1,myu,ijk )-26.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3933 acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3934 acoef(2,n,myu,ijk) = ( gcn( n+1,myu,ijk )- 2.0_rp*gcn( n,myu,ijk )+gcn( n-1,myu,ijk ) ) / 48.0_rp
3936 acoef_c(0,n,myu,ijk) = - ( gcrgn( n+1,myu,ijk )-26.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3937 acoef_c(1,n,myu,ijk) = ( gcrgn( n+1,myu,ijk ) -gcrgn( n-1,myu,ijk ) ) / 16.0_rp
3938 acoef_c(2,n,myu,ijk) = ( gcrgn( n+1,myu,ijk )- 2.0_rp*gcrgn( n,myu,ijk )+gcrgn( n-1,myu,ijk ) ) / 48.0_rp
3946 do nn = 1, loopflg(ijk)
3948 do mm = 1, iflg( myu,ijk )
3950 cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3952 aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3953 + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3954 + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3956 aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
3958 aip_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3959 + acoef_c(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3960 + acoef_c(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3962 aip_c(n,myu,ijk) = max( aip_c(n,myu,ijk), 0.0_rp )
3970 do nn = 1, loopflg(ijk)
3972 do mm = 1, iflg( myu,ijk )
3974 cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3976 aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3977 - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3978 + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3980 aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
3982 aim_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3983 - acoef_c(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3984 + acoef_c(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3986 aim_c(n,myu,ijk) = max( aim_c(n,myu,ijk), 0.0_rp )
3994 do nn = 1, loopflg(ijk)
3996 do mm = 1, iflg( myu,ijk )
3998 ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
3999 + acoef(2,n,myu,ijk) * 2.0_rp
4001 ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
4003 ai_c(n,myu,ijk) = acoef_c(0,n,myu,ijk) * 2.0_rp &
4004 + acoef_c(2,n,myu,ijk) * 2.0_rp
4006 ai_c(n,myu,ijk) = max( ai_c(n,myu,ijk), aip_c(n,myu,ijk)+aim_c(n,myu,ijk)+cldmin )
4014 do nn = 1, loopflg(ijk)
4016 do mm = 1, iflg( myu,ijk )
4018 flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
4019 - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
4021 flq_c(n,myu,ijk) = ( aip_c(n-1,myu,ijk)/ai_c(n-1,myu,ijk)*gcrgn( n-1,myu,ijk ) &
4022 - aim_c(n ,myu,ijk)/ai_c(n ,myu,ijk)*gcrgn( n ,myu,ijk ) )/dtcnd(ijk)
4030 do nn = 1, loopflg(ijk)
4032 do mm = 1, iflg( myu,ijk )
4033 dumm_regene(ijk) = dumm_regene(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
4034 * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk)+eps )
4041 do nn = 1, loopflg(ijk)
4043 do mm = 1, iflg( myu,ijk )
4045 gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
4046 gcrgn( n,myu,ijk ) = gcrgn( n,myu,ijk ) - ( flq_c(n+1,myu,ijk)-flq_c(n,myu,ijk) )*dtcnd(ijk)
4057 do nn = 1, loopflg(ijk)
4059 gclnew(ijk) = 0.0_rp
4061 gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )*dxmic
4064 gcinew(ijk) = 0.0_rp
4067 gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
4072 cndmss = gclnew(ijk) - gclold(ijk)
4073 sblmss = gcinew(ijk) - gciold(ijk)
4075 qvap(ijk) = qvap(ijk) - ( cndmss + sblmss ) / dens(ijk)
4076 temp(ijk) = temp(ijk) + ( cndmss*qlevp(ijk)+sblmss*qlsbl(ijk) ) / dens(ijk) / cv(ijk)
4080 gclold(ijk) = gclnew(ijk)
4081 gciold(ijk) = gcinew(ijk)
4092 gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
4093 gcrg( n,myu,ijk ) = gcrgn( n,myu,ijk )
4098 if ( .NOT. flg_regeneration )
then
4099 dumm_regene(:) = 0.0_rp
4105 end subroutine mixphase
4108 subroutine ice_nucleat( &
4124 atmos_saturation_pres2qsat_ice
4126 atmos_hydrometeor_lhs, &
4133 integer,
intent(in) :: ijkmax
4134 integer,
intent(in) :: num_cold
4135 integer,
intent(in) :: index_cold(ijkmax)
4136 real(
rp),
intent(in) :: dens(ijkmax)
4137 real(
rp),
intent(in) :: pres(ijkmax)
4138 real(
rp),
intent(in) :: qdry(ijkmax)
4139 real(
rp),
intent(inout) :: temp(ijkmax)
4140 real(
rp),
intent(inout) :: qvap(ijkmax)
4141 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4142 real(
rp),
intent(inout) :: cp(ijkmax)
4143 real(
rp),
intent(inout) :: cv(ijkmax)
4144 real(
dp),
intent(in) :: dtime
4147 real(
rp) :: numin, tdel, qdel
4149 real(
rp),
parameter :: acoef = -0.639_rp, bcoef = 12.96_rp
4151 real(
rp),
parameter :: tcolmu = 269.0_rp, tcolml = 265.0_rp
4152 real(
rp),
parameter :: tdendu = 257.0_rp, tdendl = 255.0_rp
4153 real(
rp),
parameter :: tplatu = 250.6_rp
4155 real(
rp) :: qsati(ijkmax), qlsbl(ijkmax)
4156 integer :: ijk, indirect
4163 call atmos_saturation_pres2qsat_ice( ijkmax, 1, ijkmax, &
4164 temp(:), pres(:), qdry(:), &
4167 call atmos_hydrometeor_lhs( ijkmax, 1, ijkmax, temp(:), qlsbl(:) )
4169 do indirect = 1, num_cold
4170 ijk = index_cold(indirect)
4173 ssice = qvap(ijk)/qsati(ijk) - 1.0_rp
4175 if( ssice <= 0.0_rp ) cycle
4180 numice = numice + gc( n,myu,ijk )*rexpxctr( n )*dxmic
4184 numin = n0_icenucl * exp( acoef + bcoef * ssice * 1.e+2_rp )
4185 if( numin > numice )
then
4187 if ( temp(ijk) <= tplatu .OR. ( temp(ijk) >= tcolml .AND. temp(ijk) < tcolmu ) )
then
4190 elseif( temp(ijk) <= tdendu .AND. temp(ijk) >= tdendl )
then
4197 numin = (numin-numice) * expxctr( 1 )
4198 numin = min( numin,qvap(ijk)*dens(ijk) )
4199 gc( 1,ispc,ijk ) = gc( 1,ispc,ijk ) + numin / dxmic
4201 tdel = numin/dens(ijk)*qlsbl(ijk)/cv(ijk)
4202 temp(ijk) = temp(ijk) + tdel
4203 qdel = numin/dens(ijk)
4204 qvap(ijk) = qvap(ijk) - qdel
4215 end subroutine ice_nucleat
4218 subroutine freezing( &
4236 atmos_hydrometeor_lhf, &
4243 integer,
intent(in) :: ijkmax
4244 integer,
intent(in) :: num_cold
4245 integer,
intent(in) :: index_cold(ijkmax)
4246 logical,
intent(in) :: flg_lt
4247 real(
rp),
intent(in) :: dens(ijkmax)
4248 real(
rp),
intent(inout) :: temp(ijkmax)
4249 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4250 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4251 real(
rp),
intent(inout) :: cp(ijkmax)
4252 real(
rp),
intent(inout) :: cv(ijkmax)
4253 real(
dp),
intent(in) :: dtime
4255 integer :: nbound, n
4256 real(
rp) :: xbound, tc, rate, dmp, frz, sumfrz, tdel
4257 real(
rp),
parameter :: coefa = 1.0e-01_rp
4258 real(
rp),
parameter :: coefb = 0.66_rp
4259 real(
rp),
parameter :: rbound = 2.0e-04_rp
4263 integer :: ijk, indirect
4264 real(
rp) :: qlmlt(ijkmax)
4269 call atmos_hydrometeor_lhf( ijkmax, 1, ijkmax, temp(:), qlmlt(:) )
4271 xbound = log( rhow * 4.0_rp*pi/3.0_rp * rbound**3 )
4272 nbound = int( ( xbound-xbnd( 1 ) )/dxmic ) + 1
4274 do indirect = 1, num_cold
4275 ijk = index_cold(indirect)
4280 rate = coefa*exp( -coefb*tc )
4285 dmp = rate*expxctr( n )
4286 frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
4287 frz = min( frz, gc( n,il,ijk ) )
4290 dcrg = frz / ( gc(n,il,ijk)+eps ) * gcrg( n,il,ijk )
4291 gcrg( n,il,ijk ) = gcrg( n,il,ijk ) - dcrg
4292 gcrg( n,ip,ijk ) = gcrg( n,ip,ijk ) + dcrg
4295 gc( n,il,ijk ) = gc( n,il,ijk ) - frz
4296 gc( n,ip,ijk ) = gc( n,ip,ijk ) + frz
4298 sumfrz = sumfrz + frz
4301 dmp = rate*expxctr( n )
4302 frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
4303 frz = min( frz, gc( n,il,ijk ) )
4306 dcrg = frz / ( gc(n,il,ijk)+eps ) * gcrg( n,il,ijk )
4307 gcrg( n,il,ijk ) = gcrg( n,il,ijk ) - dcrg
4308 gcrg( n,ih,ijk ) = gcrg( n,ih,ijk ) + dcrg
4311 gc( n,il,ijk ) = gc( n,il,ijk ) - frz
4312 gc( n,ih,ijk ) = gc( n,ih,ijk ) + frz
4314 sumfrz = sumfrz + frz
4351 sumfrz = sumfrz*dxmic
4353 tdel = sumfrz/dens(ijk)*qlmlt(ijk)/cv(ijk)
4354 temp(ijk) = temp(ijk) + tdel
4362 end subroutine freezing
4365 subroutine melting( &
4378 atmos_hydrometeor_lhf, &
4385 integer,
intent(in) :: ijkmax
4386 integer,
intent(in) :: num_warm
4387 integer,
intent(in) :: index_warm(ijkmax)
4388 logical,
intent(in) :: flg_lt
4389 real(
rp),
intent(in) :: dens(ijkmax)
4390 real(
rp),
intent(inout) :: temp(ijkmax)
4391 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4392 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4393 real(
rp),
intent(inout) :: cp(ijkmax)
4394 real(
rp),
intent(inout) :: cv(ijkmax)
4395 real(
dp),
intent(in) :: dtime
4398 real(
rp) :: summlt, sumice, tdel
4399 integer :: ijk, indirect
4400 real(
rp) :: qlmlt(ijkmax)
4405 call atmos_hydrometeor_lhf( ijkmax, 1, ijkmax, temp(:), qlmlt(:) )
4407 do indirect = 1, num_warm
4408 ijk = index_warm(indirect)
4414 sumice = sumice + gc( n,m,ijk )
4415 gc( n,m,ijk ) = 0.0_rp
4417 gc( n,il,ijk ) = gc( n,il,ijk ) + sumice
4421 dcrg = dcrg + gcrg( n,m,ijk )
4422 gcrg( n,m,ijk ) = 0.0_rp
4424 gcrg( n,il,ijk ) = gcrg( n,il,ijk ) + dcrg
4426 summlt = summlt + sumice
4428 summlt = summlt*dxmic
4430 tdel = - summlt/dens(ijk)*qlmlt(ijk)/cv(ijk)
4431 temp(ijk) = temp(ijk) + tdel
4440 end subroutine melting
4443 subroutine collmain( &
4458 integer,
intent(in) :: ka
4459 integer,
intent(in) :: ia
4460 integer,
intent(in) :: ja
4462 integer,
intent(in) :: ijkmax
4463 logical,
intent(in) :: flg_lt
4464 real(
rp),
intent(in) :: v0_crg, d0_crg
4465 real(
rp),
intent(in) :: dq(ijkmax)
4466 real(
rp),
intent(in) :: beta(ijkmax)
4467 real(
rp),
intent(in) :: temp(ijkmax)
4468 real(
rp),
intent(inout) :: ghyd(
nbin,nspc,ijkmax)
4469 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4470 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
4471 real(
dp),
intent(in) :: dt
4474 if ( flg_rndm )
then
4475 call r_collcoag( ka, ia, ja, &
4482 call collcoag( ijkmax, &
4496 end subroutine collmain
4499 subroutine collmainf( &
4514 integer,
intent(in) :: ka
4515 integer,
intent(in) :: ia
4516 integer,
intent(in) :: ja
4518 integer,
intent(in) :: ijkmax
4519 logical,
intent(in) :: flg_lt
4520 real(
rp),
intent(in) :: v0_crg, d0_crg
4521 real(
rp),
intent(in) :: dq(ijkmax)
4522 real(
rp),
intent(in) :: beta(ijkmax)
4523 real(
rp),
intent(in) :: temp(ijkmax)
4524 real(
rp),
intent(inout) :: ghyd(
nbin,nspc,ijkmax)
4525 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4526 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
4527 real(
dp),
intent(in) :: dt
4530 if ( flg_rndm )
then
4531 call r_collcoag( ka, ia, ja, &
4538 call collcoag( ijkmax, &
4552 end subroutine collmainf
4557 subroutine collcoag( &
4571 integer,
intent(in) :: ijkmax
4572 logical,
intent(in) :: flg_lt
4573 real(
rp),
intent(in) :: d0_crg, v0_crg
4574 real(
rp),
intent(in) :: dq(ijkmax)
4575 real(
rp),
intent(in) :: beta(ijkmax)
4576 real(
rp),
intent(in) :: temp(ijkmax)
4577 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
4578 real(
rp),
intent(inout) :: gcrg(
nbin,nspc,ijkmax)
4579 real(
rp),
intent(out) :: crg_sep(nspc,ijkmax)
4580 real(
dp),
intent(in) :: dtime
4582 integer :: i, j, k, l
4583 real(
rp) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
4584 real(
rp) :: gprime, gprimk, wgt, crn, sum, flux
4585 integer,
parameter :: ldeg = 2
4586 real(
rp) :: acoef( 0:ldeg )
4587 real(
rp),
parameter :: dmpmin = 1.e-01_rp
4588 real(
rp) :: suri, surj
4590 integer :: myu, n, irsl, ilrg, isml
4591 integer :: ibnd( ijkmax )
4592 integer :: iflg( nspc,ijkmax )
4593 integer :: iexst(
nbin,nspc,ijkmax )
4594 real(
rp) :: csum( nspc,ijkmax )
4595 integer :: ijk, nn, mm, pp, qq
4598 real(
rp) :: drhoi, drhoj, drhok, alpha, dgenei, dgenej
4606 csum( :,: ) = 0.0_rp
4607 crg_sep( :,: ) = 0.0_rp
4612 csum( myu,ijk ) = csum( myu,ijk ) + gc( n,myu,ijk )*dxmic
4616 if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
4619 if ( temp(ijk) < tcrit )
then
4627 if ( gc( n,myu,ijk ) > cldmin )
then
4628 iexst( n,myu,ijk ) = 1
4638 do nn = 1, iflg( isml,ijk )
4641 do mm = 1, iflg( ilrg,ijk )
4643 irsl = ifrsl( ibnd(ijk),isml,ilrg )
4645 if ( isml /= ilrg )
then
4647 else if ( isml == ilrg )
then
4652 do pp = 1, iexst( i,isml,ijk )
4655 do qq = 1, iexst( j,ilrg,ijk )
4662 dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
4663 dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
4665 if ( dmpi <= dmpmin )
then
4666 frci = gc( i,isml,ijk )*dmpi
4668 frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
4671 if ( dmpj <= dmpmin )
then
4672 frcj = gc( j,ilrg,ijk )*dmpj
4674 frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
4679 if ( gprime > 0.0_rp .AND. k <
nbin )
then
4681 suri = gc( i,isml,ijk )
4682 surj = gc( j,ilrg,ijk )
4683 gc( i,isml,ijk ) = gc( i,isml,ijk )-frci
4684 gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj
4685 gc( i,isml,ijk ) = max( gc( i,isml,ijk )-frci, 0.0_rp )
4686 gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk )-frcj, 0.0_rp )
4687 frci = suri - gc( i,isml,ijk )
4688 frcj = surj - gc( j,ilrg,ijk )
4693 drhoi = frci/suri*gcrg( i,isml,ijk )
4694 drhoj = frcj/surj*gcrg( j,ilrg,ijk )
4695 drhok = drhoi + drhoj
4698 alpha = 5.0_rp * ( 2.0_rp*radc( i )/d0_crg )**2 &
4699 * abs( vt(ilrg,j)-vt(isml,i) )/v0_crg
4700 alpha = min( 10.0_rp, alpha )
4702 dgenei = frci / xi * rcoll( isml,ilrg,i,j ) * ( -dq( ijk ) ) &
4703 * alpha * beta( ijk ) * flg_noninduct( isml,ilrg )
4705 gcrg( i,isml,ijk ) = gcrg( i,isml,ijk ) + ( dgenei-drhoi )
4706 gcrg( j,ilrg,ijk ) = gcrg( j,ilrg,ijk ) + ( dgenej-drhoj )
4707 crg_sep( isml,ijk ) = crg_sep( isml,ijk ) + dgenei
4708 crg_sep( ilrg,ijk ) = crg_sep( ilrg,ijk ) + dgenej
4711 gprimk = gc( k,irsl,ijk ) + gprime
4712 wgt = gprime / gprimk
4713 crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
4715 acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
4716 acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.50_rp
4717 acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
4721 sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
4722 *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
4726 flux = min( max( flux,0.0_rp ),gprime )
4728 gc( k,irsl,ijk ) = gprimk - flux
4729 gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
4733 if( gprime /= 0.0_rp )
then
4734 gcrg( k,irsl,ijk ) = gcrg( k,irsl,ijk ) + drhok * ( gprime-flux )/gprime
4735 gcrg( k+1,irsl,ijk ) = gcrg( k+1,irsl,ijk ) + drhok * flux/gprime
4757 end subroutine collcoag
4760 subroutine getrule( ifrsl, indx )
4762 integer,
intent(out) :: indx(
nbin,
nbin )
4763 integer,
intent(out) :: ifrsl( 2,nspc_mk,nspc_mk )
4769 xnew = log( expxctr( i )+expxctr( j ) )
4770 k = int( ( xnew-xctr( 1 ) )/dxmic ) + 1
4771 k = max( max( k,j ),i )
4778 ifrsl( 1:2,il,il ) = il
4781 ifrsl( 1:2,il,ic ) = ic
4782 ifrsl( 1,ic,il ) = ig
4783 ifrsl( 2,ic,il ) = ih
4786 ifrsl( 1:2,il,ip ) = ip
4787 ifrsl( 1,ip,il ) = ig
4788 ifrsl( 2,ip,il ) = ih
4791 ifrsl( 1:2,il,id ) = id
4792 ifrsl( 1,id,il ) = ig
4793 ifrsl( 2,id,il ) = ih
4796 ifrsl( 1:2,il,iss ) = iss
4797 ifrsl( 1,iss,il ) = ig
4798 ifrsl( 2,iss,il ) = ih
4801 ifrsl( 1:2,il,ig ) = ig
4802 ifrsl( 1,ig,il ) = ig
4803 ifrsl( 2,ig,il ) = ih
4806 ifrsl( 1:2,il,ih ) = ih
4807 ifrsl( 1,ih,il ) = ig
4808 ifrsl( 2,ih,il ) = ih
4812 ifrsl( 1:2,ic,ic ) = iss
4815 ifrsl( 1:2,ic,ip ) = iss
4816 ifrsl( 1:2,ip,ic ) = iss
4819 ifrsl( 1:2,ic,id ) = iss
4820 ifrsl( 1:2,id,ic ) = iss
4823 ifrsl( 1:2,ic,iss ) = iss
4824 ifrsl( 1:2,iss,ic ) = iss
4827 ifrsl( 1:2,ic,ig ) = ig
4828 ifrsl( 1:2,ig,ic ) = ic
4831 ifrsl( 1:2,ih,ic ) = ic
4832 ifrsl( 1,ic,ih ) = ig
4833 ifrsl( 2,ic,ih ) = ih
4837 ifrsl( 1:2,ip,ip ) = iss
4840 ifrsl( 1:2,ip,id ) = iss
4841 ifrsl( 1:2,id,ip ) = iss
4844 ifrsl( 1:2,ip,iss ) = iss
4845 ifrsl( 1:2,iss,ip ) = iss
4848 ifrsl( 1:2,ip,ig ) = ig
4849 ifrsl( 1:2,ig,ip ) = ip
4852 ifrsl( 1:2,ih,ip ) = ip
4853 ifrsl( 1,ip,ih ) = ig
4854 ifrsl( 2,ip,ih ) = ih
4858 ifrsl( 1:2,id,id ) = iss
4861 ifrsl( 1:2,id,iss ) = iss
4862 ifrsl( 1:2,iss,id ) = iss
4865 ifrsl( 1:2,id,ig ) = ig
4866 ifrsl( 1:2,ig,id ) = id
4869 ifrsl( 1:2,ih,id ) = id
4870 ifrsl( 1,id,ih ) = ig
4871 ifrsl( 2,id,ih ) = ih
4875 ifrsl( 1:2,iss,iss ) = iss
4878 ifrsl( 1:2,iss,ig ) = ig
4879 ifrsl( 1:2,ig,iss ) = iss
4882 ifrsl( 1:2,ih,iss ) = iss
4883 ifrsl( 1,iss,ih ) = ig
4884 ifrsl( 2,iss,ih ) = ih
4888 ifrsl( 1:2,ig,ig ) = ig
4891 ifrsl( 1,ig,ih ) = ig
4892 ifrsl( 1,ih,ig ) = ig
4893 ifrsl( 2,ig,ih ) = ih
4894 ifrsl( 2,ih,ig ) = ih
4897 ifrsl( 1:2,ih,ih ) = ih
4900 end subroutine getrule
4909 integer ,
intent(in) :: ijkmax
4910 real(
rp),
intent(in) :: f0(ijkmax)
4911 real(
rp),
intent(inout) :: ga(
nccn,ijkmax)
4922 ga(n,ijk) = ga(n,ijk) + f0(ijk) * marate(n) * expxactr(n) / dxaer
4929 end subroutine faero
4934 subroutine random_setup( mset )
4940 integer,
intent(in) :: mset
4944 real(
rp) :: nbinr, tmp1
4945 real(
rp) :: rans( mbin ), ranl( mbin )
4947 real(
rp),
allocatable :: ranstmp( : )
4948 real(
rp),
allocatable :: ranltmp( : )
4951 integer,
allocatable :: orderb( : )
4954 real(
rp),
allocatable :: randnum(:,:,:)
4957 allocate( blrg( mset, mbin ) )
4958 allocate( bsml( mset, mbin ) )
4959 allocate( ranstmp( pq ) )
4960 allocate( ranltmp( pq ) )
4961 allocate( orderb( pq ) )
4962 allocate( randnum(1,1,pq) )
4964 a = real(
nbin )*real(
nbin-1 )*0.50_rp
4965 if ( a < mbin )
then
4966 log_error(
"ATMOS_PHY_MP_SUZUKI10_random_setup",*)
"mbin should be smaller than {nbin}_C_{2}"
4970 wgtbin = a/real( mbin )
4971 nbinr = real(
nbin )
4978 ranstmp( (p-1)*
nbin-(p*(p-1))/2+1 : p*
nbin-(p*(p+1))/2 ) = p
4980 ranltmp( (p-1)*
nbin-(p*(p-1))/2+q ) = p+q
4985 call random_uniform( randnum )
4987 abq1 = randnum( 1,1,p )
4988 k = int( abq1*( pq-p-1 ) ) + p
4990 orderb( p ) = orderb( k )
4996 rans( p ) = ranstmp( orderb( p ) )
4997 ranl( p ) = ranltmp( orderb( p ) )
4999 rans( p ) = ranstmp( orderb( p-pq ) )
5000 ranl( p ) = ranltmp( orderb( p-pq ) )
5002 if ( rans( p ) >= ranl( p ) )
then
5004 rans( p ) = ranl( p )
5008 blrg( n,1:mbin ) = int( ranl( 1:mbin ) )
5009 bsml( n,1:mbin ) = int( rans( 1:mbin ) )
5012 deallocate( ranstmp )
5013 deallocate( ranltmp )
5014 deallocate( orderb )
5015 deallocate( randnum )
5017 end subroutine random_setup
5023 subroutine r_collcoag( &
5034 integer,
intent(in) :: ka
5035 integer,
intent(in) :: ia
5036 integer,
intent(in) :: ja
5038 integer,
intent(in) :: ijkmax
5039 real(
rp),
intent(in) :: swgt
5040 real(
rp),
intent(in) :: temp(ijkmax)
5041 real(
rp),
intent(inout) :: gc (
nbin,nspc,ijkmax)
5042 real(
dp),
intent(in) :: dtime
5044 integer :: i, j, k, l
5045 real(
rp) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
5046 real(
rp) :: gprime, gprimk, wgt, crn, sum, flux
5047 integer,
parameter :: ldeg = 2
5048 real(
rp),
parameter :: dmpmin = 1.e-01_rp, cmin = 1.e-10_rp
5049 real(
rp) :: acoef( 0:ldeg )
5052 integer :: nums( mbin ), numl( mbin )
5053 real(
rp),
parameter :: gt = 1.0_rp
5055 real(
rp) :: nbinr, mbinr
5057 real(
rp) :: tmpi, tmpj
5059 integer :: ibnd( ijkmax )
5060 integer :: iflg( nspc,ijkmax )
5061 integer :: iexst(
nbin,nspc,ijkmax )
5062 real(
rp) :: csum( nspc,ijkmax )
5063 integer :: ijk, nn, mm, pp, qq, myu, n, isml, ilrg, irsl
5070 csum( :,: ) = 0.0_rp
5074 csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
5075 csum( ic,ijk ) = csum( ic,ijk ) + gc( n,ic,ijk )*dxmic
5076 csum( ip,ijk ) = csum( ip,ijk ) + gc( n,ip,ijk )*dxmic
5077 csum( id,ijk ) = csum( id,ijk ) + gc( n,id,ijk )*dxmic
5078 csum( iss,ijk ) = csum( iss,ijk ) + gc( n,iss,ijk )*dxmic
5079 csum( ig,ijk ) = csum( ig,ijk ) + gc( n,ig,ijk )*dxmic
5080 csum( ih,ijk ) = csum( ih,ijk ) + gc( n,ih,ijk )*dxmic
5082 if ( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
5083 if ( csum( ic,ijk ) > cldmin ) iflg( ic,ijk ) = 1
5084 if ( csum( ip,ijk ) > cldmin ) iflg( ip,ijk ) = 1
5085 if ( csum( id,ijk ) > cldmin ) iflg( id,ijk ) = 1
5086 if ( csum( iss,ijk ) > cldmin ) iflg( iss,ijk ) = 1
5087 if ( csum( ig,ijk ) > cldmin ) iflg( ig,ijk ) = 1
5088 if ( csum( ih,ijk ) > cldmin ) iflg( ih,ijk ) = 1
5090 if ( temp(ijk) < tcrit )
then
5098 if ( gc( n,myu,ijk ) > cldmin )
then
5099 iexst( n,myu,ijk ) = 1
5109 do nn = 1, iflg( isml,ijk )
5112 do mm = 1, iflg( ilrg,ijk )
5114 irsl = ifrsl( ibnd(ijk),isml,ilrg )
5116 call random_uniform( rndm )
5117 det = int( rndm(1,1,1)*ia*ja*ka )
5118 nbinr = real(
nbin )
5119 mbinr = real( mbin )
5120 nums( 1:mbin ) = bsml( det,1:mbin )
5121 numl( 1:mbin ) = blrg( det,1:mbin )
5127 do pp = 1, iexst( i,isml,ijk )
5128 do qq = 1, iexst( j,ilrg,ijk )
5135 dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
5136 dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
5138 if ( dmpi <= dmpmin )
then
5139 frci = gc( i,isml,ijk )*dmpi
5141 frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
5144 if ( dmpj <= dmpmin )
then
5145 frcj = gc( j,ilrg,ijk )*dmpj
5147 frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
5149 tmpi = gc( i,isml,ijk )
5150 tmpj = gc( j,ilrg,ijk )
5152 gc( i,isml,ijk ) = gc( i,isml,ijk )-frci*swgt
5153 gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj*swgt
5156 gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk ), 0.0_rp )
5158 gc( i,isml,ijk ) = max( gc( i,isml,ijk ), 0.0_rp )
5160 frci = tmpi - gc( i,isml,ijk )
5161 frcj = tmpj - gc( j,ilrg,ijk )
5186 if ( gprime > 0.0_rp .AND. k <
nbin )
then
5188 gprimk = gc( k,irsl,ijk ) + gprime
5189 wgt = gprime / gprimk
5190 crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
5192 acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
5193 acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.5_rp
5194 acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
5198 sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
5199 *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
5203 flux = min( max( flux,0.0_rp ),gprime )
5205 gc( k,irsl,ijk ) = gprimk - flux
5206 gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
5225 end subroutine r_collcoag
5238 allocate( radc_mk(
nbin ) )
5239 allocate( xctr_mk(
nbin ) )
5240 allocate( xbnd_mk(
nbin+1 ) )
5241 allocate( cctr_mk( nspc_mk,
nbin ) )
5242 allocate( cbnd_mk( nspc_mk,
nbin+1 ) )
5243 allocate( ck_mk( nspc_mk,nspc_mk,
nbin,
nbin ) )
5244 allocate( vt_mk( nspc_mk,
nbin ) )
5245 allocate( br_mk( nspc_mk,
nbin ) )
5268 deallocate( radc_mk )
5269 deallocate( xctr_mk )
5270 deallocate( xbnd_mk )
5271 deallocate( cctr_mk )
5272 deallocate( cbnd_mk )
5277 end subroutine mkpara
5283 integer,
parameter :: il = 1, ic = 2, ip = 3, id = 4
5284 integer,
parameter :: is = 5, ig = 6, ih = 7
5285 integer,
parameter :: icemax = 3
5287 integer :: k, kk, i, j
5288 real(
dp) :: xl( ndat ), rlec( ndat ), vrl( ndat )
5289 real(
dp) :: blkradl( ndat ), blkdnsl( ndat )
5291 real(
dp) :: xi( ndat,icemax ), riec( ndat,icemax ), vri( ndat,icemax )
5292 real(
dp) :: blkradi( ndat,icemax ), blkdnsi( ndat,icemax )
5294 real(
dp) :: xs( ndat ), rsec( ndat ), vrs( ndat )
5295 real(
dp) :: blkrads( ndat ), blkdnss( ndat )
5297 real(
dp) :: xg( ndat ), rgec( ndat ), vrg( ndat )
5298 real(
dp) :: blkradg( ndat ), blkdnsg( ndat )
5300 real(
dp) :: xh( ndat ), rhec( ndat ), vrh( ndat )
5301 real(
dp) :: blkradh( ndat ), blkdnsh( ndat )
5304 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5305 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5306 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5307 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5308 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5309 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5310 data xi(1:ndat,1) / &
5311 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5312 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5313 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5314 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5315 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5316 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5317 data xi(1:ndat,2) / &
5318 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5319 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5320 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5321 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5322 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5323 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00 /
5324 data xi(1:ndat,3) / &
5325 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5326 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5327 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5328 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5329 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5330 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5332 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5333 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5334 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5335 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5336 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5337 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5339 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5340 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5341 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5342 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5343 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5344 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5346 0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
5347 0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
5348 0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
5349 0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
5350 0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
5351 0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
5352 data rlec(1:ndat) / &
5353 0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
5354 0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
5355 0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
5356 0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
5357 0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
5358 0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
5359 data riec(1:ndat,1) / &
5360 0.31936e-03_dp,0.40397e-03_dp,0.51099e-03_dp,0.64638e-03_dp,0.81764e-03_dp,0.10343e-02_dp, &
5361 0.13084e-02_dp,0.16551e-02_dp,0.20937e-02_dp,0.26486e-02_dp,0.33506e-02_dp,0.42387e-02_dp, &
5362 0.64360e-02_dp,0.81426e-02_dp,0.10302e-01_dp,0.13035e-01_dp,0.16494e-01_dp,0.20872e-01_dp, &
5363 0.26412e-01_dp,0.33426e-01_dp,0.42304e-01_dp,0.53543e-01_dp,0.67770e-01_dp,0.85783e-01_dp, &
5364 0.10859e+00_dp,0.13746e+00_dp,0.17403e+00_dp,0.22032e+00_dp,0.27895e+00_dp,0.35319e+00_dp, &
5365 0.44722e+00_dp,0.56630e+00_dp,0.71712e+00_dp /
5366 data riec(1:ndat,2) / &
5367 0.13188e-03_dp,0.16615e-03_dp,0.20953e-03_dp,0.27728e-03_dp,0.36694e-03_dp,0.48559e-03_dp, &
5368 0.64261e-03_dp,0.85040e-03_dp,0.11254e-02_dp,0.14893e-02_dp,0.19709e-02_dp,0.26082e-02_dp, &
5369 0.34515e-02_dp,0.45676e-02_dp,0.60446e-02_dp,0.79991e-02_dp,0.10586e-01_dp,0.14009e-01_dp, &
5370 0.18539e-01_dp,0.24533e-01_dp,0.32466e-01_dp,0.42964e-01_dp,0.56857e-01_dp,0.75242e-01_dp, &
5371 0.99573e-01_dp,0.13177e+00_dp,0.17438e+00_dp,0.23077e+00_dp,0.30539e+00_dp,0.40414e+00_dp, &
5372 0.53482e+00_dp,0.70775e+00_dp,0.93661e+00_dp /
5373 data riec(1:ndat,3) / 0.14998e-03_dp,0.18896e-03_dp,0.23808e-03_dp, &
5374 0.29996e-03_dp,0.37793e-03_dp,0.47616e-03_dp,0.61048e-03_dp,0.81343e-03_dp,0.10839e-02_dp, &
5375 0.14442e-02_dp,0.19243e-02_dp,0.25640e-02_dp,0.34164e-02_dp,0.45522e-02_dp,0.60656e-02_dp, &
5376 0.80820e-02_dp,0.10769e-01_dp,0.14349e-01_dp,0.19119e-01_dp,0.25475e-01_dp,0.44576e-01_dp, &
5377 0.62633e-01_dp,0.88006e-01_dp,0.12366e+00_dp,0.17375e+00_dp,0.24414e+00_dp,0.34304e+00_dp, &
5378 0.48201e+00_dp,0.67728e+00_dp,0.95164e+00_dp,0.13372e+01_dp,0.18788e+01_dp,0.26400e+01_dp /
5379 data rsec(1:ndat) / &
5380 0.92832e-03_dp,0.11696e-02_dp,0.14736e-02_dp,0.18566e-02_dp,0.23392e-02_dp,0.29472e-02_dp, &
5381 0.37133e-02_dp,0.46784e-02_dp,0.58944e-02_dp,0.74265e-02_dp,0.93569e-02_dp,0.11789e-01_dp, &
5382 0.14853e-01_dp,0.18714e-01_dp,0.23578e-01_dp,0.29706e-01_dp,0.37427e-01_dp,0.47156e-01_dp, &
5383 0.59412e-01_dp,0.74855e-01_dp,0.94311e-01_dp,0.11882e+00_dp,0.14971e+00_dp,0.18862e+00_dp, &
5384 0.23765e+00_dp,0.29942e+00_dp,0.37724e+00_dp,0.47530e+00_dp,0.59884e+00_dp,0.75449e+00_dp, &
5385 0.95060e+00_dp,0.11977e+01_dp,0.15090e+01_dp /
5386 data rgec(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
5387 0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
5388 0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
5389 0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
5390 0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
5391 0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
5392 data rhec(1:ndat) / &
5393 0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
5394 0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
5395 0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
5396 0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
5397 0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
5398 0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
5399 data vrl(1:ndat) / &
5400 0.50000e-01_dp,0.78000e-01_dp,0.12000e+00_dp,0.19000e+00_dp,0.31000e+00_dp,0.49000e+00_dp, &
5401 0.77000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp,0.48000e+01_dp,0.74000e+01_dp, &
5402 0.11000e+02_dp,0.17000e+02_dp,0.26000e+02_dp,0.37000e+02_dp,0.52000e+02_dp,0.71000e+02_dp, &
5403 0.94000e+02_dp,0.12000e+03_dp,0.16000e+03_dp,0.21000e+03_dp,0.26000e+03_dp,0.33000e+03_dp, &
5404 0.41000e+03_dp,0.48000e+03_dp,0.57000e+03_dp,0.66000e+03_dp,0.75000e+03_dp,0.82000e+03_dp, &
5405 0.88000e+03_dp,0.90000e+03_dp,0.90000e+03_dp /
5406 data vri(1:ndat,1) / 0.30000e-01_dp,0.40000e-01_dp,0.60000e-01_dp, &
5407 0.80000e-01_dp,0.11000e+00_dp,0.15000e+00_dp,0.17000e+00_dp,0.18000e+00_dp,0.20000e+00_dp, &
5408 0.25000e+00_dp,0.40000e+00_dp,0.60000e+01_dp,0.10000e+02_dp,0.15000e+02_dp,0.20000e+02_dp, &
5409 0.25000e+02_dp,0.31000e+02_dp,0.37000e+02_dp,0.41000e+02_dp,0.46000e+02_dp,0.51000e+02_dp, &
5410 0.55000e+02_dp,0.59000e+02_dp,0.62000e+02_dp,0.64000e+02_dp,0.67000e+02_dp,0.68000e+02_dp, &
5411 0.69000e+02_dp,0.70000e+02_dp,0.71000e+02_dp,0.71500e+02_dp,0.71750e+02_dp,0.72000e+02_dp /
5412 data vri(1:ndat,2) / &
5413 0.30000e-01_dp,0.40000e-01_dp,0.50000e-01_dp,0.70000e-01_dp,0.90000e-01_dp,0.12000e+00_dp, &
5414 0.50000e+00_dp,0.80000e+00_dp,0.16000e+01_dp,0.18000e+01_dp,0.20000e+01_dp,0.30000e+01_dp, &
5415 0.40000e+01_dp,0.50000e+01_dp,0.80000e+01_dp,0.13000e+02_dp,0.19000e+02_dp,0.26000e+02_dp, &
5416 0.32000e+02_dp,0.38000e+02_dp,0.47000e+02_dp,0.55000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
5417 0.77000e+02_dp,0.79000e+02_dp,0.80000e+02_dp,0.81000e+02_dp,0.81000e+02_dp,0.82000e+02_dp, &
5418 0.82000e+02_dp,0.82000e+02_dp,0.82000e+02_dp /
5419 data vri(1:ndat,3) / 0.35000e-01_dp,0.45000e-01_dp,0.55000e-01_dp, &
5420 0.75000e-01_dp,0.95000e-01_dp,0.13000e+00_dp,0.60000e+00_dp,0.90000e+00_dp,0.17000e+01_dp, &
5421 0.20000e+01_dp,0.25000e+01_dp,0.38000e+01_dp,0.50000e+01_dp,0.70000e+01_dp,0.90000e+01_dp, &
5422 0.11000e+02_dp,0.14000e+02_dp,0.17000e+02_dp,0.21000e+02_dp,0.25000e+02_dp,0.32000e+02_dp, &
5423 0.38000e+02_dp,0.44000e+02_dp,0.49000e+02_dp,0.53000e+02_dp,0.55000e+02_dp,0.58000e+02_dp, &
5424 0.59000e+02_dp,0.61000e+02_dp,0.62000e+02_dp,0.63000e+02_dp,0.64000e+02_dp,0.65000e+02_dp /
5425 data vrs(1:ndat) / &
5426 0.20000e-01_dp,0.31000e-01_dp,0.49000e-01_dp,0.77000e-01_dp,0.12000e+00_dp,0.19000e+00_dp, &
5427 0.30000e+00_dp,0.48000e+00_dp,0.76000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp, &
5428 0.48000e+01_dp,0.75000e+01_dp,0.11000e+02_dp,0.16000e+02_dp,0.21000e+02_dp,0.26000e+02_dp, &
5429 0.34000e+02_dp,0.41000e+02_dp,0.49000e+02_dp,0.57000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
5430 0.81000e+02_dp,0.87000e+02_dp,0.93000e+02_dp,0.99000e+02_dp,0.10750e+03_dp,0.11500e+03_dp, &
5431 0.12500e+03_dp,0.13500e+03_dp,0.14500e+03_dp /
5432 data vrg(1:ndat) / 0.39000e-01_dp,0.62000e-01_dp,0.97000e-01_dp, &
5433 0.15000e+00_dp,0.24000e+00_dp,0.38000e+00_dp,0.61000e+00_dp,0.96000e+00_dp,0.15000e+01_dp, &
5434 0.24000e+01_dp,0.38000e+01_dp,0.61000e+01_dp,0.96000e+01_dp,0.15000e+02_dp,0.23000e+02_dp, &
5435 0.31000e+02_dp,0.39000e+02_dp,0.49000e+02_dp,0.59000e+02_dp,0.68000e+02_dp,0.79000e+02_dp, &
5436 0.88000e+02_dp,0.10000e+03_dp,0.11000e+03_dp,0.13000e+03_dp,0.15000e+03_dp,0.17000e+03_dp, &
5437 0.20000e+03_dp,0.23000e+03_dp,0.26000e+03_dp,0.30000e+03_dp,0.35000e+03_dp,0.40000e+03_dp /
5438 data vrh(1:ndat) / &
5439 0.53000e-01_dp,0.84000e-01_dp,0.13000e+00_dp,0.21000e+00_dp,0.33000e+00_dp,0.52000e+00_dp, &
5440 0.82000e+00_dp,0.13000e+01_dp,0.21000e+01_dp,0.33000e+01_dp,0.52000e+01_dp,0.82000e+01_dp, &
5441 0.13000e+02_dp,0.20000e+02_dp,0.28000e+02_dp,0.36000e+02_dp,0.46000e+02_dp,0.56000e+02_dp, &
5442 0.67000e+02_dp,0.80000e+02_dp,0.97000e+02_dp,0.12000e+03_dp,0.14000e+03_dp,0.17000e+03_dp, &
5443 0.20000e+03_dp,0.24000e+03_dp,0.29000e+03_dp,0.35000e+03_dp,0.42000e+03_dp,0.51000e+03_dp, &
5444 0.61000e+03_dp,0.74000e+03_dp,0.89000e+03_dp /
5445 data blkradl(1:ndat) / &
5446 0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
5447 0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
5448 0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
5449 0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
5450 0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
5451 0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
5452 data blkradi(1:ndat,1) / 0.57452e-03_dp,0.72384e-03_dp,0.91199e-03_dp, &
5453 0.11490e-02_dp,0.14477e-02_dp,0.18240e-02_dp,0.22981e-02_dp,0.28954e-02_dp,0.36479e-02_dp, &
5454 0.45961e-02_dp,0.57908e-02_dp,0.72959e-02_dp,0.11572e-01_dp,0.14770e-01_dp,0.18851e-01_dp, &
5455 0.24060e-01_dp,0.30709e-01_dp,0.39194e-01_dp,0.50025e-01_dp,0.63848e-01_dp,0.81491e-01_dp, &
5456 0.10401e+00_dp,0.13275e+00_dp,0.16943e+00_dp,0.21625e+00_dp,0.27601e+00_dp,0.35228e+00_dp, &
5457 0.44962e+00_dp,0.57387e+00_dp,0.73244e+00_dp,0.93484e+00_dp,0.11932e+01_dp,0.15229e+01_dp /
5458 data blkradi(1:ndat,2) / &
5459 0.20715e-03_dp,0.26099e-03_dp,0.32912e-03_dp,0.43555e-03_dp,0.57638e-03_dp,0.76276e-03_dp, &
5460 0.10094e-02_dp,0.13358e-02_dp,0.17677e-02_dp,0.23394e-02_dp,0.30958e-02_dp,0.40969e-02_dp, &
5461 0.54216e-02_dp,0.71748e-02_dp,0.94948e-02_dp,0.12565e-01_dp,0.16628e-01_dp,0.22005e-01_dp, &
5462 0.29120e-01_dp,0.38537e-01_dp,0.50998e-01_dp,0.67488e-01_dp,0.89311e-01_dp,0.11819e+00_dp, &
5463 0.15641e+00_dp,0.20698e+00_dp,0.27391e+00_dp,0.36249e+00_dp,0.47970e+00_dp,0.63481e+00_dp, &
5464 0.84009e+00_dp,0.11117e+01_dp,0.14712e+01_dp /
5465 data blkradi(1:ndat,3) / 0.23559e-03_dp,0.29682e-03_dp,0.37397e-03_dp, &
5466 0.47118e-03_dp,0.59365e-03_dp,0.74795e-03_dp,0.95894e-03_dp,0.12777e-02_dp,0.17025e-02_dp, &
5467 0.22685e-02_dp,0.30227e-02_dp,0.40275e-02_dp,0.53665e-02_dp,0.71506e-02_dp,0.95278e-02_dp, &
5468 0.12695e-01_dp,0.16916e-01_dp,0.22539e-01_dp,0.30032e-01_dp,0.40017e-01_dp,0.70019e-01_dp, &
5469 0.98384e-01_dp,0.13824e+00_dp,0.19424e+00_dp,0.27293e+00_dp,0.38350e+00_dp,0.53885e+00_dp, &
5470 0.75714e+00_dp,0.10639e+01_dp,0.14948e+01_dp,0.21004e+01_dp,0.29513e+01_dp,0.41469e+01_dp /
5471 data blkrads(1:ndat) / &
5472 0.20715e-03_dp,0.26148e-03_dp,0.33067e-03_dp,0.41710e-03_dp,0.52691e-03_dp,0.66640e-03_dp, &
5473 0.84289e-03_dp,0.10674e-02_dp,0.13513e-02_dp,0.17129e-02_dp,0.21670e-02_dp,0.27521e-02_dp, &
5474 0.34989e-02_dp,0.44777e-02_dp,0.57347e-02_dp,0.75389e-02_dp,0.99020e-02_dp,0.13161e-01_dp, &
5475 0.17372e-01_dp,0.23337e-01_dp,0.31058e-01_dp,0.41194e-01_dp,0.55153e-01_dp,0.74854e-01_dp, &
5476 0.99806e-01_dp,0.13463e+00_dp,0.18136e+00_dp,0.24282e+00_dp,0.32955e+00_dp,0.44123e+00_dp, &
5477 0.59884e+00_dp,0.77090e+00_dp,0.99387e+00_dp /
5478 data blkradg(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
5479 0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
5480 0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
5481 0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
5482 0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
5483 0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
5484 data blkradh(1:ndat) / &
5485 0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
5486 0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
5487 0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
5488 0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
5489 0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
5490 0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
5491 data blkdnsl(1:ndat) / &
5492 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5493 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5494 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5495 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5496 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
5497 0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp /
5498 data blkdnsi(1:ndat,1) / 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5499 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5500 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.87368e+00_dp,0.87072e+00_dp,0.86777e+00_dp, &
5501 0.86483e+00_dp,0.86189e+00_dp,0.85897e+00_dp,0.85606e+00_dp,0.85316e+00_dp,0.85026e+00_dp, &
5502 0.84738e+00_dp,0.84451e+00_dp,0.84164e+00_dp,0.83879e+00_dp,0.83595e+00_dp,0.83311e+00_dp, &
5503 0.83029e+00_dp,0.82747e+00_dp,0.82467e+00_dp,0.82187e+00_dp,0.81908e+00_dp,0.81631e+00_dp /
5504 data blkdnsi(1:ndat,2) / &
5505 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5506 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5507 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5508 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5509 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5510 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
5511 data blkdnsi(1:ndat,3) / 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
5512 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
5513 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
5514 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.51790e+00_dp, &
5515 0.45557e+00_dp,0.40075e+00_dp,0.35252e+00_dp,0.31010e+00_dp,0.27278e+00_dp,0.23995e+00_dp, &
5516 0.21108e+00_dp,0.18567e+00_dp,0.16333e+00_dp,0.14367e+00_dp,0.12638e+00_dp,0.11118e+00_dp /
5517 data blkdnss(1:ndat) / &
5518 0.90000e+00_dp,0.89500e+00_dp,0.88500e+00_dp,0.88200e+00_dp,0.87500e+00_dp,0.86500e+00_dp, &
5519 0.85500e+00_dp,0.84200e+00_dp,0.83000e+00_dp,0.81500e+00_dp,0.80500e+00_dp,0.78600e+00_dp, &
5520 0.76500e+00_dp,0.73000e+00_dp,0.69500e+00_dp,0.61183e+00_dp,0.54000e+00_dp,0.46000e+00_dp, &
5521 0.40000e+00_dp,0.33000e+00_dp,0.28000e+00_dp,0.24000e+00_dp,0.20000e+00_dp,0.16000e+00_dp, &
5522 0.13500e+00_dp,0.11000e+00_dp,0.90000e-01_dp,0.75000e-01_dp,0.60000e-01_dp,0.50000e-01_dp, &
5523 0.40000e-01_dp,0.37500e-01_dp,0.35000e-01_dp /
5524 data blkdnsg(1:ndat) / &
5525 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5526 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5527 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5528 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5529 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
5530 0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp /
5531 data blkdnsh(1:ndat) / &
5532 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5533 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5534 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5535 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5536 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
5537 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
5541 xmss( il,k ) = log( dble( xl( k ) )*1.e-03_dp )
5542 xmss( ic,k ) = log( dble( xi( k,1 ) )*1.e-03_dp )
5543 xmss( ip,k ) = log( dble( xi( k,2 ) )*1.e-03_dp )
5544 xmss( id,k ) = log( dble( xi( k,3 ) )*1.e-03_dp )
5545 xmss( is,k ) = log( dble( xs( k ) )*1.e-03_dp )
5546 xmss( ig,k ) = log( dble( xg( k ) )*1.e-03_dp )
5547 xmss( ih,k ) = log( dble( xh( k ) )*1.e-03_dp )
5552 zcap( il,k ) = dble(rlec( k ))*1.e-02_dp
5553 zcap( ic,k ) = dble(riec( k,1 ))*1.e-02_dp
5554 zcap( ip,k ) = dble(riec( k,2 ))*1.e-02_dp
5555 zcap( id,k ) = dble(riec( k,3 ))*1.e-02_dp
5556 zcap( is,k ) = dble(rsec( k ))*1.e-02_dp
5557 zcap( ig,k ) = dble(rgec( k ))*1.e-02_dp
5558 zcap( ih,k ) = dble(rhec( k ))*1.e-02_dp
5563 vtrm( il,k ) = dble(vrl( k ))*1.e-02_dp
5564 vtrm( ic,k ) = dble(vri( k,1 ))*1.e-02_dp
5565 vtrm( ip,k ) = dble(vri( k,2 ))*1.e-02_dp
5566 vtrm( id,k ) = dble(vri( k,3 ))*1.e-02_dp
5567 vtrm( is,k ) = dble(vrs( k ))*1.e-02_dp
5568 vtrm( ig,k ) = dble(vrg( k ))*1.e-02_dp
5569 vtrm( ih,k ) = dble(vrh( k ))*1.e-02_dp
5574 blkr( il,k ) = dble(blkradl( k ))*10.0_dp
5575 blkr( ic,k ) = dble(blkradi( k,1 ))*10.0_dp
5576 blkr( ip,k ) = dble(blkradi( k,2 ))*10.0_dp
5577 blkr( id,k ) = dble(blkradi( k,3 ))*10.0_dp
5578 blkr( is,k ) = dble(blkrads( k ))*10.0_dp
5579 blkr( ig,k ) = dble(blkradg( k ))*10.0_dp
5580 blkr( ih,k ) = dble(blkradh( k ))*10.0_dp
5585 blkd( il,k ) = dble(blkdnsl( k ))*1000._dp
5586 blkd( ic,k ) = dble(blkdnsi( k,1 ))*1000._dp
5587 blkd( ip,k ) = dble(blkdnsi( k,2 ))*1000._dp
5588 blkd( id,k ) = dble(blkdnsi( k,3 ))*1000._dp
5589 blkd( is,k ) = dble(blkdnss( k ))*1000._dp
5590 blkd( ig,k ) = dble(blkdnsg( k ))*1000._dp
5591 blkd( ih,k ) = dble(blkdnsh( k ))*1000._dp
5600 ykrn( k,kk,i,j ) = kernels( k,kk,i,j )*1.e-06_dp
5607 end subroutine rdkdat
5611 real(
dp) :: xsta, xend
5613 real(
dp):: pi_dp = 3.1415920_dp
5614 real(
dp):: rhow_dp = 1.0e+03_dp
5616 xsta = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-06_dp )**3.0_dp )
5617 xend = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-03_dp )**3.0_dp )
5619 dxmic_mk = ( xend-xsta )/
nbin
5622 xbnd_mk( n ) = xsta + dxmic_mk*( n-1 )
5625 xctr_mk( n ) = ( xbnd_mk( n )+xbnd_mk( n+1 ) )*0.50_dp
5626 radc_mk( n ) = ( exp( xctr_mk( n ) )*3.0_dp/4.0_dp/pi_dp/rhow_dp )**( 1.0_dp/3.0_dp )
5629 end subroutine sdfgrid
5637 cctr_mk( myu,n ) = fcpc( myu,xctr_mk( n ) )
5640 cbnd_mk( myu,n ) = fcpc( myu,xbnd_mk( n ) )
5644 end subroutine getcp
5646 function fcpc( myu,x )
5648 integer,
intent(in) :: myu
5649 real(
dp),
intent(in) :: x
5652 real(
dp) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5655 ( ndat, kdeg, xmss( myu,: ), &
5659 ( ndat, kdeg, qknt, xmss( myu,: ), &
5663 ( ndat, kdeg, elm, zcap( myu,: ), &
5666 fcpc = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
5672 integer :: myu, nyu, i, j
5674 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Create micpara.dat'
5675 if( kphase == 0 )
then
5676 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Hydro-dynamic kernel'
5677 else if( kphase == 1 )
then
5678 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Long Kernel'
5679 else if( kphase == 2 )
then
5680 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
'Golovin Kernel'
5685 log_info(
"ATMOS_PHY_MP_SUZUKI10_getck",*)
' myu, nyu :', myu, nyu
5688 ck_mk( myu,nyu,i,j ) = fckrn( myu,nyu,xctr_mk( i ),xctr_mk( j ) )
5696 end subroutine getck
5698 function fckrn( myu,nyu,x,y )
5700 integer,
intent(in) :: myu, nyu
5701 real(
dp),
intent(in) :: x, y
5704 real(
dp) :: qknt( ndat+kdeg ), rknt( ndat+kdeg )
5705 real(
dp) :: coef( ndat,ndat )
5707 real(
dp) :: xlrg, xsml, vlrg, vsml, rlrg
5708 real(
dp):: pi_dp = 3.1415920_dp
5709 real(
dp):: rhow_dp = 1.0e+03_dp
5711 if( kphase == 0 )
then
5713 ( ndat, kdeg, xmss( myu,: ), &
5716 rknt( : ) = qknt( : )
5719 ( ndat, ndat, kdeg, kdeg, &
5720 xmss( myu,: ), xmss( nyu,: ), &
5722 ykrn( myu,nyu,:,: ), &
5726 ( ndat, ndat, kdeg, kdeg, &
5728 xmss( myu,: ), xmss( nyu,: ), &
5730 else if( kphase == 1 )
then
5734 vlrg = (exp( xlrg ) / rhow_dp )*1.e+06_dp
5735 vsml = (exp( xsml ) / rhow_dp )*1.e+06_dp
5737 rlrg = ( exp( xlrg )/( 4.0_dp*pi_dp*rhow_dp )*3.0_dp )**(1.0_dp/3.0_dp )*1.e+06_dp
5739 if( rlrg <=50.0_dp )
then
5740 fckrn = 9.44e+03_dp*( vlrg*vlrg + vsml*vsml )
5742 fckrn = 5.78e-03_dp*( vlrg+vsml )
5744 else if( kphase == 2 )
then
5745 fckrn = 1.5_dp*( exp(x) +exp(y) )
5758 vt_mk( myu,n ) = max( fvterm( myu,xctr_mk( n ) ), 0.0_dp )
5762 end subroutine getvt
5764 function fvterm( myu,x )
5766 integer,
intent(in) :: myu
5767 real(
dp),
intent(in) :: x
5770 real(
dp) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5773 ( ndat, kdeg, xmss( myu,: ), &
5777 ( ndat, kdeg, qknt, xmss( myu,: ), &
5781 ( ndat, kdeg, elm, vtrm( myu,: ), &
5784 fvterm = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
5794 br_mk( myu,n ) = fbulkrad( myu, xctr_mk( n ) )
5798 end subroutine getbr
5800 function fbulkrad( myu,x )
5802 integer,
intent(in) :: myu
5803 real(
dp),
intent(in) :: x
5804 real(
dp) :: fbulkrad
5806 real(
dp) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
5809 ( ndat, kdeg, xmss( myu,: ), &
5813 ( ndat, kdeg, qknt, xmss( myu,: ), &
5817 ( ndat, kdeg, elm, blkr( myu,: ), &
5820 fbulkrad = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
5822 end function fbulkrad
5826 integer :: myu, nyu, i, j, n
5828 open ( fid_micpara, file = fname_micpara, form =
'formatted', status=
'new' )
5830 write( fid_micpara,* ) nspc_mk,
nbin
5834 xctr( n ) = xctr_mk( n )
5835 radc( n ) = radc_mk( n )
5836 write( fid_micpara,* ) n, xctr( n ), radc( n )
5839 xbnd( n ) = xbnd_mk( n )
5840 write( fid_micpara,* ) n, xbnd( n )
5842 write( fid_micpara,* ) dxmic_mk
5847 cctr( n,myu ) = cctr_mk( myu,n )
5848 write( fid_micpara,* ) myu, n, cctr( n,myu )
5851 cbnd( n,myu ) = cbnd_mk( myu,n )
5852 write( fid_micpara,* ) myu, n, cbnd( n,myu )
5861 ck( myu,nyu,i,j ) = ck_mk( myu,nyu,i,j )
5862 write( fid_micpara,* ) myu, nyu, i, j, ck( myu,nyu,i,j )
5871 vt( myu,n ) = vt_mk( myu,n )
5872 write( fid_micpara,* ) myu, n, vt( myu,n )
5879 br( myu,n ) = br_mk( myu,n )
5880 write( fid_micpara,* ) myu, n, br( myu,n )
5884 close ( fid_micpara )
5886 end subroutine paraout
5890 subroutine tinvss(n,a,dt,e,nn,iw,inder)
5894 integer,
intent(in) :: n, nn
5895 integer,
intent(inout) :: inder
5896 real(
dp),
intent(inout) :: a(nn,n)
5897 real(
dp),
intent(inout) :: dt
5898 real(
dp),
intent(in) :: e
5899 integer,
intent(inout) :: iw( 2*n )
5900 integer :: i, j, k, kk, ij, nnn
5901 real(
dp) :: work, aa, az, eps
5902 integer :: ipiv, jpiv
5909 elseif( n == 1 )
then
5911 elseif( n > 1 )
then
5918 write(6,690)
"n= ", n,
"nn= ", nn, &
5919 "n should be less than or equal to nn in TINVSS"
5920 690
format(a8,i5,5x,a4,i5,a55)
5930 if( abs(a(i,j)) <= abs(piv) )
goto 110
5938 if( abs(piv) <= eps )
goto 920
5939 if( k == 1 ) eps = abs(piv)*e
5940 if( ipiv == k )
goto 130
5948 if( jpiv == k )
goto 150
5963 if( i == k )
goto 220
5965 if( az == 0.0_dp )
goto 220
5967 a(i,j) = a(i,j)-a(k,j)*az
5977 if( ij == k )
goto 420
5985 if( ij == k )
goto 400
5998 write(*,691)
"n= ", n,
"should be positive in TINVSS"
5999 691
format(a8,i5,5x,a30)
6007 write(*,692)
'given matrix A to TINVSS is ill conditioned, or sigular withrank =', nnn, &
6008 'return with no further calculation'
6009 692
format(a,1x,i4,1x,a)
6015 if( dt == 0.0_dp )
goto 920
6016 a(1,1) = 1.0_dp/a(1,1)
6019 end subroutine tinvss
6021 subroutine getknot &
6022 ( ndat, kdeg, xdat, &
6025 integer,
intent(in) :: ndat
6026 integer,
intent(in) :: kdeg
6027 real(
dp),
intent(in) :: xdat( ndat )
6029 real(
dp),
intent(out) :: qknt( ndat+kdeg )
6035 qknt( i ) = xdat( 1 )
6039 qknt( i+kdeg ) = ( xdat( i )+xdat( i+kdeg ) )*0.50_dp
6043 qknt( ndat+i ) = xdat( ndat )
6048 end subroutine getknot
6050 recursive function fbspl ( ndat, inum, kdeg, qknt, xdat, x ) &
6055 integer,
intent(in) :: ndat
6056 integer,
intent(in) :: inum
6057 integer,
intent(in) :: kdeg
6058 real(
dp),
intent(in) :: qknt( ndat+kdeg )
6059 real(
dp),
intent(in) :: xdat( ndat )
6060 real(
dp),
intent(in) :: x
6063 real(
dp) :: bsp1, bsp2
6065 if ( ( inum == 1 .AND. x == xdat( 1 ) ) .OR. &
6066 ( inum == ndat .AND. x == xdat( ndat ) ) )
then
6071 if ( kdeg == 1 )
then
6072 if ( x >= qknt( inum ) .AND. x < qknt( inum+1 ) )
then
6078 if ( qknt( inum+kdeg-1 ) /= qknt( inum ) )
then
6079 bsp1 = ( x-qknt( inum ) ) &
6080 /( qknt( inum+kdeg-1 )-qknt( inum ) ) &
6081 * fbspl( ndat, inum, kdeg-1, qknt, xdat, x )
6085 if ( qknt( inum+kdeg ) /= qknt( inum+1 ) )
then
6086 bsp2 = ( qknt( inum+kdeg )-x ) &
6087 /( qknt( inum+kdeg )-qknt( inum+1 ) ) &
6088 * fbspl( ndat, inum+1, kdeg-1, qknt, xdat, x )
6097 function fpb( ndat, i, kdeg, qknt, xdat, elm, x )
6100 integer :: ndat, i, kdeg
6101 real(
dp) :: qknt( ndat+kdeg ), xdat( ndat ), elm( ndat,ndat )
6109 sum = sum + elm( l,i )*fbspl( ndat, l, kdeg, qknt, xdat, x )
6116 subroutine getmatrx &
6117 ( ndat, kdeg, qknt, xdat, &
6122 integer,
intent(in) :: ndat
6123 integer,
intent(in) :: kdeg
6124 real(
dp),
intent(in) :: qknt( ndat+kdeg )
6125 real(
dp),
intent(in) :: xdat( ndat )
6127 real(
dp),
intent(out) :: elm( ndat,ndat )
6131 integer :: iw( 2*ndat ), i, j, inder
6132 real(
dp),
parameter :: eps = 0.
6136 elm( i,j ) = fbspl( ndat, j, kdeg, qknt, xdat, xdat( i ) )
6140 call tinvss( ndat, elm, dt, eps, ndat, iw, inder )
6144 end subroutine getmatrx
6146 subroutine getcoef &
6147 ( ndat, kdeg, elm, ydat, &
6150 integer,
intent(in) :: ndat
6151 integer,
intent(in) :: kdeg
6152 real(
dp),
intent(in) :: elm( ndat,ndat )
6153 real(
dp),
intent(in) :: ydat( ndat )
6155 real(
dp),
intent(out) :: coef( ndat )
6164 sum = sum + elm( i,j )*ydat( j )
6171 end subroutine getcoef
6173 function fspline ( ndat, kdeg, coef, qknt, xdat, x )
6175 integer,
intent(in) :: ndat
6176 integer,
intent(in) :: kdeg
6177 real(
dp),
intent(in) :: coef( ndat )
6178 real(
dp),
intent(in) :: qknt( ndat+kdeg )
6179 real(
dp),
intent(in) :: xdat( ndat )
6180 real(
dp),
intent(in) :: x
6190 sum = sum + coef( i )*fbspl( ndat, i, kdeg, qknt, xdat, x )
6197 end function fspline
6199 subroutine getcoef2 &
6200 ( mdat, ndat, kdeg, ldeg, &
6201 xdat, ydat, qknt, rknt, zdat, &
6206 integer,
intent(in) :: mdat
6207 integer,
intent(in) :: ndat
6208 integer,
intent(in) :: kdeg
6209 integer,
intent(in) :: ldeg
6210 real(
dp),
intent(in) :: xdat( mdat )
6211 real(
dp),
intent(in) :: ydat( ndat )
6212 real(
dp),
intent(in) :: qknt( mdat+kdeg )
6213 real(
dp),
intent(in) :: rknt( ndat+ldeg )
6214 real(
dp),
intent(in) :: zdat( mdat,ndat )
6216 real(
dp),
intent(out) :: coef( mdat,ndat )
6219 real(
dp) :: elmx( mdat,mdat ), elmy( ndat,ndat )
6220 integer :: iw1( 2*mdat ), iw2( 2*ndat )
6221 real(
dp) :: beta( mdat,ndat ), sum, dt
6222 real(
dp),
parameter :: eps = 0.0_dp
6223 integer :: i, j, k, l, inder
6227 elmx( i,j ) = fbspl( mdat, j, kdeg, qknt, xdat, xdat( i ) )
6230 call tinvss( mdat, elmx, dt, eps, mdat, iw1, inder )
6236 sum = sum + elmx( i,j )*zdat( j,l )
6244 elmy( i,j ) = fbspl( ndat, j, ldeg, rknt, ydat, ydat( i ) )
6247 call tinvss( ndat, elmy, dt, eps, ndat, iw2, inder )
6253 sum = sum + elmy( i,j )*beta( k,j )
6261 end subroutine getcoef2
6264 ( mdat, ndat, kdeg, ldeg, &
6265 coef, qknt, rknt, xdat, ydat, &
6268 integer,
intent(in) :: mdat
6269 integer,
intent(in) :: ndat
6270 integer,
intent(in) :: kdeg
6271 integer,
intent(in) :: ldeg
6272 real(
dp),
intent(in) :: coef( mdat,ndat )
6273 real(
dp),
intent(in) :: qknt( mdat+kdeg )
6274 real(
dp),
intent(in) :: rknt( ndat+ldeg )
6275 real(
dp),
intent(in) :: xdat( mdat ), ydat( ndat )
6276 real(
dp),
intent(in) :: x, y
6278 real(
dp) :: fspline2
6281 real(
dp) :: sum, add
6287 add = coef( i,j )*fbspl( mdat, i, kdeg, qknt, xdat, x ) &
6288 *fbspl( ndat, j, ldeg, rknt, ydat, y )
6297 end function fspline2