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