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