SCALE-RM
scale_atmos_phy_mp_suzuki10.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
29 !-------------------------------------------------------------------------------
30 #include "macro_thermodyn.h"
31 module scale_atmos_phy_mp_suzuki10
32  !-----------------------------------------------------------------------------
33  !
34  !++ Used modules
35  !
36  use scale_precision
37  use scale_stdio
38  use scale_prof
39  use scale_grid_index
40 
41  use scale_const, only: &
42  pi => const_pi, &
43  eps => const_eps, &
44  const_cpvap, &
45  const_cpdry, &
46  const_cvdry, &
47  const_cl, &
48  const_ci, &
49  const_dwatr, &
50  const_grav, &
51  const_rvap, &
52  const_rdry, &
53  const_lhv0, &
54  const_lhs0, &
55  const_emelt, &
56  const_tem00, &
57  const_tmelt, &
58  const_psat0, &
59  const_pre00, &
60  const_epsvap, &
61  const_thermodyn_type
62  use scale_const, only: &
63  cp => const_cpdry, &
64  rvap => const_rvap, &
65  esat0 => const_psat0, &
66  qlmlt => const_emelt, &
67  tmlt => const_tmelt, &
68  temp0 => const_tem00, &
69  rhow => const_dwatr
70  use scale_atmos_saturation, only: &
71  saturation_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
72  saturation_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
73  lovr_liq, &
74  lovr_ice, &
75  cpovr_liq, &
76  cpovr_ice
77  use scale_atmos_thermodyn, only: &
78  thermodyn_temp_pres => atmos_thermodyn_temp_pres, &
79  thermodyn_pott => atmos_thermodyn_pott
80  use scale_atmos_hydrometeor, only: &
81  n_hyd
82  !-----------------------------------------------------------------------------
83  implicit none
84  private
85  !-----------------------------------------------------------------------------
86  !
87  !++ Public procedure
88  !
89  public :: atmos_phy_mp_suzuki10_config
90  public :: atmos_phy_mp_suzuki10_setup
91  public :: atmos_phy_mp_suzuki10
92  public :: atmos_phy_mp_suzuki10_cloudfraction
93  public :: atmos_phy_mp_suzuki10_effectiveradius
94  public :: atmos_phy_mp_suzuki10_mixingratio
95 
96  !-----------------------------------------------------------------------------
97  !
98  !++ Public parameters & variables
99  !
100  !-----------------------------------------------------------------------------
101  integer, private :: QA_MP
102 
103  character(len=H_SHORT), public, target, allocatable :: atmos_phy_mp_suzuki10_name(:)
104  character(len=H_MID) , public, target, allocatable :: atmos_phy_mp_suzuki10_desc(:)
105  character(len=H_SHORT), public, target, allocatable :: atmos_phy_mp_suzuki10_unit(:)
106 
107  real(RP), public, target :: atmos_phy_mp_suzuki10_dens(n_hyd) ! hydrometeor density [kg/m3]=[g/L]
108 
109  integer, public :: nbin = 33
110  integer, public :: nspc = 7
111  integer, public :: nccn = 0
112  integer, public :: nccn1 = 0
113  integer, public :: kphase = 0
114  integer, public :: iceflg = 1
115 
116  character(len=3) :: namspc(8) =(/'Qcl','Qic','Qip','Qid','Qis','Qig','Qih','Qae'/)
117  character(len=27) :: lnamspc(8) = &
118  (/'Mixing ratio of cloud bin', &
119  'Mixing ratio of colum bin', &
120  'Mixing ratio of plate bin', &
121  'Mixing ratio of dendrit bin', &
122  'Mixing ratio of snow bin', &
123  'Mixing ratio of graupel bin', &
124  'Mixing ratio of hail bin', &
125  'Mixing ratio of aerosol bin' /)
126 
127 # include "kernels.h"
128  !-----------------------------------------------------------------------------
129  !
130  !++ Private procedure
131  !
132  !-----------------------------------------------------------------------------
133  !
134  !++ Private parameters
135  !
136  integer, private, parameter :: i_mp_qc = 1
137  integer, private, parameter :: i_mp_qp = 2
138  integer, private, parameter :: i_mp_qcl = 3
139  integer, private, parameter :: i_mp_qd = 4
140  integer, private, parameter :: i_mp_qs = 5
141  integer, private, parameter :: i_mp_qg = 6
142  integer, private, parameter :: i_mp_qh = 7
143  integer, private :: qs_mp
144  integer, private :: qe_mp
145  integer, private :: i_qv
146 
147  logical, private :: mp_doautoconversion = .true. ! apply collision process ?
148  logical, private :: mp_doprecipitation = .true. ! apply sedimentation of hydrometeor ?
149  logical, private :: mp_donegative_fixer = .true. ! apply negative fixer?
150  logical, private :: mp_couple_aerosol = .false. ! apply CCN effect?
151  real(RP), private :: mp_limit_negative = 1.0_rp ! Abort if abs(fixed negative vaue) > abs(MP_limit_negative)
152 
153 
154  !--- Indeces for determining species of cloud particle
155  integer, parameter :: il = 1 !--- index for liquid water
156  integer, parameter :: ic = 2 !--- index for columnar ice
157  integer, parameter :: ip = 3 !--- index for plate ice
158  integer, parameter :: id = 4 !--- index for dendrite ice
159  integer, parameter :: iss= 5 !--- index for snow
160  integer, parameter :: ig = 6 !--- index for graupel
161  integer, parameter :: ih = 7 !--- index for hail
162 
163  !--- bin information of hydrometeors
164  real(RP) :: dxmic !--- d( log(m) ) of hydrometeor bin
165  real(RP), allocatable :: xctr( : ) !--- log( m ) value of bin center [xctr = exp( 4/3*pi*DWATER*radc^3 )]
166  real(RP), allocatable :: xbnd( : ) !--- log( m ) value of bin boundary
167  real(RP), allocatable :: radc( : ) !--- radius of hydrometeor at bin center [m]
168  real(RP), allocatable :: cctr( :,: ) !--- capacitance of hydrometeor at bin center (C of equation A.17 in Suzuki 2004)
169  real(RP), allocatable :: cbnd( :,: ) !--- capacitance of hydrometeor at bin boundary (C of equation A.17 in Suzuki 2004)
170  real(RP), allocatable :: ck( :,:,:,: ) !-- collection kernel (K of equation A.20 in Suzuki 2004)
171  real(RP), allocatable :: vt( :,: ) !--- terminal velocity of hydrometeor [m/s]
172  real(RP), allocatable :: br( :,: ) !--- bulk density of hydrometeor [kg/m^3]
173  integer, allocatable :: ifrsl( :,:,: ) !--- type of species after collision
174  !--- bin information of aerosol (not supported)
175  real(RP), allocatable :: xactr( : ) !--- log( ma ) value of bin center ( ma is mass of aerosol )
176  real(RP), allocatable :: xabnd( : ) !--- log( ma ) value of bin boundary ( ma is mass of aerosol )
177  real(RP), allocatable :: rada( : ) !--- radius of aerosol at bin center [m]
178 
179  real(RP), allocatable :: expxctr( : ) !--- exp( xctr )
180  real(RP), allocatable :: expxbnd( : ) !--- exp( xbnd )
181  real(RP), allocatable :: expxactr( : ) !--- exp( xactr )
182  real(RP), allocatable :: expxabnd( : ) !--- exp( xabnd )
183  real(RP), allocatable :: rexpxctr( : ) !--- 1/exp( xctr )
184  real(RP), allocatable :: rexpxbnd( : ) !--- 1/exp( xbnd )
185  real(RP), allocatable :: rexpxactr( : ) !--- 1/exp( xactr )
186  real(RP), allocatable :: rexpxabnd( : ) !--- 1/exp( xabnd1 )
187 
188  real(RP) :: dxaer !--- d( log(ma) ) of aerosol bin
189  real(RP) :: xasta !--- exponential of mass of aerosol for smallest aerosol bin
190  real(RP) :: xaend !--- exponential of mass of aerosol for largest aerosol bin
191  real(RP), allocatable, save :: vterm( :,:,:,: ) ! terminal velocity
192  integer, private, save :: mp_nstep_sedimentation ! number of fractional step for sedimentation
193  real(RP), private, save :: mp_rnstep_sedimentation ! 1/MP_NSTEP_SEDIMENTATION
194  real(DP), private, save :: mp_dtsec_sedimentation ! DT for sedimentation
195  integer, private, save :: mp_ntmax_sedimentation= 1 ! maxinum fractional step
196  real(RP), private :: flg_thermodyn !--- flg for lhv and lhs (0 -> SIMPLE, 1 -> EXACT )
197  real(RP), private :: rtem00 !--- 1/CONST_TEM00
198 
199  !--- constant for bin
200  real(RP), parameter :: cldmin = 1.0e-10_rp !--- threshould for cloud is regarded existing
201  real(RP), parameter :: oneovthird = 1.0_rp/3.0_rp
202  real(RP), parameter :: thirdovforth = 3.0_rp/4.0_rp
203  real(RP), parameter :: twoovthird = 2.0_rp/3.0_rp
204 
205  real(RP) :: rbnd = 40.e-06_rp ! boundary radius of cloud and rian
206  integer :: nbnd ! boundary bin number corresponding to rbnd
207  !--- constant for aerosol
208  real(RP) :: rhoa = 2.25e+03_rp ! density of aerosol ( NaCl )
209  real(RP) :: emaer = 58.0_rp ! molecular weight of aerosol ( salt )
210  real(RP) :: emwtr = 18.0_rp ! molecular weight of water
211  real(RP) :: rasta = 1.e-08_rp ! minimum radius of aerosol (m)
212  real(RP) :: raend = 1.e-06_rp ! maximum radius of aerosol (m)
213  real(RP) :: r0a = 1.e-07_rp ! average radius of aerosol (m)
214  logical :: flg_regeneration=.false. ! flag regeneration of aerosol
215  logical :: flg_nucl=.false. ! flag nucleated cloud move into smallest bin
216  logical :: flg_icenucl=.false. ! flag ice nucleation
217  logical :: flg_sf_aero =.false. ! flag surface flux of aerosol
218  integer, private, save :: rndm_flgp = 0 ! flag for sthastic integration for coll.-coag.
219 
220  real(RP), allocatable :: marate( : ) ! mass rate of each aerosol bin to total aerosol mass
221  integer, allocatable, save :: ncld( : ) ! bin number of aerosol in bin of hydrometeor
222 ! integer, private, save :: K10_1, K10_2 ! scaling factor for 10m value (momentum)
223 ! real(RP), private :: R10M1, R10M2 ! scaling factor for 10m value (momentum)
224 ! real(RP), private :: R10H1, R10H2 ! scaling factor for 10m value (heat)
225 ! real(RP), private :: R10E1, R10E2 ! scaling factor for 10m value (tracer)
226 
227  character(len=11), parameter :: fname_micpara="micpara.dat" !--- file name
228  integer(4) :: fid_micpara
229 
230  !--- Use for stochastic method
231  integer, allocatable :: blrg( :,: ), bsml( :,: )
232  real(RP) :: wgtbin
233  integer :: mspc, mbin
234  real(RP), private :: rndm(1,1,1)
235 
236  !--- use for model without aerosol
237  real(RP), private :: c_ccn = 100.e+6_rp ! N0 of Nc = N0*s^kappa
238  real(RP), private :: kappa = 0.462_rp ! kappa of Nc = N0*s^kappa
239  !--- use for aerosol coupled model
240  real(RP), private :: sigma = 7.5e-02_rp ! water surface tension [ N/m2 ] (sigma in eq. (A.11) of Suzuki (2004) )
241  real(RP), private :: vhfct = 2.0_rp ! van't hoff factor (i in eq.(A.11) of Suzuki (2004))
242 
243  real(RP), parameter :: tcrit = 271.15_rp
244  integer, private, allocatable :: kindx( :,: )
245 
246  !--- for creating micpara.dat (mkpara)
247  integer, parameter :: ndat = 33, icemax = 3
248  integer, parameter :: kdeg = 4, ldeg = 4, nspc_mk = 7
249  real(DP) :: dxmic_mk
250  real(DP), allocatable :: radc_mk( : ), xctr_mk( : ), xbnd_mk( : )
251  real(DP), allocatable :: cctr_mk( :,: ), cbnd_mk( :,: )
252  real(DP), allocatable :: ck_mk( :,:,:,: )
253  real(DP), allocatable :: vt_mk( :,: )
254  real(DP), allocatable :: br_mk( :,: )
255  real(DP) :: xmss( nspc_mk,ndat ), zcap( nspc_mk,ndat ), vtrm( nspc_mk,ndat )
256  real(DP) :: blkr( nspc_mk,ndat ), blkd( nspc_mk,ndat ), ykrn( nspc_mk,nspc_mk,ndat,ndat )
257 
258  real(DP) :: ywll( ndat,ndat ), ywli( ndat,ndat,icemax ), ywls( ndat,ndat )
259  real(DP) :: ywlg( ndat,ndat ), ywlh( ndat,ndat )
260 
261  real(DP) :: ywil( ndat,ndat,icemax ), ywii( ndat,ndat,icemax,icemax )
262  real(DP) :: ywis( ndat,ndat,icemax ), ywig( ndat,ndat,icemax )
263  real(DP) :: ywih( ndat,ndat,icemax )
264 
265  real(DP) :: ywsl( ndat,ndat ), ywsi( ndat,ndat,icemax ), ywss( ndat,ndat )
266  real(DP) :: ywsg( ndat,ndat ), ywsh( ndat,ndat )
267 
268  real(DP) :: ywgl( ndat,ndat ), ywgi( ndat,ndat,icemax ), ywgs( ndat,ndat )
269  real(DP) :: ywgg( ndat,ndat ), ywgh( ndat,ndat )
270 
271  real(DP) :: ywhl( ndat,ndat ), ywhi( ndat,ndat,icemax ), ywhs( ndat,ndat )
272  real(DP) :: ywhg( ndat,ndat ), ywhh( ndat,ndat )
273 
274  !----------------------------------------------------------------------------
275 contains
276  !-----------------------------------------------------------------------------
278  subroutine atmos_phy_mp_suzuki10_config( &
279  MP_TYPE, &
280  QA, QS )
281  use scale_process, only: &
282  prc_mpistop
283  use scale_tracer, only: &
284  tracer_regist
285  use scale_atmos_hydrometeor, only: &
286  atmos_hydrometeor_regist
287  implicit none
288 
289  character(len=*), intent(in) :: mp_type
290  integer, intent(out) :: qa
291  integer, intent(out) :: qs
292 
293  namelist / param_bin / &
294  nbin, &
295  nccn, &
296  iceflg, &
297  kphase
298 
299  integer :: nl, ni
300  integer :: qs2
301  integer :: m, n, ierr
302  !---------------------------------------------------------------------------
303 
304  if( io_l ) write(io_fid_log,*)
305  if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
306  if( io_l ) write(io_fid_log,*) '*** Tracers for Suzuki (2010) Spectral BIN model'
307 
308  if ( mp_type /= 'SUZUKI10' ) then
309  write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not SUZUKI10. Check!'
310  call prc_mpistop
311  endif
312 
313 
314  if( io_l ) write(io_fid_log,*)
315  if( io_l ) write(io_fid_log,*) '+++ READ BIN NUMBER'
316 
317  rewind(io_fid_conf)
318  read(io_fid_conf,nml=param_bin,iostat=ierr)
319 
320  if( ierr < 0 ) then !--- missing
321  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
322  elseif( ierr > 0 ) then !--- fatal error
323  write(*,*) 'xxx Not appropriate names in namelist PARAM_BIN, Check!'
324  call prc_mpistop
325  end if
326 
327  if( io_nml ) write(io_fid_nml,nml=param_bin)
328 
329  if( iceflg == 0 ) then
330  nspc = 1
331  elseif( iceflg == 1 ) then
332  nspc = 7
333  else
334  write(*,*) "ICEFLG should be 0(warm rain) or 1(mixed rain) check!!"
335  call prc_mpistop
336  endif
337 
338  nccn1 = max(nccn,1)
339 
340  !-- setup QA_MP ...
341  qa_mp = 1 + nbin*nspc + nccn
342 
343  allocate( atmos_phy_mp_suzuki10_name(qa_mp) )
344  allocate( atmos_phy_mp_suzuki10_desc(qa_mp) )
345  allocate( atmos_phy_mp_suzuki10_unit(qa_mp) )
346 
347  !---------------------------------------------------------------------------
348  !
349  !++ calculate each category and aerosol
350  !
351  !---------------------------------------------------------------------------
352  do n = 1, qa_mp
353  write(atmos_phy_mp_suzuki10_unit(n),'(a)') 'kg/kg'
354  enddo
355 
356  write(atmos_phy_mp_suzuki10_name(1),'(a)') 'QV'
357 
358  do m = 1, nspc
359  do n = 1, nbin
360  write(atmos_phy_mp_suzuki10_name(1+nbin*(m-1)+n),'(a,i0)') trim(namspc(m)), n
361  enddo
362  enddo
363 
364  do n = 1, nccn
365  write(atmos_phy_mp_suzuki10_name(1+nbin*nspc+n),'(a,i0)') trim(namspc(8)), n
366  enddo
367 
368  write(atmos_phy_mp_suzuki10_desc(1),'(a)') 'Water Vapor mixing ratio'
369 
370  do m = 1, nspc
371  do n = 1, nbin
372  write(atmos_phy_mp_suzuki10_desc(1+nbin*(m-1)+n),'(a,i0)') trim(lnamspc(m)), n
373  enddo
374  enddo
375 
376  do n = 1, nccn
377  write(atmos_phy_mp_suzuki10_desc(1+nbin*nspc+n),'(a,i0)') trim(lnamspc(8)), n
378  enddo
379 
380  nl = nbin ! number of liquid water
381  ni = nbin * (nspc-1) ! number of ice water
382 
383  call atmos_hydrometeor_regist( qs, & ! [OUT]
384  1, nl, ni, & ! [IN]
385  atmos_phy_mp_suzuki10_name(1:nl+ni+1), & ! [IN]
386  atmos_phy_mp_suzuki10_desc(1:nl+ni+1), & ! [IN]
387  atmos_phy_mp_suzuki10_unit(1:nl+ni+1) ) ! [IN]
388 
389  if ( nccn > 0 ) then
390  call tracer_regist( qs2, & ! [OUT]
391  nccn, & ! [IN]
392  atmos_phy_mp_suzuki10_name(nl+ni+2:), & ! [IN]
393  atmos_phy_mp_suzuki10_desc(nl+ni+2:), & ! [IN]
394  atmos_phy_mp_suzuki10_unit(nl+ni+2:) ) ! [IN]
395  end if
396 
397  i_qv = qs
398  qa = qa_mp
399  qs_mp = qs
400  qe_mp = qs + qa_mp - 1
401 
402  return
403  end subroutine atmos_phy_mp_suzuki10_config
404 
405  !-----------------------------------------------------------------------------
407  subroutine atmos_phy_mp_suzuki10_setup
408  use scale_process, only: &
409  prc_mpistop, &
410  prc_masterrank, &
411  prc_ismaster
412  use scale_const, only: &
413  const_dwatr, &
414  const_dice
415  use scale_comm, only: &
416  comm_world, &
417  comm_datatype
418  use scale_grid, only: &
419  cdz => grid_cdz, &
420  cz => grid_cz
421  use scale_time, only: &
422  time_dtsec_atmos_phy_mp
423  use scale_tracer, only: &
424  qa
425  use scale_atmos_hydrometeor, only: &
426  i_hc, &
427  i_hr, &
428  i_hi, &
429  i_hs, &
430  i_hg, &
431  i_hh
432  implicit none
433 
434  real(RP) :: rho_aero !--- density of aerosol
435  real(RP) :: r0_aero !--- center radius of aerosol (um)
436  real(RP) :: r_min !--- minimum radius of aerosol (um)
437  real(RP) :: r_max !--- maximum radius of aerosol (um)
438  real(RP) :: s10_emaer !--- moleculer weight of aerosol
439  logical :: s10_flag_regene !--- flag of regeneration
440  logical :: s10_flag_nucleat !--- flag of regeneration
441  logical :: s10_flag_icenucleat !--- flag of regeneration
442  logical :: s10_flag_sfaero !--- flag of surface flux of aeorol
443  integer :: s10_rndm_flgp !--- flag of surface flux of aeorol
444  integer :: s10_rndm_mspc
445  integer :: s10_rndm_mbin
446 
447  namelist / param_atmos_phy_mp / &
448  mp_doprecipitation, &
449  mp_donegative_fixer, &
450  mp_limit_negative, &
451  mp_ntmax_sedimentation, &
452  mp_doautoconversion, &
453  mp_couple_aerosol
454 
455  namelist / param_atmos_phy_mp_suzuki10 / &
456  rho_aero, &
457  r_min, &
458  r_max, &
459  r0_aero, &
460  s10_emaer, &
461  s10_flag_regene, &
462  s10_flag_nucleat, &
463  s10_flag_icenucleat, &
464  s10_flag_sfaero, &
465  s10_rndm_flgp, &
466  s10_rndm_mspc, &
467  s10_rndm_mbin, &
468  c_ccn, kappa, &
469  sigma, vhfct
470 
471  real(RP), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
472  integer :: nstep_max
473  integer :: nnspc, nnbin
474  integer :: nn, mm, mmyu, nnyu
475  integer :: myu, nyu, i, j, k, n, ierr
476  !---------------------------------------------------------------------------
477 
478  if( io_l ) write(io_fid_log,*)
479  if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
480  if( io_l ) write(io_fid_log,*) '*** Suzuki (2010) Spectral BIN model'
481 
482  !--- allocation
483  allocate( xctr( nbin ) )
484  allocate( xbnd( nbin+1 ) )
485  allocate( radc( nbin ) )
486  allocate( cctr( nbin,nspc_mk ) )
487  allocate( cbnd( nbin+1,nspc_mk ) )
488  allocate( ck( nspc_mk,nspc_mk,nbin,nbin ) )
489  allocate( vt( nspc_mk,nbin ) )
490  allocate( br( nspc_mk,nbin ) )
491  allocate( ifrsl( 2,nspc_mk,nspc_mk ) )
492  allocate( expxctr( nbin ) )
493  allocate( expxbnd( nbin+1 ) )
494  allocate( rexpxctr( nbin ) )
495  allocate( rexpxbnd( nbin+1 ) )
496  if ( nccn /= 0 ) then
497  allocate( xactr( nccn ) )
498  allocate( xabnd( nccn+1 ) )
499  allocate( rada( nccn ) )
500  allocate( expxactr( nccn ) )
501  allocate( expxabnd( nccn+1 ) )
502  allocate( rexpxactr( nccn ) )
503  allocate( rexpxabnd( nccn+1 ) )
504  endif
505 
506  mbin = nbin/2
507  mspc = nspc_mk*nspc_mk
508 
509  rho_aero = rhoa
510  s10_emaer = emaer
511  r_min = rasta
512  r_max = raend
513  r0_aero = r0a
514  s10_flag_regene = flg_regeneration
515  s10_flag_nucleat = flg_nucl
516  s10_flag_icenucleat = flg_icenucl
517  s10_flag_sfaero = flg_sf_aero
518  s10_rndm_flgp = rndm_flgp
519  s10_rndm_mspc = mspc
520  s10_rndm_mbin = mbin
521 
522  rewind(io_fid_conf)
523  read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
524 
525  if ( ierr < 0 ) then !--- missing
526  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
527  elseif( ierr > 0 ) then !--- fatal error
528  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP, Check!'
529  call prc_mpistop
530  endif
531  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_mp)
532 
533  rewind(io_fid_conf)
534  read(io_fid_conf,nml=param_atmos_phy_mp_suzuki10,iostat=ierr)
535 
536  if ( ierr < 0 ) then !--- missing
537  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
538  elseif( ierr > 0 ) then !--- fatal error
539  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10, Check!'
540  call prc_mpistop
541  endif
542  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_mp_suzuki10)
543 
544  if ( nspc /= 1 .AND. nspc /= 7 ) then
545  write(*,*) 'xxx nspc should be set as 1(warm rain) or 7(mixed phase) check!'
546  call prc_mpistop
547  endif
548 
549  rhoa = rho_aero
550  emaer = s10_emaer
551  rasta = r_min
552  raend = r_max
553  r0a = r0_aero
554  flg_regeneration = s10_flag_regene
555  flg_nucl = s10_flag_nucleat
556  flg_icenucl = s10_flag_icenucleat
557  flg_sf_aero = s10_flag_sfaero
558  rndm_flgp = s10_rndm_flgp
559  mspc = s10_rndm_mspc
560  mbin = s10_rndm_mbin
561 
562  !--- read micpara.dat (microphysical parameter) and broad cast
563  if ( prc_ismaster ) then
564 
565  fid_micpara = io_get_available_fid()
566  !--- open parameter of cloud microphysics
567  open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )
568 
569  !--- micpara.dat does not exist
570  if ( ierr == 0 ) then
571 
572  read( fid_micpara,* ) nnspc, nnbin
573 
574  if ( nnbin /= nbin ) then
575  write(*,*) 'xxx nbin in inc_tracer and nbin in micpara.dat is different check!'
576  call prc_mpistop
577  endif
578 
579  ! grid parameter
580  if( io_l ) write(io_fid_log,*) '*** Radius of cloud ****'
581  do n = 1, nbin
582  read( fid_micpara,* ) nn, xctr( n ), radc( n )
583  if( io_l ) write(io_fid_log,'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
584  "Radius of ", n, "th cloud bin (bin center)= ", radc( n ) , "[m]"
585  enddo
586  do n = 1, nbin+1
587  read( fid_micpara,* ) nn, xbnd( n )
588  enddo
589  read( fid_micpara,* ) dxmic
590  if( io_l ) write(io_fid_log,*) '*** Width of Cloud SDF= ', dxmic
591 
592  ! capacity
593  do myu = 1, nspc_mk
594  do n = 1, nbin
595 ! read( fid_micpara,* ) mmyu, nn, cctr( myu,n )
596  read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
597  enddo
598  do n = 1, nbin+1
599 ! read( fid_micpara,* ) mmyu, nn, cbnd( myu,n )
600  read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
601  enddo
602  enddo
603 
604  ! collection kernel
605  do myu = 1, nspc_mk
606  do nyu = 1, nspc_mk
607  do i = 1, nbin
608  do j = 1, nbin
609  read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
610  enddo
611  enddo
612  enddo
613  enddo
614 
615  ! terminal velocity
616  do myu = 1, nspc_mk
617  do n = 1, nbin
618  read( fid_micpara,* ) mmyu, nn, vt( myu,n )
619  enddo
620  enddo
621 
622  ! bulk density
623  do myu = 1, nspc_mk
624  do n = 1, nbin
625  read( fid_micpara,* ) mmyu, nn, br( myu,n )
626  enddo
627  enddo
628 
629  close ( fid_micpara )
630 
631  !--- micpara.dat does not exist
632  else
633 
634  if( io_l ) write(io_fid_log,*) 'micpara.dat is created'
635  call mkpara
636 
637  fid_micpara = io_get_available_fid()
638  !--- open parameter of cloud microphysics
639  open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )
640 
641  read( fid_micpara,* ) nnspc, nnbin
642 
643  if ( nnbin /= nbin ) then
644  write(*,*) 'xxx nbin in inc_tracer and nbin in micpara.dat is different check!'
645  call prc_mpistop
646  endif
647 
648  ! grid parameter
649  if( io_l ) write(io_fid_log,*) '*** Radius of cloud ****'
650  do n = 1, nbin
651  read( fid_micpara,* ) nn, xctr( n ), radc( n )
652  if( io_l ) write(io_fid_log,'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
653  "Radius of ", n, "th cloud bin (bin center)= ", radc( n ) , "[m]"
654  enddo
655  do n = 1, nbin+1
656  read( fid_micpara,* ) nn, xbnd( n )
657  enddo
658  read( fid_micpara,* ) dxmic
659  if( io_l ) write(io_fid_log,*) '*** Width of Cloud SDF= ', dxmic
660 
661  ! capacity
662  do myu = 1, nspc_mk
663  do n = 1, nbin
664 ! read( fid_micpara,* ) mmyu, nn, cctr( myu,n )
665  read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
666  enddo
667  do n = 1, nbin+1
668 ! read( fid_micpara,* ) mmyu, nn, cbnd( myu,n )
669  read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
670  enddo
671  enddo
672 
673  ! collection kernel
674  do myu = 1, nspc_mk
675  do nyu = 1, nspc_mk
676  do i = 1, nbin
677  do j = 1, nbin
678  read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
679  enddo
680  enddo
681  enddo
682  enddo
683 
684  ! terminal velocity
685  do myu = 1, nspc_mk
686  do n = 1, nbin
687  read( fid_micpara,* ) mmyu, nn, vt( myu,n )
688  enddo
689  enddo
690 
691  ! bulk density
692  do myu = 1, nspc_mk
693  do n = 1, nbin
694  read( fid_micpara,* ) mmyu, nn, br( myu,n )
695  enddo
696  enddo
697 
698  close ( fid_micpara )
699 
700  endif
701 
702  endif
703 
704  call mpi_bcast( radc, nbin, comm_datatype, prc_masterrank, comm_world, ierr )
705  call mpi_bcast( xctr, nbin, comm_datatype, prc_masterrank, comm_world, ierr )
706  call mpi_bcast( dxmic,1, comm_datatype, prc_masterrank, comm_world, ierr )
707  call mpi_bcast( xbnd, nbin+1, comm_datatype, prc_masterrank, comm_world, ierr )
708  call mpi_bcast( cctr, nbin*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
709  call mpi_bcast( cbnd, (nbin+1)*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
710  call mpi_bcast( ck, nspc_mk*nspc_mk*nbin*nbin, comm_datatype, prc_masterrank, comm_world, ierr )
711  call mpi_bcast( br, nbin*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
712  call mpi_bcast( vt, nbin*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
713 
714  !--- aerosol ( CCN ) (not supported)
715  if ( nccn /= 0 ) then
716 
717  allocate ( ncld( 1:nccn ) )
718  xasta = log( rhoa*4.0_rp/3.0_rp*pi * ( rasta )**3 )
719  xaend = log( rhoa*4.0_rp/3.0_rp*pi * ( raend )**3 )
720 
721  dxaer = ( xaend-xasta )/nccn
722 
723  do n = 1, nccn+1
724  xabnd( n ) = xasta + dxaer*( n-1 )
725  enddo
726  do n = 1, nccn
727  xactr( n ) = ( xabnd( n )+xabnd( n+1 ) )*0.50_rp
728  rada( n ) = ( exp( xactr( n ) )*thirdovforth/pi/rhoa )**( oneovthird )
729  if( io_l ) write(io_fid_log,'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
730  "Radius of ", n, "th aerosol bin (bin center)= ", rada( n ) , "[m]"
731  enddo
732 
733  if ( flg_sf_aero ) then
734  write(*,*) "flg_sf_aero=true is not supported stop!! "
735  call prc_mpistop
736 ! if ( CZ(KS) >= 10.0_RP ) then
737 ! R10M1 = 10.0_RP / CZ(KS) * 0.50_RP ! scale with height
738 ! R10M2 = 10.0_RP / CZ(KS) * 0.50_RP ! scale with height
739 ! R10H1 = 1.0_RP * 0.50_RP
740 ! R10H2 = 1.0_RP * 0.0_RP
741 ! R10E1 = 1.0_RP * 0.50_RP
742 ! R10E2 = 1.0_RP * 0.50_RP
743 ! K10_1 = KS
744 ! K10_2 = KS
745 ! else
746 ! k = 1
747 ! do while ( CZ(k) < 10.0_RP )
748 ! k = k + 1
749 ! K10_1 = k
750 ! K10_2 = k + 1
751 ! R10M1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
752 ! R10M2 = ( 10.0_RP - CZ(k) ) / CDZ(k)
753 ! R10H1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
754 ! R10H2 = ( 10.0_RP - CZ(k) ) / CDZ(k)
755 ! R10E1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
756 ! R10E2 = ( 10.0_RP - CZ(k) ) / CDZ(k)
757 ! enddo
758 ! endif
759  endif
760 
761  endif
762 
763  atmos_phy_mp_suzuki10_dens(i_hc) = const_dwatr
764  atmos_phy_mp_suzuki10_dens(i_hr) = const_dice
765  atmos_phy_mp_suzuki10_dens(i_hi) = const_dice
766  atmos_phy_mp_suzuki10_dens(i_hs) = const_dice
767  atmos_phy_mp_suzuki10_dens(i_hg) = const_dice
768  atmos_phy_mp_suzuki10_dens(i_hh) = const_dice
769 
770  !--- determine nbnd
771  do n = 1, nbin
772  if( radc( n ) > rbnd ) then
773  nbnd = n
774  exit
775  endif
776  enddo
777  if( io_l ) write(io_fid_log,'(A,ES15.7,A)') '*** Radius between cloud and rain is ', radc(nbnd), '[m]'
778 
779  !--- random number setup for stochastic method
780  if ( rndm_flgp > 0 ) then
781  call random_setup( ia*ja*ka )
782  endif
783 
784  if ( mp_couple_aerosol .AND. nccn /=0 ) then
785  write(*,*) 'xxx nccn should be 0 when MP_couple_aerosol = .true. !! stop'
786  call prc_mpistop
787  endif
788 
789  if ( nccn /= 0 ) then
790  do n = 1, nccn
791  expxactr( n ) = exp( xactr( n ) )
792  rexpxactr( n ) = 1.0_rp / exp( xactr( n ) )
793  enddo
794  do n = 1, nccn+1
795  expxabnd( n ) = exp( xabnd( n ) )
796  rexpxabnd( n ) = 1.0_rp / exp( xabnd( n ) )
797  enddo
798  endif
799 
800  allocate( vterm(ka,ia,ja,qa_mp-1) )
801  vterm(:,:,:,:) = 0.0_rp
802  do myu = 1, nspc
803  do n = 1, nbin
804  vterm(:,:,:,(myu-1)*nbin+n) = -vt( myu,n )
805  enddo
806  enddo
807  do n = 1, nbin
808  expxctr( n ) = exp( xctr( n ) )
809  rexpxctr( n ) = 1.0_rp / exp( xctr( n ) )
810  enddo
811  do n = 1, nbin+1
812  expxbnd( n ) = exp( xbnd( n ) )
813  rexpxbnd( n ) = 1.0_rp / exp( xbnd( n ) )
814  enddo
815 
816  allocate( kindx(nbin,nbin) )
817  call getrule( ifrsl,kindx )
818 
819  if ( const_thermodyn_type == 'EXACT' ) then
820  flg_thermodyn = 1.0_rp
821  elseif( const_thermodyn_type == 'SIMPLE' ) then
822  flg_thermodyn = 0.0_rp
823  endif
824  rtem00 = 1.0_rp / const_tem00
825 
826  nstep_max = int( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz ) )
827  mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
828 
829  mp_nstep_sedimentation = mp_ntmax_sedimentation
830  mp_rnstep_sedimentation = 1.0_rp / real(mp_ntmax_sedimentation,kind=rp)
831  mp_dtsec_sedimentation = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation
832 
833  if( io_l ) write(io_fid_log,*)
834  if( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, ' step'
835  if( io_l ) write(io_fid_log,*) '*** DT of sedimentation is : ', mp_dtsec_sedimentation, '[s]'
836 
837  return
838 
839  end subroutine atmos_phy_mp_suzuki10_setup
840 
841  !-----------------------------------------------------------------------------
843  subroutine atmos_phy_mp_suzuki10( &
844  DENS, &
845  MOMZ, &
846  MOMX, &
847  MOMY, &
848  RHOT, &
849  QTRC, &
850  CCN, &
851  EVAPORATE, &
852  SFLX_rain, &
853  SFLX_snow )
854  use scale_grid_index
855  use scale_time, only: &
856  dt => time_dtsec_atmos_phy_mp
857  use scale_tracer, only: &
858  qa, &
859  tracer_r, &
860  tracer_cv, &
861  tracer_mass
862  use scale_atmos_thermodyn, only: &
863  thermodyn_rhoe => atmos_thermodyn_rhoe, &
864  thermodyn_rhot => atmos_thermodyn_rhot, &
865  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
866  use scale_atmos_phy_mp_common, only: &
867  mp_negative_fixer => atmos_phy_mp_negative_fixer, &
868  mp_precipitation => atmos_phy_mp_precipitation
869  use scale_history, only: &
870  hist_in
871  implicit none
872 
873  real(RP), intent(inout) :: dens (ka,ia,ja)
874  real(RP), intent(inout) :: momz (ka,ia,ja)
875  real(RP), intent(inout) :: momx (ka,ia,ja)
876  real(RP), intent(inout) :: momy (ka,ia,ja)
877  real(RP), intent(inout) :: rhot (ka,ia,ja)
878  real(RP), intent(inout) :: qtrc (ka,ia,ja,qa)
879  real(RP), intent(in) :: ccn (ka,ia,ja)
880  real(RP), intent(out) :: evaporate(ka,ia,ja) !--- number of evaporated cloud [/m3]
881  real(RP), intent(out) :: sflx_rain(ia,ja)
882  real(RP), intent(out) :: sflx_snow(ia,ja)
883 
884  real(RP) :: rhoe (ka,ia,ja)
885  real(RP) :: pott (ka,ia,ja)
886  real(RP) :: temp (ka,ia,ja)
887  real(RP) :: pres (ka,ia,ja)
888  real(RP) :: qsat (ka,ia,ja)
889  real(RP) :: ssliq(ka,ia,ja)
890  real(RP) :: ssice(ka,ia,ja)
891 
892  integer :: ijk_index (kijmax,3)
893  integer :: index_cld (kijmax)
894  integer :: index_cold(kijmax)
895  integer :: index_warm(kijmax)
896  integer :: ijkcount, ijkcount_cold, ijkcount_warm
897  integer :: ijk, indirect
898 
899  real(RP) :: dens_ijk(kijmax)
900  real(RP) :: pres_ijk(kijmax)
901  real(RP) :: temp_ijk(kijmax)
902  real(RP) :: qvap_ijk(kijmax)
903  real(RP) :: ccn_ijk(kijmax)
904  real(RP) :: evaporate_ijk(kijmax)
905  real(RP) :: ghyd_ijk(nbin,nspc,kijmax)
906  real(RP) :: gaer_ijk(nccn1 ,kijmax)
907  real(RP) :: cldsum
908  integer :: countbin
909 
910 ! logical, save :: ofirst_sdfa = .true.
911 ! real(RP) :: VELX (IA,JA)
912 ! real(RP) :: VELY (IA,JA)
913 ! real(RP) :: SFLX_AERO(IA,JA,nccn)
914 ! real(RP) :: Uabs, bparam
915 ! real(RP) :: AMR(KA,IA,JA)
916 
917  real(RP) :: pflux (ka,ia,ja,qa_mp-1) ! precipitation flux of each tracer [kg/m2/s]
918  real(RP) :: flx_hydro(ka,ia,ja,qa_mp-1)
919  real(RP) :: qhyd_out (ka,ia,ja,6)
920 
921  integer :: step
922  integer :: k, i, j, m, n, iq
923  !---------------------------------------------------------------------------
924 
925  if ( nspc == 1 ) then
926  if( io_l ) write(io_fid_log,*) '*** Atmos physics step: Cloud microphysics(SBM Liquid water only)'
927  elseif( nspc > 1 ) then
928  if( io_l ) write(io_fid_log,*) '*** Atmos physics step: Cloud microphysics(SBM Mixed phase)'
929  endif
930 
931  if ( mp_donegative_fixer ) then
932  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
933  rhot(:,:,:), & ! [INOUT]
934  qtrc(:,:,:,:), & ! [INOUT]
935  i_qv, & ! [IN]
936  mp_limit_negative ) ! [IN]
937  endif
938 
939  call prof_rapstart('MP_ijkconvert', 3)
940 
941  ! Clear EVAPORATE
942  do j = js, je
943  do i = is, ie
944  do k = ks, ke
945  evaporate(k,i,j) = 0.0_rp
946  enddo
947  enddo
948  enddo
949 
950  ijk = 0
951  do j = js, je
952  do i = is, ie
953  do k = ks, ke
954  ijk = ijk + 1
955  ijk_index(ijk,1) = i
956  ijk_index(ijk,2) = j
957  ijk_index(ijk,3) = k
958  enddo
959  enddo
960  enddo
961 
962  call thermodyn_temp_pres( temp(:,:,:), & ! [OUT]
963  pres(:,:,:), & ! [OUT]
964  dens(:,:,:), & ! [IN]
965  rhot(:,:,:), & ! [IN]
966  qtrc(:,:,:,:), & ! [IN]
967  tracer_cv(:), & ! [IN]
968  tracer_r(:), & ! [IN]
969  tracer_mass(:) ) ! [IN]
970 
971  call saturation_pres2qsat_liq( qsat(:,:,:), & ! [OUT]
972  temp(:,:,:), & ! [IN]
973  pres(:,:,:) ) ! [IN]
974 
975  do j = js, je
976  do i = is, ie
977  do k = ks, ke
978  ssliq(k,i,j) = qtrc(k,i,j,i_qv) / qsat(k,i,j) - 1.0_rp
979  enddo
980  enddo
981  enddo
982 
983  if ( nspc == 1 ) then
984  ssice(:,:,:) = 0.0_rp
985  else
986  call saturation_pres2qsat_ice( qsat(:,:,:), & ! [OUT]
987  temp(:,:,:), & ! [IN]
988  pres(:,:,:) ) ! [IN]
989 
990  do j = js, je
991  do i = is, ie
992  do k = ks, ke
993  ssice(k,i,j) = qtrc(k,i,j,i_qv) / qsat(k,i,j) - 1.0_rp
994  enddo
995  enddo
996  enddo
997  endif
998 
999 !--- store initial SDF of aerosol
1000 !--- this option is not supported
1001 ! if ( ofirst_sdfa ) then
1002 ! allocate( marate( nccn ) )
1003 ! do j = JS, JE
1004 ! do i = IS, IE
1005 ! do k = KS, KE
1006 ! sum2 = 0.0_RP
1007 ! do n = 1, nccn
1008 ! marate( n ) = gdga(k,i,j,n)*rexpxactr( n )
1009 ! sum2 = sum2 + gdga(k,i,j,n)*rexpxactr( n )
1010 ! enddo
1011 ! enddo
1012 ! enddo
1013 ! enddo
1014 ! if ( sum2 /= 0.0_RP ) then
1015 ! marate( 1:nccn ) = marate( 1:nccn )/sum2
1016 ! ofirst_sdfa = .false.
1017 ! endif
1018 ! endif
1019 
1020  !--- Arrange array for microphysics
1021  ijkcount = 0
1022  ijkcount_cold = 0
1023  ijkcount_warm = 0
1024 
1025  ijk = 0
1026  do j = js, je
1027  do i = is, ie
1028  do k = ks, ke
1029  ijk = ijk + 1
1030 
1031  cldsum = 0.0_rp
1032  countbin = qs_mp + 1
1033  do m = 1, nspc
1034  do n = 1, nbin
1035  cldsum = cldsum + qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
1036  countbin = countbin + 1
1037  enddo
1038  enddo
1039 
1040  if ( cldsum > cldmin &
1041  .OR. ssliq(k,i,j) > 0.0_rp &
1042  .OR. ssice(k,i,j) > 0.0_rp ) then
1043 
1044  ijkcount = ijkcount + 1
1045 
1046  index_cld(ijkcount) = ijk
1047 
1048  dens_ijk(ijkcount) = dens(k,i,j)
1049  pres_ijk(ijkcount) = pres(k,i,j)
1050  temp_ijk(ijkcount) = temp(k,i,j)
1051  ccn_ijk(ijkcount) = ccn(k,i,j)
1052  qvap_ijk(ijkcount) = qtrc(k,i,j,i_qv)
1053 
1054  countbin = qs_mp + 1
1055  do m = 1, nspc
1056  do n = 1, nbin
1057  ghyd_ijk(n,m,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
1058  countbin = countbin + 1
1059  enddo
1060  enddo
1061 
1062  do n = 1, nccn
1063  gaer_ijk(n,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxaer
1064  countbin = countbin + 1
1065  enddo
1066 
1067  if ( temp(k,i,j) < const_tem00 .AND. nspc > 1 ) then ! cold
1068  ijkcount_cold = ijkcount_cold + 1
1069  index_cold(ijkcount_cold) = ijkcount
1070  else ! warm
1071  ijkcount_warm = ijkcount_warm + 1
1072  index_warm(ijkcount_warm) = ijkcount
1073  endif
1074  endif
1075 
1076  enddo
1077  enddo
1078  enddo
1079 
1080  call prof_rapend ('MP_ijkconvert', 3)
1081 
1082  ! tentative timername registration
1083  call prof_rapstart('MP_suzuki10', 3)
1084  call prof_rapend ('MP_suzuki10', 3)
1085  call prof_rapstart('_SBM_Nucleat', 3)
1086  call prof_rapend ('_SBM_Nucleat', 3)
1087 ! call PROF_rapstart('_SBM_NucleatA', 3)
1088 ! call PROF_rapend ('_SBM_NucleatA', 3)
1089  call prof_rapstart('_SBM_Liqphase', 3)
1090  call prof_rapend ('_SBM_Liqphase', 3)
1091  call prof_rapstart('_SBM_Icephase', 3)
1092  call prof_rapend ('_SBM_Icephase', 3)
1093  call prof_rapstart('_SBM_Mixphase', 3)
1094  call prof_rapend ('_SBM_Mixphase', 3)
1095  call prof_rapstart('_SBM_AdvLiq', 3)
1096  call prof_rapend ('_SBM_AdvLiq', 3)
1097  call prof_rapstart('_SBM_AdvIce', 3)
1098  call prof_rapend ('_SBM_AdvIce', 3)
1099  call prof_rapstart('_SBM_AdvMix', 3)
1100  call prof_rapend ('_SBM_AdvMix', 3)
1101 ! call PROF_rapstart('_SBM_FAero', 3)
1102 ! call PROF_rapend ('_SBM_FAero', 3)
1103  call prof_rapstart('_SBM_Freezing', 3)
1104  call prof_rapend ('_SBM_Freezing', 3)
1105  call prof_rapstart('_SBM_IceNucleat', 3)
1106  call prof_rapend ('_SBM_IceNucleat', 3)
1107  call prof_rapstart('_SBM_Melting', 3)
1108  call prof_rapend ('_SBM_Melting', 3)
1109  call prof_rapstart('_SBM_CollCoag', 3)
1110  call prof_rapend ('_SBM_CollCoag', 3)
1111 ! call PROF_rapstart('_SBM_CollCoagR', 3)
1112 ! call PROF_rapend ('_SBM_CollCoagR', 3)
1113 
1114  if ( ijkcount > 0 ) then
1115 
1116  call prof_rapstart('MP_suzuki10', 3)
1117 
1118  call mp_suzuki10( ijkcount, & ! [IN]
1119  ijkcount_cold, & ! [IN]
1120  ijkcount_warm, & ! [IN]
1121  index_cold( 1:ijkcount), & ! [IN]
1122  index_warm( 1:ijkcount), & ! [IN]
1123  dens_ijk( 1:ijkcount), & ! [IN]
1124  pres_ijk( 1:ijkcount), & ! [IN]
1125  ccn_ijk( 1:ijkcount), & ! [IN]
1126  temp_ijk( 1:ijkcount), & ! [INOUT]
1127  qvap_ijk( 1:ijkcount), & ! [INOUT]
1128  ghyd_ijk(:,:,1:ijkcount), & ! [INOUT]
1129  gaer_ijk(:, 1:ijkcount), & ! [INOUT]
1130  evaporate_ijk( 1:ijkcount), & ! [OUT]
1131  dt ) ! [IN]
1132 
1133  call prof_rapend ('MP_suzuki10', 3)
1134 
1135 ! if ( flg_sf_aero ) then
1136 ! do j = JS-2, JE+2
1137 ! do i = IS-2, IE+1
1138 ! VELX(i,j) = MOMX(K10_1,i,j) / ( DENS(K10_1,i+1,j)+DENS(K10_1,i,j) ) * R10M1 &
1139 ! + MOMX(K10_2,i,j) / ( DENS(K10_2,i+1,j)+DENS(K10_2,i,j) ) * R10M2
1140 ! enddo
1141 ! enddo
1142 !
1143 ! do j = JS-2, JE+1
1144 ! do i = IS-2, IE+2
1145 ! VELY(i,j) = MOMY(K10_1,i,j) / ( DENS(K10_1,i,j+1)+DENS(K10_1,i,j) ) * R10M1 &
1146 ! + MOMY(K10_2,i,j) / ( DENS(K10_2,i,j+1)+DENS(K10_2,i,j) ) * R10M2
1147 ! enddo
1148 ! enddo
1149 ! endif
1150 !
1151 ! !--- SURFACE FLUX by Monahan et al. (1986)
1152 ! if ( flg_sf_aero .AND. nccn /= 0 ) then
1153 ! do j = JS, JE
1154 ! do i = IS, IE
1155 ! ijk = ( j - JS ) * KMAX * IMAX &
1156 ! + ( i - IS ) * KMAX &
1157 ! + ( KS - KS ) &
1158 ! + 1
1159 ! Uabs = sqrt( ( ( VELX(i,j) + VELX(i-1,j ) ) * 0.50_RP )**2 &
1160 ! + ( ( VELY(i,j) + VELY(i ,j-1) ) * 0.50_RP )**2 )
1161 ! do n = 1, nccn
1162 ! if ( rada( n ) <= 2.0E-5_RP .AND. rada( n ) >= 3.0E-7_RP ) then
1163 ! bparam = ( 0.38_RP - log( rada( n ) ) )/0.65_RP
1164 ! SFLX_AERO(i,j,n) = 1.373_RP * Uabs**( 3.41_RP ) * rada( n )**( -3.0_RP ) &
1165 ! * ( 1.0_RP + 0.057_RP * rada( n )**( 1.05_RP ) ) &
1166 ! * 10.0_RP**( 1.19_RP * exp( -bparam*bparam ) )
1167 ! ! convert from [#/m^2/um/s] -> [kg/m^3/unit log (m)]
1168 ! SFLX_AERO(i,j,n) = SFLX_AERO(i,j,n) / DENS(KS,i,j) &
1169 ! / GRID_CDZ(KS) * rada( n ) / 3.0_RP * dt * expxactr( n )
1170 ! Gaer_ijk(n,ijk) = Gaer_ijk(n,ijk) + SFLX_AERO(i,j,n)/dxaer
1171 ! endif
1172 ! enddo
1173 ! enddo
1174 ! enddo
1175 ! endif
1176 
1177  call prof_rapstart('MP_ijkconvert', 3)
1178 
1179  !---- return original array
1180  do ijk = 1, ijkcount
1181  indirect = index_cld(ijk)
1182  i = ijk_index(indirect,1)
1183  j = ijk_index(indirect,2)
1184  k = ijk_index(indirect,3)
1185 
1186  temp(k,i,j) = temp_ijk(ijk)
1187  qtrc(k,i,j,i_qv) = qvap_ijk(ijk)
1188  evaporate(k,i,j) = evaporate_ijk(ijk) / dt ! [#/m3/s]
1189 
1190  countbin = qs_mp + 1
1191  do m = 1, nspc
1192  do n = 1, nbin
1193  qtrc(k,i,j,countbin) = ghyd_ijk(n,m,ijk) / dens(k,i,j) * dxmic
1194  countbin = countbin + 1
1195  enddo
1196  enddo
1197 
1198  do n = 1, nccn
1199  qtrc(k,i,j,countbin) = gaer_ijk(n,ijk) / dens(k,i,j) * dxaer
1200  countbin = countbin + 1
1201  enddo
1202  enddo
1203 
1204  do j = js, je
1205  do i = is, ie
1206  do k = ks, ke
1207  countbin = qs_mp + 1
1208  do m = 1, nspc
1209  do n = 1, nbin
1210  if ( qtrc(k,i,j,countbin) < eps ) then
1211  qtrc(k,i,j,countbin) = 0.0_rp
1212  endif
1213  countbin = countbin + 1
1214  enddo
1215  enddo
1216  enddo
1217  enddo
1218  enddo
1219 
1220  call thermodyn_pott( pott(:,:,:), & ! [OUT]
1221  temp(:,:,:), & ! [IN]
1222  pres(:,:,:), & ! [IN]
1223  qtrc(:,:,:,:), & ! [IN]
1224  tracer_cv(:), & ! [IN]
1225  tracer_r(:), & ! [IN]
1226  tracer_mass(:) ) ! [IN]
1227 
1228  do j = js, je
1229  do i = is, ie
1230  do k = ks, ke
1231  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
1232  enddo
1233  enddo
1234  enddo
1235 
1236 ! if ( nccn /= 0 ) then
1237 ! AMR(:,:,:) = 0.0_RP
1238 ! do j = JS, JE
1239 ! do i = IS, IE
1240 ! do k = KS, KE
1241 ! do n = 1, nccn
1242 ! AMR(k,i,j) = AMR(k,i,j) + QTRC(k,i,j,QQE-1+n)
1243 ! enddo
1244 ! enddo
1245 ! enddo
1246 ! enddo
1247 ! endif
1248 
1249  call prof_rapend ('MP_ijkconvert', 3)
1250 
1251  endif
1252 
1253  flx_hydro(:,:,:,:) = 0.0_rp
1254 
1255  if ( mp_doprecipitation ) then
1256 
1257  call thermodyn_rhoe( rhoe(:,:,:), & ! [OUT]
1258  rhot(:,:,:), & ! [IN]
1259  qtrc(:,:,:,:), & ! [IN]
1260  tracer_cv(:), & ! [IN]
1261  tracer_r(:), & ! [IN]
1262  tracer_mass(:) ) ! [IN]
1263 
1264  do step = 1, mp_nstep_sedimentation
1265 
1266  call thermodyn_temp_pres_e( temp(:,:,:), & ! [OUT]
1267  pres(:,:,:), & ! [OUT]
1268  dens(:,:,:), & ! [IN]
1269  rhoe(:,:,:), & ! [IN]
1270  qtrc(:,:,:,:), & ! [IN]
1271  tracer_cv(:), & ! [IN]
1272  tracer_r(:), & ! [IN]
1273  tracer_mass(:) ) ! [IN]
1274 
1275  call mp_precipitation( qa_mp, & ! [IN]
1276  qs_mp, & ! [IN]
1277  pflux(:,:,:,:), & ! [OUT]
1278  vterm(:,:,:,:), & ! [INOUT]
1279  dens(:,:,:), & ! [INOUT]
1280  momz(:,:,:), & ! [INOUT]
1281  momx(:,:,:), & ! [INOUT]
1282  momy(:,:,:), & ! [INOUT]
1283  rhoe(:,:,:), & ! [INOUT]
1284  qtrc(:,:,:,:), & ! [INOUT]
1285  temp(:,:,:), & ! [IN]
1286  tracer_cv(:), & ! [IN]
1287  mp_dtsec_sedimentation ) ! [IN]
1288 
1289  do iq = 1, qa_mp-1
1290  do j = js, je
1291  do i = is, ie
1292  do k = ks-1, ke-1
1293  flx_hydro(k,i,j,iq) = flx_hydro(k,i,j,iq) + pflux(k,i,j,iq) * mp_rnstep_sedimentation
1294  enddo
1295  enddo
1296  enddo
1297  enddo
1298 
1299  enddo
1300 
1301  call thermodyn_rhot( rhot(:,:,:), & ! [OUT]
1302  rhoe(:,:,:), & ! [IN]
1303  qtrc(:,:,:,:), & ! [IN]
1304  tracer_cv(:), & ! [IN]
1305  tracer_r(:), & ! [IN]
1306  tracer_mass(:) ) ! [IN]
1307  endif
1308 
1309  sflx_rain(:,:) = 0.0_rp
1310  sflx_snow(:,:) = 0.0_rp
1311 
1312  !--- lowermost flux is saved for land process
1313  do n = 1, nbin
1314  iq = n
1315 
1316  do j = js, je
1317  do i = is, ie
1318  sflx_rain(i,j) = sflx_rain(i,j) - flx_hydro(ks-1,i,j,iq)
1319  enddo
1320  enddo
1321  enddo
1322 
1323  if ( nspc > 1 ) then
1324  do m = ic, ih
1325  do n = 1, nbin
1326  iq = (m-1)*nbin + n
1327 
1328  do j = js, je
1329  do i = is, ie
1330  sflx_snow(i,j) = sflx_snow(i,j) - flx_hydro(ks-1,i,j,iq)
1331  enddo
1332  enddo
1333  enddo
1334  enddo
1335  endif
1336 
1337  !##### END MP Main #####
1338 
1339  if ( mp_donegative_fixer ) then
1340  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
1341  rhot(:,:,:), & ! [INOUT]
1342  qtrc(:,:,:,:), & ! [INOUT]
1343  i_qv, & ! [IN]
1344  mp_limit_negative ) ! [IN]
1345  endif
1346 
1347  qhyd_out(:,:,:,:) = 0.0_rp
1348 
1349  do n = 1, nbnd
1350  iq = qs_mp + n
1351 
1352  do j = js, je
1353  do i = is, ie
1354  do k = ks, ke
1355  qhyd_out(k,i,j,1) = qhyd_out(k,i,j,1) + qtrc(k,i,j,iq)
1356  enddo
1357  enddo
1358  enddo
1359  enddo
1360 
1361  do n = nbnd+1, nbin
1362  iq = qs_mp + n
1363 
1364  do j = js, je
1365  do i = is, ie
1366  do k = ks, ke
1367  qhyd_out(k,i,j,2) = qhyd_out(k,i,j,2) + qtrc(k,i,j,iq)
1368  enddo
1369  enddo
1370  enddo
1371  enddo
1372 
1373  call hist_in( qhyd_out(:,:,:,1), 'QC', 'Mixing ratio of QC', 'kg/kg' )
1374  call hist_in( qhyd_out(:,:,:,2), 'QR', 'Mixing ratio of QR', 'kg/kg' )
1375 
1376  if ( nspc > 1 ) then
1377  do m = ic, id ! columnar,plate,dendrite = ice
1378  do n = 1, nbin
1379  iq = qs_mp + (m-1)*nbin + n
1380 
1381  do j = js, je
1382  do i = is, ie
1383  do k = ks, ke
1384  qhyd_out(k,i,j,3) = qhyd_out(k,i,j,3) + qtrc(k,i,j,iq)
1385  enddo
1386  enddo
1387  enddo
1388  enddo
1389  enddo
1390 
1391  do n = 1, nbin
1392  iq = qs_mp + (iss-1)*nbin + n
1393 
1394  do j = js, je
1395  do i = is, ie
1396  do k = ks, ke
1397  qhyd_out(k,i,j,4) = qhyd_out(k,i,j,4) + qtrc(k,i,j,iq)
1398  enddo
1399  enddo
1400  enddo
1401  enddo
1402 
1403  do n = 1, nbin
1404  iq = qs_mp + (ig-1)*nbin + n
1405 
1406  do j = js, je
1407  do i = is, ie
1408  do k = ks, ke
1409  qhyd_out(k,i,j,5) = qhyd_out(k,i,j,5) + qtrc(k,i,j,iq)
1410  enddo
1411  enddo
1412  enddo
1413  enddo
1414 
1415  do n = 1, nbin
1416  iq = qs_mp + (ih-1)*nbin + n
1417 
1418  do j = js, je
1419  do i = is, ie
1420  do k = ks, ke
1421  qhyd_out(k,i,j,6) = qhyd_out(k,i,j,6) + qtrc(k,i,j,iq)
1422  enddo
1423  enddo
1424  enddo
1425  enddo
1426 
1427  call hist_in( qhyd_out(:,:,:,3), 'QI', 'Mixing ratio of QI', 'kg/kg' )
1428  call hist_in( qhyd_out(:,:,:,4), 'QS', 'Mixing ratio of QS', 'kg/kg' )
1429  call hist_in( qhyd_out(:,:,:,5), 'QG', 'Mixing ratio of QG', 'kg/kg' )
1430  call hist_in( qhyd_out(:,:,:,6), 'QH', 'Mixing ratio of QH', 'kg/kg' )
1431  endif
1432 
1433  return
1434  end subroutine atmos_phy_mp_suzuki10
1435 
1436  !-----------------------------------------------------------------------------
1437  subroutine mp_suzuki10( &
1438  ijkmax, &
1439  ijkmax_cold, &
1440  ijkmax_warm, &
1441  index_cold, &
1442  index_warm, &
1443  dens, &
1444  pres, &
1445  ccn, &
1446  temp, &
1447  qvap, &
1448  ghyd, &
1449  gaer, &
1450  evaporate, &
1451  dt )
1452  implicit none
1453 
1454  integer, intent(in) :: ijkmax
1455  integer, intent(in) :: ijkmax_cold
1456  integer, intent(in) :: ijkmax_warm
1457  integer , intent(in) :: index_cold(ijkmax)
1458  integer , intent(in) :: index_warm(ijkmax)
1459  real(RP), intent(in) :: dens (ijkmax) ! Density [kg/m3]
1460  real(RP), intent(in) :: pres (ijkmax) ! Pressure [Pa]
1461  real(RP), intent(in) :: ccn (ijkmax) ! Number concentration of CCN [#/m3]
1462  real(RP), intent(inout) :: temp (ijkmax) ! Temperature [K]
1463  real(RP), intent(inout) :: qvap (ijkmax) ! Specific humidity [kg/kg]
1464  real(RP), intent(inout) :: ghyd (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
1465  real(RP), intent(inout) :: gaer (nccn1, ijkmax) ! Mass size distribution function of aerosol
1466  real(RP), intent(out) :: evaporate (ijkmax) ! Number concentration of evaporated cloud [/m3/s]
1467  real(DP), intent(in) :: dt ! Time step interval
1468  !---------------------------------------------------------------------------
1469 
1470  if ( nccn /= 0 ) then
1471  if ( nspc == 1 ) then
1472  !---< warm rain only with aerosol tracer >---
1473 
1474  ! nucleation from aerosol
1475  call nucleata( ijkmax, & ! [IN]
1476  dens(:), & ! [IN]
1477  pres(:), & ! [IN]
1478  temp(:), & ! [INOUT]
1479  qvap(:), & ! [INOUT]
1480  ghyd(:,:,:), & ! [INOUT]
1481  gaer(:,:), & ! [INOUT]
1482  dt ) ! [IN]
1483 
1484  ! condensation / evaporation
1485  call cndevpsbla( ijkmax, & ! [IN]
1486  dens(:), & ! [IN]
1487  pres(:), & ! [IN]
1488  temp(:), & ! [INOUT]
1489  qvap(:), & ! [INOUT]
1490  ghyd(:,:,:), & ! [INOUT]
1491  gaer(:,:), & ! [INOUT]
1492  evaporate(:),& ! [OUT]
1493  dt ) ! [IN]
1494 
1495  if ( mp_doautoconversion ) then
1496  ! collision-coagulation
1497  call collmain( ijkmax, & ! [IN]
1498  temp(:), & ! [IN]
1499  ghyd(:,:,:), & ! [INOUT]
1500  dt ) ! [IN]
1501  endif
1502 
1503  elseif( nspc > 1 ) then
1504  !---< mixed phase rain only with aerosol tracer >---
1505 
1506  ! nucleation from aerosol
1507  call nucleata( ijkmax, & ! [IN]
1508  dens(:), & ! [IN]
1509  pres(:), & ! [IN]
1510  temp(:), & ! [INOUT]
1511  qvap(:), & ! [INOUT]
1512  ghyd(:,:,:), & ! [INOUT]
1513  gaer(:,:), & ! [INOUT]
1514  dt ) ! [IN]
1515 
1516  ! freezing / melting
1517  call freezing( ijkmax, & ! [IN]
1518  ijkmax_cold, & ! [IN]
1519  index_cold(:), & ! [IN]
1520  dens(:), & ! [IN]
1521  temp(:), & ! [INOUT]
1522  ghyd(:,:,:), & ! [INOUT]
1523  dt ) ! [IN]
1524 
1525 ! call ice_nucleat( ijkmax, & ! [IN]
1526 ! ijkmax_cold, & ! [IN]
1527 ! index_cold(:), & ! [IN]
1528 ! dens(:), & ! [IN]
1529 ! pres(:), & ! [IN]
1530 ! temp(:), & ! [INOUT]
1531 ! qvap(:), & ! [INOUT]
1532 ! ghyd(:,:,:), & ! [INOUT]
1533 ! dt ) ! [IN]
1534 
1535  call melting( ijkmax, & ! [IN]
1536  ijkmax_warm, & ! [IN]
1537  index_warm(:), & ! [IN]
1538  dens(:), & ! [IN]
1539  temp(:), & ! [INOUT]
1540  ghyd(:,:,:), & ! [INOUT]
1541  dt ) ! [IN]
1542 
1543  ! condensation / evaporation
1544  call cndevpsbla( ijkmax, & ! [IN]
1545  dens(:), & ! [IN]
1546  pres(:), & ! [IN]
1547  temp(:), & ! [INOUT]
1548  qvap(:), & ! [INOUT]
1549  ghyd(:,:,:), & ! [INOUT]
1550  gaer(:,:), & ! [INOUT]
1551  evaporate(:),& ! [OUT]
1552  dt ) ! [IN]
1553 
1554  if ( mp_doautoconversion ) then
1555  ! collision-coagulation
1556  call collmainf( ijkmax, & ! [IN]
1557  temp(:), & ! [IN]
1558  ghyd(:,:,:), & ! [INOUT]
1559  dt ) ! [IN]
1560  endif
1561 
1562  endif
1563 
1564  elseif( nccn == 0 ) then
1565 
1566  if ( nspc == 1 ) then
1567  !---< warm rain only without aerosol tracer >---
1568 
1569  ! nucleation
1570  call nucleat( ijkmax, & ! [IN]
1571  dens(:), & ! [IN]
1572  pres(:), & ! [IN]
1573  ccn(:), & ! [IN]
1574  temp(:), & ! [INOUT]
1575  qvap(:), & ! [INOUT]
1576  ghyd(:,:,:), & ! [INOUT]
1577  dt ) ! [IN]
1578 
1579  ! condensation / evaporation
1580  call cndevpsbl( ijkmax, & ! [IN]
1581  dens(:), & ! [IN]
1582  pres(:), & ! [IN]
1583  temp(:), & ! [INOUT]
1584  qvap(:), & ! [INOUT]
1585  ghyd(:,:,:), & ! [INOUT]
1586  evaporate(:),& ! [OUT]
1587  dt ) ! [IN]
1588 
1589  if ( mp_doautoconversion ) then
1590  ! collision-coagulation
1591  call collmain( ijkmax, & ! [IN]
1592  temp(:), & ! [IN]
1593  ghyd(:,:,:), & ! [INOUT]
1594  dt ) ! [IN]
1595  endif
1596 
1597  elseif( nspc > 1 ) then
1598  !---< mixed phase rain only without aerosol tracer >---
1599 
1600  ! nucleation
1601  call nucleat( ijkmax, & ! [IN]
1602  dens(:), & ! [IN]
1603  pres(:), & ! [IN]
1604  ccn(:), & ! [IN]
1605  temp(:), & ! [INOUT]
1606  qvap(:), & ! [INOUT]
1607  ghyd(:,:,:), & ! [INOUT]
1608  dt ) ! [IN]
1609 
1610  ! freezing / melting
1611  call freezing( ijkmax, & ! [IN]
1612  ijkmax_cold, & ! [IN]
1613  index_cold(:), & ! [IN]
1614  dens(:), & ! [IN]
1615  temp(:), & ! [INOUT]
1616  ghyd(:,:,:), & ! [INOUT]
1617  dt ) ! [IN]
1618 
1619  call ice_nucleat( ijkmax, & ! [IN]
1620  ijkmax_cold, & ! [IN]
1621  index_cold(:), & ! [IN]
1622  dens(:), & ! [IN]
1623  pres(:), & ! [IN]
1624  temp(:), & ! [INOUT]
1625  qvap(:), & ! [INOUT]
1626  ghyd(:,:,:), & ! [INOUT]
1627  dt ) ! [IN]
1628 
1629  call melting( ijkmax, & ! [IN]
1630  ijkmax_warm, & ! [IN]
1631  index_warm(:), & ! [IN]
1632  dens(:), & ! [IN]
1633  temp(:), & ! [INOUT]
1634  ghyd(:,:,:), & ! [INOUT]
1635  dt ) ! [IN]
1636 
1637  ! condensation / evaporation
1638  call cndevpsbl( ijkmax, & ! [IN]
1639  dens(:), & ! [IN]
1640  pres(:), & ! [IN]
1641  temp(:), & ! [INOUT]
1642  qvap(:), & ! [INOUT]
1643  ghyd(:,:,:), & ! [INOUT]
1644  evaporate(:),& ! [OUT]
1645  dt ) ! [IN]
1646 
1647  if ( mp_doautoconversion ) then
1648  ! collision-coagulation
1649  call collmainf( ijkmax, & ! [IN]
1650  temp(:), & ! [IN]
1651  ghyd(:,:,:), & ! [INOUT]
1652  dt ) ! [IN]
1653  endif
1654 
1655  endif
1656 
1657  endif
1658 
1659  return
1660  end subroutine mp_suzuki10
1661 
1662  !-----------------------------------------------------------------------------
1663  subroutine nucleat( &
1664  ijkmax, &
1665  dens, &
1666  pres, &
1667  ccn, &
1668  temp, &
1669  qvap, &
1670  gc, &
1671  dtime )
1672  implicit none
1673 
1674  integer, intent(in) :: ijkmax
1675  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
1676  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
1677  real(RP), intent(in) :: ccn(ijkmax) ! CCN number concentration [#/m3]
1678  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
1679  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
1680  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
1681  real(DP), intent(in) :: dtime ! Time step interval
1682 
1683  real(RP) :: ssliq(ijkmax)
1684  real(RP) :: qlevp(ijkmax) ! LH
1685  real(RP) :: dmp
1686  integer :: n
1687  !
1688  real(RP) :: n_c
1689  real(RP) :: sumnum(ijkmax)
1690  real(RP) :: gcn( nbin,ijkmax ) ! number of cloud particles in each bin (=gc/exp(xctr))
1691  real(RP) :: psat
1692  integer :: ijk
1693 
1694  call prof_rapstart('_SBM_Nucleat', 3)
1695 
1696  if( mp_couple_aerosol ) then
1697 
1698  do ijk = 1, ijkmax
1699  dmp = ccn(ijk) * expxctr( 1 )
1700  dmp = min( dmp,qvap(ijk)*dens(ijk) )
1701  gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
1702  qvap(ijk) = qvap(ijk) - dmp/dens(ijk)
1703  qvap(ijk) = max( qvap(ijk),0.0_rp )
1704  temp(ijk) = temp(ijk) + dmp/dens(ijk)*qlevp(ijk)/cp
1705  enddo
1706 
1707  else
1708 
1709  do ijk = 1, ijkmax
1710 
1711  !--- lhv
1712  qlevp(ijk) = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
1713  !--- supersaturation
1714  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
1715  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
1716  ssliq(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
1717  ssliq(ijk) = qvap(ijk)/ssliq(ijk)-1.0_rp
1718 
1719  enddo
1720 
1721  sumnum(:) = 0.0_rp
1722  do ijk = 1, ijkmax
1723 ! if ( ssliq <= 0.0_RP ) cycle
1724  if ( ssliq(ijk) > 0.0_rp ) then
1725 
1726  !--- use for aerosol coupled model
1727  !--- mass -> number
1728  do n = 1, nbin
1729  gcn( n,ijk ) = gc( n,il,ijk )*rexpxctr( n )
1730  enddo
1731 
1732  do n = 1, nbin
1733  sumnum(ijk) = sumnum(ijk) + gcn( n,ijk )*dxmic
1734  enddo
1735  n_c = c_ccn * ( ssliq(ijk) * 1.e+2_rp )**( kappa )
1736  if ( n_c > sumnum(ijk) ) then
1737  dmp = ( n_c - sumnum(ijk) ) * expxctr( 1 )
1738  dmp = min( dmp,qvap(ijk)*dens(ijk) )
1739  gc( 1,il,ijk ) = gc( 1,il,ijk ) + dmp/dxmic
1740  qvap(ijk) = qvap(ijk) - dmp/dens(ijk)
1741  qvap(ijk) = max( qvap(ijk),0.0_rp )
1742  temp(ijk) = temp(ijk) + dmp/dens(ijk)*qlevp(ijk)/cp
1743  endif
1744  endif
1745 
1746  enddo
1747 
1748  endif
1749 
1750  call prof_rapend ('_SBM_Nucleat', 3)
1751 
1752  return
1753  end subroutine nucleat
1754 
1755  !-----------------------------------------------------------------------------
1756  subroutine nucleata( &
1757  ijkmax, &
1758  dens, &
1759  pres, &
1760  temp, &
1761  qvap, &
1762  gc, &
1763  ga, &
1764  dtime )
1765  implicit none
1766 
1767  integer, intent(in) :: ijkmax
1768  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
1769  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
1770  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
1771  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
1772  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
1773  real(RP), intent(inout) :: ga (nccn ,ijkmax) ! Mass size distribution function of aerosol
1774  real(DP), intent(in) :: dtime ! Time step interval
1775 
1776  real(RP) :: gan( nccn ) ! size distribution function ( aerosol ) : number ( gan = ga/exp( xactr ) )
1777  real(RP) :: ssliq, ssice, qlevp ! supersaturatioin of liq. and ice, and LH
1778  real(RP) :: acoef, bcoef ! A and B in eq. (A.11) of Suzuki (2004)
1779  real(RP) :: rcrit ! critical radius (rcrit, r_N,crit of (A.11) of Suzuki (2004))
1780  real(RP) :: xcrit ! exp of hydrometeror whose radi is corresponding to rcrit (xcrit)
1781  real(RP) :: ractr, rcld, xcld, part, dmp
1782  integer :: n, nc, ncrit
1783 ! integer, allocatable, save :: ncld( : )
1784 ! integer, save :: ncld( 1:nccn )
1785 ! logical, save :: ofirst(1:ijkmax) = .true.
1786  !
1787  real(RP) :: psat
1788  integer :: ijk
1789 
1790  call prof_rapstart('_SBM_NucleatA', 3)
1791 
1792  do ijk = 1, ijkmax
1793  !
1794  !--- lhv
1795  qlevp = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
1796  !--- supersaturation
1797  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
1798  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
1799  ssliq = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
1800  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
1801  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
1802  ssice = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
1803  ssliq = qvap(ijk)/ssliq-1.0_rp
1804  ssice = qvap(ijk)/ssice-1.0_rp
1805 
1806  if ( ssliq <= 0.0_rp ) cycle
1807  !--- use for aerosol coupled model
1808  !--- mass -> number
1809  do n = 1, nccn
1810  gan( n ) = ga( n,ijk )*rexpxactr( n )
1811  enddo
1812 
1813  acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk) ! A in (A.11) of Suzuki (2004)
1814  bcoef = vhfct* rhoa/rhow * emwtr/emaer ! B in (A.11) of Suzuki (2004)
1815 
1816  !--- relationship of bin number
1817  do n = 1, nccn
1818  ractr = ( expxactr( n )*thirdovforth/pi/rhoa )**( oneovthird )
1819  rcld = sqrt( 3.0_rp*bcoef*ractr*ractr*ractr / acoef )
1820  xcld = log( rhow * 4.0_rp*pi*oneovthird*rcld*rcld*rcld )
1821  if ( flg_nucl ) then
1822  ncld( n ) = 1
1823  else
1824  ncld( n ) = int( ( xcld-xctr( 1 ) )/dxmic ) + 1
1825  ncld( n ) = min( max( ncld( n ),1 ),nbin )
1826  endif
1827  enddo
1828 
1829  !--- nucleation
1830  do n = nccn, 1, -1
1831  !--- super saturation
1832  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
1833  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
1834  ssliq = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
1835  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
1836  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
1837  ssice = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
1838  ssliq = qvap(ijk)/ssliq-1.0_rp
1839  ssice = qvap(ijk)/ssice-1.0_rp
1840 
1841  if ( ssliq <= 0.0_rp ) exit
1842  !--- use for aerosol coupled model
1843  acoef = 2.0_rp*sigma/rvap/rhow/temp(ijk) ! A in (A.11) of Suzuki (2004)
1844  rcrit = acoef*oneovthird * ( 4.0_rp/bcoef )**( oneovthird ) / ssliq**( twoovthird ) ! r_{N,crit} in (A.11) of Suzuki (2004)
1845  xcrit = log( rhoa * 4.0_rp*pi*oneovthird * rcrit*rcrit*rcrit )
1846  ncrit = int( ( xcrit-xabnd( 1 ) )/dxaer ) + 1
1847 
1848  if ( n == ncrit ) then
1849  part = ( xabnd( ncrit+1 )-xcrit )/dxaer
1850  elseif ( n > ncrit ) then
1851  part = 1.0_rp
1852  else
1853  exit
1854  endif
1855 
1856  !--- calculate mass change
1857  nc = ncld( n )
1858  dmp = part*gan( n )*dxaer*expxctr( nc )
1859  dmp = min( dmp,qvap(ijk)*dens(ijk) )
1860  gc( nc,il,ijk ) = gc( nc,il,ijk ) + dmp/dxmic
1861  gan( n ) = gan( n ) - dmp/dxaer*rexpxctr( nc )
1862  gan( n ) = max( gan( n ), 0.0_rp )
1863  qvap(ijk) = qvap(ijk) - dmp/dens(ijk)
1864  qvap(ijk) = max( qvap(ijk),0.0_rp )
1865  temp(ijk) = temp(ijk) + dmp/dens(ijk)*qlevp/cp
1866  enddo
1867 
1868  !--- number -> mass
1869  do n = 1, nccn
1870  ga( n,ijk ) = gan( n )*expxactr( n )
1871  enddo
1872 
1873  enddo
1874 
1875  call prof_rapend ('_SBM_NucleatA', 3)
1876 
1877  return
1878  end subroutine nucleata
1879 
1880  !-----------------------------------------------------------------------------
1881  subroutine cndevpsbl( &
1882  ijkmax, &
1883  dens, &
1884  pres, &
1885  temp, &
1886  qvap, &
1887  gc, &
1888  evaporate, & ! [OUT]
1889  dtime )
1890  implicit none
1891 
1892  integer, intent(in) :: ijkmax
1893  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
1894  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
1895  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
1896  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
1897  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
1898  real(RP), intent(out) :: evaporate(ijkmax) ! Number concentration of evaporated cloud [/m3]
1899  real(DP), intent(in) :: dtime ! Time step interval
1900  !---------------------------------------------------------------------------
1901 
1902  call liqphase( ijkmax, & ! [IN]
1903  dens(:), & ! [IN]
1904  pres(:), & ! [IN]
1905  temp(:), & ! [INOUT]
1906  qvap(:), & ! [INOUT]
1907  gc(:,:,:), & ! [INOUT]
1908  evaporate(:), & ! [OUT]
1909  dtime ) ! [IN]
1910 
1911  if( nspc > 1 ) then
1912  call icephase( ijkmax, & ! [IN]
1913  dens(:), & ! [IN]
1914  pres(:), & ! [IN]
1915  temp(:), & ! [INOUT]
1916  qvap(:), & ! [INOUT]
1917  gc(:,:,:), & ! [INOUT]
1918  dtime ) ! [IN]
1919 
1920  call mixphase( ijkmax, & ! [IN]
1921  dens(:), & ! [IN]
1922  pres(:), & ! [IN]
1923  temp(:), & ! [INOUT]
1924  qvap(:), & ! [INOUT]
1925  gc(:,:,:), & ! [INOUT]
1926  dtime ) ! [IN]
1927  endif
1928 
1929  return
1930  end subroutine cndevpsbl
1931 
1932  !-----------------------------------------------------------------------------
1933  subroutine cndevpsbla( &
1934  ijkmax, &
1935  dens, &
1936  pres, &
1937  temp, &
1938  qvap, &
1939  gc, &
1940  ga, &
1941  evaporate, &
1942  dtime )
1943  implicit none
1944 
1945  integer, intent(in) :: ijkmax
1946  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
1947  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
1948  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
1949  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
1950  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
1951  real(RP), intent(inout) :: ga (nccn ,ijkmax) ! Mass size distribution function of aerosol
1952  real(RP), intent(out) :: evaporate(ijkmax) ! Number concentration of evaporated cloud [/m3]
1953  real(DP), intent(in) :: dtime ! Time step interval
1954  !---------------------------------------------------------------------------
1955 
1956  call liqphase( ijkmax, & ! [IN]
1957  dens(:), & ! [IN]
1958  pres(:), & ! [IN]
1959  temp(:), & ! [INOUT]
1960  qvap(:), & ! [INOUT]
1961  gc(:,:,:), & ! [INOUT]
1962  evaporate(:), & ! [OUT]
1963  dtime ) ! [IN]
1964 
1965  ! regeneration of aerosol
1966  if ( flg_regeneration ) then
1967  call faero( ijkmax, & ! [IN]
1968 ! regene_gcn(:), & ! [IN]
1969  evaporate(:), & ! [IN]
1970  ga(:,:) ) ! [INOUT]
1971  endif
1972 
1973  if( nspc > 1 ) then
1974  call icephase( ijkmax, & ! [IN]
1975  dens(:), & ! [IN]
1976  pres(:), & ! [IN]
1977  temp(:), & ! [INOUT]
1978  qvap(:), & ! [INOUT]
1979  gc(:,:,:), & ! [INOUT]
1980  dtime ) ! [IN]
1981 
1982  call mixphase( ijkmax, & ! [IN]
1983  dens(:), & ! [IN]
1984  pres(:), & ! [IN]
1985  temp(:), & ! [INOUT]
1986  qvap(:), & ! [INOUT]
1987  gc(:,:,:), & ! [INOUT]
1988  dtime ) ! [IN]
1989  endif
1990 
1991  return
1992  end subroutine cndevpsbla
1993 
1994  !-----------------------------------------------------------------------------
1995  subroutine liqphase( &
1996  ijkmax, &
1997  dens, &
1998  pres, &
1999  temp, &
2000  qvap, &
2001  gc, &
2002  regene_gcn, &
2003  dtime )
2004  implicit none
2005 
2006  integer, intent(in) :: ijkmax
2007  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
2008  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
2009  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
2010  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
2011  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
2012  real(RP), intent(out) :: regene_gcn(ijkmax) ! mass of regenerated aerosol
2013  real(DP), intent(in) :: dtime ! Time step interval
2014 
2015  integer :: n, myu, ncount
2016  integer :: nloop(ijkmax) ! number of fractional step for condensation
2017  real(RP) :: gclold(ijkmax)
2018  real(RP) :: gclnew(ijkmax)
2019  real(RP) :: cndmss(ijkmax)
2020  real(RP) :: dtcnd(ijkmax) ! dt for condensation with considering CFL condition
2021  real(RP) :: gtliq(ijkmax) ! G of eq. (A.17) of Suzuki (2004)
2022  real(RP) :: qlevp(ijkmax) ! LH
2023  real(RP) :: cefliq, a, sliqtnd
2024  real(RP) :: sumliq(ijkmax), umax(ijkmax)
2025  real(RP) :: ssliq(ijkmax), ssice(ijkmax) ! super saturation
2026  real(RP) :: gcn( -1:nbin+2,nspc,ijkmax ) ! size distribution function(hydrometeor): number
2027 ! real(RP) :: gcnold( nbin,ijkmax )
2028  real(RP), parameter :: cflfct = 0.50_rp ! CFL limiter
2029 ! real(RP) :: old_sum_gcn, new_sum_gcn
2030  integer :: iflg( nspc,ijkmax ) ! flag whether calculation is conduct or not
2031  real(RP) :: csum( nspc,ijkmax )
2032  real(RP) :: f1, f2, emu, cefd, cefk, festl, psat
2033  real(RP), parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2034  real(RP), parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2035  real(RP) :: zerosw, qvtmp
2036  integer :: ijk, nn, mm, loopflg(ijkmax)
2037 
2038  !--- local for advection
2039  real(RP) :: uadv ( 0:nbin+2,nspc,ijkmax )
2040  real(RP) :: flq ( 1:nbin+1,nspc,ijkmax )
2041  real(RP) :: acoef( 0:2,0:nbin+1,nspc,ijkmax )
2042  real(RP) :: crn ( 0:nbin+2,nspc,ijkmax )
2043  real(RP) :: aip ( 0:nbin+1,nspc,ijkmax )
2044  real(RP) :: aim ( 0:nbin+1,nspc,ijkmax )
2045  real(RP) :: ai ( 0:nbin+1,nspc,ijkmax )
2046  real(RP) :: cmins, cplus
2047  integer :: nloopmax
2048 
2049  call prof_rapstart('_SBM_Liqphase', 3)
2050 
2051 
2052  iflg(:,:) = 0
2053  csum(:,:) = 0.0_rp
2054  do ijk = 1, ijkmax
2055 
2056  do n = 1, nbin
2057  csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
2058  enddo
2059 
2060  if( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
2061 
2062  enddo
2063 
2064  nloop(:) = 0
2065  gclold(:) = 0.0_rp
2066  do ijk = 1, ijkmax
2067 
2068  do n = 1, nbin
2069  gclold(ijk) = gclold(ijk) + gc( n,il,ijk ) * dxmic
2070  enddo
2071  !
2072  !------- mass -> number
2073  do n = 1, nbin
2074  gcn( n,il,ijk ) = gc( n,il,ijk ) * rexpxctr( n )
2075  enddo
2076  gcn( -1,il,ijk ) = 0.0_rp
2077  gcn( 0,il,ijk ) = 0.0_rp
2078  gcn( nbin+1,il,ijk ) = 0.0_rp
2079  gcn( nbin+2,il,ijk ) = 0.0_rp
2080  !
2081  !--- lhv
2082  qlevp(ijk) = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2083  !--- super saturation
2084  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
2085  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
2086  ssliq(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2087  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
2088  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
2089  ssice(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2090  ssliq(ijk) = qvap(ijk)/ssliq(ijk)-1.0_rp
2091  ssice(ijk) = qvap(ijk)/ssice(ijk)-1.0_rp
2092 
2093  emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2094  cefd = emu/dens(ijk)
2095  cefk = fct*emu
2096 
2097  festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2098  f1 = rvap*temp(ijk)/festl/cefd
2099  f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2100  gtliq(ijk) = 4.0_rp*pi/( f1+f2 ) !--- G of eq. (A.17) of Suzuki (2004)
2101  !------- CFL condition
2102  umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq(ijk) )
2103  dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2104  nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2105  dtcnd(ijk) = dtime / nloop(ijk)
2106 
2107  nloop(ijk) = nloop(ijk) * iflg( il,ijk ) !--- for determing trivial loop
2108  enddo
2109  nloopmax = maxval(nloop,1)
2110 
2111  !
2112  regene_gcn(:) = 0.0_rp
2113 !OCL LOOP_NOFISSION
2114 !OCL LOOP_NOINTERCHANGE
2115  do ncount = 1, nloopmax
2116 
2117  do ijk = 1, ijkmax
2118  loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) ) ! 0 or 1
2119  enddo
2120 
2121 !OCL LOOP_NOFUSION
2122  do ijk = 1, ijkmax
2123  do nn = 1, loopflg(ijk)
2124 
2125  !--- lhv
2126  qlevp(ijk) = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2127  !----- matrix for supersaturation tendency
2128  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
2129  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
2130  ssliq(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2131  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
2132  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
2133  ssice(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2134  ssliq(ijk) = qvap(ijk)/ssliq(ijk)-1.0_rp
2135  ssice(ijk) = qvap(ijk)/ssice(ijk)-1.0_rp
2136 
2137  emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2138  cefd = emu/dens(ijk)
2139  cefk = fct*emu
2140  festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2141  f1 = rvap*temp(ijk)/festl/cefd
2142  f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2143  gtliq(ijk) = 4.0_rp*pi/( f1+f2 ) !--- G of eq. (A.17) of Suzuki (2004)
2144 
2145  sumliq(ijk) = 0.0_rp
2146 ! old_sum_gcn = 0.0_RP
2147  do n = 1, nbin
2148  sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
2149  enddo
2150 
2151  enddo
2152  enddo
2153 
2154 !OCL LOOP_NOFUSION
2155  do ijk = 1, ijkmax
2156  do nn = 1, loopflg(ijk)
2157 
2158  !----- supersaturation tendency
2159  zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps ) !--- zerosw = 1 (qv>0), zerosw=0 (qv=0)
2160  qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
2161  cefliq = ( ssliq(ijk)+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)*qlevp(ijk)/cp/rvap/temp(ijk)/temp(ijk) )
2162  a = - cefliq*sumliq(ijk)*gtliq(ijk)/dens(ijk) ! a of eq. (A.19) of Suzuki (2004)
2163  a = a + eps * ( 1.0_rp - zerosw ) !--- avoiding division by zero when qv = 0
2164 
2165  sliqtnd = zerosw * &
2166  ( ssliq(ijk) * ( exp( a*dtcnd(ijk) )-1.0_rp )/( a*dtcnd(ijk) ) &
2167  * ( 0.5_rp + sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
2168  + ssliq(ijk) &
2169  * ( 0.5_rp - sign( 0.5_rp,abs(a*dtcnd(ijk)-0.1_rp) ) ) &
2170  ) &
2171  + ssliq(ijk) * ( 1.0_rp - zerosw )
2172  !
2173  !----- change of SDF
2174  do mm = 1, iflg( il,ijk )
2175  !--- advection speed
2176  do n = 1, nbin+1
2177  uadv( n,il,ijk ) = cbnd( n,il )*rexpxbnd( n )*gtliq(ijk)*sliqtnd ! U of eq. (A.18) of Suzuki (2004)
2178  enddo
2179  uadv( 0 ,il,ijk ) = 0.0_rp
2180  uadv( nbin+2,il,ijk ) = 0.0_rp
2181 
2182 ! do n = 1, nbin
2183 ! gcnold( n ) = gcn( n,il,ijk )
2184 ! enddo
2185  enddo
2186 
2187  enddo
2188  enddo
2189 
2190  call prof_rapstart('_SBM_AdvLiq', 3)
2191 
2192  do ijk = 1, ijkmax
2193  do nn = 1, loopflg(ijk)
2194  do myu = 1, il
2195  do mm = 1, iflg( myu,ijk )
2196  do n = 0, nbin+2
2197  crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
2198  enddo
2199  enddo
2200  enddo
2201  enddo
2202  enddo
2203 
2204  do ijk = 1, ijkmax
2205  do nn = 1, loopflg(ijk)
2206  do myu = 1, il
2207  do mm = 1, iflg( myu,ijk )
2208  do n = 0, nbin+1
2209  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
2210  acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
2211  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
2212  enddo
2213  enddo
2214  enddo
2215  enddo
2216  enddo
2217 
2218  do ijk = 1, ijkmax
2219  do nn = 1, loopflg(ijk)
2220  do myu = 1, il
2221  do mm = 1, iflg( myu,ijk )
2222  do n = 0, nbin+1
2223  cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
2224 
2225  aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
2226  + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
2227  + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
2228 
2229  aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
2230  enddo
2231  enddo
2232  enddo
2233  enddo
2234  enddo
2235 
2236  do ijk = 1, ijkmax
2237  do nn = 1, loopflg(ijk)
2238  do myu = 1, il
2239  do mm = 1, iflg( myu,ijk )
2240  do n = 0, nbin+1
2241  cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
2242 
2243  aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
2244  - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
2245  + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
2246 
2247  aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
2248  enddo
2249  enddo
2250  enddo
2251  enddo
2252  enddo
2253 
2254  do ijk = 1, ijkmax
2255  do nn = 1, loopflg(ijk)
2256  do myu = 1, il
2257  do mm = 1, iflg( myu,ijk )
2258  do n = 0, nbin+1
2259  ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
2260  + acoef(2,n,myu,ijk) * 2.0_rp
2261 
2262  ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
2263  enddo
2264  enddo
2265  enddo
2266  enddo
2267  enddo
2268 
2269  do ijk = 1, ijkmax
2270  do nn = 1, loopflg(ijk)
2271  do myu = 1, il
2272  do mm = 1, iflg( myu,ijk )
2273  do n = 1, nbin+1
2274  flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
2275  - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
2276  enddo
2277  enddo
2278  enddo
2279  enddo
2280  enddo
2281 
2282  do ijk = 1, ijkmax
2283  do nn = 1, loopflg(ijk)
2284  do myu = 1, il
2285  do mm = 1, iflg( myu,ijk )
2286  regene_gcn(ijk) = regene_gcn(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
2287  * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
2288  enddo
2289  enddo
2290  enddo
2291  enddo
2292 
2293  do ijk = 1, ijkmax
2294  do nn = 1, loopflg(ijk)
2295  do myu = 1, il
2296  do mm = 1, iflg( myu,ijk )
2297  do n = 1, nbin
2298  gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
2299  enddo
2300  enddo
2301  enddo
2302  enddo
2303  enddo
2304 
2305  call prof_rapend ('_SBM_AdvLiq', 3)
2306 
2307 !OCL LOOP_NOFUSION
2308  do ijk = 1, ijkmax
2309  do nn = 1, loopflg(ijk)
2310 
2311  !----- new mass
2312  gclnew(ijk) = 0.0_rp
2313 ! new_sum_gcn = 0.0_RP
2314  do n = 1, nbin
2315  gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )
2316 ! old_sum_gcn = old_sum_gcn + gcnold( n )*dxmic
2317 ! new_sum_gcn = new_sum_gcn + gcn( n,ijk )*dxmic
2318  enddo
2319 
2320  gclnew(ijk) = gclnew(ijk)*dxmic
2321  !
2322  !----- change of humidity and temperature
2323  cndmss(ijk) = gclnew(ijk) - gclold(ijk)
2324  qvap(ijk) = qvap(ijk) - cndmss(ijk)/dens(ijk)
2325  temp(ijk) = temp(ijk) + cndmss(ijk)/dens(ijk)*qlevp(ijk)/cp
2326  !
2327  gclold(ijk) = gclnew(ijk)
2328  !
2329  !----- continue/end
2330  enddo
2331  enddo
2332 
2333  enddo !nloop
2334  !
2335 !OCL NORECURRENCE(gc)
2336  do ijk = 1, ijkmax
2337  !------- number -> mass
2338  do n = 1 , nbin
2339  gc( n,il,ijk ) = gcn( n,il,ijk )*expxctr( n )
2340  enddo
2341  enddo
2342 
2343 !OCL NORECURRENCE(gc)
2344  do ijk = 1, ijkmax
2345  !--- lhv
2346  qlevp(ijk) = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2347  do n = 1 , nbin
2348  if ( gc( n,il,ijk ) < 0.0_rp ) then
2349  cndmss(ijk) = -gc( n,il,ijk )
2350  gc( n,il,ijk ) = 0.0_rp
2351  qvap(ijk) = qvap(ijk) + cndmss(ijk)/dens(ijk)
2352  temp(ijk) = temp(ijk) - cndmss(ijk)/dens(ijk)*qlevp(ijk)/cp
2353  endif
2354  enddo
2355  enddo
2356 
2357  call prof_rapend ('_SBM_Liqphase', 3)
2358 
2359  return
2360  end subroutine liqphase
2361 
2362  !-----------------------------------------------------------------------------
2363  subroutine icephase( &
2364  ijkmax, &
2365  dens, &
2366  pres, &
2367  temp, &
2368  qvap, &
2369  gc, &
2370  dtime )
2371  implicit none
2372 
2373  integer, intent(in) :: ijkmax
2374  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
2375  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
2376  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
2377  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
2378  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
2379  real(DP), intent(in) :: dtime ! Time step interval
2380 
2381  integer :: myu, n, ncount
2382  integer :: nloop(ijkmax) !number of fractional step for condensation
2383  real(RP) :: gciold(ijkmax), gcinew(ijkmax)
2384  real(RP) :: dtcnd(ijkmax) ! dt for condensation with considering CFL condition
2385  real(RP) :: sblmss(ijkmax)
2386  real(RP) :: gtice(ijkmax), umax(ijkmax)
2387  real(RP) :: qlsbl(ijkmax)
2388  real(RP) :: cefice, d, uval, sicetnd
2389  real(RP) :: sumice(ijkmax)
2390  real(RP) :: ssliq(ijkmax), ssice(ijkmax)
2391  real(RP) :: gcn( -1:nbin+2,nspc,ijkmax ) ! size distribution function (Hydrometeor): number
2392  real(RP), parameter :: cflfct = 0.50_rp
2393  real(RP) :: dumm_regene(ijkmax)
2394  integer :: iflg( nspc,ijkmax )
2395  real(RP) :: csum( nspc,ijkmax )
2396  real(RP) :: f1, f2, emu, cefd, cefk, festi, psat
2397  real(RP), parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2398  real(RP), parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2399  integer :: ijk, mm, nn, loopflg(ijkmax)
2400  real(RP) :: zerosw, qvtmp
2401 
2402  !--- local for advection
2403 ! real(RP) :: qadv( -1:nbin+2 ), uadv( 0:nbin+2 )
2404  real(RP) :: uadv ( 0:nbin+2,nspc,ijkmax )
2405  real(RP) :: flq ( 1:nbin+1,nspc,ijkmax )
2406  real(RP) :: acoef( 0:2,0:nbin+1,nspc,ijkmax )
2407  real(RP) :: crn ( 0:nbin+2,nspc,ijkmax )
2408  real(RP) :: aip ( 0:nbin+1,nspc,ijkmax )
2409  real(RP) :: aim ( 0:nbin+1,nspc,ijkmax )
2410  real(RP) :: ai ( 0:nbin+1,nspc,ijkmax )
2411  real(RP) :: cmins, cplus
2412  integer :: nloopmax
2413  !---------------------------------------------------------------------------
2414 
2415  call prof_rapstart('_SBM_Icephase', 3)
2416 
2417  iflg(:,:) = 0
2418  csum( :,: ) = 0.0_rp
2419  do ijk = 1, ijkmax
2420 
2421  do myu = 2, nspc
2422  do n = 1, nbin
2423  csum( myu,ijk ) = csum( myu,ijk ) + gc( n,myu,ijk )*dxmic
2424  enddo
2425  enddo
2426 
2427  do myu = 2, nspc
2428  if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
2429  enddo
2430 
2431  enddo
2432 
2433  gciold(:) = 0.0_rp
2434  nloop(:) = 0
2435  do ijk = 1, ijkmax
2436 
2437  !----- old mass
2438  do myu = 2, nspc
2439  do n = 1, nbin
2440  gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
2441  enddo
2442  enddo
2443 
2444  !----- mass -> number
2445  do myu = 2, nspc
2446  do n = 1, nbin
2447  gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
2448  enddo
2449  gcn( -1,myu,ijk ) = 0.0_rp
2450  gcn( 0,myu,ijk ) = 0.0_rp
2451  gcn( nbin+1,myu,ijk ) = 0.0_rp
2452  gcn( nbin+2,myu,ijk ) = 0.0_rp
2453  enddo
2454 
2455  !--- lhv
2456  qlsbl(ijk) = const_lhs0 + ( const_cpvap - const_ci ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2457  !--- supersaturation
2458  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
2459  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
2460  ssliq(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2461  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
2462  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
2463  ssice(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2464  ssliq(ijk) = qvap(ijk)/ssliq(ijk)-1.0_rp
2465  ssice(ijk) = qvap(ijk)/ssice(ijk)-1.0_rp
2466 
2467  emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2468  cefd = emu/dens(ijk)
2469  cefk = fct*emu
2470  festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2471  f1 = rvap*temp(ijk)/festi/cefd
2472  f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
2473  gtice(ijk) = 4.0_rp*pi/( f1+f2 )
2474  !----- CFL condition
2475  umax(ijk) = 0.0_rp
2476  do myu = 2, nspc
2477  uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice(ijk) )
2478  umax(ijk) = max( umax(ijk),uval )
2479  enddo
2480 
2481  dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2482  nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2483  dtcnd(ijk) = dtime/nloop(ijk)
2484 
2485  nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) ) !--- for determing trivial loop
2486 
2487  enddo
2488  nloopmax = maxval(nloop,1)
2489 
2490 !OCL LOOP_NOFISSION
2491 !OCL LOOP_NOINTERCHANGE
2492  do ncount = 1, nloopmax
2493 
2494  do ijk = 1, ijkmax
2495  loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) ) ! 0 or 1
2496  enddo
2497 
2498 !OCL LOOP_NOFUSION
2499  do ijk = 1, ijkmax
2500  do nn = 1, loopflg(ijk)
2501  !--- lhv
2502  qlsbl(ijk) = const_lhs0 + ( const_cpvap - const_ci ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2503  !----- matrix for supersaturation tendency
2504  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
2505  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
2506  ssliq(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2507  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
2508  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
2509  ssice(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2510  ssliq(ijk) = qvap(ijk)/ssliq(ijk)-1.0_rp
2511  ssice(ijk) = qvap(ijk)/ssice(ijk)-1.0_rp
2512 
2513  emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2514  cefd = emu/dens(ijk)
2515  cefk = fct*emu
2516  festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2517  f1 = rvap*temp(ijk)/festi/cefd
2518  f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
2519  gtice(ijk) = 4.0_rp*pi/( f1+f2 ) ! G of (A.17) of Suzuki (2004)
2520 
2521  sumice(ijk) = 0.0_rp
2522  do myu = 2, nspc
2523  do n = 1, nbin
2524  sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
2525  enddo
2526  enddo
2527 
2528  enddo
2529  enddo
2530 
2531 !OCL LOOP_NOFUSION
2532  do ijk = 1, ijkmax
2533  do nn = 1, loopflg(ijk)
2534 
2535  !----- supersaturation tendency
2536  zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps ) !--- zerosw = 1 (qv>0), zerosw=0 (qv=0)
2537  qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
2538 
2539  cefice = ( ssice(ijk)+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)*qlsbl(ijk)/cp/rvap/temp(ijk)/temp(ijk) )
2540  d = - cefice*sumice(ijk)*gtice(ijk)/dens(ijk) ! d of (A.19) of Suzuki (2004)
2541  d = d + eps * ( 1.0_rp - zerosw ) !--- avoiding division by zero when qv = 0
2542 
2543  sicetnd = zerosw * &
2544  ( ssice(ijk) * ( exp( d*dtcnd(ijk) )-1.0_rp )/( d*dtcnd(ijk) ) &
2545  * ( 0.5_rp + sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
2546  + ssice(ijk) &
2547  * ( 0.5_rp - sign( 0.5_rp,abs(d*dtcnd(ijk)-0.1_rp) ) ) &
2548  ) &
2549  + ssice(ijk) * ( 1.0_rp - zerosw )
2550  !
2551  !----- change of SDF
2552  do myu = 2, nspc
2553  do mm = 1, iflg( myu,ijk )
2554  !--- advection speed
2555  do n = 1, nbin+1
2556  uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd ! U of eq. (A.18) of Suzuki (2004)
2557  enddo
2558  uadv( 0, myu,ijk ) = 0.0_rp
2559  uadv( nbin+2,myu,ijk ) = 0.0_rp
2560  enddo
2561  enddo
2562 
2563  enddo
2564  enddo
2565 
2566  call prof_rapstart('_SBM_AdvIce', 3)
2567 
2568  do ijk = 1, ijkmax
2569  do nn = 1, loopflg(ijk)
2570  do myu = 2, nspc
2571  do mm = 1, iflg( myu,ijk )
2572  do n = 0, nbin+2
2573  crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
2574  enddo
2575  enddo
2576  enddo
2577  enddo
2578  enddo
2579 
2580  do ijk = 1, ijkmax
2581  do nn = 1, loopflg(ijk)
2582  do myu = 2, nspc
2583  do mm = 1, iflg( myu,ijk )
2584  do n = 0, nbin+1
2585  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
2586  acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
2587  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
2588  enddo
2589  enddo
2590  enddo
2591  enddo
2592  enddo
2593 
2594  do ijk = 1, ijkmax
2595  do nn = 1, loopflg(ijk)
2596  do myu = 2, nspc
2597  do mm = 1, iflg( myu,ijk )
2598  do n = 0, nbin+1
2599  cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
2600 
2601  aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
2602  + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
2603  + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
2604 
2605  aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
2606  enddo
2607  enddo
2608  enddo
2609  enddo
2610  enddo
2611 
2612  do ijk = 1, ijkmax
2613  do nn = 1, loopflg(ijk)
2614  do myu = 2, nspc
2615  do mm = 1, iflg( myu,ijk )
2616  do n = 0, nbin+1
2617  cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
2618 
2619  aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
2620  - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
2621  + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
2622 
2623  aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
2624  enddo
2625  enddo
2626  enddo
2627  enddo
2628  enddo
2629 
2630  do ijk = 1, ijkmax
2631  do nn = 1, loopflg(ijk)
2632  do myu = 2, nspc
2633  do mm = 1, iflg( myu,ijk )
2634  do n = 0, nbin+1
2635  ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
2636  + acoef(2,n,myu,ijk) * 2.0_rp
2637 
2638  ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
2639  enddo
2640  enddo
2641  enddo
2642  enddo
2643  enddo
2644 
2645  do ijk = 1, ijkmax
2646  do nn = 1, loopflg(ijk)
2647  do myu = 2, nspc
2648  do mm = 1, iflg( myu,ijk )
2649  do n = 1, nbin+1
2650  flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
2651  - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
2652  enddo
2653  enddo
2654  enddo
2655  enddo
2656  enddo
2657 
2658  do ijk = 1, ijkmax
2659  do nn = 1, loopflg(ijk)
2660  do myu = 2, nspc
2661  do mm = 1, iflg( myu,ijk )
2662  dumm_regene(ijk) = dumm_regene(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
2663  * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk) + eps )
2664  enddo
2665  enddo
2666  enddo
2667  enddo
2668 
2669  do ijk = 1, ijkmax
2670  do nn = 1, loopflg(ijk)
2671  do myu = 2, nspc
2672  do mm = 1, iflg( myu,ijk )
2673  do n = 1, nbin
2674  gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
2675  enddo
2676  enddo
2677  enddo
2678  enddo
2679  enddo
2680 
2681  call prof_rapend ('_SBM_AdvIce', 3)
2682 
2683 !OCL LOOP_NOFUSION
2684  do ijk = 1, ijkmax
2685  do nn = 1, loopflg(ijk)
2686 
2687  !----- new mass
2688  gcinew(ijk) = 0.0_rp
2689  do n = 1, nbin
2690  do myu = 2, nspc
2691  gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
2692  enddo
2693  enddo
2694  !
2695  !----- change of humidity and temperature
2696  sblmss(ijk) = gcinew(ijk) - gciold(ijk)
2697  qvap(ijk) = qvap(ijk) - sblmss(ijk)/dens(ijk)
2698  temp(ijk) = temp(ijk) + sblmss(ijk)/dens(ijk)*qlsbl(ijk)/cp
2699 
2700  gciold(ijk) = gcinew(ijk)
2701  !
2702  !----- continue / end
2703  enddo
2704  enddo
2705 
2706  enddo ! nloop
2707  !
2708 !OCL NORECURRENCE(gc)
2709  do ijk = 1, ijkmax
2710  !------- number -> mass
2711  do myu = 2, nspc
2712  do n = 1, nbin
2713  gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
2714  enddo
2715  enddo
2716  enddo
2717 
2718  call prof_rapend ('_SBM_Icephase', 3)
2719 
2720  return
2721  end subroutine icephase
2722 
2723  !-----------------------------------------------------------------------------
2724  subroutine mixphase( &
2725  ijkmax, &
2726  dens, &
2727  pres, &
2728  temp, &
2729  qvap, &
2730  gc, &
2731  dtime )
2732  implicit none
2733 
2734  integer, intent(in) :: ijkmax
2735  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
2736  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
2737  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
2738  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
2739  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
2740  real(DP), intent(in) :: dtime ! Time step interval
2741 
2742  integer :: n, myu, mm, ncount
2743  integer :: nloop(ijkmax)
2744  real(RP) :: gclold(ijkmax), gclnew(ijkmax)
2745  real(RP) :: gciold(ijkmax), gcinew(ijkmax)
2746  real(RP) :: cndmss(ijkmax), sblmss(ijkmax)
2747  real(RP) :: gtliq(ijkmax), gtice(ijkmax)
2748  real(RP) :: umax(ijkmax), uval, dtcnd(ijkmax)
2749  real(RP) :: qlevp(ijkmax), qlsbl(ijkmax)
2750  real(RP) :: cef1, cef2, cef3, cef4, a, b, c, d
2751  real(RP) :: rmdplus, rmdmins, ssplus, ssmins, tplus, tmins
2752  real(RP) :: sliqtnd, sicetnd
2753  real(RP) :: ssliq(ijkmax), ssice(ijkmax)
2754  real(RP) :: sumliq(ijkmax), sumice(ijkmax)
2755  real(RP) :: gcn( -1:nbin+2,nspc,ijkmax )
2756  real(RP), parameter :: cflfct = 0.50_rp
2757  real(RP) :: dumm_regene(ijkmax)
2758  real(RP) :: csum( nspc,ijkmax )
2759  integer :: iflg( nspc,ijkmax )
2760  real(RP) :: f1, f2, emu, cefd, cefk, festl, festi, psat
2761  real(RP), parameter :: afmyu = 1.72e-05_rp, bfmyu = 3.93e+2_rp
2762  real(RP), parameter :: cfmyu = 1.2e+02_rp, fct = 1.4e+3_rp
2763  real(RP) :: qvtmp, zerosw
2764  integer :: ijk, nn, loopflg(ijkmax)
2765 
2766  !--- local for advection
2767 ! real(RP) :: qadv( -1:nbin+2 ), uadv( 0:nbin+2 )
2768  real(RP) :: uadv ( 0:nbin+2,nspc,ijkmax )
2769  real(RP) :: flq ( 1:nbin+1,nspc,ijkmax )
2770  real(RP) :: acoef( 0:2,0:nbin+1,nspc,ijkmax )
2771  real(RP) :: crn ( 0:nbin+2,nspc,ijkmax )
2772  real(RP) :: aip ( 0:nbin+1,nspc,ijkmax )
2773  real(RP) :: aim ( 0:nbin+1,nspc,ijkmax )
2774  real(RP) :: ai ( 0:nbin+1,nspc,ijkmax )
2775  real(RP) :: cmins, cplus
2776  integer :: nloopmax
2777 
2778  call prof_rapstart('_SBM_Mixphase', 3)
2779 
2780  dumm_regene( : ) = 0.0_rp
2781  iflg( :,: ) = 0
2782  csum( :,: ) = 0.0_rp
2783  do ijk = 1, ijkmax
2784 
2785  do myu = 1, nspc
2786  do n = 1, nbin
2787  csum( myu,ijk ) = csum( myu,ijk )+gc( n,myu,ijk )*dxmic
2788  enddo
2789  enddo
2790 
2791  do myu = 1, nspc
2792  if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
2793  enddo
2794  enddo
2795 
2796  gclold(:) = 0.0_rp
2797  gciold(:) = 0.0_rp
2798  nloop(:) = 0
2799  do ijk = 1, ijkmax
2800  !----- old mass
2801  do n = 1, nbin
2802  gclold(ijk) = gclold(ijk) + gc( n,il,ijk )*dxmic
2803  enddo
2804 
2805  do myu = 2, nspc
2806  do n = 1, nbin
2807  gciold(ijk) = gciold(ijk) + gc( n,myu,ijk )*dxmic
2808  enddo
2809  enddo
2810 
2811  !----- mass -> number
2812  do myu = 1, nspc
2813  do n = 1, nbin
2814  gcn( n,myu,ijk ) = gc( n,myu,ijk ) * rexpxctr( n )
2815  enddo
2816  gcn( -1,myu,ijk ) = 0.0_rp
2817  gcn( 0,myu,ijk ) = 0.0_rp
2818  gcn( nbin+1,myu,ijk ) = 0.0_rp
2819  gcn( nbin+2,myu,ijk ) = 0.0_rp
2820  enddo
2821 
2822  !-- thermodyn
2823  qlevp(ijk) = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2824  qlsbl(ijk) = const_lhs0 + ( const_cpvap - const_ci ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2825  !-- supersaturation
2826  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
2827  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
2828  ssliq(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2829  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
2830  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
2831  ssice(ijk) = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
2832  ssliq(ijk) = qvap(ijk)/ssliq(ijk)-1.0_rp
2833  ssice(ijk) = qvap(ijk)/ssice(ijk)-1.0_rp
2834 
2835  emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2836  cefd = emu/dens(ijk)
2837  cefk = fct*emu
2838 
2839  festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2840  f1 = rvap*temp(ijk)/festl/cefd
2841  f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2842  gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2843 
2844  festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2845  f1 = rvap*temp(ijk)/festi/cefd
2846  f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
2847  gtice(ijk) = 4.0_rp*pi/( f1+f2 )
2848 
2849  !----- CFL condition
2850  umax(ijk) = cbnd( 1,il )*rexpxbnd( 1 )*gtliq(ijk)*abs( ssliq(ijk) )
2851  do myu = 2, nspc
2852  uval = cbnd( 1,myu )*rexpxbnd( 1 )*gtice(ijk)*abs( ssice(ijk) )
2853  umax(ijk) = max( umax(ijk),uval )
2854  enddo
2855 
2856  dtcnd(ijk) = cflfct*dxmic/umax(ijk)
2857  nloop(ijk) = int( dtime/dtcnd(ijk) ) + 1
2858  dtcnd(ijk) = dtime/nloop(ijk)
2859 
2860  nloop(ijk) = nloop(ijk) * iflg( il,ijk ) !--- for determing trivial loop
2861  nloop(ijk) = nloop(ijk) * maxval( iflg( 2:nspc,ijk ) ) !--- for determing trivial loop
2862 ! nloop(ijk) = nloop(ijk) * maxval( iflg( 1:nspc,ijk ) ) !--- for determing trivial loop
2863  enddo
2864  nloopmax = maxval(nloop,1)
2865 
2866 
2867 !OCL LOOP_NOFISSION
2868 !OCL LOOP_NOINTERCHANGE
2869  do ncount = 1, nloopmax
2870 
2871 
2872  do ijk = 1, ijkmax
2873  loopflg(ijk) = min( 1, int(nloop(ijk)/ncount) ) ! 0 or 1
2874  enddo
2875 
2876 !OCL LOOP_NOFUSION
2877  do ijk = 1, ijkmax
2878  do nn = 1, loopflg(ijk)
2879  !-- thermodyn
2880  qlevp(ijk) = const_lhv0 + ( const_cpvap - const_cl ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2881  qlsbl(ijk) = const_lhs0 + ( const_cpvap - const_ci ) * ( temp(ijk) - const_tem00 ) * flg_thermodyn
2882  !-- matrix for supersaturation tendency
2883  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
2884  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
2885  ssliq(ijk) = qvap(ijk) / ( const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat ) )-1.0_rp
2886  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
2887  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
2888  ssice(ijk) = qvap(ijk) / ( const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat ) )-1.0_rp
2889 
2890  emu = afmyu*( bfmyu/( temp(ijk)+cfmyu ) )*( temp(ijk)/tmlt )**1.50_rp
2891  cefd = emu/dens(ijk)
2892  cefk = fct*emu
2893  festl = esat0*exp( qlevp(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2894  f1 = rvap*temp(ijk)/festl/cefd
2895  f2 = qlevp(ijk)/cefk/temp(ijk)*( qlevp(ijk)/rvap/temp(ijk) - 1.0_rp )
2896  gtliq(ijk) = 4.0_rp*pi/( f1+f2 )
2897  festi = esat0*exp( qlsbl(ijk)/rvap*( 1.0_rp/temp0 - 1.0_rp/temp(ijk) ) )
2898  f1 = rvap*temp(ijk)/festi/cefd
2899  f2 = qlsbl(ijk)/cefk/temp(ijk)*( qlsbl(ijk)/rvap/temp(ijk) - 1.0_rp )
2900  gtice(ijk) = 4.0_rp*pi/( f1+f2 )
2901 
2902  sumliq(ijk) = 0.0_rp
2903  do n = 1, nbin
2904  sumliq(ijk) = sumliq(ijk) + gcn( n,il,ijk )*cctr( n,il )*dxmic
2905  enddo
2906 
2907  sumice(ijk) = 0.0_rp
2908  do myu = 2, nspc
2909  do n = 1, nbin
2910  sumice(ijk) = sumice(ijk) + gcn( n,myu,ijk )*cctr( n,myu )*dxmic
2911  enddo
2912  enddo
2913  enddo
2914  enddo
2915 
2916 
2917 !OCL LOOP_NOFUSION
2918  do ijk = 1, ijkmax
2919  do nn = 1, loopflg(ijk)
2920  zerosw = 0.5_rp + sign( 0.5_rp,qvap(ijk)-eps ) !--- zerosw = 1 (qv>0), zerosw=0 (qv=0)
2921  qvtmp = qvap(ijk) * zerosw + ( qvap(ijk)+eps ) * ( 1.0_rp-zerosw )
2922  cef1 = ( ssliq(ijk)+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cp )
2923  cef2 = ( ssliq(ijk)+1.0_rp )*( 1.0_rp/qvtmp + qlevp(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cp )
2924  cef3 = ( ssice(ijk)+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlevp(ijk)/cp )
2925  cef4 = ( ssice(ijk)+1.0_rp )*( 1.0_rp/qvtmp + qlsbl(ijk)/rvap/temp(ijk)/temp(ijk)*qlsbl(ijk)/cp )
2926 
2927  a = - cef1*sumliq(ijk)*gtliq(ijk)/dens(ijk) ! a of (A.19) of Suzuki (2004)
2928  b = - cef2*sumice(ijk)*gtice(ijk)/dens(ijk) ! b of (A.19) of Suzuki (2004)
2929  c = - cef3*sumliq(ijk)*gtliq(ijk)/dens(ijk) ! c of (A.19) of Suzuki (2004)
2930  d = - cef4*sumice(ijk)*gtice(ijk)/dens(ijk) ! d of (A.19) of Suzuki (2004)
2931 
2932  b = b + eps * ( 1.0_rp - zerosw ) !--- avoiding division by zero when qv = 0
2933  !--- eigenvalues
2934  rmdplus = ( ( a+d ) + sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
2935  rmdmins = ( ( a+d ) - sqrt( ( a-d )**2 + 4.0_rp*b*c ) ) * 0.50_rp
2936 
2937  rmdplus = rmdplus + eps * ( 1.0_rp - zerosw ) !--- avoiding division by zero when qv = 0
2938  rmdmins = rmdmins + eps * ( 1.0_rp - zerosw ) !--- avoiding division by zero when qv = 0
2939 
2940  !--- supersaturation tendency
2941  ssplus = ( ( rmdmins-a )*ssliq(ijk) - b*ssice(ijk) )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
2942  ssmins = ( ( a-rmdplus )*ssliq(ijk) + b*ssice(ijk) )/b/( rmdmins-rmdplus + eps * ( 1.0_rp - zerosw ) )
2943 
2944  tplus = ( exp( rmdplus*dtcnd(ijk) )-1.0_rp )/( rmdplus*dtcnd(ijk) ) &
2945  * ( 0.5_rp + sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) ) &
2946  + 1.0_rp &
2947  * ( 0.5_rp - sign( 0.5_rp,abs(rmdplus*dtcnd(ijk)-0.1_rp) ) )
2948  tmins = ( exp( rmdmins*dtcnd(ijk) )-1.0_rp )/( rmdmins*dtcnd(ijk) ) &
2949  * ( 0.5_rp + sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) ) &
2950  + 1.0_rp &
2951  * ( 0.5_rp - sign( 0.5_rp,abs(rmdmins*dtcnd(ijk)-0.1_rp) ) )
2952 
2953  sliqtnd = ssliq(ijk) * ( 1.0_rp - zerosw ) &
2954  + zerosw * &
2955  ( b*tplus*ssplus + b*tmins*ssmins ) ! sliwtnd in page 116 of Suzuki (2004)
2956  sicetnd = ssice(ijk) * ( 1.0_rp - zerosw ) &
2957  + zerosw * &
2958  ( ( rmdplus-a )*tplus*ssplus & ! sicetnd in page 116 of Suzuki (2004)
2959  + ( rmdmins-a )*tmins*ssmins )
2960 
2961  !--- change of SDF
2962  do myu = 1, nspc
2963  do mm = 1, iflg( myu,ijk )
2964  !--- advection speed
2965  do n = 1, nbin+1
2966  !--- myu = 1 -> ssliq, myu > 1 -> ssice
2967  uadv( n,myu,ijk ) = cbnd( n,myu )*rexpxbnd( n )*gtliq(ijk)*sliqtnd & ! U of eq. (A.18) of Suzuki (2004)
2968  * ( 0.5_rp - sign( 0.5_rp,real(myu)-1.5_RP) ) &
2969  + cbnd( n,myu )*rexpxbnd( n )*gtice(ijk)*sicetnd & ! U of eq. (A.18) of Suzuki (2004)
2970  * ( 0.5_RP + sign( 0.5_RP,real(myu)-1.5_RP) )
2971  enddo
2972  uadv( 0, myu,ijk ) = 0.0_rp
2973  uadv( nbin+2,myu,ijk ) = 0.0_rp
2974  enddo
2975  enddo
2976  enddo
2977  enddo
2978 
2979  call prof_rapstart('_SBM_AdvMix', 3)
2980 
2981  do ijk = 1, ijkmax
2982  do nn = 1, loopflg(ijk)
2983  do myu = 1, nspc
2984  do mm = 1, iflg( myu,ijk )
2985  do n = 0, nbin+2
2986  crn( n,myu,ijk ) = uadv( n,myu,ijk )*dtcnd(ijk)/dxmic
2987  enddo
2988  enddo
2989  enddo
2990  enddo
2991  enddo
2992 
2993  do ijk = 1, ijkmax
2994  do nn = 1, loopflg(ijk)
2995  do myu = 1, nspc
2996  do mm = 1, iflg( myu,ijk )
2997  do n = 0, nbin+1
2998  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
2999  acoef(1,n,myu,ijk) = ( gcn( n+1,myu,ijk ) -gcn( n-1,myu,ijk ) ) / 16.0_rp
3000  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
3001  enddo
3002  enddo
3003  enddo
3004  enddo
3005  enddo
3006 
3007  do ijk = 1, ijkmax
3008  do nn = 1, loopflg(ijk)
3009  do myu = 1, nspc
3010  do mm = 1, iflg( myu,ijk )
3011  do n = 0, nbin+1
3012  cplus = 1.0_rp - ( crn(n+1,myu,ijk) + abs(crn(n+1,myu,ijk)) )
3013 
3014  aip(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cplus**1 ) &
3015  + acoef(1,n,myu,ijk) * ( 1.0_rp-cplus**2 ) &
3016  + acoef(2,n,myu,ijk) * ( 1.0_rp-cplus**3 )
3017 
3018  aip(n,myu,ijk) = max( aip(n,myu,ijk), 0.0_rp )
3019  enddo
3020  enddo
3021  enddo
3022  enddo
3023  enddo
3024 
3025  do ijk = 1, ijkmax
3026  do nn = 1, loopflg(ijk)
3027  do myu = 1, nspc
3028  do mm = 1, iflg( myu,ijk )
3029  do n = 0, nbin+1
3030  cmins = 1.0_rp - ( abs(crn(n,myu,ijk)) - crn(n,myu,ijk) )
3031 
3032  aim(n,myu,ijk) = acoef(0,n,myu,ijk) * ( 1.0_rp-cmins**1 ) &
3033  - acoef(1,n,myu,ijk) * ( 1.0_rp-cmins**2 ) &
3034  + acoef(2,n,myu,ijk) * ( 1.0_rp-cmins**3 )
3035 
3036  aim(n,myu,ijk) = max( aim(n,myu,ijk), 0.0_rp )
3037  enddo
3038  enddo
3039  enddo
3040  enddo
3041  enddo
3042 
3043  do ijk = 1, ijkmax
3044  do nn = 1, loopflg(ijk)
3045  do myu = 1, nspc
3046  do mm = 1, iflg( myu,ijk )
3047  do n = 0, nbin+1
3048  ai(n,myu,ijk) = acoef(0,n,myu,ijk) * 2.0_rp &
3049  + acoef(2,n,myu,ijk) * 2.0_rp
3050 
3051  ai(n,myu,ijk) = max( ai(n,myu,ijk), aip(n,myu,ijk)+aim(n,myu,ijk)+cldmin )
3052  enddo
3053  enddo
3054  enddo
3055  enddo
3056  enddo
3057 
3058  do ijk = 1, ijkmax
3059  do nn = 1, loopflg(ijk)
3060  do myu = 1, nspc
3061  do mm = 1, iflg( myu,ijk )
3062  do n = 1, nbin+1
3063  flq(n,myu,ijk) = ( aip(n-1,myu,ijk)/ai(n-1,myu,ijk)*gcn( n-1,myu,ijk ) &
3064  - aim(n ,myu,ijk)/ai(n ,myu,ijk)*gcn( n ,myu,ijk ) )*dxmic/dtcnd(ijk)
3065  enddo
3066  enddo
3067  enddo
3068  enddo
3069  enddo
3070 
3071  do ijk = 1, ijkmax
3072  do nn = 1, loopflg(ijk)
3073  do myu = 1, nspc
3074  do mm = 1, iflg( myu,ijk )
3075  dumm_regene(ijk) = dumm_regene(ijk)+( -flq(1,myu,ijk)*dtcnd(ijk)/dxmic ) &
3076  * min( uadv(1,myu,ijk),0.0_rp )/( uadv(1,myu,ijk)+eps )
3077  enddo
3078  enddo
3079  enddo
3080  enddo
3081 
3082  do ijk = 1, ijkmax
3083  do nn = 1, loopflg(ijk)
3084  do myu = 1, nspc
3085  do mm = 1, iflg( myu,ijk )
3086  do n = 1, nbin
3087  gcn( n,myu,ijk ) = gcn( n,myu,ijk ) - ( flq(n+1,myu,ijk)-flq(n,myu,ijk) )*dtcnd(ijk)/dxmic
3088  enddo
3089  enddo
3090  enddo
3091  enddo
3092  enddo
3093 
3094  call prof_rapend ('_SBM_AdvMix', 3)
3095 
3096 !OCL LOOP_NOFUSION
3097  do ijk = 1, ijkmax
3098  do nn = 1, loopflg(ijk)
3099  !--- new mass
3100  gclnew(ijk) = 0.0_rp
3101  do n = 1, nbin
3102  gclnew(ijk) = gclnew(ijk) + gcn( n,il,ijk )*expxctr( n )*dxmic
3103  enddo
3104 
3105  gcinew(ijk) = 0.0_rp
3106  do myu = 2, nspc
3107  do n = 1, nbin
3108  gcinew(ijk) = gcinew(ijk) + gcn( n,myu,ijk )*expxctr( n )*dxmic
3109  enddo
3110  enddo
3111 
3112  !--- change of humidity and temperature
3113  cndmss(ijk) = gclnew(ijk) - gclold(ijk)
3114  sblmss(ijk) = gcinew(ijk) - gciold(ijk)
3115 
3116  qvap(ijk) = qvap(ijk) - ( cndmss(ijk)+sblmss(ijk) )/dens(ijk)
3117  temp(ijk) = temp(ijk) + ( cndmss(ijk)*qlevp(ijk)+sblmss(ijk)*qlsbl(ijk) )/dens(ijk)/cp
3118 
3119  gclold(ijk) = gclnew(ijk)
3120  gciold(ijk) = gcinew(ijk)
3121  enddo
3122  enddo
3123 
3124  enddo ! ncount
3125 
3126 !OCL NORECURRENCE(gc)
3127  do ijk = 1, ijkmax
3128  !----- number -> mass
3129  do myu = 1, nspc
3130  do n = 1, nbin
3131  gc( n,myu,ijk ) = gcn( n,myu,ijk )*expxctr( n )
3132  enddo
3133  enddo
3134  enddo
3135 
3136  if ( .NOT. flg_regeneration ) then
3137  dumm_regene(:) = 0.0_rp
3138  endif
3139 
3140  call prof_rapend ('_SBM_Mixphase', 3)
3141 
3142  return
3143  end subroutine mixphase
3144 
3145  !-----------------------------------------------------------------------------
3146  subroutine ice_nucleat( &
3147  ijkmax, &
3148  num_cold, &
3149  index_cold, &
3150  dens, &
3151  pres, &
3152  temp, &
3153  qvap, &
3154  gc, &
3155  dtime )
3156  implicit none
3157 
3158  integer, intent(in) :: ijkmax
3159  integer, intent(in) :: num_cold
3160  integer, intent(in) :: index_cold(ijkmax)
3161  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
3162  real(RP), intent(in) :: pres(ijkmax) ! Pressure [Pa]
3163  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
3164  real(RP), intent(inout) :: qvap(ijkmax) ! Specific humidity [kg/kg]
3165  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3166  real(DP), intent(in) :: dtime ! Time step interval
3167 
3168  real(RP) :: ssliq, ssice
3169  real(RP) :: numin, tdel, qdel
3170 ! real(RP), parameter :: n0 = 1.E+3_RP ! N_{IN0} in page 112 of Suzuki (2004)
3171  real(RP), parameter :: acoef = -0.639_rp, bcoef = 12.96_rp ! A and B in paeg 112 of Suzuki (2004)
3172  ! threshould to determine the type of freezed hydrometeor (the detail is described in page 113 of Suzuki(2004))
3173  real(RP), parameter :: tcolmu = 269.0_rp, tcolml = 265.0_rp! -4[degC], -8[degC]
3174  real(RP), parameter :: tdendu = 257.0_rp, tdendl = 255.0_rp! -14[degC], -18[degC]
3175  real(RP), parameter :: tplatu = 250.6_rp ! -22.4[degC]
3176  !
3177  real(RP) :: psat
3178  integer :: ijk, indirect
3179 
3180 
3181  call prof_rapstart('_SBM_IceNucleat', 3)
3182 
3183  do indirect = 1, num_cold
3184  ijk = index_cold(indirect)
3185 
3186  !--- supersaturation
3187  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_liq &
3188  * exp( lovr_liq * ( rtem00 - 1.0_rp/temp(ijk) ) )
3189  ssliq = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
3190  psat = const_psat0 * ( temp(ijk) * rtem00 )**cpovr_ice &
3191  * exp( lovr_ice * ( rtem00 - 1.0_rp/temp(ijk) ) )
3192  ssice = const_epsvap * psat / ( pres(ijk) - ( 1.0_rp-const_epsvap ) * psat )
3193  ssliq = qvap(ijk)/ssliq-1.0_rp
3194  ssice = qvap(ijk)/ssice-1.0_rp
3195 
3196  if( ssice <= 0.0_rp ) cycle
3197 
3198  numin = bcoef * exp( acoef + bcoef * ssice )
3199  numin = numin * expxctr( 1 )/dxmic
3200  numin = min( numin,qvap(ijk)*dens(ijk) )
3201  !--- -4 [deg] > T >= -8 [deg] and T < -22.4 [deg] -> column
3202  if ( temp(ijk) <= tplatu .OR. ( temp(ijk) >= tcolml .AND. temp(ijk) < tcolmu ) ) then
3203  gc( 1,ic,ijk ) = gc( 1,ic,ijk ) + numin
3204  !--- -14 [deg] > T >= -18 [deg] -> dendrite
3205  elseif( temp(ijk) <= tdendu .AND. temp(ijk) >= tdendl ) then
3206  gc( 1,id,ijk ) = gc( 1,id,ijk ) + numin
3207  !--- else -> plate
3208  else
3209  gc( 1,ip,ijk ) = gc( 1,ip,ijk ) + numin
3210  endif
3211 
3212  tdel = numin/dens(ijk)*qlmlt/cp
3213  temp(ijk) = temp(ijk) + tdel
3214  qdel = numin/dens(ijk)
3215  qvap(ijk) = qvap(ijk) - qdel
3216 
3217  enddo
3218 
3219 
3220  call prof_rapend ('_SBM_IceNucleat', 3)
3221 
3222  return
3223  end subroutine ice_nucleat
3224 
3225  !-----------------------------------------------------------------------------
3226  subroutine freezing( &
3227  ijkmax, &
3228  num_cold, &
3229  index_cold, &
3230  dens, &
3231  temp, &
3232  gc, &
3233  dtime )
3234  implicit none
3235 
3236  integer, intent(in) :: ijkmax
3237  integer, intent(in) :: num_cold
3238  integer, intent(in) :: index_cold(ijkmax)
3239  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
3240  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
3241  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3242  real(DP), intent(in) :: dtime ! Time step interval
3243 
3244  integer :: nbound, n
3245  real(RP) :: xbound, tc, rate, dmp, frz, sumfrz, tdel
3246  real(RP), parameter :: coefa = 1.0e-01_rp ! a_{fr} of eq.(3.19) of Suzuki (2004)
3247  real(RP), parameter :: coefb = 0.66_rp ! b_{fr} of eq.(3.19) of Suzuki (2004)
3248  real(RP), parameter :: rbound = 2.0e-04_rp ! 200 um
3249 ! real(RP), parameter :: tthreth = 235.0_RP ! -38 [degC] threshold for using Bigg's parameterization
3250 ! real(RP), parameter :: ncoefim = 1.0E+7_RP ! N_{im0} of eq.(3.18) of Suzuki (2004)
3251 ! real(RP), parameter :: gamm = 3.3_RP ! gamma of eq.(3.18) of Suzuki (2004)
3252  integer :: ijk, indirect
3253 
3254  call prof_rapstart('_SBM_Freezing', 3)
3255 
3256  do indirect = 1, num_cold
3257  ijk = index_cold(indirect)
3258 
3259 ! if ( temp <= tthreth ) then !--- Bigg (1975)
3260  xbound = log( rhow * 4.0_rp*pi/3.0_rp * rbound**3 )
3261  nbound = int( ( xbound-xbnd( 1 ) )/dxmic ) + 1
3262 
3263  tc = temp(ijk)-tmlt
3264  rate = coefa*exp( -coefb*tc )
3265 
3266  sumfrz = 0.0_rp
3267  do n = 1, nbound-1
3268  dmp = rate*expxctr( n )
3269  frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
3270 
3271  gc( n,il,ijk ) = gc( n,il,ijk ) - frz
3272  gc( n,il,ijk ) = max( gc( n,il,ijk ),0.0_rp )
3273 
3274  gc( n,ip,ijk ) = gc( n,ip,ijk ) + frz
3275 
3276  sumfrz = sumfrz + frz
3277  enddo
3278  do n = nbound, nbin
3279  dmp = rate*expxctr( n )
3280  frz = gc( n,il,ijk )*( 1.0_rp-exp( -dmp*dtime ) )
3281 
3282  gc( n,il,ijk ) = gc( n,il,ijk ) - frz
3283  gc( n,il,ijk ) = max( gc( n,il,ijk ),0.0_rp )
3284 
3285  gc( n,ih,ijk ) = gc( n,ih,ijk ) + frz
3286 
3287  sumfrz = sumfrz + frz
3288  enddo
3289 ! elseif( temp > tthreth ) then !--- Vali (1975)
3290 !
3291 ! tc = temp-tmlt
3292 ! dmp = ncoefim * ( 0.1_RP * ( tmlt-temp )**gamm )
3293 ! sumfrz = 0.0_RP
3294 ! do n = 1, nbin
3295 ! frz = gc( (il-1)*nbin+n )*dmp*xctr( n )/dens
3296 ! frz = max( frz,gc( (il-1)*nbin+n )*dmp*xctr( n )/dens )
3297 ! gc( (il-1)*nbin+n ) = gc( (il-1)*nbin+n ) - frz
3298 ! if ( n >= nbound ) then
3299 ! gc( (ih-1)*nbin+n ) = gc( (ih-1)*nbin+n ) + frz
3300 ! else
3301 ! gc( (ip-1)*nbin+n ) = gc( (ip-1)*nbin+n ) + frz
3302 ! endif
3303 !
3304 ! sumfrz = sumfrz + frz
3305 ! enddo
3306 ! endif
3307  sumfrz = sumfrz*dxmic
3308 
3309  tdel = sumfrz/dens(ijk)*qlmlt/cp
3310  temp(ijk) = temp(ijk) + tdel
3311  enddo
3312 
3313  call prof_rapend ('_SBM_Freezing', 3)
3314 
3315  return
3316  end subroutine freezing
3317 
3318  !-----------------------------------------------------------------------------
3319  subroutine melting( &
3320  ijkmax, &
3321  num_warm, &
3322  index_warm, &
3323  dens, &
3324  temp, &
3325  gc, &
3326  dtime )
3327  implicit none
3328 
3329  integer, intent(in) :: ijkmax
3330  integer, intent(in) :: num_warm
3331  integer, intent(in) :: index_warm(ijkmax)
3332  real(RP), intent(in) :: dens(ijkmax) ! Density [kg/m3]
3333  real(RP), intent(inout) :: temp(ijkmax) ! Temperature [K]
3334  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3335  real(DP), intent(in) :: dtime ! Time step interval
3336 
3337  integer :: n, m
3338  real(RP) :: summlt, sumice, tdel
3339  integer :: ijk, indirect
3340 
3341  call prof_rapstart('_SBM_Melting', 3)
3342 
3343  do indirect = 1, num_warm
3344  ijk = index_warm(indirect)
3345 
3346  summlt = 0.0_rp
3347  do n = 1, nbin
3348  sumice = 0.0_rp
3349  do m = ic, ih
3350  sumice = sumice + gc( n,m,ijk )
3351  gc( n,m,ijk ) = 0.0_rp
3352  enddo
3353  gc( n,il,ijk ) = gc( n,il,ijk ) + sumice
3354  summlt = summlt + sumice !--- All freezed particle melt instantaneously
3355  enddo
3356  summlt = summlt*dxmic
3357 
3358  tdel = - summlt/dens(ijk)*qlmlt/cp
3359  temp(ijk) = temp(ijk) + tdel
3360  !
3361  enddo
3362 
3363  call prof_rapend ('_SBM_Melting', 3)
3364 
3365  return
3366  end subroutine melting
3367 
3368  !-----------------------------------------------------------------------------
3369  subroutine collmain( &
3370  ijkmax, &
3371  temp, &
3372  ghyd, &
3373  dt )
3374  implicit none
3375 
3376  integer, intent(in) :: ijkmax
3377  real(RP), intent(in) :: temp(ijkmax) ! Temperature [K]
3378  real(RP), intent(inout) :: ghyd(nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3379  real(DP), intent(in) :: dt ! Time step interval
3380  !---------------------------------------------------------------------------
3381 
3382  if ( rndm_flgp == 1 ) then ! stochastic method
3383  call r_collcoag( ijkmax, & ! [IN]
3384  wgtbin, & ! [IN]
3385  temp(:), & ! [IN]
3386  ghyd(:,:,:), & ! [INOUT]
3387  dt ) ! [IN]
3388  else ! default
3389  call collcoag( ijkmax, & ! [IN]
3390  temp(:), & ! [IN]
3391  ghyd(:,:,:), & ! [INOUT]
3392  dt ) ! [IN]
3393  endif
3394 
3395  return
3396  end subroutine collmain
3397 
3398  !-----------------------------------------------------------------------------
3399  subroutine collmainf( &
3400  ijkmax, &
3401  temp, &
3402  ghyd, &
3403  dt )
3404  implicit none
3405 
3406  integer, intent(in) :: ijkmax
3407  real(RP), intent(in) :: temp(ijkmax) ! Temperature [K]
3408  real(RP), intent(inout) :: ghyd(nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3409  real(DP), intent(in) :: dt ! Time step interval
3410  !---------------------------------------------------------------------------
3411 
3412  if ( rndm_flgp == 1 ) then ! stochastic method
3413  call r_collcoag( ijkmax, & ! [IN]
3414  wgtbin, & ! [IN]
3415  temp(:), & ! [IN]
3416  ghyd(:,:,:), & ! [INOUT]
3417  dt ) ! [IN]
3418  else ! default
3419  call collcoag( ijkmax, & ! [IN]
3420  temp(:), & ! [IN]
3421  ghyd(:,:,:), & ! [INOUT]
3422  dt ) ! [IN]
3423  endif
3424 
3425  return
3426  end subroutine collmainf
3427 
3428  !-----------------------------------------------------------------------------
3429  !--- reference paper
3430  ! Bott et al. (1998) J. Atmos. Sci. vol.55, pp. 2284-
3431  subroutine collcoag( &
3432  ijkmax, &
3433  temp, &
3434  gc, &
3435  dtime )
3436  implicit none
3437 
3438  integer, intent(in) :: ijkmax
3439  real(RP), intent(in) :: temp(ijkmax) ! Temperature [K]
3440  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3441  real(DP), intent(in) :: dtime ! Time step interval
3442 
3443  integer :: i, j, k, l
3444  real(RP) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
3445  real(RP) :: gprime, gprimk, wgt, crn, sum, flux
3446  integer, parameter :: ldeg = 2
3447  real(RP) :: acoef( 0:ldeg )
3448  real(RP), parameter :: dmpmin = 1.e-01_rp
3449  real(RP) :: suri, surj
3450 
3451  integer :: myu, n, irsl, ilrg, isml
3452  integer :: ibnd( ijkmax )
3453  integer :: iflg( nspc,ijkmax )
3454  integer :: iexst( nbin,nspc,ijkmax )
3455  real(RP) :: csum( nspc,ijkmax )
3456  integer :: ijk, nn, mm, pp, qq
3457  !---------------------------------------------------------------------------
3458 
3459  call prof_rapstart('_SBM_CollCoag', 3)
3460 
3461  iflg( :,: ) = 0
3462  iexst( :,:,: ) = 0
3463  csum( :,: ) = 0.0_rp
3464  do ijk = 1, ijkmax
3465  !--- judgement of particle existence
3466  do myu = 1, nspc
3467  do n = 1, nbin
3468  csum( myu,ijk ) = csum( myu,ijk ) + gc( n,myu,ijk )*dxmic
3469  enddo
3470  enddo
3471  do myu = 1, nspc
3472  if ( csum( myu,ijk ) > cldmin ) iflg( myu,ijk ) = 1
3473  enddo
3474 
3475  if ( temp(ijk) < tcrit ) then
3476  ibnd(ijk) = 1
3477  else
3478  ibnd(ijk) = 2
3479  endif
3480 
3481  do myu = 1, nspc
3482  do n = 1, nbin
3483  if ( gc( n,myu,ijk ) > cldmin ) then
3484  iexst( n,myu,ijk ) = 1
3485  endif
3486  enddo
3487  enddo
3488 
3489  enddo
3490 
3491 !OCL PARALLEL
3492  do ijk = 1, ijkmax
3493  do isml = 1, nspc
3494  do nn = 1, iflg( isml,ijk )
3495 
3496  do ilrg = 1, nspc
3497  do mm = 1, iflg( ilrg,ijk )
3498  !--- rule of interaction
3499  irsl = ifrsl( ibnd(ijk),isml,ilrg )
3500 
3501  do i = 1, nbin-1 ! small
3502  do pp = 1, iexst( i,isml,ijk )
3503 
3504  do j = i+1, nbin ! large
3505  do qq = 1, iexst( j,ilrg,ijk )
3506 
3507  k = kindx( i,j )
3508  xi = expxctr( i )
3509  xj = expxctr( j )
3510  xnew = log( xi+xj )
3511 
3512  dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime ! dg_{i}/dt*dt in page 119 of Suzuki (2004)
3513  dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime ! dg_{j}/dt*dt in page 119 of Suzuki (2004)
3514 
3515  if ( dmpi <= dmpmin ) then
3516  frci = gc( i,isml,ijk )*dmpi ! Dg_{i} in page 119 of Suzuki (2004)
3517  else
3518  frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) ) ! Dg_{i} in page 119 of Suzuki (2004)
3519  endif
3520 
3521  if ( dmpj <= dmpmin ) then
3522  frcj = gc( j,ilrg,ijk )*dmpj ! Dg_{j} in page 119 of Suzuki (2004)
3523  else
3524  frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) ) ! Dg_{j} in page 119 of Suzuki (2004)
3525  endif
3526 
3527  gprime = frci+frcj
3528 ! if ( gprime <= 0.0_RP ) cycle large
3529  if ( gprime > 0.0_rp .AND. k < nbin ) then
3530 
3531  suri = gc( i,isml,ijk )
3532  surj = gc( j,ilrg,ijk )
3533  gc( i,isml,ijk ) = gc( i,isml,ijk )-frci
3534  gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj
3535  gc( i,isml,ijk ) = max( gc( i,isml,ijk )-frci, 0.0_rp )
3536  gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk )-frcj, 0.0_rp )
3537  frci = suri - gc( i,isml,ijk )
3538  frcj = surj - gc( j,ilrg,ijk )
3539  gprime = frci+frcj ! g' in page 119 of Suzuki (2004)
3540 
3541  gprimk = gc( k,irsl,ijk ) + gprime ! g'_{k} in page 119 of Suzuki (2004)
3542  wgt = gprime / gprimk ! w in page 119 of Suzuki (2004)
3543  crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) ) ! c_{k} in page 119 of Suzuki (2004)
3544 
3545  acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp ! a_{k,0} in page 119 of Suzuki (2004)
3546  acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.50_rp ! a_{k,1} in page 119 of Suzuki (2004)
3547  acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp ! a_{k,2} in page 119 of Suzuki (2004)
3548 
3549  sum = 0.0_rp
3550  do l = 0, ldeg
3551  sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
3552  *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
3553  enddo
3554 
3555  flux = wgt*sum ! f_{k+1/2} in page 119 of Suzuki (2004)
3556  flux = min( max( flux,0.0_rp ),gprime )
3557 
3558  gc( k,irsl,ijk ) = gprimk - flux ! tilda{g_{k}} in page 119 of Suzuki (2004)
3559  gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux ! tilda{g_{k+1}} in page 119 of Suzuki (2004)
3560 
3561  endif
3562 
3563  enddo
3564  enddo !large
3565  enddo
3566  enddo !small
3567  !
3568  enddo
3569  enddo
3570 
3571  enddo
3572  enddo
3573 
3574  enddo
3575 
3576  call prof_rapend ('_SBM_CollCoag', 3)
3577 
3578  return
3579  end subroutine collcoag
3580 
3581  !-------------------------------------------------------------------
3582  subroutine getrule &
3583  ( ifrsl,indx ) !--- out
3584  ! subroutine for creating lookup table (Table A.2 of Suzuki 2004)
3585  integer, intent(out) :: indx( nbin,nbin )
3586  integer, intent(out) :: ifrsl( 2,nspc_mk,nspc_mk )
3587  integer :: i, j, k
3588  real(RP) :: xnew
3589  !
3590  do i = 1, nbin
3591  do j = 1, nbin
3592  xnew = log( expxctr( i )+expxctr( j ) )
3593  k = int( ( xnew-xctr( 1 ) )/dxmic ) + 1
3594  k = max( max( k,j ),i )
3595  indx( i,j ) = k
3596  enddo
3597  enddo
3598  !
3599  !
3600  !--- liquid + liquid -> liquid
3601  ifrsl( 1:2,il,il ) = il
3602  !
3603  !--- liquid + column -> ( graupel, hail ) + column
3604  ifrsl( 1:2,il,ic ) = ic
3605  ifrsl( 1,ic,il ) = ig
3606  ifrsl( 2,ic,il ) = ih
3607  !
3608  !--- liquid + plate -> ( graupel, hail ) + plate
3609  ifrsl( 1:2,il,ip ) = ip
3610  ifrsl( 1,ip,il ) = ig
3611  ifrsl( 2,ip,il ) = ih
3612  !
3613  !--- liquid + dendrite -> ( graupel, hail ) + dendrite
3614  ifrsl( 1:2,il,id ) = id
3615  ifrsl( 1,id,il ) = ig
3616  ifrsl( 2,id,il ) = ih
3617  !
3618  !--- liquid + snowflake -> ( graupel, hail ) + snowflake
3619  ifrsl( 1:2,il,iss ) = iss
3620  ifrsl( 1,iss,il ) = ig
3621  ifrsl( 2,iss,il ) = ih
3622  !
3623  !--- liquid + graupel -> ( graupel, hail )
3624  ifrsl( 1:2,il,ig ) = ig
3625  ifrsl( 1,ig,il ) = ig
3626  ifrsl( 2,ig,il ) = ih
3627  !
3628  !--- liquid + hail -> ( graupel, hail )
3629  ifrsl( 1:2,il,ih ) = ih
3630  ifrsl( 1,ih,il ) = ig
3631  ifrsl( 2,ih,il ) = ih
3632  !
3633  !
3634  !--- column + column -> snowflake
3635  ifrsl( 1:2,ic,ic ) = iss
3636  !
3637  !--- column + plate -> snowflake
3638  ifrsl( 1:2,ic,ip ) = iss
3639  ifrsl( 1:2,ip,ic ) = iss
3640  !
3641  !--- column + dendrite -> snowflake
3642  ifrsl( 1:2,ic,id ) = iss
3643  ifrsl( 1:2,id,ic ) = iss
3644  !
3645  !--- column + snowflake -> snowflake
3646  ifrsl( 1:2,ic,iss ) = iss
3647  ifrsl( 1:2,iss,ic ) = iss
3648  !
3649  !--- column + graupel -> column + graupel
3650  ifrsl( 1:2,ic,ig ) = ig
3651  ifrsl( 1:2,ig,ic ) = ic
3652  !
3653  !--- column + hail -> column + ( graupel, hail )
3654  ifrsl( 1:2,ih,ic ) = ic
3655  ifrsl( 1,ic,ih ) = ig
3656  ifrsl( 2,ic,ih ) = ih
3657  !
3658  !
3659  !--- plate + plate -> snowflake
3660  ifrsl( 1:2,ip,ip ) = iss
3661  !
3662  !--- plate + dendrite -> snowflake
3663  ifrsl( 1:2,ip,id ) = iss
3664  ifrsl( 1:2,id,ip ) = iss
3665  !
3666  !--- plate + snowflake -> snowflake
3667  ifrsl( 1:2,ip,iss ) = iss
3668  ifrsl( 1:2,iss,ip ) = iss
3669  !
3670  !--- plate + graupel -> plate + graupel
3671  ifrsl( 1:2,ip,ig ) = ig
3672  ifrsl( 1:2,ig,ip ) = ip
3673  !
3674  !--- plate + hail -> plate + ( graupel, hail )
3675  ifrsl( 1:2,ih,ip ) = ip
3676  ifrsl( 1,ip,ih ) = ig
3677  ifrsl( 2,ip,ih ) = ih
3678  !
3679  !
3680  !--- dendrite + dendrite -> snowflake
3681  ifrsl( 1:2,id,id ) = iss
3682  !
3683  !--- dendrite + snowflake -> snowflake
3684  ifrsl( 1:2,id,iss ) = iss
3685  ifrsl( 1:2,iss,id ) = iss
3686  !
3687  !--- dendrite + graupel -> dendrite + graupel
3688  ifrsl( 1:2,id,ig ) = ig
3689  ifrsl( 1:2,ig,id ) = id
3690  !
3691  !--- dendrite + hail -> dendrite + ( graupel, hail )
3692  ifrsl( 1:2,ih,id ) = id
3693  ifrsl( 1,id,ih ) = ig
3694  ifrsl( 2,id,ih ) = ih
3695  !
3696  !
3697  !--- snowflake + snowflake -> snowflake
3698  ifrsl( 1:2,iss,iss ) = iss
3699  !
3700  !--- snowflake + graupel -> snowflake + graupel
3701  ifrsl( 1:2,iss,ig ) = ig
3702  ifrsl( 1:2,ig,iss ) = iss
3703  !
3704  !--- snowflake + hail -> snowflake + ( graupel, hail )
3705  ifrsl( 1:2,ih,iss ) = iss
3706  ifrsl( 1,iss,ih ) = ig
3707  ifrsl( 2,iss,ih ) = ih
3708  !
3709  !
3710  !--- graupel + graupel -> graupel
3711  ifrsl( 1:2,ig,ig ) = ig
3712  !
3713  !--- graupel + hail -> ( graupel, hail )
3714  ifrsl( 1,ig,ih ) = ig
3715  ifrsl( 1,ih,ig ) = ig
3716  ifrsl( 2,ig,ih ) = ih
3717  ifrsl( 2,ih,ig ) = ih
3718  !
3719  !--- hail + hail -> hail
3720  ifrsl( 1:2,ih,ih ) = ih
3721  !
3722  !
3723  return
3724  !
3725  end subroutine getrule
3726 
3727  !-----------------------------------------------------------------------------
3728  subroutine faero( &
3729  ijkmax, &
3730  f0, &
3731  ga )
3732  implicit none
3733 
3734  integer , intent(in) :: ijkmax
3735  real(RP), intent(in) :: f0(ijkmax)
3736  real(RP), intent(inout) :: ga(nccn,ijkmax)
3737 
3738 ! real(RP), parameter :: alpha = 3.0_RP
3739 
3740  integer :: ijk, n
3741  !---------------------------------------------------------------------------
3742 
3743  call prof_rapstart('_SBM_FAero', 3)
3744 
3745  do ijk = 1, ijkmax
3746  do n = 1, nccn
3747  ga(n,ijk) = ga(n,ijk) + f0(ijk) * marate(n) * expxactr(n) / dxaer
3748  enddo
3749  enddo
3750 
3751  call prof_rapend ('_SBM_FAero', 3)
3752 
3753  return
3754  end subroutine faero
3755 
3756  !-----------------------------------------------------------------------------
3757  ! + Y. Sato added for stochastic method
3758  ! + Reference Sato et al. (2009) JGR, doi:10.1029/2008JD011247
3759  subroutine random_setup( mset ) !--- in
3760 
3761  use scale_random, only: &
3762  random_get
3763  use scale_process, only: &
3764  prc_mpistop
3765 
3766  integer, intent(in) :: mset
3767 
3768  !--- local ----
3769  integer :: n
3770  real(RP) :: nbinr, tmp1
3771  real(RP) :: rans( mbin ), ranl( mbin )
3772  integer :: pq
3773  real(RP), allocatable :: ranstmp( : )
3774  real(RP), allocatable :: ranltmp( : )
3775  integer :: p, q
3776  integer :: k, temp
3777  integer, allocatable :: orderb( : )
3778  real(RP) :: abq1
3779  real(RP) :: a
3780  real(RP), allocatable :: randnum(:,:,:)
3781  !-------------------------------------------------------
3782  pq = nbin*(nbin-1)/2
3783  allocate( blrg( mset, mbin ) )
3784  allocate( bsml( mset, mbin ) )
3785  allocate( ranstmp( pq ) )
3786  allocate( ranltmp( pq ) )
3787  allocate( orderb( pq ) )
3788  allocate( randnum(1,1,pq) )
3789 
3790  a = real( nbin )*real( nbin-1 )*0.50_RP
3791  if ( a < mbin ) then
3792  write(*,*) "mbin should be smaller than {nbin}_C_{2}"
3793  call prc_mpistop
3794  endif
3795 
3796  wgtbin = a/real( mbin )
3797  nbinr = real( nbin )
3798 
3799  do p = 1, pq
3800  orderb( p ) = p
3801  enddo
3802 
3803  do p = 1, nbin-1
3804  ranstmp( (p-1)*nbin-(p*(p-1))/2+1 : p*nbin-(p*(p+1))/2 ) = p
3805  do q = 1, nbin-p
3806  ranltmp( (p-1)*nbin-(p*(p-1))/2+q ) = p+q
3807  enddo
3808  enddo
3809 
3810  do n = 1, mset
3811  call random_get( randnum )
3812  do p = 1, pq
3813  abq1 = randnum( 1,1,p )
3814  k = int( abq1*( pq-p-1 ) ) + p
3815  temp = orderb( p )
3816  orderb( p ) = orderb( k )
3817  orderb( k ) = temp
3818  enddo
3819 
3820  do p = 1, mbin
3821  if ( p <= pq ) then
3822  rans( p ) = ranstmp( orderb( p ) )
3823  ranl( p ) = ranltmp( orderb( p ) )
3824  else
3825  rans( p ) = ranstmp( orderb( p-pq ) )
3826  ranl( p ) = ranltmp( orderb( p-pq ) )
3827  endif
3828  if ( rans( p ) >= ranl( p ) ) then
3829  tmp1 = rans( p )
3830  rans( p ) = ranl( p )
3831  ranl( p ) = tmp1
3832  endif
3833  enddo
3834  blrg( n,1:mbin ) = int( ranl( 1:mbin ) )
3835  bsml( n,1:mbin ) = int( rans( 1:mbin ) )
3836  enddo
3837 
3838  deallocate( ranstmp )
3839  deallocate( ranltmp )
3840  deallocate( orderb )
3841  deallocate( randnum )
3842 
3843  end subroutine random_setup
3844 
3845  !-----------------------------------------------------------------------------
3846  !--- reference paper
3847  ! Bott et al. (1998) J. Atmos. Sci. vol.55, pp. 2284-
3848  ! Bott et al. (2000) J. Atmos. Sci. Vol.57, pp. 284-
3849  subroutine r_collcoag( &
3850  ijkmax, &
3851  swgt, &
3852  temp, &
3853  gc, &
3854  dtime )
3855  use scale_random, only: &
3856  random_get
3857  implicit none
3858 
3859  integer, intent(in) :: ijkmax
3860  real(RP), intent(in) :: swgt
3861  real(RP), intent(in) :: temp(ijkmax) ! Temperature [K]
3862  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3863  real(DP), intent(in) :: dtime ! Time step interval
3864 
3865  integer :: i, j, k, l
3866  real(RP) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
3867  real(RP) :: gprime, gprimk, wgt, crn, sum, flux
3868  integer, parameter :: ldeg = 2
3869  real(RP), parameter :: dmpmin = 1.e-01_rp, cmin = 1.e-10_rp
3870  real(RP) :: acoef( 0:ldeg )
3871  !
3872  !--- Y.sato added to use code6
3873  integer :: nums( mbin ), numl( mbin )
3874  real(RP), parameter :: gt = 1.0_rp
3875  integer :: s, det
3876  real(RP) :: nbinr, mbinr ! use to weight
3877 ! real(RP) :: beta
3878  real(RP) :: tmpi, tmpj
3879 
3880  integer :: ibnd( ijkmax )
3881  integer :: iflg( nspc,ijkmax )
3882  integer :: iexst( nbin,nspc,ijkmax )
3883  real(RP) :: csum( nspc,ijkmax )
3884  integer :: ijk, nn, mm, pp, qq, myu, n, isml, ilrg, irsl
3885  !---------------------------------------------------------------------------
3886 
3887  call prof_rapstart('_SBM_CollCoagR', 3)
3888 
3889  iflg( :,: ) = 0
3890  iexst( :,:,: ) = 0
3891  csum( :,: ) = 0.0_rp
3892  do ijk = 1, ijkmax
3893  !--- judgement of particle existence
3894  do n = 1, nbin
3895  csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
3896  csum( ic,ijk ) = csum( ic,ijk ) + gc( n,ic,ijk )*dxmic
3897  csum( ip,ijk ) = csum( ip,ijk ) + gc( n,ip,ijk )*dxmic
3898  csum( id,ijk ) = csum( id,ijk ) + gc( n,id,ijk )*dxmic
3899  csum( iss,ijk ) = csum( iss,ijk ) + gc( n,iss,ijk )*dxmic
3900  csum( ig,ijk ) = csum( ig,ijk ) + gc( n,ig,ijk )*dxmic
3901  csum( ih,ijk ) = csum( ih,ijk ) + gc( n,ih,ijk )*dxmic
3902  enddo
3903  if ( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
3904  if ( csum( ic,ijk ) > cldmin ) iflg( ic,ijk ) = 1
3905  if ( csum( ip,ijk ) > cldmin ) iflg( ip,ijk ) = 1
3906  if ( csum( id,ijk ) > cldmin ) iflg( id,ijk ) = 1
3907  if ( csum( iss,ijk ) > cldmin ) iflg( iss,ijk ) = 1
3908  if ( csum( ig,ijk ) > cldmin ) iflg( ig,ijk ) = 1
3909  if ( csum( ih,ijk ) > cldmin ) iflg( ih,ijk ) = 1
3910 
3911  if ( temp(ijk) < tcrit ) then
3912  ibnd(ijk) = 1
3913  else
3914  ibnd(ijk) = 2
3915  endif
3916 
3917  do myu = 1, nspc
3918  do n = 1, nbin
3919  if ( gc( n,myu,ijk ) > cldmin ) then
3920  iexst( n,myu,ijk ) = 1
3921  endif
3922  enddo
3923  enddo
3924 
3925  enddo
3926 
3927 !OCL PARALLEL
3928  do ijk = 1, ijkmax
3929  do isml = 1, nspc
3930  do nn = 1, iflg( isml,ijk )
3931 
3932  do ilrg = 1, nspc
3933  do mm = 1, iflg( ilrg,ijk )
3934  !--- rule of interaction
3935  irsl = ifrsl( ibnd(ijk),isml,ilrg )
3936 
3937  call random_get( rndm )
3938  det = int( rndm(1,1,1)*ia*ja*ka )
3939  nbinr = real( nbin )
3940  mbinr = real( mbin )
3941  nums( 1:mbin ) = bsml( det,1:mbin )
3942  numl( 1:mbin ) = blrg( det,1:mbin )
3943 
3944  do s = 1, mbin
3945  i = nums( s )
3946  j = numl( s )
3947 
3948  do pp = 1, iexst( i,isml,ijk )
3949  do qq = 1, iexst( j,ilrg,ijk )
3950 
3951  k = kindx( i,j )
3952  xi = expxctr( i )
3953  xj = expxctr( j )
3954  xnew = log( xi+xj )
3955 
3956  dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
3957  dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
3958 
3959  if ( dmpi <= dmpmin ) then
3960  frci = gc( i,isml,ijk )*dmpi
3961  else
3962  frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
3963  endif
3964 
3965  if ( dmpj <= dmpmin ) then
3966  frcj = gc( j,ilrg,ijk )*dmpj
3967  else
3968  frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
3969  endif
3970  tmpi = gc( i,isml,ijk )
3971  tmpj = gc( j,ilrg,ijk )
3972 
3973  gc( i,isml,ijk ) = gc( i,isml,ijk )-frci*swgt
3974  gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj*swgt
3975 
3976  if ( j /= k ) then
3977  gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk ), 0.0_rp )
3978  endif
3979  gc( i,isml,ijk ) = max( gc( i,isml,ijk ), 0.0_rp )
3980 
3981  frci = tmpi - gc( i,isml,ijk )
3982  frcj = tmpj - gc( j,ilrg,ijk )
3983 
3984  gprime = frci+frcj
3985 
3986  !-----------------------------------------------
3987  !--- Exponential Flux Method (Bott, 2000, JAS)
3988  !-----------------------------------------------
3989  ! if ( gprime <= 0.0_RP ) cycle !large
3990  ! if ( gprime > 0.0_RP .AND. k < nbin ) then
3991  ! gprimk = gc( (irsl-1)*nbin+k ) + gprime
3992  !
3993  ! beta = log( gc( (irsl-1)*nbin+k+1 )/gprimk+1.E-60_RP )
3994  ! crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
3995  !
3996  ! flux = ( gprime/beta )*( exp( beta*0.50_RP ) -exp( beta*( 0.50_RP-crn ) ) )
3997  ! flux = min( gprimk ,gprime )
3998  !
3999  ! gc( (irsl-1)*nbin+k ) = gprimk - flux
4000  ! gc( (irsl-1)*nbin+k+1 ) = gc( (irsl-1)*nbin+k+1 ) + flux
4001  ! endif
4002 
4003  !-----------------------------------------------
4004  !--- Flux Method (Bott, 1998, JAS)
4005  !-----------------------------------------------
4006 ! if ( gprime <= 0.0_RP ) cycle !large
4007  if ( gprime > 0.0_rp .AND. k < nbin ) then
4008 
4009  gprimk = gc( k,irsl,ijk ) + gprime
4010  wgt = gprime / gprimk
4011  crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
4012 
4013  acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
4014  acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.5_rp
4015  acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
4016 
4017  sum = 0.0_rp
4018  do l = 0, ldeg
4019  sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
4020  *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
4021  enddo
4022 
4023  flux = wgt*sum
4024  flux = min( max( flux,0.0_rp ),gprime )
4025 
4026  gc( k,irsl,ijk ) = gprimk - flux
4027  gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
4028  endif
4029 
4030  enddo
4031  enddo
4032 
4033  enddo ! bin
4034  !
4035  enddo
4036  enddo
4037 
4038  enddo
4039  enddo
4040 
4041  enddo
4042 
4043  call prof_rapend ('_SBM_CollCoagR', 3)
4044 
4045  return
4046  end subroutine r_collcoag
4047 
4048  !-----------------------------------------------------------------------------
4050  subroutine atmos_phy_mp_suzuki10_cloudfraction( &
4051  cldfrac, &
4052  QTRC, &
4053  mask_criterion )
4054  use scale_grid_index
4055  use scale_tracer, only: &
4056  qa
4057  implicit none
4058 
4059  real(RP), intent(out) :: cldfrac(ka,ia,ja)
4060  real(RP), intent(in) :: qtrc (ka,ia,ja,qa)
4061  real(RP), intent(in) :: mask_criterion
4062 
4063  real(RP) :: qhydro
4064  integer :: k, i, j, iq, ihydro
4065  !---------------------------------------------------------------------------
4066 
4067  if( nspc > 1 ) then
4068  do j = js, je
4069  do i = is, ie
4070  do k = ks, ke
4071  qhydro = 0.0_rp
4072  do ihydro = 1, nspc
4073  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4074  qhydro = qhydro + qtrc(k,i,j,iq)
4075  enddo
4076  enddo
4077  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
4078  enddo
4079  enddo
4080  enddo
4081  elseif( nspc == 1 ) then
4082  do j = js, je
4083  do i = is, ie
4084  do k = ks, ke
4085  qhydro = 0.0_rp
4086  do ihydro = 1, i_mp_qc
4087  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4088  qhydro = qhydro + qtrc(k,i,j,iq)
4089  enddo
4090  enddo
4091  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
4092  enddo
4093  enddo
4094  enddo
4095  endif
4096 
4097  return
4098  end subroutine atmos_phy_mp_suzuki10_cloudfraction
4099 
4100  !-----------------------------------------------------------------------------
4102  subroutine atmos_phy_mp_suzuki10_effectiveradius( &
4103  Re, &
4104  QTRC0, &
4105  DENS0, &
4106  TEMP0 )
4107  use scale_grid_index
4108  use scale_const, only: &
4109  eps => const_eps
4110  use scale_tracer, only: &
4111  qa
4112  use scale_atmos_hydrometeor, only: &
4113  n_hyd, &
4114  i_hc, &
4115  i_hr, &
4116  i_hi, &
4117  i_hs, &
4118  i_hg, &
4119  i_hh
4120  implicit none
4121 
4122  real(RP), intent(out) :: re (ka,ia,ja,n_hyd) ! effective radius [cm]
4123  real(RP), intent(in) :: qtrc0(ka,ia,ja,qa) ! tracer mass concentration [kg/kg]
4124  real(RP), intent(in) :: dens0(ka,ia,ja) ! density [kg/m3]
4125  real(RP), intent(in) :: temp0(ka,ia,ja) ! temperature [K]
4126 
4127  real(RP), parameter :: um2cm = 100.0_rp
4128 
4129  real(RP) :: sum0(nspc), sum2, sum3, re_tmp(nspc)
4130  integer :: i, j, k, iq, ihydro
4131  !---------------------------------------------------------------------------
4132 
4133  do k = ks, ke
4134  do j = js, je
4135  do i = is, ie
4136  re(k,i,j,:) = 0.0_rp
4137 
4138  ! HC
4139  sum3 = 0.0_rp
4140  sum2 = 0.0_rp
4141  ihydro = i_mp_qc
4142  do iq = qs_mp+1, qs_mp+nbnd
4143  sum3 = sum3 &
4144  + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
4145  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
4146  * radc( iq-(i_qv+nbin*(ihydro-1)) )**3.0_rp )
4147  sum2 = sum2 &
4148  + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
4149  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
4150  * radc( iq-(i_qv+nbin*(ihydro-1)) )**2.0_rp )
4151  enddo
4152  sum3 = max( sum3, 0.0_rp )
4153  sum2 = max( sum2, 0.0_rp )
4154  if ( sum2 /= 0.0_rp ) then
4155  re(k,i,j,i_hc) = sum3 / sum2 * um2cm
4156  else
4157  re(k,i,j,i_hc) = 0.0_rp
4158  endif
4159 
4160  ! HR
4161  sum3 = 0.0_rp
4162  sum2 = 0.0_rp
4163  ihydro = i_mp_qc
4164  do iq = qs_mp+nbnd+1, qs_mp+nbin
4165  sum3 = sum3 &
4166  + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
4167  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
4168  * radc( iq-(i_qv+nbin*(ihydro-1)) )**3.0_rp )
4169  sum2 = sum2 &
4170  + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
4171  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
4172  * radc( iq-(i_qv+nbin*(ihydro-1)) )**2.0_rp )
4173  enddo
4174  sum3 = max( sum3, 0.0_rp )
4175  sum2 = max( sum2, 0.0_rp )
4176  if ( sum2 /= 0.0_rp ) then
4177  re(k,i,j,i_hr) = sum3 / sum2 * um2cm
4178  else
4179  re(k,i,j,i_hr) = 0.0_rp
4180  endif
4181 
4182  enddo
4183  enddo
4184  enddo
4185 
4186  ! other hydrometeors
4187  if ( nspc > 1 ) then
4188  do k = ks, ke
4189  do j = js, je
4190  do i = is, ie
4191  do ihydro = 2, nspc
4192  sum0(ihydro) = 0.0_rp
4193  sum2 = 0.0_rp
4194  sum3 = 0.0_rp
4195  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4196  sum0(ihydro) = sum0(ihydro) &
4197  + ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) !--- [kg/kg] -> [kg/m3]
4198  sum3 = sum3 &
4199  + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
4200  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
4201  * radc( iq-(i_qv+nbin*(ihydro-1)) )**3.0_rp )
4202  sum2 = sum2 &
4203  + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
4204  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
4205  * radc( iq-(i_qv+nbin*(ihydro-1)) )**2.0_rp )
4206  enddo
4207  sum3 = max( sum3, 0.0_rp )
4208  sum2 = max( sum2, 0.0_rp )
4209  if ( sum2 == 0.0_rp ) then
4210  re_tmp(ihydro) = 0.0_rp
4211  else
4212  re_tmp(ihydro) = sum3 / sum2 * um2cm
4213  end if
4214  end do
4215 
4216  re(k,i,j,i_hi) = ( re_tmp(i_mp_qp ) * sum0(i_mp_qp ) &
4217  + re_tmp(i_mp_qcl) * sum0(i_mp_qcl) &
4218  + re_tmp(i_mp_qd ) * sum0(i_mp_qd ) ) &
4219  / ( sum0(i_mp_qp) + sum0(i_mp_qcl) + sum0(i_mp_qd) + eps )
4220  re(k,i,j,i_hs) = re_tmp(i_mp_qs)
4221  re(k,i,j,i_hg) = re_tmp(i_mp_qg)
4222  re(k,i,j,i_hh) = re_tmp(i_mp_qh)
4223 
4224  enddo
4225  enddo
4226  enddo
4227 
4228  end if
4229 
4230  return
4231  end subroutine atmos_phy_mp_suzuki10_effectiveradius
4232 
4233  !-----------------------------------------------------------------------------
4235  subroutine atmos_phy_mp_suzuki10_mixingratio( &
4236  Qe, &
4237  QTRC0 )
4238  use scale_grid_index
4239  use scale_const, only: &
4240  eps => const_eps
4241  use scale_tracer, only: &
4242  qa
4243  use scale_atmos_hydrometeor, only: &
4244  n_hyd, &
4245  i_hc, &
4246  i_hr, &
4247  i_hi, &
4248  i_hs, &
4249  i_hg, &
4250  i_hh
4251  implicit none
4252 
4253  real(RP), intent(out) :: qe (ka,ia,ja,n_hyd) ! mixing ratio of each cateory [kg/kg]
4254  real(RP), intent(in) :: qtrc0(ka,ia,ja,qa) ! tracer mass concentration [kg/kg]
4255 
4256  real(RP) :: sum2
4257  integer :: i, j, k, iq, ihydro
4258  !---------------------------------------------------------------------------
4259 
4260 
4261  do k = ks, ke
4262  do j = js, je
4263  do i = is, ie
4264 
4265  qe(k,i,j,:) = 0.0_rp
4266 
4267  ! HC
4268  sum2 = 0.0_rp
4269  do iq = qs_mp+1, qs_mp+nbnd
4270  sum2 = sum2 + qtrc0(k,i,j,iq)
4271  enddo
4272  qe(k,i,j,i_hc) = sum2
4273 
4274  ! HR
4275  sum2 = 0.0_rp
4276  do iq = qs_mp+nbnd+1, qs_mp+nbin
4277  sum2 = sum2 + qtrc0(k,i,j,iq)
4278  enddo
4279  qe(k,i,j,i_hr) = sum2
4280 
4281  enddo
4282  enddo
4283  enddo
4284 
4285  ! other hydrometeors
4286  if ( nspc > 1 ) then
4287 
4288  do k = ks, ke
4289  do j = js, je
4290  do i = is, ie
4291  ! HI
4292  sum2 = 0.0_rp
4293  do ihydro = i_mp_qp, i_mp_qd
4294  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4295  sum2 = sum2 + qtrc0(k,i,j,iq)
4296  enddo
4297  enddo
4298  qe(k,i,j,i_hi) = sum2
4299 
4300  ! HS
4301  sum2 = 0.0_rp
4302  ihydro = i_mp_qs
4303  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4304  sum2 = sum2 + qtrc0(k,i,j,iq)
4305  enddo
4306  qe(k,i,j,i_hs) = sum2
4307 
4308  ! HG
4309  sum2 = 0.0_rp
4310  ihydro = i_mp_qg
4311  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4312  sum2 = sum2 + qtrc0(k,i,j,iq)
4313  enddo
4314  qe(k,i,j,i_hg) = sum2
4315 
4316  ! HS
4317  sum2 = 0.0_rp
4318  ihydro = i_mp_qh
4319  do iq = qs_mp+nbin*(ihydro-1)+1, qs_mp+nbin*ihydro
4320  sum2 = sum2 + qtrc0(k,i,j,iq)
4321  enddo
4322  qe(k,i,j,i_hh) = sum2
4323 
4324  enddo
4325  enddo
4326  enddo
4327  endif
4328 
4329  return
4330  end subroutine atmos_phy_mp_suzuki10_mixingratio
4331 
4332  !-----------------------------------------------------------------------------
4333  !----- mkpara is module to create micpara.dat, which is parameter file of
4334  !----- micrphyisical proprties of hydrometeors (collision kernel, radius...).
4335  !----- Imported from preprocess/mk_para2 at 2013/12/26 (Y.Sato)
4336  subroutine mkpara
4337 
4338  implicit none
4339 
4340  integer :: i, j
4341  !-----------------------------------------------------------------------------
4342 
4343  allocate( radc_mk( nbin ) )
4344  allocate( xctr_mk( nbin ) )
4345  allocate( xbnd_mk( nbin+1 ) )
4346  allocate( cctr_mk( nspc_mk,nbin ) )
4347  allocate( cbnd_mk( nspc_mk,nbin+1 ) )
4348  allocate( ck_mk( nspc_mk,nspc_mk,nbin,nbin ) )
4349  allocate( vt_mk( nspc_mk,nbin ) )
4350  allocate( br_mk( nspc_mk,nbin ) )
4351 
4352  !--- file reading
4353  call rdkdat
4354 
4355  !--- grid setting
4356  call sdfgrid
4357 
4358  !--- capacity
4359  call getcp
4360 
4361  !--- collection kernel
4362  call getck
4363 
4364  !--- terminal velocity
4365  call getvt
4366 
4367  !--- bulk radius
4368  call getbr
4369 
4370  !--- output
4371  call paraout
4372 
4373  deallocate( radc_mk )
4374  deallocate( xctr_mk )
4375  deallocate( xbnd_mk )
4376  deallocate( cctr_mk )
4377  deallocate( cbnd_mk )
4378  deallocate( ck_mk )
4379  deallocate( vt_mk )
4380  deallocate( br_mk )
4381 
4382  end subroutine mkpara
4383 
4384  !---------------------------------------------------------------------------------------
4385  subroutine rdkdat
4386 
4387  implicit none
4388  integer, parameter :: il = 1, ic = 2, ip = 3, id = 4
4389  integer, parameter :: is = 5, ig = 6, ih = 7
4390  integer, parameter :: icemax = 3
4391 
4392  integer :: k, kk, i, j
4393  real(DP) :: xl( ndat ), rlec( ndat ), vrl( ndat )
4394  real(DP) :: blkradl( ndat ), blkdnsl( ndat )
4395 
4396  real(DP) :: xi( ndat,icemax ), riec( ndat,icemax ), vri( ndat,icemax )
4397  real(DP) :: blkradi( ndat,icemax ), blkdnsi( ndat,icemax )
4398 
4399  real(DP) :: xs( ndat ), rsec( ndat ), vrs( ndat )
4400  real(DP) :: blkrads( ndat ), blkdnss( ndat )
4401 
4402  real(DP) :: xg( ndat ), rgec( ndat ), vrg( ndat )
4403  real(DP) :: blkradg( ndat ), blkdnsg( ndat )
4404 
4405  real(DP) :: xh( ndat ), rhec( ndat ), vrh( ndat )
4406  real(DP) :: blkradh( ndat ), blkdnsh( ndat )
4407 
4408  data xl(1:ndat) / &
4409  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4410  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4411  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4412  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4413  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4414  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4415  data xi(1:ndat,1) / &
4416  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4417  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4418  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4419  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4420  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4421  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4422  data xi(1:ndat,2) / &
4423  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4424  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4425  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4426  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4427  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4428  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00 /
4429  data xi(1:ndat,3) / &
4430  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4431  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4432  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4433  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4434  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4435  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4436  data xs(1:ndat) / &
4437  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4438  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4439  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4440  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4441  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4442  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4443  data xg(1:ndat) / &
4444  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4445  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4446  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4447  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4448  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4449  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4450  data xh(1:ndat) / &
4451  0.33510e-10_dp,0.67021e-10_dp,0.13404e-09_dp,0.26808e-09_dp,0.53617e-09_dp,0.10723e-08_dp, &
4452  0.21447e-08_dp,0.42893e-08_dp,0.85786e-08_dp,0.17157e-07_dp,0.34315e-07_dp,0.68629e-07_dp, &
4453  0.13726e-06_dp,0.27452e-06_dp,0.54903e-06_dp,0.10981e-05_dp,0.21961e-05_dp,0.43923e-05_dp, &
4454  0.87845e-05_dp,0.17569e-04_dp,0.35138e-04_dp,0.70276e-04_dp,0.14055e-03_dp,0.28110e-03_dp, &
4455  0.56221e-03_dp,0.11244e-02_dp,0.22488e-02_dp,0.44977e-02_dp,0.89954e-02_dp,0.17991e-01_dp, &
4456  0.35981e-01_dp,0.71963e-01_dp,0.14393e+00_dp /
4457  data rlec(1:ndat) / &
4458  0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
4459  0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
4460  0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
4461  0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
4462  0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
4463  0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
4464  data riec(1:ndat,1) / &
4465  0.31936e-03_dp,0.40397e-03_dp,0.51099e-03_dp,0.64638e-03_dp,0.81764e-03_dp,0.10343e-02_dp, &
4466  0.13084e-02_dp,0.16551e-02_dp,0.20937e-02_dp,0.26486e-02_dp,0.33506e-02_dp,0.42387e-02_dp, &
4467  0.64360e-02_dp,0.81426e-02_dp,0.10302e-01_dp,0.13035e-01_dp,0.16494e-01_dp,0.20872e-01_dp, &
4468  0.26412e-01_dp,0.33426e-01_dp,0.42304e-01_dp,0.53543e-01_dp,0.67770e-01_dp,0.85783e-01_dp, &
4469  0.10859e+00_dp,0.13746e+00_dp,0.17403e+00_dp,0.22032e+00_dp,0.27895e+00_dp,0.35319e+00_dp, &
4470  0.44722e+00_dp,0.56630e+00_dp,0.71712e+00_dp /
4471  data riec(1:ndat,2) / &
4472  0.13188e-03_dp,0.16615e-03_dp,0.20953e-03_dp,0.27728e-03_dp,0.36694e-03_dp,0.48559e-03_dp, &
4473  0.64261e-03_dp,0.85040e-03_dp,0.11254e-02_dp,0.14893e-02_dp,0.19709e-02_dp,0.26082e-02_dp, &
4474  0.34515e-02_dp,0.45676e-02_dp,0.60446e-02_dp,0.79991e-02_dp,0.10586e-01_dp,0.14009e-01_dp, &
4475  0.18539e-01_dp,0.24533e-01_dp,0.32466e-01_dp,0.42964e-01_dp,0.56857e-01_dp,0.75242e-01_dp, &
4476  0.99573e-01_dp,0.13177e+00_dp,0.17438e+00_dp,0.23077e+00_dp,0.30539e+00_dp,0.40414e+00_dp, &
4477  0.53482e+00_dp,0.70775e+00_dp,0.93661e+00_dp /
4478  data riec(1:ndat,3) / 0.14998e-03_dp,0.18896e-03_dp,0.23808e-03_dp, &
4479  0.29996e-03_dp,0.37793e-03_dp,0.47616e-03_dp,0.61048e-03_dp,0.81343e-03_dp,0.10839e-02_dp, &
4480  0.14442e-02_dp,0.19243e-02_dp,0.25640e-02_dp,0.34164e-02_dp,0.45522e-02_dp,0.60656e-02_dp, &
4481  0.80820e-02_dp,0.10769e-01_dp,0.14349e-01_dp,0.19119e-01_dp,0.25475e-01_dp,0.44576e-01_dp, &
4482  0.62633e-01_dp,0.88006e-01_dp,0.12366e+00_dp,0.17375e+00_dp,0.24414e+00_dp,0.34304e+00_dp, &
4483  0.48201e+00_dp,0.67728e+00_dp,0.95164e+00_dp,0.13372e+01_dp,0.18788e+01_dp,0.26400e+01_dp /
4484  data rsec(1:ndat) / &
4485  0.92832e-03_dp,0.11696e-02_dp,0.14736e-02_dp,0.18566e-02_dp,0.23392e-02_dp,0.29472e-02_dp, &
4486  0.37133e-02_dp,0.46784e-02_dp,0.58944e-02_dp,0.74265e-02_dp,0.93569e-02_dp,0.11789e-01_dp, &
4487  0.14853e-01_dp,0.18714e-01_dp,0.23578e-01_dp,0.29706e-01_dp,0.37427e-01_dp,0.47156e-01_dp, &
4488  0.59412e-01_dp,0.74855e-01_dp,0.94311e-01_dp,0.11882e+00_dp,0.14971e+00_dp,0.18862e+00_dp, &
4489  0.23765e+00_dp,0.29942e+00_dp,0.37724e+00_dp,0.47530e+00_dp,0.59884e+00_dp,0.75449e+00_dp, &
4490  0.95060e+00_dp,0.11977e+01_dp,0.15090e+01_dp /
4491  data rgec(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
4492  0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
4493  0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
4494  0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
4495  0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
4496  0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
4497  data rhec(1:ndat) / &
4498  0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
4499  0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
4500  0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
4501  0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
4502  0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
4503  0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
4504  data vrl(1:ndat) / &
4505  0.50000e-01_dp,0.78000e-01_dp,0.12000e+00_dp,0.19000e+00_dp,0.31000e+00_dp,0.49000e+00_dp, &
4506  0.77000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp,0.48000e+01_dp,0.74000e+01_dp, &
4507  0.11000e+02_dp,0.17000e+02_dp,0.26000e+02_dp,0.37000e+02_dp,0.52000e+02_dp,0.71000e+02_dp, &
4508  0.94000e+02_dp,0.12000e+03_dp,0.16000e+03_dp,0.21000e+03_dp,0.26000e+03_dp,0.33000e+03_dp, &
4509  0.41000e+03_dp,0.48000e+03_dp,0.57000e+03_dp,0.66000e+03_dp,0.75000e+03_dp,0.82000e+03_dp, &
4510  0.88000e+03_dp,0.90000e+03_dp,0.90000e+03_dp /
4511  data vri(1:ndat,1) / 0.30000e-01_dp,0.40000e-01_dp,0.60000e-01_dp, &
4512  0.80000e-01_dp,0.11000e+00_dp,0.15000e+00_dp,0.17000e+00_dp,0.18000e+00_dp,0.20000e+00_dp, &
4513  0.25000e+00_dp,0.40000e+00_dp,0.60000e+01_dp,0.10000e+02_dp,0.15000e+02_dp,0.20000e+02_dp, &
4514  0.25000e+02_dp,0.31000e+02_dp,0.37000e+02_dp,0.41000e+02_dp,0.46000e+02_dp,0.51000e+02_dp, &
4515  0.55000e+02_dp,0.59000e+02_dp,0.62000e+02_dp,0.64000e+02_dp,0.67000e+02_dp,0.68000e+02_dp, &
4516  0.69000e+02_dp,0.70000e+02_dp,0.71000e+02_dp,0.71500e+02_dp,0.71750e+02_dp,0.72000e+02_dp /
4517  data vri(1:ndat,2) / &
4518  0.30000e-01_dp,0.40000e-01_dp,0.50000e-01_dp,0.70000e-01_dp,0.90000e-01_dp,0.12000e+00_dp, &
4519  0.50000e+00_dp,0.80000e+00_dp,0.16000e+01_dp,0.18000e+01_dp,0.20000e+01_dp,0.30000e+01_dp, &
4520  0.40000e+01_dp,0.50000e+01_dp,0.80000e+01_dp,0.13000e+02_dp,0.19000e+02_dp,0.26000e+02_dp, &
4521  0.32000e+02_dp,0.38000e+02_dp,0.47000e+02_dp,0.55000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
4522  0.77000e+02_dp,0.79000e+02_dp,0.80000e+02_dp,0.81000e+02_dp,0.81000e+02_dp,0.82000e+02_dp, &
4523  0.82000e+02_dp,0.82000e+02_dp,0.82000e+02_dp /
4524  data vri(1:ndat,3) / 0.35000e-01_dp,0.45000e-01_dp,0.55000e-01_dp, &
4525  0.75000e-01_dp,0.95000e-01_dp,0.13000e+00_dp,0.60000e+00_dp,0.90000e+00_dp,0.17000e+01_dp, &
4526  0.20000e+01_dp,0.25000e+01_dp,0.38000e+01_dp,0.50000e+01_dp,0.70000e+01_dp,0.90000e+01_dp, &
4527  0.11000e+02_dp,0.14000e+02_dp,0.17000e+02_dp,0.21000e+02_dp,0.25000e+02_dp,0.32000e+02_dp, &
4528  0.38000e+02_dp,0.44000e+02_dp,0.49000e+02_dp,0.53000e+02_dp,0.55000e+02_dp,0.58000e+02_dp, &
4529  0.59000e+02_dp,0.61000e+02_dp,0.62000e+02_dp,0.63000e+02_dp,0.64000e+02_dp,0.65000e+02_dp /
4530  data vrs(1:ndat) / &
4531  0.20000e-01_dp,0.31000e-01_dp,0.49000e-01_dp,0.77000e-01_dp,0.12000e+00_dp,0.19000e+00_dp, &
4532  0.30000e+00_dp,0.48000e+00_dp,0.76000e+00_dp,0.12000e+01_dp,0.19000e+01_dp,0.30000e+01_dp, &
4533  0.48000e+01_dp,0.75000e+01_dp,0.11000e+02_dp,0.16000e+02_dp,0.21000e+02_dp,0.26000e+02_dp, &
4534  0.34000e+02_dp,0.41000e+02_dp,0.49000e+02_dp,0.57000e+02_dp,0.65000e+02_dp,0.73000e+02_dp, &
4535  0.81000e+02_dp,0.87000e+02_dp,0.93000e+02_dp,0.99000e+02_dp,0.10750e+03_dp,0.11500e+03_dp, &
4536  0.12500e+03_dp,0.13500e+03_dp,0.14500e+03_dp /
4537  data vrg(1:ndat) / 0.39000e-01_dp,0.62000e-01_dp,0.97000e-01_dp, &
4538  0.15000e+00_dp,0.24000e+00_dp,0.38000e+00_dp,0.61000e+00_dp,0.96000e+00_dp,0.15000e+01_dp, &
4539  0.24000e+01_dp,0.38000e+01_dp,0.61000e+01_dp,0.96000e+01_dp,0.15000e+02_dp,0.23000e+02_dp, &
4540  0.31000e+02_dp,0.39000e+02_dp,0.49000e+02_dp,0.59000e+02_dp,0.68000e+02_dp,0.79000e+02_dp, &
4541  0.88000e+02_dp,0.10000e+03_dp,0.11000e+03_dp,0.13000e+03_dp,0.15000e+03_dp,0.17000e+03_dp, &
4542  0.20000e+03_dp,0.23000e+03_dp,0.26000e+03_dp,0.30000e+03_dp,0.35000e+03_dp,0.40000e+03_dp /
4543  data vrh(1:ndat) / &
4544  0.53000e-01_dp,0.84000e-01_dp,0.13000e+00_dp,0.21000e+00_dp,0.33000e+00_dp,0.52000e+00_dp, &
4545  0.82000e+00_dp,0.13000e+01_dp,0.21000e+01_dp,0.33000e+01_dp,0.52000e+01_dp,0.82000e+01_dp, &
4546  0.13000e+02_dp,0.20000e+02_dp,0.28000e+02_dp,0.36000e+02_dp,0.46000e+02_dp,0.56000e+02_dp, &
4547  0.67000e+02_dp,0.80000e+02_dp,0.97000e+02_dp,0.12000e+03_dp,0.14000e+03_dp,0.17000e+03_dp, &
4548  0.20000e+03_dp,0.24000e+03_dp,0.29000e+03_dp,0.35000e+03_dp,0.42000e+03_dp,0.51000e+03_dp, &
4549  0.61000e+03_dp,0.74000e+03_dp,0.89000e+03_dp /
4550  data blkradl(1:ndat) / &
4551  0.20000e-03_dp,0.25198e-03_dp,0.31748e-03_dp,0.40000e-03_dp,0.50397e-03_dp,0.63496e-03_dp, &
4552  0.80000e-03_dp,0.10079e-02_dp,0.12699e-02_dp,0.16000e-02_dp,0.20159e-02_dp,0.25398e-02_dp, &
4553  0.32000e-02_dp,0.40317e-02_dp,0.50797e-02_dp,0.64000e-02_dp,0.80635e-02_dp,0.10159e-01_dp, &
4554  0.12800e-01_dp,0.16127e-01_dp,0.20319e-01_dp,0.25600e-01_dp,0.32254e-01_dp,0.40637e-01_dp, &
4555  0.51200e-01_dp,0.64508e-01_dp,0.81275e-01_dp,0.10240e+00_dp,0.12902e+00_dp,0.16255e+00_dp, &
4556  0.20480e+00_dp,0.25803e+00_dp,0.32510e+00_dp /
4557  data blkradi(1:ndat,1) / 0.57452e-03_dp,0.72384e-03_dp,0.91199e-03_dp, &
4558  0.11490e-02_dp,0.14477e-02_dp,0.18240e-02_dp,0.22981e-02_dp,0.28954e-02_dp,0.36479e-02_dp, &
4559  0.45961e-02_dp,0.57908e-02_dp,0.72959e-02_dp,0.11572e-01_dp,0.14770e-01_dp,0.18851e-01_dp, &
4560  0.24060e-01_dp,0.30709e-01_dp,0.39194e-01_dp,0.50025e-01_dp,0.63848e-01_dp,0.81491e-01_dp, &
4561  0.10401e+00_dp,0.13275e+00_dp,0.16943e+00_dp,0.21625e+00_dp,0.27601e+00_dp,0.35228e+00_dp, &
4562  0.44962e+00_dp,0.57387e+00_dp,0.73244e+00_dp,0.93484e+00_dp,0.11932e+01_dp,0.15229e+01_dp /
4563  data blkradi(1:ndat,2) / &
4564  0.20715e-03_dp,0.26099e-03_dp,0.32912e-03_dp,0.43555e-03_dp,0.57638e-03_dp,0.76276e-03_dp, &
4565  0.10094e-02_dp,0.13358e-02_dp,0.17677e-02_dp,0.23394e-02_dp,0.30958e-02_dp,0.40969e-02_dp, &
4566  0.54216e-02_dp,0.71748e-02_dp,0.94948e-02_dp,0.12565e-01_dp,0.16628e-01_dp,0.22005e-01_dp, &
4567  0.29120e-01_dp,0.38537e-01_dp,0.50998e-01_dp,0.67488e-01_dp,0.89311e-01_dp,0.11819e+00_dp, &
4568  0.15641e+00_dp,0.20698e+00_dp,0.27391e+00_dp,0.36249e+00_dp,0.47970e+00_dp,0.63481e+00_dp, &
4569  0.84009e+00_dp,0.11117e+01_dp,0.14712e+01_dp /
4570  data blkradi(1:ndat,3) / 0.23559e-03_dp,0.29682e-03_dp,0.37397e-03_dp, &
4571  0.47118e-03_dp,0.59365e-03_dp,0.74795e-03_dp,0.95894e-03_dp,0.12777e-02_dp,0.17025e-02_dp, &
4572  0.22685e-02_dp,0.30227e-02_dp,0.40275e-02_dp,0.53665e-02_dp,0.71506e-02_dp,0.95278e-02_dp, &
4573  0.12695e-01_dp,0.16916e-01_dp,0.22539e-01_dp,0.30032e-01_dp,0.40017e-01_dp,0.70019e-01_dp, &
4574  0.98384e-01_dp,0.13824e+00_dp,0.19424e+00_dp,0.27293e+00_dp,0.38350e+00_dp,0.53885e+00_dp, &
4575  0.75714e+00_dp,0.10639e+01_dp,0.14948e+01_dp,0.21004e+01_dp,0.29513e+01_dp,0.41469e+01_dp /
4576  data blkrads(1:ndat) / &
4577  0.20715e-03_dp,0.26148e-03_dp,0.33067e-03_dp,0.41710e-03_dp,0.52691e-03_dp,0.66640e-03_dp, &
4578  0.84289e-03_dp,0.10674e-02_dp,0.13513e-02_dp,0.17129e-02_dp,0.21670e-02_dp,0.27521e-02_dp, &
4579  0.34989e-02_dp,0.44777e-02_dp,0.57347e-02_dp,0.75389e-02_dp,0.99020e-02_dp,0.13161e-01_dp, &
4580  0.17372e-01_dp,0.23337e-01_dp,0.31058e-01_dp,0.41194e-01_dp,0.55153e-01_dp,0.74854e-01_dp, &
4581  0.99806e-01_dp,0.13463e+00_dp,0.18136e+00_dp,0.24282e+00_dp,0.32955e+00_dp,0.44123e+00_dp, &
4582  0.59884e+00_dp,0.77090e+00_dp,0.99387e+00_dp /
4583  data blkradg(1:ndat) / 0.27144e-03_dp,0.34200e-03_dp,0.43089e-03_dp, &
4584  0.54288e-03_dp,0.68399e-03_dp,0.86177e-03_dp,0.10858e-02_dp,0.13680e-02_dp,0.17235e-02_dp, &
4585  0.21715e-02_dp,0.27360e-02_dp,0.34471e-02_dp,0.43431e-02_dp,0.54719e-02_dp,0.68942e-02_dp, &
4586  0.86861e-02_dp,0.10944e-01_dp,0.13788e-01_dp,0.17372e-01_dp,0.21888e-01_dp,0.27577e-01_dp, &
4587  0.34745e-01_dp,0.43775e-01_dp,0.55154e-01_dp,0.69489e-01_dp,0.87551e-01_dp,0.11031e+00_dp, &
4588  0.13898e+00_dp,0.17510e+00_dp,0.22061e+00_dp,0.27796e+00_dp,0.35020e+00_dp,0.44123e+00_dp /
4589  data blkradh(1:ndat) / &
4590  0.20715e-03_dp,0.26099e-03_dp,0.32883e-03_dp,0.41430e-03_dp,0.52198e-03_dp,0.65766e-03_dp, &
4591  0.82860e-03_dp,0.10440e-02_dp,0.13153e-02_dp,0.16572e-02_dp,0.20879e-02_dp,0.26306e-02_dp, &
4592  0.33144e-02_dp,0.41759e-02_dp,0.52613e-02_dp,0.66288e-02_dp,0.83517e-02_dp,0.10523e-01_dp, &
4593  0.13258e-01_dp,0.16703e-01_dp,0.21045e-01_dp,0.26515e-01_dp,0.33407e-01_dp,0.42090e-01_dp, &
4594  0.53030e-01_dp,0.66814e-01_dp,0.84180e-01_dp,0.10606e+00_dp,0.13363e+00_dp,0.16836e+00_dp, &
4595  0.21212e+00_dp,0.26725e+00_dp,0.33672e+00_dp /
4596  data blkdnsl(1:ndat) / &
4597  0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4598  0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4599  0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4600  0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4601  0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp, &
4602  0.10000e+01_dp,0.10000e+01_dp,0.10000e+01_dp /
4603  data blkdnsi(1:ndat,1) / 0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4604  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4605  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.87368e+00_dp,0.87072e+00_dp,0.86777e+00_dp, &
4606  0.86483e+00_dp,0.86189e+00_dp,0.85897e+00_dp,0.85606e+00_dp,0.85316e+00_dp,0.85026e+00_dp, &
4607  0.84738e+00_dp,0.84451e+00_dp,0.84164e+00_dp,0.83879e+00_dp,0.83595e+00_dp,0.83311e+00_dp, &
4608  0.83029e+00_dp,0.82747e+00_dp,0.82467e+00_dp,0.82187e+00_dp,0.81908e+00_dp,0.81631e+00_dp /
4609  data blkdnsi(1:ndat,2) / &
4610  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4611  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4612  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4613  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4614  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4615  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
4616  data blkdnsi(1:ndat,3) / 0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
4617  0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
4618  0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp, &
4619  0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.61183e+00_dp,0.51790e+00_dp, &
4620  0.45557e+00_dp,0.40075e+00_dp,0.35252e+00_dp,0.31010e+00_dp,0.27278e+00_dp,0.23995e+00_dp, &
4621  0.21108e+00_dp,0.18567e+00_dp,0.16333e+00_dp,0.14367e+00_dp,0.12638e+00_dp,0.11118e+00_dp /
4622  data blkdnss(1:ndat) / &
4623  0.90000e+00_dp,0.89500e+00_dp,0.88500e+00_dp,0.88200e+00_dp,0.87500e+00_dp,0.86500e+00_dp, &
4624  0.85500e+00_dp,0.84200e+00_dp,0.83000e+00_dp,0.81500e+00_dp,0.80500e+00_dp,0.78600e+00_dp, &
4625  0.76500e+00_dp,0.73000e+00_dp,0.69500e+00_dp,0.61183e+00_dp,0.54000e+00_dp,0.46000e+00_dp, &
4626  0.40000e+00_dp,0.33000e+00_dp,0.28000e+00_dp,0.24000e+00_dp,0.20000e+00_dp,0.16000e+00_dp, &
4627  0.13500e+00_dp,0.11000e+00_dp,0.90000e-01_dp,0.75000e-01_dp,0.60000e-01_dp,0.50000e-01_dp, &
4628  0.40000e-01_dp,0.37500e-01_dp,0.35000e-01_dp /
4629  data blkdnsg(1:ndat) / &
4630  0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4631  0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4632  0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4633  0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4634  0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp, &
4635  0.40000e+00_dp,0.40000e+00_dp,0.40000e+00_dp /
4636  data blkdnsh(1:ndat) / &
4637  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4638  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4639  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4640  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4641  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp, &
4642  0.90000e+00_dp,0.90000e+00_dp,0.90000e+00_dp /
4643 
4644  !--- mass
4645  do k = 1, ndat
4646  xmss( il,k ) = log( dble( xl( k ) )*1.e-03_dp )
4647  xmss( ic,k ) = log( dble( xi( k,1 ) )*1.e-03_dp )
4648  xmss( ip,k ) = log( dble( xi( k,2 ) )*1.e-03_dp )
4649  xmss( id,k ) = log( dble( xi( k,3 ) )*1.e-03_dp )
4650  xmss( is,k ) = log( dble( xs( k ) )*1.e-03_dp )
4651  xmss( ig,k ) = log( dble( xg( k ) )*1.e-03_dp )
4652  xmss( ih,k ) = log( dble( xh( k ) )*1.e-03_dp )
4653  end do
4654 
4655  !--- capacity
4656  do k = 1, ndat ! cm -> m
4657  zcap( il,k ) = dble(rlec( k ))*1.e-02_dp
4658  zcap( ic,k ) = dble(riec( k,1 ))*1.e-02_dp
4659  zcap( ip,k ) = dble(riec( k,2 ))*1.e-02_dp
4660  zcap( id,k ) = dble(riec( k,3 ))*1.e-02_dp
4661  zcap( is,k ) = dble(rsec( k ))*1.e-02_dp
4662  zcap( ig,k ) = dble(rgec( k ))*1.e-02_dp
4663  zcap( ih,k ) = dble(rhec( k ))*1.e-02_dp
4664  end do
4665 
4666  !--- terminal velocity
4667  do k = 1, ndat ! cm/s -> m/s
4668  vtrm( il,k ) = dble(vrl( k ))*1.e-02_dp
4669  vtrm( ic,k ) = dble(vri( k,1 ))*1.e-02_dp
4670  vtrm( ip,k ) = dble(vri( k,2 ))*1.e-02_dp
4671  vtrm( id,k ) = dble(vri( k,3 ))*1.e-02_dp
4672  vtrm( is,k ) = dble(vrs( k ))*1.e-02_dp
4673  vtrm( ig,k ) = dble(vrg( k ))*1.e-02_dp
4674  vtrm( ih,k ) = dble(vrh( k ))*1.e-02_dp
4675  end do
4676 
4677  !--- bulk radii
4678  do k = 1, ndat ! cm -> mm
4679  blkr( il,k ) = dble(blkradl( k ))*10.0_dp
4680  blkr( ic,k ) = dble(blkradi( k,1 ))*10.0_dp
4681  blkr( ip,k ) = dble(blkradi( k,2 ))*10.0_dp
4682  blkr( id,k ) = dble(blkradi( k,3 ))*10.0_dp
4683  blkr( is,k ) = dble(blkrads( k ))*10.0_dp
4684  blkr( ig,k ) = dble(blkradg( k ))*10.0_dp
4685  blkr( ih,k ) = dble(blkradh( k ))*10.0_dp
4686  end do
4687 
4688  !--- bulk density
4689  do k = 1, ndat ! g/cm^3 -> kg/m^3
4690  blkd( il,k ) = dble(blkdnsl( k ))*1000._dp
4691  blkd( ic,k ) = dble(blkdnsi( k,1 ))*1000._dp
4692  blkd( ip,k ) = dble(blkdnsi( k,2 ))*1000._dp
4693  blkd( id,k ) = dble(blkdnsi( k,3 ))*1000._dp
4694  blkd( is,k ) = dble(blkdnss( k ))*1000._dp
4695  blkd( ig,k ) = dble(blkdnsg( k ))*1000._dp
4696  blkd( ih,k ) = dble(blkdnsh( k ))*1000._dp
4697  end do
4698 
4699  !--- collection kernel
4700  ! cm**3/s -> m**3/s
4701  do i = 1, ndat
4702  do j = 1, ndat
4703  do k = 1, 7
4704  do kk = 1, 7
4705  ykrn( k,kk,i,j ) = kernels( k,kk,i,j )*1.e-06_dp
4706  enddo
4707  enddo
4708 
4709  end do
4710  end do
4711 
4712  end subroutine rdkdat
4713  !---------------------------------------------------------------------------------------
4714  subroutine sdfgrid
4715 
4716  real(DP) :: xsta, xend
4717  integer :: n
4718  real(DP):: pi_dp = 3.1415920_dp
4719  real(DP):: rhow_dp = 1.0e+03_dp
4720 
4721  xsta = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-06_dp )**3.0_dp )
4722  xend = log( rhow_dp * 4.0_dp*pi_dp/3.0_dp * ( 3.e-03_dp )**3.0_dp )
4723 
4724  dxmic_mk = ( xend-xsta )/nbin
4725  dxmic = dxmic_mk
4726  do n = 1, nbin+1
4727  xbnd_mk( n ) = xsta + dxmic_mk*( n-1 )
4728  end do
4729  do n = 1, nbin
4730  xctr_mk( n ) = ( xbnd_mk( n )+xbnd_mk( n+1 ) )*0.50_dp
4731  radc_mk( n ) = ( exp( xctr_mk( n ) )*3.0_dp/4.0_dp/pi_dp/rhow_dp )**( 1.0_dp/3.0_dp )
4732  end do
4733 
4734  end subroutine sdfgrid
4735  !---------------------------------------------------------------------------------------
4736  subroutine getcp
4737 
4738  integer :: myu, n
4739 
4740  do myu = 1, nspc_mk
4741  do n = 1, nbin
4742  cctr_mk( myu,n ) = fcpc( myu,xctr_mk( n ) )
4743  end do
4744  do n = 1, nbin+1
4745  cbnd_mk( myu,n ) = fcpc( myu,xbnd_mk( n ) )
4746  end do
4747  end do
4748 
4749  end subroutine getcp
4750  !---------------------------------------------------------------------------------------
4751  function fcpc( myu,x )
4752 
4753  integer, intent(in) :: myu
4754  real(DP), intent(in) :: x
4755  real(DP) :: fcpc
4756 
4757  real(DP) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
4758 
4759  call getknot &
4760  ( ndat, kdeg, xmss( myu,: ), & !--- in
4761  qknt ) !--- out
4762 
4763  call getmatrx &
4764  ( ndat, kdeg, qknt, xmss( myu,: ), & !--- in
4765  elm ) !--- out
4766 
4767  call getcoef &
4768  ( ndat, kdeg, elm, zcap( myu,: ), & !--- in
4769  coef ) !--- out
4770 
4771  fcpc = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
4772 
4773  end function fcpc
4774  !---------------------------------------------------------------------------------------
4775  subroutine getck
4776 
4777  integer :: myu, nyu, i, j
4778 
4779  if( io_l ) write(io_fid_log,*) 'Create micpara.dat'
4780  if( kphase == 0 ) then
4781  if( io_l ) write(io_fid_log,*) 'Hydro-dynamic kernel'
4782  else if( kphase == 1 ) then
4783  if( io_l ) write(io_fid_log,*) 'Long Kernel'
4784  else if( kphase == 2 ) then
4785  if( io_l ) write(io_fid_log,*) 'Golovin Kernel'
4786  endif
4787 
4788  do myu = 1, nspc_mk
4789  do nyu = 1, nspc_mk
4790  if( io_l ) write(io_fid_log,*) ' myu, nyu :', myu, nyu
4791  do i = 1, nbin
4792  do j = 1, nbin
4793  ck_mk( myu,nyu,i,j ) = fckrn( myu,nyu,xctr_mk( i ),xctr_mk( j ) )
4794  end do
4795  end do
4796  end do
4797  end do
4798 
4799  return
4800 
4801  end subroutine getck
4802  !---------------------------------------------------------------------------------------
4803  function fckrn( myu,nyu,x,y )
4804 
4805  integer, intent(in) :: myu, nyu
4806  real(DP), intent(in) :: x, y
4807  real(DP) :: fckrn
4808 
4809  real(DP) :: qknt( ndat+kdeg ), rknt( ndat+kdeg )
4810  real(DP) :: coef( ndat,ndat )
4811 
4812  real(DP) :: xlrg, xsml, vlrg, vsml, rlrg
4813  real(DP):: pi_dp = 3.1415920_dp
4814  real(DP):: rhow_dp = 1.0e+03_dp
4815 
4816  if( kphase == 0 ) then
4817  call getknot &
4818  ( ndat, kdeg, xmss( myu,: ), & !--- in
4819  qknt ) !--- out
4820 
4821  rknt( : ) = qknt( : )
4822 
4823  call getcoef2 &
4824  ( ndat, ndat, kdeg, kdeg, & !--- in
4825  xmss( myu,: ), xmss( nyu,: ), & !--- in
4826  qknt, rknt, & !--- in
4827  ykrn( myu,nyu,:,: ), & !--- in
4828  coef ) !--- out
4829 
4830  fckrn = fspline2 &
4831  ( ndat, ndat, kdeg, kdeg, & !--- in
4832  coef, qknt, rknt, & !--- in
4833  xmss( myu,: ), xmss( nyu,: ), & !--- in
4834  x, y ) !--- in
4835  else if( kphase == 1 ) then
4836  xlrg = max( x, y )
4837  xsml = min( x, y )
4838 
4839  vlrg = (exp( xlrg ) / rhow_dp )*1.e+06_dp
4840  vsml = (exp( xsml ) / rhow_dp )*1.e+06_dp
4841 
4842  rlrg = ( exp( xlrg )/( 4.0_dp*pi_dp*rhow_dp )*3.0_dp )**(1.0_dp/3.0_dp )*1.e+06_dp
4843 
4844  if( rlrg <=50.0_dp ) then
4845  fckrn = 9.44e+03_dp*( vlrg*vlrg + vsml*vsml )
4846  else
4847  fckrn = 5.78e-03_dp*( vlrg+vsml )
4848  end if
4849  else if( kphase == 2 ) then
4850  fckrn = 1.5_dp*( exp(x) +exp(y) )
4851  end if
4852 
4853  return
4854 
4855  end function fckrn
4856  !---------------------------------------------------------------------------------------
4857  subroutine getvt
4858 
4859  integer :: myu, n
4860 
4861  do myu = 1, nspc_mk
4862  do n = 1, nbin
4863  vt_mk( myu,n ) = max( fvterm( myu,xctr_mk( n ) ), 0.0_dp )
4864  end do
4865  end do
4866 
4867  end subroutine getvt
4868  !---------------------------------------------------------------------------------------
4869  function fvterm( myu,x )
4870 
4871  integer, intent(in) :: myu
4872  real(DP), intent(in) :: x
4873  real(DP) :: fvterm
4874 
4875  real(DP) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
4876 
4877  call getknot &
4878  ( ndat, kdeg, xmss( myu,: ), & !--- in
4879  qknt ) !--- out
4880 
4881  call getmatrx &
4882  ( ndat, kdeg, qknt, xmss( myu,: ), & !--- in
4883  elm ) !--- out
4884 
4885  call getcoef &
4886  ( ndat, kdeg, elm, vtrm( myu,: ), & !--- in
4887  coef ) !--- out
4888 
4889  fvterm = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
4890 
4891  end function fvterm
4892  !---------------------------------------------------------------------------------------
4893  subroutine getbr
4894 
4895  integer :: myu, n
4896 
4897  do myu = 1, nspc_mk
4898  do n = 1, nbin
4899  br_mk( myu,n ) = fbulkrad( myu, xctr_mk( n ) )
4900  end do
4901  end do
4902 
4903  end subroutine getbr
4904  !---------------------------------------------------------------------------------------
4905  function fbulkrad( myu,x )
4906 
4907  integer, intent(in) :: myu
4908  real(DP), intent(in) :: x
4909  real(DP) :: fbulkrad
4910 
4911  real(DP) :: qknt( ndat+kdeg ), elm( ndat,ndat ), coef( ndat )
4912 
4913  call getknot &
4914  ( ndat, kdeg, xmss( myu,: ), & !--- in
4915  qknt ) !--- out
4916 
4917  call getmatrx &
4918  ( ndat, kdeg, qknt, xmss( myu,: ), & !--- in
4919  elm ) !--- out
4920 
4921  call getcoef &
4922  ( ndat, kdeg, elm, blkr( myu,: ), & !--- in
4923  coef ) !--- out
4924 
4925  fbulkrad = fspline( ndat, kdeg, coef, qknt, xmss( myu,: ), x )
4926 
4927  end function fbulkrad
4928  !---------------------------------------------------------------------------------------
4929  subroutine paraout
4930 
4931  integer :: myu, nyu, i, j, n
4932 
4933  open ( fid_micpara, file = fname_micpara, form = 'formatted', status='new' )
4934 
4935  write( fid_micpara,* ) nspc_mk, nbin
4936 
4937  ! grid parameter
4938  do n = 1, nbin
4939  xctr( n ) = xctr_mk( n )
4940  radc( n ) = radc_mk( n )
4941  write( fid_micpara,* ) n, xctr( n ), radc( n )
4942  end do
4943  do n = 1, nbin+1
4944  xbnd( n ) = xbnd_mk( n )
4945  write( fid_micpara,* ) n, xbnd( n )
4946  end do
4947  write( fid_micpara,* ) dxmic_mk
4948 
4949  ! capacity
4950  do myu = 1, nspc_mk
4951  do n = 1, nbin
4952  cctr( n,myu ) = cctr_mk( myu,n )
4953  write( fid_micpara,* ) myu, n, cctr( n,myu )
4954  end do
4955  do n = 1, nbin+1
4956  cbnd( n,myu ) = cbnd_mk( myu,n )
4957  write( fid_micpara,* ) myu, n, cbnd( n,myu )
4958  end do
4959  end do
4960 
4961  ! collection kernel
4962  do myu = 1, nspc_mk
4963  do nyu = 1, nspc_mk
4964  do i = 1, nbin
4965  do j = 1, nbin
4966  ck( myu,nyu,i,j ) = ck_mk( myu,nyu,i,j )
4967  write( fid_micpara,* ) myu, nyu, i, j, ck( myu,nyu,i,j )
4968  end do
4969  end do
4970  end do
4971  end do
4972 
4973  ! falling velocity
4974  do myu = 1, nspc_mk
4975  do n = 1, nbin
4976  vt( myu,n ) = vt_mk( myu,n )
4977  write( fid_micpara,* ) myu, n, vt( myu,n )
4978  end do
4979  end do
4980 
4981  ! bulk radius
4982  do myu = 1, nspc_mk
4983  do n = 1, nbin
4984  br( myu,n ) = br_mk( myu,n )
4985  write( fid_micpara,* ) myu, n, br( myu,n )
4986  end do
4987  end do
4988 
4989  close ( fid_micpara )
4990 
4991  end subroutine paraout
4992  !---------------------------------------------------------------------------------------
4993  !---- unify from other files
4994  !---------------------------------------------------------------------------------------
4995  subroutine tinvss(n,a,dt,e,nn,iw,inder)
4996 
4997  implicit none
4998 
4999  integer, intent(in) :: n, nn
5000  integer, intent(inout) :: inder
5001  real(DP), intent(inout) :: a(nn,n)
5002  real(DP), intent(inout) :: dt
5003  real(DP), intent(in) :: e
5004  integer, intent(inout) :: iw( 2*n )
5005  integer :: i, j, k, kk, ij, nnn
5006  real(DP) :: work, aa, az, eps
5007  integer :: ipiv, jpiv
5008  real(DP) :: piv
5009  !-----------------------------------------
5010 
5011  inder = 0
5012  if( n < 1 ) then
5013  goto 910
5014  elseif( n == 1 ) then
5015  goto 930
5016  elseif( n > 1 ) then
5017  goto 101
5018  endif
5019 
5020 101 continue
5021  if( n > nn ) then
5022  inder = -1
5023  write(6,690) "n= ", n, "nn= ", nn, &
5024  "n should be less than or equal to nn in TINVSS"
5025 690 format(a8,i5,5x,a4,i5,a55)
5026  return
5027  endif
5028 
5029  eps = 0.0_dp
5030  dt = 1.0_dp
5031  do k = 1, n
5032  piv = 0.0_dp
5033  do i = k, n
5034  do j = k, n
5035  if( abs(a(i,j)) <= abs(piv) ) goto 110 !
5036  ipiv = i
5037  jpiv = j
5038  piv = a(i,j)
5039 110 continue
5040  enddo
5041  enddo
5042  dt = dt * piv
5043  if( abs(piv) <= eps ) goto 920
5044  if( k == 1 ) eps = abs(piv)*e
5045  if( ipiv == k ) goto 130
5046  dt = -dt
5047  do j = 1, n
5048  work = a(ipiv,j)
5049  a(ipiv,j) = a(k,j)
5050  a(k,j) = work
5051  enddo
5052 130 continue
5053  if( jpiv == k ) goto 150
5054  dt = -dt
5055  do i = 1, n
5056  work = a(i,jpiv)
5057  a(i,jpiv) = a(i,k)
5058  a(i,k) = work
5059  enddo
5060 150 continue
5061  iw(2*k-1) = ipiv
5062  aa=1.0_dp/piv
5063  iw(2*k) = jpiv
5064  do j = 1, n
5065  a(k,j) = a(k,j)*aa
5066  enddo
5067  do i = 1, n
5068  if( i == k ) goto 220
5069  az = a(i,k)
5070  if( az == 0.0_dp ) goto 220
5071  do j = 1, n
5072  a(i,j) = a(i,j)-a(k,j)*az
5073  enddo
5074  a(i,k) = -aa*az
5075 220 continue
5076  enddo
5077  a(k,k) = aa
5078  enddo
5079  do kk = 2, n
5080  k=n+1-kk
5081  ij=iw(2*k)
5082  if( ij == k ) goto 420
5083  do j = 1, n
5084  work=a(ij,j)
5085  a(ij,j) = a(k,j)
5086  a(k,j) = work
5087  enddo
5088 420 continue
5089  ij = iw(2*k-1)
5090  if( ij == k ) goto 400
5091  do i = 1, n
5092  work=a(i,ij)
5093  a(i,ij)=a(i,k)
5094  a(i,k)=work
5095  enddo
5096 400 continue
5097  enddo
5098 
5099  return
5100 
5101 910 continue
5102  inder = -1
5103  write(*,691) "n= ", n, "should be positive in TINVSS"
5104 691 format(a8,i5,5x,a30)
5105  return
5106 
5107 
5108 920 continue
5109  dt = 0.0_dp
5110  inder = n-k+1
5111  nnn = k-1
5112  write(*,692) 'given matrix A to TINVSS is ill conditioned, or sigular withrank =', nnn, &
5113  'return with no further calculation'
5114 692 format(a,1x,i4,1x,a)
5115  return
5116 
5117 930 continue
5118  dt=a(1,1)
5119  k=1
5120  if( dt == 0.0_dp ) goto 920
5121  a(1,1) = 1.0_dp/a(1,1)
5122  return
5123 
5124  end subroutine tinvss
5125  !---------------------------------------------------------------
5126  subroutine getknot &
5127  ( ndat, kdeg, xdat, & !--- in
5128  qknt ) !--- out
5129 
5130  integer, intent(in) :: ndat ! number of data
5131  integer, intent(in) :: kdeg ! degree of Spline + 1
5132  real(DP), intent(in) :: xdat( ndat ) ! data of independent var.
5133 
5134  real(DP), intent(out) :: qknt( ndat+kdeg ) ! knots for B-Spline
5135 
5136  !--- local
5137  integer :: i
5138 
5139  do i = 1, kdeg
5140  qknt( i ) = xdat( 1 )
5141  end do
5142 
5143  do i = 1, ndat-kdeg
5144  qknt( i+kdeg ) = ( xdat( i )+xdat( i+kdeg ) )*0.50_dp
5145  end do
5146 
5147  do i = 1, kdeg
5148  qknt( ndat+i ) = xdat( ndat )
5149  end do
5150 
5151  return
5152 
5153  end subroutine getknot
5154  !---------------------------------------------------------------
5155  recursive function fbspl ( ndat, inum, kdeg, qknt, xdat, x ) &
5156  result(bspl)
5157 
5158  real(DP) :: bspl
5159 
5160  integer, intent(in) :: ndat ! number of data
5161  integer, intent(in) :: inum ! index of B-Spline
5162  integer, intent(in) :: kdeg ! degree of B-Spline + 1
5163  real(DP), intent(in) :: qknt( ndat+kdeg ) ! knot of B-Spline
5164  real(DP), intent(in) :: xdat( ndat ) ! data of independent variable
5165  real(DP), intent(in) :: x ! interpolation point
5166 
5167  !--- local
5168  real(DP) :: bsp1, bsp2
5169 
5170  if ( ( inum == 1 .AND. x == xdat( 1 ) ) .OR. &
5171  ( inum == ndat .AND. x == xdat( ndat ) ) ) then
5172  bspl = 1.
5173  return
5174  end if
5175 
5176  if ( kdeg == 1 ) then
5177  if ( x >= qknt( inum ) .AND. x < qknt( inum+1 ) ) then
5178  bspl = 1.0_dp
5179  else
5180  bspl = 0.0_dp
5181  end if
5182  else
5183  if ( qknt( inum+kdeg-1 ) /= qknt( inum ) ) then
5184  bsp1 = ( x-qknt( inum ) ) &
5185  /( qknt( inum+kdeg-1 )-qknt( inum ) ) &
5186  * fbspl( ndat, inum, kdeg-1, qknt, xdat, x )
5187  else
5188  bsp1 = 0.0_dp
5189  end if
5190  if ( qknt( inum+kdeg ) /= qknt( inum+1 ) ) then
5191  bsp2 = ( qknt( inum+kdeg )-x ) &
5192  /( qknt( inum+kdeg )-qknt( inum+1 ) ) &
5193  * fbspl( ndat, inum+1, kdeg-1, qknt, xdat, x )
5194  else
5195  bsp2 = 0.0_dp
5196  end if
5197  bspl = bsp1 + bsp2
5198  end if
5199 
5200  end function fbspl
5201  !---------------------------------------------------------------
5202  function fpb( ndat, i, kdeg, qknt, xdat, elm, x )
5203 
5204  real :: fpb
5205  integer :: ndat, i, kdeg
5206  real(DP) :: qknt( ndat+kdeg ), xdat( ndat ), elm( ndat,ndat )
5207  real(DP) :: x
5208 
5209  integer :: l
5210  real(DP) :: sum
5211 
5212  sum = 0.0_dp
5213  do l = 1, ndat
5214  sum = sum + elm( l,i )*fbspl( ndat, l, kdeg, qknt, xdat, x )
5215  end do
5216 
5217  fpb = sum
5218 
5219  end function fpb
5220  !---------------------------------------------------------------
5221  subroutine getmatrx &
5222  ( ndat, kdeg, qknt, xdat, & !--- in
5223  elm ) !--- out
5224 
5225 ! use scale_tinvss, only: TINVSS
5226 
5227  integer, intent(in) :: ndat
5228  integer, intent(in) :: kdeg
5229  real(DP), intent(in) :: qknt( ndat+kdeg )
5230  real(DP), intent(in) :: xdat( ndat )
5231 
5232  real(DP), intent(out) :: elm( ndat,ndat )
5233 
5234  !--- local
5235  real(DP) :: dt
5236  integer :: iw( 2*ndat ), i, j, inder
5237  real(DP), parameter :: eps = 0.
5238 
5239  do i = 1, ndat
5240  do j = 1, ndat
5241  elm( i,j ) = fbspl( ndat, j, kdeg, qknt, xdat, xdat( i ) )
5242  end do
5243  end do
5244 
5245  call tinvss( ndat, elm, dt, eps, ndat, iw, inder )
5246 
5247  return
5248 
5249  end subroutine getmatrx
5250  !---------------------------------------------------------------
5251  subroutine getcoef &
5252  ( ndat, kdeg, elm, ydat, & !--- in
5253  coef ) !--- out
5254 
5255  integer, intent(in) :: ndat ! number of data
5256  integer, intent(in) :: kdeg ! degree of Spline + 1
5257  real(DP), intent(in) :: elm( ndat,ndat ) ! matrix ( inverse )
5258  real(DP), intent(in) :: ydat( ndat ) ! data of dependent var.
5259 
5260  real(DP), intent(out) :: coef( ndat ) ! expansion coefficient
5261 
5262  !--- local
5263  integer :: i, j
5264  real(DP) :: sum
5265 
5266  do i = 1, ndat
5267  sum = 0.0_dp
5268  do j = 1, ndat
5269  sum = sum + elm( i,j )*ydat( j )
5270  end do
5271  coef( i ) = sum
5272  end do
5273 
5274  return
5275 
5276  end subroutine getcoef
5277  !---------------------------------------------------------------
5278  function fspline ( ndat, kdeg, coef, qknt, xdat, x )
5279 
5280  integer, intent(in) :: ndat
5281  integer, intent(in) :: kdeg
5282  real(DP), intent(in) :: coef( ndat )
5283  real(DP), intent(in) :: qknt( ndat+kdeg )
5284  real(DP), intent(in) :: xdat( ndat )
5285  real(DP), intent(in) :: x
5286 
5287  real(DP) :: fspline
5288 
5289  !--- local
5290  real(DP) :: sum
5291  integer :: i
5292 
5293  sum = 0.0_dp
5294  do i = 1, ndat
5295  sum = sum + coef( i )*fbspl( ndat, i, kdeg, qknt, xdat, x )
5296  end do
5297 
5298  fspline = sum
5299 
5300  return
5301 
5302  end function fspline
5303  !---------------------------------------------------------------
5304  subroutine getcoef2 &
5305  ( mdat, ndat, kdeg, ldeg, & !--- in
5306  xdat, ydat, qknt, rknt, zdat, & !--- in
5307  coef ) !--- out
5308 
5309 ! use scale_tinvss, only: TINVSS
5310 
5311  integer, intent(in) :: mdat ! number of data (x-direction)
5312  integer, intent(in) :: ndat ! number of data (y-direction)
5313  integer, intent(in) :: kdeg ! degree of Spline + 1 (x)
5314  integer, intent(in) :: ldeg ! degree of Spline + 1 (y)
5315  real(DP), intent(in) :: xdat( mdat ) ! data of independent var. (x)
5316  real(DP), intent(in) :: ydat( ndat ) ! data of independent var. (y)
5317  real(DP), intent(in) :: qknt( mdat+kdeg ) ! knots of B-Spline (x)
5318  real(DP), intent(in) :: rknt( ndat+ldeg ) ! knots of B-Spline (y)
5319  real(DP), intent(in) :: zdat( mdat,ndat ) ! data of dependent var.
5320 
5321  real(DP), intent(out) :: coef( mdat,ndat ) ! expansion coefficient
5322 
5323  !--- local
5324  real(DP) :: elmx( mdat,mdat ), elmy( ndat,ndat )
5325  integer :: iw1( 2*mdat ), iw2( 2*ndat )
5326  real(DP) :: beta( mdat,ndat ), sum, dt
5327  real(DP), parameter :: eps = 0.0_dp
5328  integer :: i, j, k, l, inder
5329 
5330  do i = 1, mdat
5331  do j = 1, mdat
5332  elmx( i,j ) = fbspl( mdat, j, kdeg, qknt, xdat, xdat( i ) )
5333  end do
5334  end do
5335  call tinvss( mdat, elmx, dt, eps, mdat, iw1, inder )
5336 
5337  do l = 1, ndat
5338  do i = 1, mdat
5339  sum = 0.0_dp
5340  do j = 1, mdat
5341  sum = sum + elmx( i,j )*zdat( j,l )
5342  end do
5343  beta( i,l ) = sum
5344  end do
5345  end do
5346 
5347  do i = 1, ndat
5348  do j = 1, ndat
5349  elmy( i,j ) = fbspl( ndat, j, ldeg, rknt, ydat, ydat( i ) )
5350  end do
5351  end do
5352  call tinvss( ndat, elmy, dt, eps, ndat, iw2, inder )
5353 
5354  do k = 1, mdat
5355  do i = 1, ndat
5356  sum = 0.0_dp
5357  do j = 1, ndat
5358  sum = sum + elmy( i,j )*beta( k,j )
5359  end do
5360  coef( k,i ) = sum
5361  end do
5362  end do
5363 
5364  return
5365 
5366  end subroutine getcoef2
5367  !---------------------------------------------------------------
5368  function fspline2 &
5369  ( mdat, ndat, kdeg, ldeg, & !--- in
5370  coef, qknt, rknt, xdat, ydat, & !--- in
5371  x, y ) !--- in
5372 
5373  integer, intent(in) :: mdat
5374  integer, intent(in) :: ndat
5375  integer, intent(in) :: kdeg
5376  integer, intent(in) :: ldeg
5377  real(DP), intent(in) :: coef( mdat,ndat )
5378  real(DP), intent(in) :: qknt( mdat+kdeg )
5379  real(DP), intent(in) :: rknt( ndat+ldeg )
5380  real(DP), intent(in) :: xdat( mdat ), ydat( ndat )
5381  real(DP), intent(in) :: x, y
5382 
5383  real(DP) :: fspline2
5384 
5385  !--- local
5386  real(DP) :: sum, add
5387  integer :: i, j
5388 
5389  sum = 0.0_dp
5390  do i = 1, mdat
5391  do j = 1, ndat
5392  add = coef( i,j )*fbspl( mdat, i, kdeg, qknt, xdat, x ) &
5393  *fbspl( ndat, j, ldeg, rknt, ydat, y )
5394  sum = sum + add
5395  end do
5396  end do
5397 
5398  fspline2 = sum
5399 
5400  return
5401 
5402  end function fspline2
5403 
5404 end module scale_atmos_phy_mp_suzuki10
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_comm::comm_datatype
integer, public comm_datatype
datatype of variable
Definition: scale_comm.F90:117
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:59
scale_const::const_psat0
real(rp), parameter, public const_psat0
saturate pressure of water vapor at 0C [Pa]
Definition: scale_const.F90:83
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_process::prc_ismaster
logical, public prc_ismaster
master process in local communicator?
Definition: scale_process.F90:80
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_config
subroutine, public atmos_phy_mp_suzuki10_config(MP_TYPE, QA, QS)
Config.
Definition: scale_atmos_phy_mp_suzuki10.F90:281
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_atmos_phy_mp_suzuki10::mp_suzuki10
subroutine mp_suzuki10(ijkmax, ijkmax_cold, ijkmax_warm, index_cold, index_warm, dens, pres, ccn, temp, qvap, ghyd, gaer, evaporate, dt)
Definition: scale_atmos_phy_mp_suzuki10.F90:1452
scale_const::const_ci
real(rp), parameter, public const_ci
specific heat (ice) [J/kg/K]
Definition: scale_const.F90:69
scale_atmos_phy_mp_suzuki10::nbin
integer, public nbin
Definition: scale_atmos_phy_mp_suzuki10.F90:109
scale_const::const_dwatr
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:84
scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer
subroutine, public atmos_phy_mp_negative_fixer(DENS, RHOT, QTRC, I_QV, limit_negative)
Negative fixer.
Definition: scale_atmos_phy_mp_common.F90:64
scale_tracer::tracer_r
real(rp), dimension(qa_max), public tracer_r
Definition: scale_tracer.F90:41
scale_time::time_dtsec_atmos_phy_mp
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
scale_atmos_hydrometeor::i_hs
integer, parameter, public i_hs
snow
Definition: scale_atmos_sub_hydrometeor.F90:75
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
scale_grid::grid_cz
real(rp), dimension(:), allocatable, public grid_cz
center coordinate [m]: z, local=global
Definition: scale_grid_cartesian.F90:56
scale_atmos_saturation::cpovr_liq
real(rp), public cpovr_liq
Definition: scale_atmos_sub_saturation.F90:128
scale_const::const_cl
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
Definition: scale_const.F90:68
scale_atmos_hydrometeor::i_hr
integer, parameter, public i_hr
liquid water rain
Definition: scale_atmos_sub_hydrometeor.F90:73
scale_atmos_hydrometeor::i_hi
integer, parameter, public i_hi
ice water cloud
Definition: scale_atmos_sub_hydrometeor.F90:74
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_unit
character(len=h_short), dimension(:), allocatable, target, public atmos_phy_mp_suzuki10_unit
Definition: scale_atmos_phy_mp_suzuki10.F90:105
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:92
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:35
scale_atmos_phy_mp_suzuki10::iceflg
integer, public iceflg
Definition: scale_atmos_phy_mp_suzuki10.F90:114
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_setup
subroutine, public atmos_phy_mp_suzuki10_setup
Setup.
Definition: scale_atmos_phy_mp_suzuki10.F90:408
scale_atmos_hydrometeor::i_hh
integer, parameter, public i_hh
hail
Definition: scale_atmos_sub_hydrometeor.F90:77
scale_const::const_tmelt
real(rp), parameter, public const_tmelt
Definition: scale_const.F90:75
scale_const::const_dice
real(rp), parameter, public const_dice
density of ice [kg/m3]
Definition: scale_const.F90:85
scale_tracer::tracer_cv
real(rp), dimension(qa_max), public tracer_cv
Definition: scale_tracer.F90:40
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
scale_const::const_lhs0
real(rp), parameter, public const_lhs0
latent heat of sublimation at 0C [J/kg]
Definition: scale_const.F90:79
scale_atmos_phy_mp_suzuki10::nspc
integer, public nspc
Definition: scale_atmos_phy_mp_suzuki10.F90:110
scale_atmos_phy_mp_common
module ATMOSPHERE / Physics Cloud Microphysics - Common
Definition: scale_atmos_phy_mp_common.F90:16
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_mixingratio
subroutine, public atmos_phy_mp_suzuki10_mixingratio(Qe, QTRC0)
Calculate mixing ratio of each category.
Definition: scale_atmos_phy_mp_suzuki10.F90:4238
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_stdio::io_nml
logical, public io_nml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:62
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_stdio::io_get_available_fid
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_stdio.F90:313
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_effectiveradius
subroutine, public atmos_phy_mp_suzuki10_effectiveradius(Re, QTRC0, DENS0, TEMP0)
Calculate Effective Radius.
Definition: scale_atmos_phy_mp_suzuki10.F90:4107
scale_grid_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_grid_index.F90:54
scale_const::const_lhv0
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
Definition: scale_const.F90:77
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:90
scale_comm::comm_world
integer, public comm_world
communication world ID
Definition: scale_comm.F90:118
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_name
character(len=h_short), dimension(:), allocatable, target, public atmos_phy_mp_suzuki10_name
Definition: scale_atmos_phy_mp_suzuki10.F90:103
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_cloudfraction
subroutine, public atmos_phy_mp_suzuki10_cloudfraction(cldfrac, QTRC, mask_criterion)
Calculate Cloud Fraction.
Definition: scale_atmos_phy_mp_suzuki10.F90:4054
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_time
module TIME
Definition: scale_time.F90:15
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_dens
real(rp), dimension(n_hyd), target, public atmos_phy_mp_suzuki10_dens
Definition: scale_atmos_phy_mp_suzuki10.F90:107
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:71
scale_atmos_saturation::lovr_ice
real(rp), public lovr_ice
Definition: scale_atmos_sub_saturation.F90:133
scale_atmos_phy_mp_suzuki10::mkpara
subroutine mkpara
Definition: scale_atmos_phy_mp_suzuki10.F90:4337
scale_atmos_phy_mp_suzuki10::nccn
integer, public nccn
Definition: scale_atmos_phy_mp_suzuki10.F90:111
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:65
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_atmos_hydrometeor::i_hc
integer, parameter, public i_hc
liquid water cloud
Definition: scale_atmos_sub_hydrometeor.F90:72
scale_process::prc_masterrank
integer, parameter, public prc_masterrank
master process in each communicator
Definition: scale_process.F90:57
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10
subroutine, public atmos_phy_mp_suzuki10(DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
Cloud Microphysics.
Definition: scale_atmos_phy_mp_suzuki10.F90:854
scale_const::const_emelt
real(rp), parameter, public const_emelt
Definition: scale_const.F90:74
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:156
scale_random::random_get
subroutine, public random_get(var)
Get random number.
Definition: scale_random.F90:136
scale_grid_index::kijmax
integer, public kijmax
of computational cells: z*x*y
Definition: scale_grid_index.F90:65
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_atmos_saturation::lovr_liq
real(rp), public lovr_liq
Definition: scale_atmos_sub_saturation.F90:132
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_grid::grid_cdz
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
Definition: scale_grid_cartesian.F90:59
scale_atmos_hydrometeor::atmos_hydrometeor_regist
subroutine, public atmos_hydrometeor_regist(Q0, NV, NL, NI, NAME, DESC, UNIT, ADVC)
Regist tracer.
Definition: scale_atmos_sub_hydrometeor.F90:196
scale_atmos_phy_mp_suzuki10::r_collcoag
subroutine r_collcoag(ijkmax, swgt, temp, gc, dtime)
Definition: scale_atmos_phy_mp_suzuki10.F90:3855
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_atmos_phy_mp_common::atmos_phy_mp_precipitation
subroutine, public atmos_phy_mp_precipitation(QA_MP, QS_MP, qflx, vterm, DENS, MOMZ, MOMX, MOMY, RHOE, QTRC, temp, CVq, dt, vt_fixed)
precipitation transport
Definition: scale_atmos_phy_mp_common.F90:763
scale_atmos_phy_mp_suzuki10::nccn1
integer, public nccn1
Definition: scale_atmos_phy_mp_suzuki10.F90:112
scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_desc
character(len=h_mid), dimension(:), allocatable, target, public atmos_phy_mp_suzuki10_desc
Definition: scale_atmos_phy_mp_suzuki10.F90:104
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_tracer::tracer_regist
subroutine, public tracer_regist(QS, NQ, NAME, DESC, UNIT, CV, CP, R, ADVC, MASS)
Regist tracer.
Definition: scale_tracer.F90:64
scale_random
module RANDOM
Definition: scale_random.F90:13
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_atmos_hydrometeor::n_hyd
integer, parameter, public n_hyd
Definition: scale_atmos_sub_hydrometeor.F90:70
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:66
scale_const::const_thermodyn_type
character(len=h_short), public const_thermodyn_type
internal energy type
Definition: scale_const.F90:98
scale_atmos_phy_mp_suzuki10
module Spectran Bin Microphysics
Definition: scale_atmos_phy_mp_suzuki10.F90:31
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:204
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
scale_atmos_hydrometeor::i_hg
integer, parameter, public i_hg
graupel
Definition: scale_atmos_sub_hydrometeor.F90:76
scale_stdio::io_fid_nml
integer, public io_fid_nml
Log file ID (only for output namelist)
Definition: scale_stdio.F90:57
scale_tracer::tracer_mass
real(rp), dimension(qa_max), public tracer_mass
Definition: scale_tracer.F90:42
scale_atmos_saturation::cpovr_ice
real(rp), public cpovr_ice
Definition: scale_atmos_sub_saturation.F90:129
scale_atmos_phy_mp_suzuki10::kphase
integer, public kphase
Definition: scale_atmos_phy_mp_suzuki10.F90:113
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.14