SCALE-RM
scale_atmos_phy_ae_kajino13.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_atmos_phy_ae_kajino13
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  use scale_const, only: &
21  mwair => const_mdry, & ! molecular weight for dry air
22  mwwat => const_mvap, & ! mean molecular weight for water vapor [g/mol]
23  dnwat => const_dwatr, & ! water density [kg/m3]
24  rgas => const_r, & ! universal gas constant [J/mol/K]
25  stdatmpa => const_pstd, & ! standard pressure [Pa]
26  stdtemp => const_tem00, & ! standard temperature [K]
27  pi => const_pi ! pi
28  use scale_atmos_aerosol, only: &
29  n_ae
30  !-----------------------------------------------------------------------------
31  implicit none
32  private
33  !-----------------------------------------------------------------------------
34  !
35  !++ Public procedure
36  !
37  public :: atmos_phy_ae_kajino13_tracer_setup
38  public :: atmos_phy_ae_kajino13_setup
39  public :: atmos_phy_ae_kajino13_tendency
40  public :: atmos_phy_ae_kajino13_effective_radius
41  public :: atmos_phy_ae_kajino13_mkinit
42  public :: atmos_phy_ae_kajino13_negative_fixer
43 
44  !-----------------------------------------------------------------------------
45  !
46  !++ Public parameters & variables
47  !
48 
49  integer, private :: QA_AE = 0
50 
51  character(len=H_SHORT), public, allocatable :: atmos_phy_ae_kajino13_name(:)
52  character(len=H_MID) , public, allocatable :: atmos_phy_ae_kajino13_desc(:)
53  character(len=H_SHORT), public, allocatable :: atmos_phy_ae_kajino13_unit(:)
54 
55  real(rp), public :: atmos_phy_ae_kajino13_h2so4dt = 5.e-6_rp ! h2so4 production rate (Temporal) [ug/m3/s]
56  real(rp), public :: atmos_phy_ae_kajino13_ocgasdt = 8.e-5_rp ! other condensational bas production rate 16*h2so4dt (see Kajino et al. 2013)
57 ! real(RP), public :: ATMOS_PHY_AE_KAJINO13_c_ratio = 16.0_RP ! ratio of condensable mass to h2so4 (after NPF) [-]
58  real(rp), public :: atmos_phy_ae_kajino13_c_kappa = 0.3_rp ! hygroscopicity of condensable mass [-]
59  logical, public :: atmos_phy_ae_kajino13_flag_npf = .false.
60  logical, public :: atmos_phy_ae_kajino13_flag_cond = .true.
61  logical, public :: atmos_phy_ae_kajino13_flag_coag = .true.
62  logical, public :: atmos_phy_ae_kajino13_flag_ccn_interactive = .true.
63  logical, public :: atmos_phy_ae_kajino13_flag_regeneration = .true.
64  real(rp), public :: atmos_phy_ae_kajino13_dg_reg = 5.e-7_rp ! dg of regenerated aerosol (droplet mode 500 nm) [m]
65  real(rp), public :: atmos_phy_ae_kajino13_sg_reg = 1.6_rp ! sg of regenerated aerosol [-]
66  real(rp), public :: atmos_phy_ae_kajino13_logk_aenucl = -12.4_rp !constant coefficient for kinetic nucleation [-]
67 
68  !-----------------------------------------------------------------------------
69  !
70  !++ Private procedure
71  !
72  !-----------------------------------------------------------------------------
73  !
74  !++ Private parameters & variables
75  !
76  integer, parameter :: gas_ctg = 2
77  integer, parameter :: n_atr = 5
78 
79  integer, parameter :: ic_mix = 1
80  integer, parameter :: ic_sea = 2
81  integer, parameter :: ic_dus = 3
82 
83  integer, parameter :: ig_h2so4 = 1
84  integer, parameter :: ig_cgas = 2
85 
86  integer :: ae_ctg = 1
87  integer, allocatable :: nkap(:)
88  integer, allocatable :: nsiz(:)
89 
90  integer :: atmos_phy_ae_kajino13_nbins_out = 1024
91 
92  ! physical parameters
93 ! real(RP), parameter :: mwair = 28.9628_RP ! molecular weight for dry air
94 ! real(RP), parameter :: mwwat = 18.0153_RP ! mean molecular weight for water vapor [g/mol]
95 ! real(RP), parameter :: dnwat = 1.E3_RP ! water density [kg/m3]
96 ! real(RP), parameter :: rgas = 8.31447_RP ! universal gas constant [J/mol/K]
97 ! real(RP), parameter :: stdatmpa = 101325._RP ! standard pressure [Pa]
98 ! real(RP), parameter :: stdtemp = 273.15_RP ! standard temperature [K]
99 
100  real(rp), parameter :: avo = 6.0221367e23_rp ! avogadro number [#/mol]
101  real(rp), parameter :: boltz = 1.38048e-23_rp ! boltzmann constant [J K-1]
102  real(rp), parameter :: conv_n2m = 1.e6_rp/avo*1.e6_rp &!
103  *98._rp ! ug/m3 = conv_n2m * #/cm3
104  real(rp), parameter :: conv_m2n = 1._rp/conv_n2m ! #/cm3 = conv_m2n * ug/m3
105  real(rp), parameter :: rhod_ae = 1.83_rp ! particle density [g/cm3] sulfate assumed
106  real(rp), parameter :: rho_kg = rhod_ae*1.e3_rp ! particle density [kg/m3] sulfate assumed
107  real(rp), parameter :: grav = 9.80622_rp ! mean gravitational acceleration [m/s2]
108  real(rp), parameter :: conv_ms_vl = 1.e-12_rp/rhod_ae ! mass[ug/m3] to M3(volume)[m3/m3] rhod
109  real(rp), parameter :: conv_vl_ms = rhod_ae/1.e-12_rp ! M3(volume)[m3/m3] to mass[ug/m3]
110  real(rp), parameter :: mwrat_s6 = 96._rp/98._rp ! molecular weight ratio (part/gas) of sulfate
111  real(rp), parameter :: mwrat_s6_i = 1._rp/mwrat_s6 ! molecular weight ratio (gas/part) of sulfate
112  real(rp), parameter :: diffsulf = 9.36e-6_rp ! std. molecular diffusivity of sulfuric acid [m2/s]
113  real(rp), parameter :: mwh2so4 = 98._rp ! molecular weight for h2so4 gas [g/mol]
114  real(rp) :: pi6 != pi / 6._RP ! pi/6
115  real(rp) :: sixpi != 6._RP / pi ! 6/pi
116  real(rp) :: forpi != 4._RP / pi ! 4/pi
117  ! parameters for condensation/coagulation
118  real(rp), parameter :: c1=1.425e-6_rp, c2=0.5039_rp, c3=108.3_rp ! Perry's [Tableo 2-312]
119  real(rp) :: c4 ! Perry's [Tableo 2-312]
120 
121  !(attributes: zero chemistry version)
122 ! integer, parameter :: n_atr = 5 !number of attributes
123  integer, parameter :: ia_m0 = 1 !1. number conc [#/m3]
124  integer, parameter :: ia_m2 = 2 !2. 2nd mom conc [m2/m3]
125  integer, parameter :: ia_m3 = 3 !3. 3rd mom conc [m3/m3]
126  integer, parameter :: ia_ms = 4 !4. mass conc [ug/m3]
127  integer, parameter :: ia_kp = 5 !5. mean kappa * volume [-] ( ia_kp*ia_m3 nisuru)
128  integer, parameter :: ik_out = 1
129  !(set in subroutine aerosol_settings)
130 ! integer :: n_ctg !number of category
131  integer, allocatable :: n_siz(:) !number of size bins (n_ctg)
132  real(rp),allocatable :: d_min(:) !lower bound of 1st size bin (n_ctg)
133  real(rp),allocatable :: d_max(:) !upper bound of last size bin (n_ctg)
134  integer, allocatable :: n_kap(:) !number of kappa bins (n_ctg)
135  real(rp),allocatable :: k_min(:) !lower bound of 1st kappa bin (n_ctg)
136  real(rp),allocatable :: k_max(:) !upper bound of last kappa bin (n_ctg)
137  !(diagnosed in subroutine aerosol_settings)
138 
139  !--- bin settings (lower bound, center, upper bound)
140  real(rp),allocatable :: d_lw(:,:), d_ct(:,:), d_up(:,:) !diameter [m]
141  real(rp),allocatable :: k_lw(:,:), k_ct(:,:), k_up(:,:) !kappa [-]
142  real(rp) :: dlogd, dk !delta log(D), delta K
143 
144  !--- coagulation rule (i+j=k)
145  integer, allocatable :: is_i(:), is_j(:), is_k(:) !(mcomb)
146  integer, allocatable :: ik_i(:), ik_j(:), ik_k(:) !(mcomb)
147  integer, allocatable :: ic_i(:), ic_j(:), ic_k(:) !(mcomb)
148  integer :: mcomb !combinations of sections
149  integer :: is1, is2, mc, ik, is0, ic, ia0
150 
151 ! real(RP) :: m0_init = 0.E0 ! initial total num. conc. of modes (Atk,Acm,Cor) [#/m3]
152 ! real(RP) :: dg_init = 80.E-9 ! initial number equivalen diameters of modes [m]
153 ! real(RP) :: sg_init = 1.6 ! initial standard deviation [-]
154 
155  integer :: is0_reg ! size bin of regenerated aerosol
156 
157  integer :: n_ctg ! number of category
158  integer :: n_trans ! number of total transport variables
159  integer :: n_siz_max ! maximum number of size bins
160  integer :: n_kap_max ! maximum number of kappa bins
161  integer, allocatable :: it_procs2trans(:,:,:,:) !procs to trans conversion
162  integer, allocatable :: ia_trans2procs(:) !trans to procs conversion
163  integer, allocatable :: is_trans2procs(:) !trans to procs conversion
164  integer, allocatable :: ik_trans2procs(:) !trans to procs conversion
165  integer, allocatable :: ic_trans2procs(:) !trans to procs conversion
166  real(rp), allocatable :: rnum_out(:)
167 
168  character(len=H_SHORT),allocatable :: ctg_name(:)
169  real(rp), parameter :: cleannumber = 1.e-3 ! tiny number of aerosol per bin for neglectable mass,
170  ! D =10 um resulted in 1.e-6 ug/m3
171 
172 
173  real(rp), allocatable :: aerosol_procs(:,:,:,:) ! (n_atr,n_siz_max,n_kap_max,n_ctg)
174  real(rp), allocatable :: aerosol_activ(:,:,:,:) ! (n_atr,n_siz_max,n_kap_max,n_ctg)
175  real(rp), allocatable :: emis_procs (:,:,:,:) ! (n_atr,n_siz_max,n_kap_max,n_ctg)
176  real(rp), allocatable :: emis_gas (:) ! emission of gas
177 
178  !-----------------------------------------------------------------------------
179 contains
180  !-----------------------------------------------------------------------------
182  subroutine atmos_phy_ae_kajino13_tracer_setup( QA_AE )
183  use scale_prc, only: &
184  prc_abort
185 
186  integer, intent(out) :: qa_ae
187 
188  integer, allocatable :: aero_idx(:,:,:,:)
189  integer :: ncat_max
190  integer :: nasiz(3), nakap(3)
191  character(len=H_SHORT) :: attribute, catego, aunit
192 
193  namelist / param_atmos_phy_ae_kajino13_tracer / &
194  ae_ctg, &
195  nasiz, &
196  nakap
197 
198  integer :: m, ierr, ik, ic, ia0, is0
199 
200  !---------------------------------------------------------------------------
201 
202  log_newline
203  log_info("ATMOS_PHY_AE_kajino13_tracer_setup",*) 'Setup'
204 
205  ! note: tentatively, aerosol module should be called at all time. we need dummy subprogram.
206  !if ( ATMOS_sw_phy_ae ) then
207 
208  log_warn("ATMOS_PHY_AE_kajino13_tracer_setup",*) '### Kajino13 scheme is still experimental'
209 
210  ncat_max = max( ic_mix, ic_sea, ic_dus )
211 
212  nasiz(:) = 64
213  nakap(:) = 1
214 
215  rewind(io_fid_conf)
216  read(io_fid_conf,nml=param_atmos_phy_ae_kajino13_tracer,iostat=ierr)
217 
218  if( ierr < 0 ) then !--- missing
219  log_info("ATMOS_PHY_AE_kajino13_tracer_setup",*) 'Not found namelist. Default used.'
220  elseif( ierr > 0 ) then !--- fatal error
221  log_error("ATMOS_PHY_AE_kajino13_tracer_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_AE_KAJINO13_TRACER, Check!'
222  call prc_abort
223  end if
224 
225  log_nml(param_atmos_phy_ae_kajino13_tracer)
226 
227  if( ae_ctg > ncat_max ) then
228  log_error("ATMOS_PHY_AE_kajino13_tracer_setup",*) 'AE_CTG should be smaller than', ncat_max+1, 'stop'
229  call prc_abort
230  endif
231 
232  allocate( nsiz(ae_ctg) )
233  allocate( nkap(ae_ctg) )
234 
235  nkap(1:ae_ctg) = nakap(1:ae_ctg)
236  nsiz(1:ae_ctg) = nasiz(1:ae_ctg)
237 
238  if( maxval( nkap ) /= 1 .OR. minval( nkap ) /= 1 ) then
239  log_error("ATMOS_PHY_AE_kajino13_tracer_setup",*) 'NKAP(:) /= 1 is not supported now, Stop!'
240  call prc_abort
241  end if
242 
243  ! do ia0 = 1, N_ATR
244  do ic = 1, ae_ctg
245  do ik = 1, nkap(ic)
246  do is0 = 1, nsiz(ic)
247  qa_ae = qa_ae + n_atr
248  enddo
249  enddo
250  enddo
251  qa_ae = qa_ae + gas_ctg
252 
253  allocate( atmos_phy_ae_kajino13_name(qa_ae) )
254  allocate( atmos_phy_ae_kajino13_desc(qa_ae) )
255  allocate( atmos_phy_ae_kajino13_unit(qa_ae) )
256 
257  n_siz_max = 0
258  n_kap_max = 0
259  do ic = 1, ae_ctg
260  n_siz_max = max(n_siz_max, nsiz(ic))
261  n_kap_max = max(n_kap_max, nkap(ic))
262  enddo
263 
264  allocate( aero_idx(n_atr,ae_ctg,n_kap_max,n_siz_max) )
265  m = 0
266  ! do ia0 = 1, N_ATR
267  do ic = 1, ae_ctg
268  do ik = 1, nkap(ic)
269  do is0 = 1, nsiz(ic)
270  do ia0 = 1, n_atr
271  m = m+1
272  aero_idx(ia0,ic,ik,is0) = m
273  enddo
274  enddo
275  enddo
276  enddo
277 
278  !-----------------------------------------------------------------------------
279  !
280  !++ calculate each category and aerosol
281  !
282  !-----------------------------------------------------------------------------
283  ic = qa_ae-gas_ctg+ig_h2so4
284  write(atmos_phy_ae_kajino13_unit(ic),'(a)') 'kg/kg'
285  ic = qa_ae-gas_ctg+ig_cgas
286  write(atmos_phy_ae_kajino13_unit(ic),'(a)') 'kg/kg'
287 
288  !do ia0 = 1, N_ATR
289  ! if( ia0 == 1 ) then
290  ! write(attribute,'(a)') "Number"
291  ! write(aunit,'(a)') "num/kg"
292  ! elseif( ia0 == 2 ) then
293  ! write(attribute,'(a)') "Section"
294  ! write(aunit,'(a)') "m2/kg"
295  ! elseif( ia0 == 3 ) then
296  ! write(attribute,'(a)') "Volume"
297  ! write(aunit,'(a)') "m3/kg"
298  ! elseif( ia0 == 4 ) then
299  ! write(attribute,'(a)') "Mass"
300  ! write(aunit,'(a)') "kg/kg"
301  ! elseif( ia0 == 5 ) then
302  ! write(attribute,'(a)') "kpXmass"
303  ! write(aunit,'(a)') "kg/kg"
304  ! endif
305  do ic = 1, ae_ctg !aerosol category
306  do ik = 1, nkap(ic) !kappa bin
307  do is0 = 1, nsiz(ic)
308  do ia0 = 1, n_atr
309  if( ia0 == 1 ) then
310  write(attribute,'(a)') "Number"
311  write(aunit,'(a)') "num/kg"
312  elseif( ia0 == 2 ) then
313  write(attribute,'(a)') "Section"
314  write(aunit,'(a)') "m2/kg"
315  elseif( ia0 == 3 ) then
316  write(attribute,'(a)') "Volume"
317  write(aunit,'(a)') "m3/kg"
318  elseif( ia0 == 4 ) then
319  write(attribute,'(a)') "Mass"
320  write(aunit,'(a)') "kg/kg"
321  elseif( ia0 == 5 ) then
322  write(attribute,'(a)') "kXm"
323  write(aunit,'(a)') "kg/kg"
324  endif
325  if( ic == ic_mix ) then
326  write(catego,'(a)') "Sulf_"
327  elseif( ic == ic_sea ) then
328  write(catego,'(a)') "Salt_"
329  elseif( ic == ic_dus ) then
330  write(catego,'(a)') "Dust_"
331  endif
332  write(atmos_phy_ae_kajino13_unit(aero_idx(ia0,ic,ik,is0)),'(a)') trim(aunit)
333  write(atmos_phy_ae_kajino13_name(aero_idx(ia0,ic,ik,is0)),'(a,a,i0)') trim(catego), trim(attribute), is0
334  write(atmos_phy_ae_kajino13_desc(aero_idx(ia0,ic,ik,is0)),'(a,a,a,i0)') trim(attribute), ' mixing radio of ', trim(catego), is0
335  enddo
336  enddo
337  enddo
338  enddo
339  ic = qa_ae-gas_ctg+ig_h2so4
340  write(atmos_phy_ae_kajino13_name(ic),'(a)') 'H2SO4_Gas'
341  ic = qa_ae-gas_ctg+ig_cgas
342  write(atmos_phy_ae_kajino13_name(ic),'(a)') 'Condensable_GAS'
343 
344  ic = qa_ae-gas_ctg+ig_h2so4
345  write(atmos_phy_ae_kajino13_desc(ic),'(a)') 'Mixing ratio of H2SO4 Gas'
346  ic = qa_ae-gas_ctg+ig_cgas
347  write(atmos_phy_ae_kajino13_desc(ic),'(a)') 'Mixing ratio of Condensable GAS'
348 
349  deallocate(aero_idx)
350 
351  return
352  end subroutine atmos_phy_ae_kajino13_tracer_setup
353 
354  !-----------------------------------------------------------------------------
356  subroutine atmos_phy_ae_kajino13_setup
357  use scale_prc, only: &
358  prc_abort
359  implicit none
360 
361  real(rp), allocatable :: d_min_inp(:)
362  real(rp), allocatable :: d_max_inp(:)
363  real(rp), allocatable :: k_min_inp(:)
364  real(rp), allocatable :: k_max_inp(:)
365  integer , allocatable :: n_kap_inp(:)
366 
367  real(rp), parameter :: d_min_def = 1.e-9_rp ! default lower bound of 1st size bin
368  real(rp), parameter :: d_max_def = 1.e-5_rp ! upper bound of last size bin
369  integer , parameter :: n_kap_def = 1 ! number of kappa bins
370  real(rp), parameter :: k_min_def = 0.e0_rp ! lower bound of 1st kappa bin
371  real(rp), parameter :: k_max_def = 1.e0_rp ! upper bound of last kappa bin
372 
373  namelist / param_atmos_phy_ae_kajino13 / &
374  atmos_phy_ae_kajino13_h2so4dt, &
375  atmos_phy_ae_kajino13_ocgasdt, &
376 ! ATMOS_PHY_AE_KAJINO13_c_ratio, &
377  atmos_phy_ae_kajino13_c_kappa, &
378 ! d_min_inp, &
379 ! d_max_inp, &
380 ! k_min_inp, &
381 ! k_max_inp, &
382 ! n_kap_inp, &
383  atmos_phy_ae_kajino13_flag_npf, &
384  atmos_phy_ae_kajino13_flag_cond, &
385  atmos_phy_ae_kajino13_flag_coag, &
386  atmos_phy_ae_kajino13_flag_ccn_interactive, &
387  atmos_phy_ae_kajino13_flag_regeneration, &
388  atmos_phy_ae_kajino13_dg_reg, &
389  atmos_phy_ae_kajino13_sg_reg, &
390  atmos_phy_ae_kajino13_logk_aenucl, &
391  atmos_phy_ae_kajino13_nbins_out
392 
393  integer :: it, ierr
394  !---------------------------------------------------------------------------
395 
396  log_newline
397  log_info("ATMOS_PHY_AE_kajino13_setup",*) 'Setup'
398  log_info("ATMOS_PHY_AE_kajino13_setup",*) 'Kajino(2013) scheme'
399 
400  !--- setup parameter
401  pi6 = pi / 6._rp ! pi/6
402  sixpi = 6._rp / pi ! 6/pi
403  forpi = 4._rp / pi ! 4/pi
404 
405 
406  n_ctg = ae_ctg
407  allocate( rnum_out(atmos_phy_ae_kajino13_nbins_out) )
408  allocate( n_siz(n_ctg) )
409  allocate( d_min(n_ctg) )
410  allocate( d_max(n_ctg) )
411  allocate( n_kap(n_ctg) )
412  allocate( k_min(n_ctg) )
413  allocate( k_max(n_ctg) )
414  allocate( d_min_inp(n_ctg) )
415  allocate( d_max_inp(n_ctg) )
416  allocate( n_kap_inp(n_ctg) )
417  allocate( k_min_inp(n_ctg) )
418  allocate( k_max_inp(n_ctg) )
419  allocate( ctg_name(n_ctg) )
420 
421  n_siz(1:n_ctg) = nsiz(1:n_ctg) ! number of size bins
422  d_min(1:n_ctg) = d_min_def ! lower bound of 1st size bin
423  d_max(1:n_ctg) = d_max_def ! upper bound of last size bin
424  n_kap(1:n_ctg) = n_kap_def ! number of kappa bins
425  k_min(1:n_ctg) = k_min_def ! lower bound of 1st kappa bin
426  k_max(1:n_ctg) = k_max_def ! upper bound of last kappa bin
427 
428  do it = 1, n_ctg
429  if( n_ctg == 1 ) then
430  write(ctg_name(it),'(a)') "Sulfate"
431  elseif( n_ctg == 2 ) then
432  write(ctg_name(it),'(a)') "Seasalt"
433  elseif( n_ctg == 3 ) then
434  write(ctg_name(it),'(a)') "Dust"
435  endif
436  enddo
437 
438  !--- read namelist
439  rewind(io_fid_conf)
440  read(io_fid_conf,nml=param_atmos_phy_ae_kajino13,iostat=ierr)
441  if( ierr < 0 ) then !--- missing
442  log_info("ATMOS_PHY_AE_kajino13_setup",*) 'Not found namelist. Default used.'
443  elseif( ierr > 0 ) then !--- fatal error
444  log_error("ATMOS_PHY_AE_kajino13_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_AE_KAJINO13. Check!'
445  call prc_abort
446  endif
447  log_nml(param_atmos_phy_ae_kajino13)
448 
449  !--- now only the default setting is supported
450 ! n_siz(1:n_ctg) = NSIZ(1:n_ctg) ! number of size bins
451 ! d_min(1:n_ctg) = d_min_inp(1:n_ctg) ! lower bound of 1st size bin
452 ! d_max(1:n_ctg) = d_max_inp(1:n_ctg) ! upper bound of last size bin
453 ! n_kap(1:n_ctg) = n_kap_inp(1:n_ctg) ! number of kappa bins
454 ! k_min(1:n_ctg) = k_min_inp(1:n_ctg) ! lower bound of 1st kappa bin
455 ! k_max(1:n_ctg) = k_max_inp(1:n_ctg) ! upper bound of last kappa bin
456 
457  !--- diagnose parameters (n_trans, n_siz_max, n_kap_max)
458  n_trans = 0
459  n_siz_max = 0
460  n_kap_max = 0
461  do ic = 1, n_ctg
462  n_trans = n_trans + n_siz(ic) * n_kap(ic) * n_atr
463  n_siz_max = max(n_siz_max, n_siz(ic))
464  n_kap_max = max(n_kap_max, n_kap(ic))
465  enddo
466 
467  !--- bin settings
468  allocate(d_lw(n_siz_max,n_ctg))
469  allocate(d_ct(n_siz_max,n_ctg))
470  allocate(d_up(n_siz_max,n_ctg))
471  allocate(k_lw(n_kap_max,n_ctg))
472  allocate(k_ct(n_kap_max,n_ctg))
473  allocate(k_up(n_kap_max,n_ctg))
474  d_lw(:,:) = 0.0_rp
475  d_ct(:,:) = 0.0_rp
476  d_up(:,:) = 0.0_rp
477  k_lw(:,:) = 0.0_rp
478  k_ct(:,:) = 0.0_rp
479  k_up(:,:) = 0.0_rp
480 
481  do ic = 1, n_ctg
482  dlogd = (log(d_max(ic)) - log(d_min(ic)))/float(n_siz(ic))
483  do is0 = 1, n_siz(ic) !size bin
484  d_lw(is0,ic) = exp(log(d_min(ic))+dlogd* float(is0-1) )
485  d_ct(is0,ic) = exp(log(d_min(ic))+dlogd*(float(is0)-0.5_rp))
486  d_up(is0,ic) = exp(log(d_min(ic))+dlogd* float(is0) )
487  enddo !is (1:n_siz(ic))
488 
489  dk = (k_max(ic) - k_min(ic))/float(n_kap(ic))
490  do ik = 1, n_kap(ic) !size bin
491  k_lw(ik,ic) = k_min(ic) + dk * float(ik-1)
492  k_ct(ik,ic) = k_min(ic) + dk *(float(ik)-0.5_rp)
493  k_up(ik,ic) = k_min(ic) + dk * float(ik)
494  enddo !ik (1:n_kap(ic))
495 
496  enddo !ic (1:n_ctg)
497 
498 !find size bin of regenerated aerosols
499  do is0 = 1, n_siz(ic_mix)
500  if ( atmos_phy_ae_kajino13_dg_reg >= d_lw(is0,ic_mix) .AND. &
501  atmos_phy_ae_kajino13_dg_reg < d_up(is0,ic_mix) ) then
502  is0_reg = is0
503  endif !d_lw < ATMOS_PHY_AE_KAJINO13_dg_reg < d_up
504  enddo
505 
506  !--- coagulation rule
507  ! [ NOTE: current version has one category and single
508  ! hygroscopicity bins and thus does not consider inter-category
509  ! nor inter-hygroscopicity-section coagulation ]
510  mcomb = 0
511  do ic = 1, n_ctg !=1
512  do ik = 1, n_kap(ic) !=1
513  do is2 = 1, n_siz(ic)
514  do is1 = 1, n_siz(ic)
515  if ( d_ct(is2,ic) >= d_ct(is1,ic) ) then
516  mcomb = mcomb + 1
517  endif !d_ct(is2) >= d_ct(is1)
518  enddo !is1 (1:n_siz(ic) )
519  enddo !is2 (1:n_siz(ic) )
520  enddo !ik(1:n_kap(ic))
521  enddo !ic(1:n_ctg)
522 
523  allocate(is_i(mcomb))
524  allocate(is_j(mcomb))
525  allocate(is_k(mcomb))
526  allocate(ik_i(mcomb))
527  allocate(ik_j(mcomb))
528  allocate(ik_k(mcomb))
529  allocate(ic_i(mcomb))
530  allocate(ic_j(mcomb))
531  allocate(ic_k(mcomb))
532 
533  mc = 0
534  do ic = 1, n_ctg !=1
535  do ik = 1, n_kap(ic) !=1
536  do is2 = 1, n_siz(ic)
537  do is1 = 1, n_siz(ic)
538  if ( d_ct(is2,ic) >= d_ct(is1,ic) ) then
539  mc = mc + 1
540  is_i(mc) = is1
541  ik_i(mc) = ik
542  ic_i(mc) = ic
543  is_j(mc) = is2
544  ik_j(mc) = ik
545  ic_j(mc) = ic
546  is_k(mc) = is2
547  ik_k(mc) = ik
548  ic_k(mc) = ic
549  endif !d_ct(is2) >= d_ct(is1)
550  enddo !is1 (1:n_siz(ic) )
551  enddo !is2 (1:n_siz(ic) )
552  enddo !ik(1:n_kap(ic))
553  enddo !ic(1:n_ctg)
554 
555  !--- gas concentration
556 ! conc_h2so4 = 0.0_RP
557 ! conc_cgas = 0.0_RP
558 
559  allocate( it_procs2trans(n_atr,n_siz_max,n_kap_max,n_ctg) )
560  allocate( ia_trans2procs(n_trans) )
561  allocate( is_trans2procs(n_trans) )
562  allocate( ik_trans2procs(n_trans) )
563  allocate( ic_trans2procs(n_trans) )
564 
565  it_procs2trans(:,:,:,:)= -999
566  ia_trans2procs(:) = 0
567  is_trans2procs(:) = 0
568  ik_trans2procs(:) = 0
569  ic_trans2procs(:) = 0
570 
571  !--- get pointer for trans2procs, procs2trans
572  it = 0
573  do ic = 1, n_ctg !aerosol category
574  do ik = 1, n_kap(ic) !kappa bin
575  do is0 = 1, n_siz(ic) !size bin
576  do ia0 = 1, n_atr !attributes
577  it = it + 1
578  it_procs2trans(ia0,is0,ik,ic)= it
579  ia_trans2procs(it) = ia0
580  is_trans2procs(it) = is0
581  ik_trans2procs(it) = ik
582  ic_trans2procs(it) = ic
583  enddo !ia (1:n_atr_prog )
584  enddo !is (1:n_siz(ic) )
585  enddo !ik (1:n_kap(ic) )
586  enddo !ic (1:n_ctg )
587 
588  allocate( aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg) )
589  allocate( aerosol_activ(n_atr,n_siz_max,n_kap_max,n_ctg) )
590  allocate( emis_procs(n_atr,n_siz_max,n_kap_max,n_ctg) )
591  allocate( emis_gas(gas_ctg) )
592  aerosol_procs(:,:,:,:) = 0.0_rp
593  aerosol_activ(:,:,:,:) = 0.0_rp
594  emis_procs(:,:,:,:) = 0.0_rp
595  emis_gas(:) = 0.0_rp
596 
597  return
598  end subroutine atmos_phy_ae_kajino13_setup
599 
600  !-----------------------------------------------------------------------------
602  subroutine atmos_phy_ae_kajino13_tendency( &
603  KA, KS, KE, &
604  IA, IS, IE, &
605  JA, JS, JE, &
606  QA_AE, &
607  TEMP, &
608  PRES, &
609  QDRY, &
610  NREG, &
611  DENS, &
612  QV, &
613  QTRC, &
614  EMIT, &
615  dt, &
616  RHOQ_t_AE, &
617  CN, &
618  CCN )
619  use scale_atmos_saturation, only: &
620  pres2qsat_liq => atmos_saturation_pres2qsat_liq
621  use scale_file_history, only: &
622  file_history_in
623  implicit none
624  integer, intent(in) :: ka, ks, ke
625  integer, intent(in) :: ia, is, ie
626  integer, intent(in) :: ja, js, je
627  integer, intent(in) :: qa_ae
628  real(rp), intent(in) :: temp(ka,ia,ja)
629  real(rp), intent(in) :: pres(ka,ia,ja)
630  real(rp), intent(in) :: qdry(ka,ia,ja)
631  real(rp), intent(in) :: nreg(ka,ia,ja)
632  real(rp), intent(in) :: dens(ka,ia,ja)
633  real(rp), intent(in) :: qv (ka,ia,ja)
634  real(rp), intent(in) :: qtrc(ka,ia,ja,qa_ae)
635  real(rp), intent(in) :: emit(ka,ia,ja,qa_ae)
636  real(dp), intent(in) :: dt
637 
638  real(rp), intent(out) :: rhoq_t_ae(ka,ia,ja,qa_ae)
639  real(rp), intent(out) :: cn(ka,ia,ja)
640  real(rp), intent(out) :: ccn(ka,ia,ja)
641 
642  real(rp) :: qtrc0(ka,ia,ja,qa_ae)
643  real(rp) :: qtrc1(ka,ia,ja,qa_ae)
644 
645  !--- local
646  real(rp) :: ssliq_ae(ka,ia,ja)
647  real(rp) :: rrhog(ka,ia,ja)
648 ! real(RP) :: t_ccn, t_cn
649  real(rp) :: qsat_tmp
650  real(rp) :: m0_reg, m2_reg, m3_reg !regenerated aerosols [m^k/m3]
651  real(rp) :: ms_reg !regenerated aerosol mass [ug/m3]
652  real(rp) :: reg_factor_m2,reg_factor_m3 !to save cpu time for moment conversion
653  !--- aerosol variables
654  real(rp) :: total_aerosol_mass(ka,ia,ja,n_ctg)
655  real(rp) :: total_aerosol_number(ka,ia,ja,n_ctg)
656  real(rp) :: total_emit_aerosol_mass(ka,ia,ja,n_ctg)
657  real(rp) :: total_emit_aerosol_number(ka,ia,ja,n_ctg)
658  !--- gas
659 ! real(RP),allocatable :: conc_h2so4(:,:,:,:) !concentration [ug/m3]
660 ! real(RP) :: conc_gas(KA,IA,JA,GAS_CTG) !concentration [ug/m3]
661  real(rp) :: conc_gas(gas_ctg) !concentration [ug/m3]
662  integer :: i, j, k, iq, it
663 
664  log_progress(*) 'atmosphere / physics / aerosol / kajino13'
665 
666  aerosol_procs(:,:,:,:) = 0.0_rp
667  aerosol_activ(:,:,:,:) = 0.0_rp
668  emis_procs(:,:,:,:) = 0.0_rp
669  emis_gas(:) = 0.0_rp
670 
671  do iq = 1, qa_ae
672  do j = js, je
673  do i = is, ie
674  do k = ks, ke
675  qtrc0(k,i,j,iq) = qtrc(k,i,j,iq) ! save
676  enddo
677  enddo
678  enddo
679  enddo
680 
681  reg_factor_m2 = atmos_phy_ae_kajino13_dg_reg**2._rp * exp( 2.0_rp *(log(atmos_phy_ae_kajino13_sg_reg)**2._rp)) !m0_reg to m2_reg
682  reg_factor_m3 = atmos_phy_ae_kajino13_dg_reg**3._rp * exp( 4.5_rp *(log(atmos_phy_ae_kajino13_sg_reg)**2._rp)) !m0_reg to m3_reg
683 
684  !--- convert SCALE variable to zerochem variable
685 
686  do j = js, je
687  do i = is, ie
688 
689  do k = ks, ke
690  rrhog(k,i,j) = 1.0_rp / dens(k,i,j)
691  enddo
692  do k = ks, ke
693  !--- calculate super saturation of water
694  call pres2qsat_liq( temp(k,i,j), pres(k,i,j), qdry(k,i,j), & ! [IN]
695  qsat_tmp ) ! [OUT]
696  ssliq_ae(k,i,j) = qv(k,i,j)/qsat_tmp - 1.0_rp
697  enddo
698 
699  enddo
700  enddo
701 
702  ! tiny number, tiny mass
703  do j = js, je
704  do i = is, ie
705  do k = ks, ke
706  do ic = 1, n_ctg !aerosol category
707  do ik = 1, n_kap(ic) !kappa bin
708  do is0 = 1, n_siz(ic) !size bin
709  if (qtrc0(k,i,j,it_procs2trans(ia_m0,is0,ik,ic))*dens(k,i,j) < cleannumber) then
710  do ia0 = 1, n_atr
711  qtrc0(k,i,j,it_procs2trans(ia0,is0,ik,ic)) = 0._rp !to save cpu time and avoid underflow
712  enddo !ia0 (1:n_atr )
713  endif
714  enddo !is (1:n_siz(ic) )
715  enddo !ik (1:n_kap(ic) )
716  enddo !ic (1:n_ctg )
717  enddo
718  enddo
719  enddo
720 
721  do iq = 1, qa_ae
722  do j = js, je
723  do i = is, ie
724  do k = ks, ke
725  qtrc1(k,i,j,iq) = qtrc0(k,i,j,iq) ! save
726  enddo
727  enddo
728  enddo
729  enddo
730 
731  !---- Calculate aerosol processs
732  cn(:,:,:) = 0.0_rp
733  ccn(:,:,:) = 0.0_rp
734  do k = ks, ke
735  do j = js, je
736  do i = is, ie
737 
738  !--- aerosol_trans at initial time
739  ! [xx/kg] -> [xx/m3]
740  do it = 1, n_trans
741  aerosol_procs(ia_trans2procs(it), &
742  is_trans2procs(it), &
743  ik_trans2procs(it), &
744  ic_trans2procs(it)) = qtrc1(k,i,j,it)*dens(k,i,j)
745  emis_procs(ia_trans2procs(it), &
746  is_trans2procs(it), &
747  ik_trans2procs(it), &
748  ic_trans2procs(it)) = emit(k,i,j,it)*dens(k,i,j)
749  enddo !it(1:n_trans)
750  ! mixing ratio [kg/kg] -> concentration [ug/m3]
751  conc_gas(1:gas_ctg) &
752  = qtrc1(k,i,j,qa_ae-gas_ctg+1:qa_ae-gas_ctg+gas_ctg)*dens(k,i,j)*1.e+9_rp
753 
754  emis_gas(1:gas_ctg) = emit(k,i,j,qa_ae-gas_ctg+ig_h2so4:qa_ae-gas_ctg+ig_cgas)
755 
756  call aerosol_zerochem( &
757  dt, & !--- in
758  temp(k,i,j), & !--- in
759  pres(k,i,j), & !--- in
760  ssliq_ae(k,i,j), & !--- in
761  atmos_phy_ae_kajino13_flag_npf, & !--- in
762  atmos_phy_ae_kajino13_flag_cond, & !--- in
763  atmos_phy_ae_kajino13_flag_coag, & !--- in
764  aerosol_procs, & !--- inout
765  conc_gas, & !--- inout
766  emis_procs, & !--- out
767  emis_gas, & !--- out
768  aerosol_activ ) !--- out
769 
770 ! aerosol loss due to activation to cloud droplets
771  if (atmos_phy_ae_kajino13_flag_ccn_interactive) then
772  do is0 = 1, n_siz(ic_mix)
773  do ia0 = 1, n_atr !attributes
774  aerosol_activ(ia0,is0,ik_out,ic_mix) = min(max(0._rp, aerosol_activ(ia0,is0,ik_out,ic_mix)), &
775  aerosol_procs(ia0,is0,ik_out,ic_mix) )
776 
777  aerosol_procs(ia0,is0,ik_out,ic_mix) = &
778  aerosol_procs(ia0,is0,ik_out,ic_mix) - aerosol_activ(ia0,is0,ik_out,ic_mix)
779  enddo
780  enddo
781  endif !flag_ccn_interactive
782 
783 ! aerosol regeneration due to evaporation of cloud droplets
784 ! using prescribed size parameters and to internal mixture category (ic_mix)
785  if (atmos_phy_ae_kajino13_flag_regeneration) then
786  m0_reg = nreg(k,i,j) !#/m3
787 ! m2_reg = m0_reg * ATMOS_PHY_AE_KAJINO13_dg_reg**2._RP * exp( 2.0_RP *(log(ATMOS_PHY_AE_KAJINO13_sg_reg)**2._RP)) !m2/m3
788 ! m3_reg = m0_reg * ATMOS_PHY_AE_KAJINO13_dg_reg**3._RP * exp( 4.5_RP *(log(ATMOS_PHY_AE_KAJINO13_sg_reg)**2._RP)) !m3/m3
789  m2_reg = m0_reg * reg_factor_m2 !m2/m3
790  m3_reg = m0_reg * reg_factor_m3 !m3/m3
791  ms_reg = m3_reg * pi6 * conv_vl_ms !ug/m3
792  aerosol_procs(ia_m0,is0_reg,ik_out,ic_mix) = &
793  aerosol_procs(ia_m0,is0_reg,ik_out,ic_mix) + m0_reg !#/m3
794  aerosol_procs(ia_m2,is0_reg,ik_out,ic_mix) = &
795  aerosol_procs(ia_m2,is0_reg,ik_out,ic_mix) + m2_reg !m2/m3
796  aerosol_procs(ia_m3,is0_reg,ik_out,ic_mix) = &
797  aerosol_procs(ia_m3,is0_reg,ik_out,ic_mix) + m3_reg !m3/m3
798  aerosol_procs(ia_ms,is0_reg,ik_out,ic_mix) = &
799  aerosol_procs(ia_ms,is0_reg,ik_out,ic_mix) + ms_reg !ug/m3
800 ! additional attirbute to be added (ia_kp)
801  endif !flag_regeneration
802 
803 ! diagnosed variables
804  do is0 = 1, n_siz(ic_mix)
805  ccn(k,i,j) = ccn(k,i,j) + aerosol_activ(ia_m0,is0,ik_out,ic_mix)
806  cn(k,i,j) = cn(k,i,j) + aerosol_procs(ia_m0,is0,ik_out,ic_mix)
807  enddo
808 
809 ! call trans_ccn(aerosol_procs, aerosol_activ, t_ccn, t_cn, &
810 ! n_ctg, n_kap_max, n_siz_max, N_ATR, &
811 ! ic_mix, ia_m0, ia_m2, ia_m3, ik_out, n_siz, &
812 ! rnum_out, ATMOS_PHY_AE_KAJINO13_nbins_out)
813 
814 ! CN(k,i,j) = t_cn
815 ! CCN(k,i,j) = t_ccn
816 
817  ! [xx/m3] -> [xx/kg]
818  do ic = 1, n_ctg !category
819  do ik = 1, n_kap(ic) !kappa bin
820  do is0 = 1, n_siz(ic) !size bin
821  do ia0 = 1, n_atr !attributes
822  qtrc1(k,i,j,it_procs2trans(ia0,is0,ik,ic)) = aerosol_procs(ia0,is0,ik,ic) / dens(k,i,j)
823  enddo !ia (1:N_ATR )
824  enddo !is (1:n_siz(ic) )
825  enddo !ik (1:n_kap(ic) )
826  enddo !ic (1:n_ctg )
827  ! [ug/m3] -> mixing ratio [kg/kg]
828  qtrc1(k,i,j,qa_ae-gas_ctg+1:qa_ae-gas_ctg+gas_ctg) = conc_gas(1:gas_ctg) / dens(k,i,j)*1.e-9_rp
829 
830  ! tiny number, tiny mass
831  do ic = 1, n_ctg !aerosol category
832  do ik = 1, n_kap(ic) !kappa bin
833  do is0 = 1, n_siz(ic) !size bin
834  if (qtrc1(k,i,j,it_procs2trans(ia_m0,is0,ik,ic))*dens(k,i,j) < cleannumber) then
835  do ia0 = 1, n_atr
836  qtrc1(k,i,j,it_procs2trans(ia0,is0,ik,ic)) = 0._rp !to save cpu time and avoid underflow
837  enddo !ia0 (1:n_atr )
838  endif
839  enddo !is (1:n_siz(ic) )
840  enddo !ik (1:n_kap(ic) )
841  enddo !ic (1:n_ctg )
842 
843  ! for history
844  total_aerosol_mass(k,i,j,:) = 0.0_rp
845  total_aerosol_number(k,i,j,:) = 0.0_rp
846  total_emit_aerosol_mass(k,i,j,:) = 0.0_rp
847  total_emit_aerosol_number(k,i,j,:) = 0.0_rp
848  do ic = 1, n_ctg
849  do ik = 1, n_kap(ic)
850  do is0 = 1, n_siz(ic)
851  total_aerosol_mass(k,i,j,ic) = total_aerosol_mass(k,i,j,ic) &
852  + qtrc1(k,i,j,it_procs2trans(ia_ms,is0,ik,ic))
853  total_aerosol_number(k,i,j,ic) = total_aerosol_number(k,i,j,ic) &
854  + qtrc1(k,i,j,it_procs2trans(ia_m0,is0,ik,ic))
855  total_emit_aerosol_mass(k,i,j,ic) = total_emit_aerosol_mass(k,i,j,ic) &
856  + emit(k,i,j,it_procs2trans(ia_ms,is0,ik,ic))
857  total_emit_aerosol_number(k,i,j,ic) = total_emit_aerosol_number(k,i,j,ic) &
858  + emit(k,i,j,it_procs2trans(ia_m0,is0,ik,ic))
859  enddo
860  enddo
861  enddo
862 
863  enddo
864  enddo
865  enddo
866 
867  do ic = 1, n_ctg
868  call file_history_in( total_aerosol_mass(:,:,:,ic), trim(ctg_name(ic))//'_mass', 'Total mass mixing ratio of aerosol', 'kg/kg' )
869  call file_history_in( total_aerosol_number(:,:,:,ic), trim(ctg_name(ic))//'_number', 'Total number mixing ratio of aerosol', 'num/kg' )
870  call file_history_in( total_emit_aerosol_mass(:,:,:,ic), trim(ctg_name(ic))//'_mass_emit', 'Total mass mixing ratio of emitted aerosol', 'kg/kg' )
871  call file_history_in( total_emit_aerosol_number(:,:,:,ic), trim(ctg_name(ic))//'_number_emit', 'Total number mixing ratio of emitted aerosol', 'num/kg' )
872  enddo
873 
874  call file_history_in( emit(:,:,:,qa_ae-gas_ctg+ig_h2so4), 'H2SO4_emit', 'Emission ratio of H2SO4 gas', 'ug/m3/s' )
875  call file_history_in( emit(:,:,:,qa_ae-gas_ctg+ig_cgas), 'CGAS_emit', 'Emission ratio of Condensabule gas', 'ug/m3/s' )
876 
877 
878  do iq = 1, qa_ae
879  do j = js, je
880  do i = is, ie
881  do k = ks, ke
882  rhoq_t_ae(k,i,j,iq) = ( qtrc1(k,i,j,iq) - qtrc0(k,i,j,iq) ) * dens(k,i,j) / dt
883  enddo
884  enddo
885  enddo
886  enddo
887 
888  return
889  end subroutine atmos_phy_ae_kajino13_tendency
890 
891  !-----------------------------------------------------------------------------
893  subroutine atmos_phy_ae_kajino13_effective_radius( &
894  KA, IA, JA, QA_AE, &
895  QTRC, RH, &
896  Re )
897  use scale_atmos_aerosol, only: &
898  n_ae
899  implicit none
900  integer, intent(in) :: ka, ia, ja, qa_ae
901 
902  real(rp), intent(in) :: qtrc(ka,ia,ja,qa_ae) ! tracer mass concentration [kg/kg]
903  real(rp), intent(in) :: rh (ka,ia,ja) ! relative humidity (0-1)
904 
905  real(rp), intent(out) :: re (ka,ia,ja,n_ae) ! effective radius
906  !---------------------------------------------------------------------------
907 
908  re(:,:,:,:) = 0.0_rp
909 
910 ! Re(:,:,:,I_ae_seasalt) = 2.E-4_RP
911 ! Re(:,:,:,I_ae_dust ) = 4.E-6_RP
912 ! Re(:,:,:,I_ae_bc ) = 4.E-8_RP
913 ! Re(:,:,:,I_ae_oc ) = RH(:,:,:)
914 ! Re(:,:,:,I_ae_sulfate) = RH(:,:,:)
915 
916  return
917  end subroutine atmos_phy_ae_kajino13_effective_radius
918 
919  !-----------------------------------------------------------------------------
920  subroutine atmos_phy_ae_kajino13_mkinit( &
921  KA, KS, KE, &
922  IA, IS, IE, &
923  JA, JS, JE, &
924  QA_AE, &
925  DENS, &
926  TEMP, &
927  PRES, &
928  QDRY, &
929  QV, &
930  m0_init, &
931  dg_init, &
932  sg_init, &
933  d_min_inp, &
934  d_max_inp, &
935  k_min_inp, &
936  k_max_inp, &
937  n_kap_inp, &
938  QTRC, &
939  CCN )
940  use scale_const, only: &
941  pi => const_pi
942  use scale_atmos_saturation, only: &
943  saturation_pres2qsat_liq => atmos_saturation_pres2qsat_liq
944  implicit none
945 
946  integer, intent(in) :: ka, ks, ke
947  integer, intent(in) :: ia, is, ie
948  integer, intent(in) :: ja, js, je
949  integer, intent(in) :: qa_ae
950 
951  real(rp), intent(in) :: dens(ka,ia,ja)
952  real(rp), intent(in) :: temp(ka,ia,ja)
953  real(rp), intent(in) :: pres(ka,ia,ja)
954  real(rp), intent(in) :: qdry(ka,ia,ja)
955  real(rp), intent(in) :: qv (ka,ia,ja)
956  real(rp), intent(in) :: m0_init ! initial total num. conc. of modes (Atk,Acm,Cor) [#/m3]
957  real(rp), intent(in) :: dg_init ! initial number equivalen diameters of modes [m]
958  real(rp), intent(in) :: sg_init ! initial standard deviation [-]
959  real(rp), intent(in) :: d_min_inp(3)
960  real(rp), intent(in) :: d_max_inp(3)
961  real(rp), intent(in) :: k_min_inp(3)
962  real(rp), intent(in) :: k_max_inp(3)
963  integer, intent(in) :: n_kap_inp(3)
964 
965  real(rp), intent(out) :: qtrc(ka,ia,ja,qa_ae)
966  real(rp), intent(out) :: ccn (ka,ia,ja)
967 
968  integer, parameter :: ia_m0 = 1 ! 1. number conc [#/m3]
969  integer, parameter :: ia_m2 = 2 ! 2. 2nd mom conc [m2/m3]
970  integer, parameter :: ia_m3 = 3 ! 3. 3rd mom conc [m3/m3]
971  integer, parameter :: ia_ms = 4 ! 4. mass conc [ug/m3]
972  integer, parameter :: ia_kp = 5
973  integer, parameter :: ik_out = 1
974 
975  real(rp) :: c_kappa = 0.3_rp ! hygroscopicity of condensable mass [-]
976  real(rp), parameter :: cleannumber = 1.e-3_rp
977  ! local variables
978  real(rp), parameter :: rhod_ae = 1.83_rp ! particle density [g/cm3] sulfate assumed
979  real(rp), parameter :: conv_vl_ms = rhod_ae/1.e-12_rp ! M3(volume)[m3/m3] to mass[m3/m3]
980 
981  real(rp) :: m0t, dgt, sgt, m2t, m3t, mst
982  !--- gas
983  real(rp) :: conc_gas(gas_ctg) !concentration [ug/m3]
984  !--- bin settings (lower bound, center, upper bound)
985  real(rp), allocatable :: d_lw(:,:), d_ct(:,:), d_up(:,:) !diameter [m]
986  real(rp), allocatable :: k_lw(:,:), k_ct(:,:), k_up(:,:) !kappa [-]
987  real(rp) :: dlogd, dk !delta log(D), delta K
988  real(rp), allocatable :: d_min(:) !lower bound of 1st size bin (n_ctg)
989  real(rp), allocatable :: d_max(:) !upper bound of last size bin (n_ctg)
990  integer, allocatable :: n_kap(:) !number of kappa bins (n_ctg)
991  real(rp), allocatable :: k_min(:) !lower bound of 1st kappa bin (n_ctg)
992  real(rp), allocatable :: k_max(:)
993 
994  real(rp) :: qsat_tmp, ssliq
995  real(rp) :: pi6
996 
997  integer :: ia0, ik, is0, ic, k, i, j, it
998  !---------------------------------------------------------------------------
999 
1000 
1001  call atmos_phy_ae_kajino13_setup
1002 
1003 
1004  pi6 = pi / 6._rp
1005  n_ctg = ae_ctg
1006 
1007  allocate( d_min(n_ctg) )
1008  allocate( d_max(n_ctg) )
1009  allocate( n_kap(n_ctg) )
1010  allocate( k_min(n_ctg) )
1011  allocate( k_max(n_ctg) )
1012 
1013 
1014  d_min(1:n_ctg) = d_min_inp(1:n_ctg) ! lower bound of 1st size bin
1015  d_max(1:n_ctg) = d_max_inp(1:n_ctg) ! upper bound of last size bin
1016  n_kap(1:n_ctg) = n_kap_inp(1:n_ctg) ! number of kappa bins
1017  k_min(1:n_ctg) = k_min_inp(1:n_ctg) ! lower bound of 1st kappa bin
1018  k_max(1:n_ctg) = k_max_inp(1:n_ctg) ! upper bound of last kappa bin
1019 
1020 
1021  !bin setting
1022  allocate(d_lw(n_siz_max,n_ctg))
1023  allocate(d_ct(n_siz_max,n_ctg))
1024  allocate(d_up(n_siz_max,n_ctg))
1025  allocate(k_lw(n_kap_max,n_ctg))
1026  allocate(k_ct(n_kap_max,n_ctg))
1027  allocate(k_up(n_kap_max,n_ctg))
1028  d_lw(:,:) = 0._rp
1029  d_ct(:,:) = 0._rp
1030  d_up(:,:) = 0._rp
1031  k_lw(:,:) = 0._rp
1032  k_ct(:,:) = 0._rp
1033  k_up(:,:) = 0._rp
1034 
1035  do ic = 1, n_ctg
1036 
1037  dlogd = (log(d_max(ic)) - log(d_min(ic)))/real(nsiz(ic),kind=rp)
1038 
1039  do is0 = 1, nsiz(ic) !size bin
1040  d_lw(is0,ic) = exp(log(d_min(ic))+dlogd* real(is0-1,kind=rp) )
1041  d_ct(is0,ic) = exp(log(d_min(ic))+dlogd*(real(is0 ,kind=rp)-0.5_rp))
1042  d_up(is0,ic) = exp(log(d_min(ic))+dlogd* real(is0 ,kind=rp) )
1043  enddo !is (1:n_siz(ic))
1044  dk = (k_max(ic) - k_min(ic))/real(n_kap(ic),kind=rp)
1045  do ik = 1, n_kap(ic) !size bin
1046  k_lw(ik,ic) = k_min(ic) + dk * real(ik-1,kind=rp)
1047  k_ct(ik,ic) = k_min(ic) + dk *(real(ik ,kind=rp)-0.5_rp)
1048  k_up(ik,ic) = k_min(ic) + dk * real(ik ,kind=rp)
1049  enddo !ik (1:n_kap(ic))
1050 
1051  enddo !ic (1:n_ctg)
1052 ! ik = 1 !only one kappa bin
1053 
1054  m0t = m0_init !total M0 [#/m3]
1055  dgt = dg_init ![m]
1056  sgt = sg_init ![-]
1057 
1058  if ( m0t <= cleannumber ) then
1059  m0t = cleannumber
1060  dgt = 0.1e-6_rp
1061  sgt = 1.3_rp
1062  log_warn("ATMOS_PHY_AE_kajino13_mkinit",*) 'Initial aerosol number is set as ', cleannumber, '[#/m3]'
1063  endif
1064 
1065  m2t = m0t*dgt**(2.d0) *dexp(2.0d0 *(dlog(real(sgt,kind=dp))**2.d0)) !total M2 [m2/m3]
1066  m3t = m0t*dgt**(3.d0) *dexp(4.5d0 *(dlog(real(sgt,kind=dp))**2.d0)) !total M3 [m3/m3]
1067  mst = m3t*pi6*conv_vl_ms !total Ms [ug/m3]
1068 
1069  do ic = 1, n_ctg
1070  !aerosol_procs initial condition
1071  do ik = 1, n_kap(ic) !kappa bin
1072  do is0 = 1, nsiz(ic)
1073  if (dgt >= d_lw(is0,ic) .and. dgt < d_up(is0,ic)) then
1074  aerosol_procs(ia_m0,is0,ik,ic) = aerosol_procs(ia_m0,is0,ik,ic) + m0t ![#/m3]
1075  aerosol_procs(ia_m2,is0,ik,ic) = aerosol_procs(ia_m2,is0,ik,ic) + m2t ![m2/m3]
1076  aerosol_procs(ia_m3,is0,ik,ic) = aerosol_procs(ia_m3,is0,ik,ic) + m3t ![m3/m3]
1077  aerosol_procs(ia_ms,is0,ik,ic) = aerosol_procs(ia_ms,is0,ik,ic) + mst*1.e-9_rp ![kg/m3]
1078  elseif (dgt < d_lw(1,ic)) then
1079  aerosol_procs(ia_m0,1 ,ik,ic) = aerosol_procs(ia_m0,1 ,ik,ic) + m0t ![#/m3]
1080  aerosol_procs(ia_m2,1 ,ik,ic) = aerosol_procs(ia_m2,1 ,ik,ic) + m2t ![m2/m3]
1081  aerosol_procs(ia_m3,1 ,ik,ic) = aerosol_procs(ia_m3,1 ,ik,ic) + m3t ![m3/m3]
1082  aerosol_procs(ia_ms,1 ,ik,ic) = aerosol_procs(ia_ms,1 ,ik,ic) + mst*1.e-9_rp ![kg/m3]
1083  elseif (dgt >= d_up(nsiz(ic),ic)) then
1084  aerosol_procs(ia_m0,nsiz(ic),ik,ic) = aerosol_procs(ia_m0,nsiz(ic),ik,ic) + m0t ![#/m3]
1085  aerosol_procs(ia_m2,nsiz(ic),ik,ic) = aerosol_procs(ia_m2,nsiz(ic),ik,ic) + m2t ![m2/m3]
1086  aerosol_procs(ia_m3,nsiz(ic),ik,ic) = aerosol_procs(ia_m3,nsiz(ic),ik,ic) + m3t ![m3/m3]
1087  aerosol_procs(ia_ms,nsiz(ic),ik,ic) = aerosol_procs(ia_ms,nsiz(ic),ik,ic) + mst*1.e-9_rp ![kg/m3]
1088  endif
1089  enddo
1090  enddo
1091  enddo
1092 
1093  conc_gas(:) = 0.0_rp
1094  do j = 1, ja
1095  do i = 1, ia
1096  do k = 1, ka
1097  it = 1
1098  do ic = 1, n_ctg ! category
1099  do ik = 1, nkap(ic) ! kappa bin
1100  do is0 = 1, nsiz(ic) ! size bin
1101  do ia0 = 1, n_atr ! attributes
1102  qtrc(k,i,j,it) = aerosol_procs(ia0,is0,ik,ic) / dens(k,i,j) !#,m2,m3,kg/m3 -> #,m2,m3,kg/kg
1103  it = it + 1
1104  enddo !ia0 (1:N_ATR )
1105  enddo !is (1:n_siz(ic) )
1106  enddo !ik (1:n_kap(ic) )
1107  enddo !ic (1:n_ctg )
1108 
1109  do ic = 1, gas_ctg ! GAS category
1110  qtrc(k,i,j,it) = conc_gas(ic) / dens(k,i,j) * 1.e-9_rp !mixing ratio [kg/kg]
1111  it = it + 1
1112  enddo
1113  enddo
1114  enddo
1115  enddo
1116 
1117  ccn(:,:,:) = 0.0_rp
1118  do j = js, je
1119  do i = is, ie
1120  do k = ks, ke
1121  !--- calculate super saturation of water
1122  call saturation_pres2qsat_liq( temp(k,i,j), pres(k,i,j), qdry(k,i,j), & ! [IN]
1123  qsat_tmp ) ! [OUT]
1124  ssliq = qv(k,i,j) / qsat_tmp - 1.0_rp
1125 
1126  call aerosol_activation( c_kappa, ssliq, temp(k,i,j), ia_m0, ia_m2, ia_m3, &
1127  n_atr, n_siz_max, n_kap_max, n_ctg, nsiz, n_kap, &
1128  d_ct, aerosol_procs, aerosol_activ )
1129 
1130  do is0 = 1, nsiz(ic_mix)
1131  ccn(k,i,j) = ccn(k,i,j) + aerosol_activ(ia_m0,is0,ik_out,ic_mix)
1132  enddo
1133  enddo
1134  enddo
1135  enddo
1136 
1137  return
1138  end subroutine atmos_phy_ae_kajino13_mkinit
1139 
1140  subroutine atmos_phy_ae_kajino13_negative_fixer( &
1141  KA, KS, KE, IA, IS, IE, JA, JS, JE, QA_AE, &
1142  QTRC )
1143  integer, intent(in) :: ka, ks, ke
1144  integer, intent(in) :: ia, is, ie
1145  integer, intent(in) :: ja, js, je
1146  integer, intent(in) :: qa_ae
1147 
1148  real(rp), intent(inout) :: qtrc(ka,ia,ja,qa_ae)
1149 
1150  integer :: k ,i, j
1151  integer :: ic, ik, is0
1152 
1153  !--- Negative fixer
1154  do j = js, je
1155  do i = is, ie
1156  do k = ks, ke
1157  do ic = 1, n_ctg !aerosol category
1158  do ik = 1, n_kap(ic) !kappa bin
1159  do is0 = 1, n_siz(ic) !size bin
1160  if (qtrc(k,i,j,it_procs2trans(ia_m0,is0,ik,ic)) < 0.0_rp .or. &
1161  qtrc(k,i,j,it_procs2trans(ia_m2,is0,ik,ic)) < 0.0_rp .or. &
1162  qtrc(k,i,j,it_procs2trans(ia_m3,is0,ik,ic)) < 0.0_rp .or. &
1163  qtrc(k,i,j,it_procs2trans(ia_ms,is0,ik,ic)) < 0.0_rp .or. &
1164  qtrc(k,i,j,it_procs2trans(ia_kp,is0,ik,ic)) < 0.0_rp ) then
1165  qtrc(k,i,j,it_procs2trans(1:n_atr,is0,ik,ic)) = 0.0_rp
1166  endif
1167  enddo
1168  enddo
1169  enddo
1170  enddo
1171  enddo
1172  enddo
1173 
1174  return
1175  end subroutine atmos_phy_ae_kajino13_negative_fixer
1176 
1177  ! private subroutines
1178 
1179  !-----------------------------------------------------------------------------
1180  ! subroutine 4. aerosol_zerochem
1181  !-----------------------------------------------------------------------------
1182  subroutine aerosol_zerochem ( &
1183  deltt, & !--- in
1184  temp_k, pres_pa, super, & !--- in
1185 ! h2so4dt, c_ratio, c_kappa, & !--- in
1186  flag_npf, flag_cond, flag_coag,& !--- in
1187  aerosol_procs, & !--- inout
1188  conc_gas, & !--- inout
1189  emis_procs, & !--- out
1190  emis_gas, & !--- out
1191  aerosol_activ )
1192  implicit none
1193  ! i/o variables
1194  real(dp),intent(in) :: deltt ! delta t [sec]
1195  real(rp),intent(in) :: temp_k ! temperature [K]
1196  real(rp),intent(in) :: pres_pa ! pressure [Pa]
1197  real(rp),intent(in) :: super ! supersaturation [-]
1198 ! real(RP),intent(in) :: h2so4dt ! h2so4 production rate [ug/m3/s]
1199 ! real(RP),intent(in) :: c_ratio ! ratio of condensable mass to h2so4 (after NPF) [-]
1200 ! real(RP),intent(in) :: c_kappa ! hygroscopicity of condensable mass [-]
1201  logical, intent(in) :: flag_npf ! (on/off) new particle formation
1202  logical, intent(in) :: flag_cond ! (on/off) condensation
1203  logical, intent(in) :: flag_coag ! (on/off) coagulation
1204  real(rp),intent(inout) :: aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
1205  real(rp),intent(inout) :: conc_gas(gas_ctg)
1206  real(rp),intent(out) :: aerosol_activ(n_atr,n_siz_max,n_kap_max,n_ctg)
1207  real(rp),intent(in) :: emis_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
1208  real(rp),intent(in) :: emis_gas(gas_ctg)
1209  ! local variables
1210  real(rp) :: j_1nm ! nucleation rate of 1nm particles [#/cm3/s]
1211  integer :: ic_nuc ! category for 1nm new particles
1212  integer :: ik_nuc ! kappa bin for 1nm new particles
1213  integer :: is_nuc ! size bin for 1nm new particles
1214  real(rp) :: c_ratio
1215  real(rp) :: chem_gas(gas_ctg)
1216 
1217  chem_gas(ig_h2so4) = atmos_phy_ae_kajino13_h2so4dt
1218  chem_gas(ig_cgas) = atmos_phy_ae_kajino13_ocgasdt
1219 
1220  !--- convert unit of aerosol mass [ia=ia_ms]
1221  do ic = 1, n_ctg !aerosol category
1222  do ik = 1, n_kap(ic) !kappa bin
1223  do is0 = 1, n_siz(ic) !size bin
1224  aerosol_procs(ia_ms,is0,ik,ic) = aerosol_procs(ia_ms,is0,ik,ic) * 1.e+9_rp ! [kg/m3] -> [ug/m3]
1225  enddo !is (1:n_siz(ic) )
1226  enddo !ik (1:n_kap(ic) )
1227  enddo !ic (1:n_ctg )
1228 
1229  ! emission
1230  do ic = 1, n_ctg !aerosol category
1231  do ik = 1, n_kap(ic) !kappa bin
1232  do is0 = 1, n_siz(ic) !size bin
1233  do ia0 = 1, n_atr !attributes (prognostic)
1234  aerosol_procs(ia0,is0,ik,ic) = aerosol_procs(ia0,is0,ik,ic) &
1235  + emis_procs(ia0,is0,ik,ic) * deltt
1236  enddo !ia0 (1:n_atr )
1237  enddo !is (1:n_siz(ic) )
1238  enddo !ik (1:n_kap(ic) )
1239  enddo !ic (1:n_ctg )
1240 
1241  ! tiny number, tiny mass
1242  do ic = 1, n_ctg !aerosol category
1243  do ik = 1, n_kap(ic) !kappa bin
1244  do is0 = 1, n_siz(ic) !size bin
1245  if (aerosol_procs(ia_m0,is0,ik,ic) < cleannumber) then
1246  do ia0 = 1, n_atr
1247  aerosol_procs(ia0,is0,ik,ic) = 0._rp !to save cpu time and avoid underflow
1248  enddo !ia0 (1:n_atr )
1249  endif
1250  enddo !is (1:n_siz(ic) )
1251  enddo !ik (1:n_kap(ic) )
1252  enddo !ic (1:n_ctg )
1253 
1254  ! update conc_h2so4
1255 ! conc_h2so4 = conc_h2so4 + h2so4dt * deltt ! [ug/m3]
1256  conc_gas(ig_h2so4) = conc_gas(ig_h2so4) + chem_gas(ig_h2so4) * deltt ! [ug/m3]
1257  conc_gas(ig_cgas) = conc_gas(ig_cgas) + chem_gas(ig_cgas) * deltt ! [ug/m3]
1258 
1259  conc_gas(ig_h2so4) = conc_gas(ig_h2so4) + emis_gas(ig_h2so4) * deltt ! [ug/m3]
1260  conc_gas(ig_cgas) = conc_gas(ig_cgas) + emis_gas(ig_cgas) * deltt ! [ug/m3]
1261 
1262  if( chem_gas(ig_h2so4)+emis_gas(ig_h2so4) /= 0.0_rp ) then
1263  c_ratio = ( chem_gas(ig_h2so4)+emis_gas(ig_h2so4)+chem_gas(ig_cgas)+emis_gas(ig_cgas) ) &
1264  / ( chem_gas(ig_h2so4)+emis_gas(ig_h2so4) )
1265  else
1266  c_ratio = 0.0_rp
1267  endif
1268 
1269  ! new particle formation
1270  ic_nuc = ic_mix
1271  ik_nuc = 1
1272  is_nuc = 1
1273  if (flag_npf) then
1274  call prof_rapstart('ATM_Aerosol_NPF',1)
1275  call aerosol_nucleation(conc_gas(ig_h2so4), j_1nm) !(i) conc_h2so4 (o) J_1nm
1276  call prof_rapend('ATM_Aerosol_NPF',1)
1277  else
1278  j_1nm = 0._rp !( new particle formation does not occur )
1279  endif !if flag_npf=.true.
1280 
1281  ! condensation
1282  if (flag_cond) then
1283  call prof_rapstart('ATM_Aerosol_cond',1)
1284  call aerosol_condensation(j_1nm, temp_k, pres_pa, deltt, &
1285  ic_nuc, ik_nuc, is_nuc, n_ctg, n_kap, n_siz, n_atr, &
1286  n_siz_max, n_kap_max, ia_m0, ia_m2, ia_m3, ia_ms, &
1287  ! d_lw, d_ct, d_up, k_lw, k_ct, k_up, &
1288  d_lw, d_ct, d_up, &
1289  conc_gas(ig_h2so4), c_ratio, aerosol_procs )
1290  call prof_rapend('ATM_Aerosol_cond',1)
1291  endif !if flag_cond=.true.
1292 
1293  ! coagulation
1294  if (flag_coag) then
1295  call prof_rapstart('ATM_Aerosol_coag',1)
1296  call aerosol_coagulation(deltt, temp_k, pres_pa, &
1297  mcomb,is_i,is_j,is_k,ik_i,ik_j,ik_k,ic_i,ic_j,ic_k, &
1298  n_atr,n_siz_max,n_kap_max,n_ctg,ia_m0,ia_m2,ia_m3,ia_ms, &
1299  n_siz,n_kap,d_lw,d_ct,d_up,aerosol_procs)
1300  call prof_rapend('ATM_Aerosol_coag',1)
1301  endif !if flag_coag=.true.
1302 
1303  call aerosol_activation(atmos_phy_ae_kajino13_c_kappa, &
1304  super, temp_k, ia_m0, ia_m2, ia_m3, &
1305  n_atr,n_siz_max,n_kap_max,n_ctg,n_siz,n_kap, &
1306  d_ct,aerosol_procs, aerosol_activ)
1307 
1308 ! conc_gas(IG_CGAS) = conc_gas(IG_H2SO4)*( c_ratio-1.0_RP )
1309 
1310  !--- convert unit of aerosol mass [ia=ia_ms]
1311  do ic = 1, n_ctg !aerosol category
1312  do ik = 1, n_kap(ic) !kappa bin
1313  do is0 = 1, n_siz(ic) !size bin
1314  aerosol_procs(ia_ms,is0,ik,ic) = aerosol_procs(ia_ms,is0,ik,ic) * 1.e-9_rp ! [ug/m3] -> [kg/m3]
1315  enddo !is (1:n_siz(ic) )
1316  enddo !ik (1:n_kap(ic) )
1317  enddo !ic (1:n_ctg )
1318 
1319  return
1320  end subroutine aerosol_zerochem
1321  !-----------------------------------------------------------------------------
1322  ! subroutine 1. aerosol_nucleation
1323  !----------------------------------------------------------------------------------------
1324  subroutine aerosol_nucleation(conc_h2so4, J_1nm) ! (i) conc_h2so4 (i/o) J_1nm
1325  implicit none
1326  !i/o variables
1327  real(rp), intent(in) :: conc_h2so4 !H2SO4 concentration [ug/m3]
1328  real(rp), intent(inout) :: j_1nm !nucleation rate of 1nm particles [#/cm3/s]
1329  !local variables
1330  real(rp) :: conc_num_h2so4 !H2SO4 number concentration [#/cm3]
1331  ! real(RP), parameter :: logk = -12.4_RP !constant coefficient for kinetic nucleation [-]
1332  ! real(RP), parameter :: logk = -11.4_RP !constant coefficient for kinetic nucleation [-]
1333 
1334  conc_num_h2so4 = conc_h2so4 * conv_m2n ![ug/m3] -> [#/cm3]
1335 
1336  !emperical formula of new particle formation rate [Kuang et al., 2008]
1337  j_1nm = (10._rp**(atmos_phy_ae_kajino13_logk_aenucl)) * conc_num_h2so4 ** 2._rp
1338 
1339  return
1340  end subroutine aerosol_nucleation
1341  !----------------------------------------------------------------------------------------
1342  ! subroutine 2. aerosol_condensation
1343  !----------------------------------------------------------------------------------------
1344  subroutine aerosol_condensation(J_1nm, temp_k, pres_pa, deltt, &
1345  ic_nuc, ik_nuc, is_nuc, n_ctg, n_kap, n_siz, n_atr, &
1346  n_siz_max, n_kap_max, ia_m0, ia_m2, ia_m3, ia_ms, &
1347  ! d_lw, d_ct, d_up, k_lw, k_ct, k_up, &
1348  d_lw, d_ct, d_up, &
1349  conc_h2so4, c_ratio, aerosol_procs )
1350 
1351  implicit none
1352  !i/o variables
1353  real(rp),intent(in) :: j_1nm ! nucleation rate of 1nm particles [#/cm3/s]
1354  real(rp),intent(in) :: temp_k ! temperature [K]
1355  real(rp),intent(in) :: pres_pa ! pressure [Pa]
1356  real(dp),intent(in) :: deltt ! delta t [sec]
1357  integer, intent(in) :: ic_nuc ! category for 1nm new particles
1358  integer, intent(in) :: ik_nuc ! kappa bin for 1nm new particles
1359  integer, intent(in) :: is_nuc ! size bin for 1nm new particles
1360  integer, intent(in) :: n_ctg ! number of aerosol category
1361  integer, intent(in) :: n_kap(n_ctg) ! number of kappa bins
1362  integer, intent(in) :: n_siz(n_ctg) ! number of size bins
1363  integer, intent(in) :: n_atr ! number of aerosol category
1364  integer, intent(in) :: n_siz_max ! max of n_siz
1365  integer, intent(in) :: n_kap_max ! max of n_kap
1366  integer, intent(in) :: ia_m0 !
1367  integer, intent(in) :: ia_m2 !
1368  integer, intent(in) :: ia_m3 !
1369  integer, intent(in) :: ia_ms !
1370  real(rp),intent(in) :: c_ratio ! ratio of condensable mass to h2so4 (after NPF) [-]
1371  real(rp),intent(in) :: d_lw(n_siz_max,n_ctg), d_ct(n_siz_max,n_ctg), d_up(n_siz_max,n_ctg)
1372  ! real(RP), dimension(n_kap_max,n_ctg), intent(in) :: k_lw, k_ct, k_up
1373  real(rp),intent(inout) :: conc_h2so4 !concentration [ug/m3]
1374  real(rp),intent(inout) :: aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
1375  !local variables
1376  real(rp), parameter :: alpha = 0.1_rp ! accomodation coefficient
1377  real(rp), parameter :: cour = 0.5_rp ! courant number for condensation
1378  integer :: isplt ! number of time splitted
1379  real(rp) :: dtsplt ! splitted time step
1380  real(rp) :: sq_cbar_h2so4 ! square of molecular speed of H2SO4 gas [m2/s2]
1381  real(rp) :: cbar_h2so4 ! molecular speed of H2SO4 gas [m/s]
1382  real(rp) :: dv_h2so4 ! diffusivity of H2SO4 gas [m2/s]
1383  real(rp) :: drive ! driving force of H2SO4 gas [m3SO4/m3]
1384  real(rp) :: dm0dt_npf, dm2dt_npf, dm3dt_npf, dmsdt_npf ! dMk/dt for new particle formation
1385  real(rp) :: dm0dt_cnd(n_siz_max,n_kap_max,n_ctg) ! dM0/dt for condensation/evaporation
1386  real(rp) :: dm2dt_cnd(n_siz_max,n_kap_max,n_ctg) ! dM2/dt for condensation/evaporation
1387  real(rp) :: dm3dt_cnd(n_siz_max,n_kap_max,n_ctg) ! dM3/dt for condensation/evaporation
1388  real(rp) :: dmsdt_cnd(n_siz_max,n_kap_max,n_ctg) ! dMs/dt for condensation/evaporation
1389  real(rp) :: rm0_pls(n_siz_max), rm0_mns(n_siz_max) ! used for moving center calc.
1390  real(rp) :: rm2_pls(n_siz_max), rm2_mns(n_siz_max) ! used for moving center calc.
1391  real(rp) :: rm3_pls(n_siz_max), rm3_mns(n_siz_max) ! used for moving center calc.
1392  real(rp) :: rms_pls(n_siz_max), rms_mns(n_siz_max) ! used for moving center calc.
1393  real(rp) :: m0t,m1t,m2t,m3t,dgt,sgt,dm2
1394  real(rp) :: gnc2,gnc3,gfm2,gfm3,harm2,harm3
1395  real(rp) :: lossrate,tmps6
1396  integer :: ic, ik, is0, i, is1, is2
1397 
1398 ! nothing happens --> in future, c_ratio is removed and condensable gas should be used
1399  if (conc_h2so4 <= 0.0_rp) return
1400 ! nothing happens --> in future, c_ratio is removed and condensable gas should be used
1401 
1402  drive = conc_h2so4 * mwrat_s6 * conv_ms_vl ![ugH2SO4/m3]=>volume [m3SO4/m3]
1403  sq_cbar_h2so4 = 8.0_rp*rgas*temp_k/(pi*mwh2so4*1.e-3_rp)
1404  cbar_h2so4 = sqrt( sq_cbar_h2so4 )
1405  dv_h2so4 = diffsulf*(stdatmpa/pres_pa)*(temp_k/stdtemp)**1.75_rp
1406 
1407  dm0dt_npf = 0.0_rp
1408  dm2dt_npf = 0.0_rp
1409  dm3dt_npf = 0.0_rp
1410  dmsdt_npf = 0.0_rp
1411  dm0dt_cnd = 0.0_rp
1412  dm2dt_cnd = 0.0_rp
1413  dm3dt_cnd = 0.0_rp
1414  dmsdt_cnd = 0.0_rp
1415 
1416  !(npf rate)
1417  dm0dt_npf = j_1nm * 1.e6_rp ! [#/m3/s]
1418  dm2dt_npf = dm0dt_npf * 1.e-18_rp ! [m2/m3/s]
1419  dm3dt_npf = dm0dt_npf * 1.e-27_rp ! [m3/m3/s]
1420  dmsdt_npf = dm3dt_npf * pi6 * conv_vl_ms ! [ug/m3/s]
1421 
1422  !(condensation rate)
1423  do ic = 1, n_ctg !aerosol category
1424  do ik = 1, n_kap(ic) !kappa bin
1425  do is0 = 1, n_siz(ic) !size bin
1426 
1427  dm2dt_cnd(is0,ik,ic) = 0._rp
1428  dm3dt_cnd(is0,ik,ic) = 0._rp
1429 
1430  if (aerosol_procs(ia_m0,is0,ik,ic) &
1431  *aerosol_procs(ia_m2,is0,ik,ic) &
1432  *aerosol_procs(ia_m3,is0,ik,ic) > 0._rp) then
1433  m0t = aerosol_procs(ia_m0,is0,ik,ic)
1434  m2t = aerosol_procs(ia_m2,is0,ik,ic)
1435  m3t = aerosol_procs(ia_m3,is0,ik,ic)
1436  call diag_ds(m0t,m2t,m3t,dgt,sgt,dm2)
1437  if (dgt <= 0._rp) dgt = d_ct(is0,ic)
1438  if (sgt <= 0._rp) sgt = 1.3_rp
1439  m1t = m0t*dgt*exp(0.5_rp*(log(sgt)**2._rp)) ![m/m3]
1440  !=4kDvM_k-2
1441  gnc2 = 8._rp * dv_h2so4 * m0t !m2/s*#/m3 = m2/m3/s
1442  gnc3 = 12._rp * dv_h2so4 * m1t !m2/s*m/m3 = m3/m3/s
1443  != k/2acM_k-1
1444  gfm2 = alpha * cbar_h2so4 * m1t !m/s *m/m3 = m2/m3/s
1445  gfm3 = 1.5_rp * alpha * cbar_h2so4 * m2t !m/s *m2/m3 = m3/m3/s
1446  != harmonic mean approach
1447  harm2 = gnc2 * gfm2 / ( gnc2 + gfm2 ) !m2/m3/s
1448  harm3 = gnc3 * gfm3 / ( gnc3 + gfm3 ) !m3/m3/s
1449 
1450  dm2dt_cnd(is0,ik,ic) = harm2*drive !m2/m3/s
1451  dm3dt_cnd(is0,ik,ic) = harm3*drive !m3/m3/s
1452  dmsdt_cnd(is0,ik,ic) = dm3dt_cnd(is0,ik,ic)*pi6*conv_vl_ms !ug/m3/s
1453 
1454  endif !aerosol number > 0
1455 
1456  enddo
1457  enddo
1458  enddo
1459 
1460  ! coagulation
1461  !=time split
1462  lossrate = dmsdt_npf
1463 
1464  do ic = 1, n_ctg !aerosol category
1465  do ik = 1, n_kap(ic) !kappa bin
1466  do is0 = 1, n_siz(ic) !size bin
1467  lossrate = lossrate + dmsdt_cnd(is0,ik,ic)
1468  enddo
1469  enddo
1470  enddo
1471 
1472  if (lossrate<=0._rp) return !nothing happens
1473 
1474  isplt = int( (lossrate*deltt) / (cour*conc_h2so4) )+1
1475  dtsplt = deltt/float(isplt)
1476 
1477  !== NPF & condensation calculation
1478  loop_split: do i = 1, isplt
1479  tmps6 = conc_h2so4
1480  conc_h2so4 = conc_h2so4 - (lossrate * dtsplt) * mwrat_s6_i !ugH2SO4/m3
1481 
1482  if (conc_h2so4 < 0._rp) then
1483  dtsplt = tmps6 * mwrat_s6 / lossrate
1484  conc_h2so4 = 0._rp
1485  endif !conc_h2so4 < 0
1486 
1487  !(npf)
1488  aerosol_procs(ia_m0,is_nuc,ik_nuc,ic_nuc) = & ! #/m3
1489  aerosol_procs(ia_m0,is_nuc,ik_nuc,ic_nuc) + dm0dt_npf * dtsplt
1490  aerosol_procs(ia_m2,is_nuc,ik_nuc,ic_nuc) = & ! m2/m3
1491  aerosol_procs(ia_m2,is_nuc,ik_nuc,ic_nuc) + dm2dt_npf * dtsplt
1492  aerosol_procs(ia_m3,is_nuc,ik_nuc,ic_nuc) = & ! m3/m3
1493  aerosol_procs(ia_m3,is_nuc,ik_nuc,ic_nuc) + dm3dt_npf * dtsplt
1494  aerosol_procs(ia_ms,is_nuc,ik_nuc,ic_nuc) = & ! ug/m3
1495  aerosol_procs(ia_ms,is_nuc,ik_nuc,ic_nuc) + dmsdt_npf * dtsplt
1496 
1497  !(condensation)
1498  do ic = 1, n_ctg !aerosol category
1499  do ik = 1, n_kap(ic) !kappa bin
1500  do is0 = 1, n_siz(ic) !size bin
1501 
1502  aerosol_procs(ia_m2,is0,ik,ic) = & ! m2/m3
1503  aerosol_procs(ia_m2,is0,ik,ic) + dm2dt_cnd(is0,ik,ic) * dtsplt & !m2/m3
1504  * c_ratio !additional mass for condensation
1505  aerosol_procs(ia_m3,is0,ik,ic) = & ! m3/m3
1506  aerosol_procs(ia_m3,is0,ik,ic) + dm3dt_cnd(is0,ik,ic) * dtsplt & !m3/m3
1507  * c_ratio !additional mass for condensation
1508  aerosol_procs(ia_ms,is0,ik,ic) = & ! ug/m3
1509  aerosol_procs(ia_ms,is0,ik,ic) + dmsdt_cnd(is0,ik,ic) * dtsplt & !ug/m3
1510  * c_ratio !additional mass for condensation
1511 
1512  enddo
1513  enddo
1514  enddo
1515 
1516 ! if (conc_h2so4 <= 0._RP) goto 777
1517  if (conc_h2so4 <= 0.0_rp) exit loop_split
1518 
1519  enddo loop_split
1520 
1521 ! 777 continue
1522 
1523  conc_h2so4 = max(conc_h2so4, 0._rp)
1524 
1525  !== moving center
1526  do ic = 1, n_ctg !aerosol category
1527  do ik = 1, n_kap(ic) !kappa bin
1528 
1529  rm0_pls = 0._rp
1530  rm2_pls = 0._rp
1531  rm3_pls = 0._rp
1532  rms_pls = 0._rp
1533  rm0_mns = 0._rp
1534  rm2_mns = 0._rp
1535  rm3_mns = 0._rp
1536  rms_mns = 0._rp
1537 
1538  do is1 = 1, n_siz(ic) !size bin
1539 
1540  if (aerosol_procs(ia_m0,is1,ik,ic) &
1541  *aerosol_procs(ia_m2,is1,ik,ic) &
1542  *aerosol_procs(ia_m3,is1,ik,ic) > 0._rp) then
1543  m0t = aerosol_procs(ia_m0,is1,ik,ic)
1544  m2t = aerosol_procs(ia_m2,is1,ik,ic)
1545  m3t = aerosol_procs(ia_m3,is1,ik,ic)
1546  call diag_ds(m0t,m2t,m3t,dgt,sgt,dm2)
1547  if (dgt <= 0._rp) dgt = d_ct(is1,ic)
1548  if (sgt <= 0._rp) sgt = 1.3_rp
1549 
1550  do is2 = 1, n_siz(ic)
1551  if (dgt >= d_lw(is2,ic) .AND. dgt < d_up(is2,ic)) then !moving center
1552  rm0_pls(is2) = rm0_pls(is2) + aerosol_procs(ia_m0,is1,ik,ic) !is2 <=
1553  rm0_mns(is1) = aerosol_procs(ia_m0,is1,ik,ic) !is1 =>
1554  rm2_pls(is2) = rm2_pls(is2) + aerosol_procs(ia_m2,is1,ik,ic) !is2 <=
1555  rm2_mns(is1) = aerosol_procs(ia_m2,is1,ik,ic) !is1 =>
1556  rm3_pls(is2) = rm3_pls(is2) + aerosol_procs(ia_m3,is1,ik,ic) !is2 <=
1557  rm3_mns(is1) = aerosol_procs(ia_m3,is1,ik,ic) !is1 =>
1558  rms_pls(is2) = rms_pls(is2) + aerosol_procs(ia_ms,is1,ik,ic) !is2 <=
1559  rms_mns(is1) = aerosol_procs(ia_ms,is1,ik,ic) !is1 =>
1560  exit
1561  endif !d_lw(is2) < dgt < d_up(is2)
1562  enddo
1563 
1564  endif !aerosol number > 0
1565 
1566  enddo
1567 
1568  do is0 = 1, n_siz(ic)
1569  aerosol_procs(ia_m0,is0,ik,ic) = aerosol_procs(ia_m0,is0,ik,ic) &
1570  + rm0_pls(is0) - rm0_mns(is0)
1571  aerosol_procs(ia_m2,is0,ik,ic) = aerosol_procs(ia_m2,is0,ik,ic) &
1572  + rm2_pls(is0) - rm2_mns(is0)
1573  aerosol_procs(ia_m3,is0,ik,ic) = aerosol_procs(ia_m3,is0,ik,ic) &
1574  + rm3_pls(is0) - rm3_mns(is0)
1575  aerosol_procs(ia_ms,is0,ik,ic) = aerosol_procs(ia_ms,is0,ik,ic) &
1576  + rms_pls(is0) - rms_mns(is0)
1577  enddo !is(1:n_siz(ic))
1578 
1579  enddo
1580  enddo
1581 
1582  return
1583  end subroutine aerosol_condensation
1584  !----------------------------------------------------------------------------------------
1585  ! subroutine 3. aerosol_coagulation
1586  !----------------------------------------------------------------------------------------
1587  subroutine aerosol_coagulation(deltt, temp_k, pres_pa, &
1588  mcomb,is_i,is_j,is_k,ik_i,ik_j,ik_k,ic_i,ic_j,ic_k, &
1589  n_atr,n_siz_max,n_kap_max,n_ctg,ia_m0,ia_m2,ia_m3,ia_ms, &
1590  n_siz,n_kap,d_lw,d_ct,d_up,aerosol_procs)
1591  implicit none
1592  !i/o variables
1593  real(dp),intent(in) :: deltt ! delta t [sec]
1594  real(rp),intent(in) :: temp_k ! temperature [K]
1595  real(rp),intent(in) :: pres_pa ! pressure [Pa]
1596  integer, intent(in) :: mcomb !combinations of sections
1597  integer, intent(in) :: is_i(mcomb), is_j(mcomb), is_k(mcomb)
1598  integer, intent(in) :: ik_i(mcomb), ik_j(mcomb), ik_k(mcomb)
1599  integer, intent(in) :: ic_i(mcomb), ic_j(mcomb), ic_k(mcomb)
1600  integer, intent(in) :: n_ctg ! number of aerosol category
1601  integer, intent(in) :: n_atr ! number of aerosol category
1602  integer, intent(in) :: n_siz_max ! max of n_siz
1603  integer, intent(in) :: n_kap_max ! max of n_kap
1604  integer, intent(in) :: ia_m0 !
1605  integer, intent(in) :: ia_m2 !
1606  integer, intent(in) :: ia_m3 !
1607  integer, intent(in) :: ia_ms !
1608  integer, intent(in) :: n_siz(n_ctg) !
1609  integer, intent(in) :: n_kap(n_ctg) !
1610  real(rp),intent(in):: d_lw(n_siz_max,n_ctg), d_ct(n_siz_max,n_ctg), d_up(n_siz_max,n_ctg)
1611  real(rp),intent(inout) :: aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
1612  !local variables
1613  real(rp) :: dt_m0(n_siz_max,n_kap_max,n_ctg)
1614  real(rp) :: dt_m3(n_siz_max,n_kap_max,n_ctg)
1615  real(rp) :: dt_m6(n_siz_max,n_kap_max,n_ctg)
1616  real(rp) :: dt_ms(n_siz_max,n_kap_max,n_ctg)
1617  real(rp) :: sixth(n_siz_max,n_kap_max,n_ctg)
1618  real(rp) :: rm0_pls(n_siz_max), rm0_mns(n_siz_max) ! used for moving center calc.
1619  real(rp) :: rm2_pls(n_siz_max), rm2_mns(n_siz_max) ! used for moving center calc.
1620  real(rp) :: rm3_pls(n_siz_max), rm3_mns(n_siz_max) ! used for moving center calc.
1621  real(rp) :: rms_pls(n_siz_max), rms_mns(n_siz_max) ! used for moving center calc.
1622  real(rp) :: visair ! viscosity of air [Pa s]=[kg/m/s]
1623  real(rp) :: lambda ! mean free path of air molecules [cm]
1624  real(rp) :: rkfm,rknc
1625  real(rp) :: m0i,m0j,m2i,m2j,m3i,m3j,msi,msj,dgi,dgj,sgi,sgj,dm2,m6i,m6j
1626  real(rp) :: dm0i,dm0j,dm0k,dm3i,dm3j,dm3k,dm6i,dm6j,dm6k,dmsi,dmsj,dmsk
1627  real(rp) :: m0t,m3t,m6t,dgt,sgt,m2t
1628  integer :: mc, ic, ik, is0, is1, is2
1629 
1630  dt_m0(:,:,:) = 0._rp
1631  dt_m3(:,:,:) = 0._rp
1632  dt_m6(:,:,:) = 0._rp
1633  dt_ms(:,:,:) = 0._rp
1634  sixth(:,:,:) = 0._rp
1635 
1636  visair=c1*temp_k**c2/(1._rp+c3/temp_k) ! viscosity of air [Pa s]=[kg/m/s]
1637  rkfm=(3._rp*boltz*temp_k /rho_kg)**0.5_rp ![J/K*K*m3/kg]=[kg*m2/s2*m3/kg]=[m2.5/s]
1638  rknc= 2._rp*boltz*temp_k /(3._rp*visair) ![J*m*s/kg]=[kg*m2/s2*m*s/kg]=[m3/s]
1639  c4 = pi*boltz*temp_k/(8._rp*mwair/avo*1.e-3_rp) !
1640  lambda = visair/(0.499_rp*pres_pa)*c4**0.5_rp*1.e2_rp ! mean free path of air molecules [cm]
1641 
1642 
1643  !--- 555 coagulation combination rule loop
1644  do mc = 1, mcomb
1645  m0i = aerosol_procs(ia_m0,is_i(mc),ik_i(mc),ic_i(mc))
1646  m0j = aerosol_procs(ia_m0,is_j(mc),ik_j(mc),ic_j(mc))
1647  m2i = aerosol_procs(ia_m2,is_i(mc),ik_i(mc),ic_i(mc))
1648  m2j = aerosol_procs(ia_m2,is_j(mc),ik_j(mc),ic_j(mc))
1649  m3i = aerosol_procs(ia_m3,is_i(mc),ik_i(mc),ic_i(mc))
1650  m3j = aerosol_procs(ia_m3,is_j(mc),ik_j(mc),ic_j(mc))
1651  msi = aerosol_procs(ia_ms,is_i(mc),ik_i(mc),ic_i(mc))
1652  msj = aerosol_procs(ia_ms,is_j(mc),ik_j(mc),ic_j(mc))
1653 
1654  if (m0i*m2i*m3i <= 0._rp .OR. m0j*m2j*m3j <= 0._rp) cycle
1655 
1656  call diag_ds(m0i,m2i,m3i,dgi,sgi,dm2)
1657  if (dgi <= 0._rp) dgi = d_ct(is_i(mc),ic_i(mc))
1658  if (sgi <= 0._rp) sgi = 1.3_rp
1659  call diag_ds(m0j,m2j,m3j,dgj,sgj,dm2)
1660  if (dgj <= 0._rp) dgj = d_ct(is_j(mc),ic_j(mc))
1661  if (sgj <= 0._rp) sgj = 1.3_rp
1662 
1663  m6i = m0i*dgi**6._rp*exp(18._rp*(log(sgi)**2._rp))
1664  m6j = m0j*dgj**6._rp*exp(18._rp*(log(sgj)**2._rp))
1665  sixth(is_i(mc),ik_i(mc),ic_i(mc)) = m6i
1666  sixth(is_j(mc),ik_j(mc),ic_j(mc)) = m6j
1667 
1668  !intra sectional coagulation
1669  if ( is_i(mc) == is_j(mc) .AND. &
1670  ik_i(mc) == ik_j(mc) .AND. &
1671  ic_i(mc) == ic_j(mc) ) then
1672  call aero_intra(m0i, m3i, & !i
1673  dm0i, dm3i, dm6i, dmsi, & !o
1674  dgi, sgi, lambda, & !i
1675  rknc, rkfm, deltt ) !i
1676  dm0j = dm0i
1677  dm3j = dm3i
1678  dm6j = dm6i
1679  dmsj = dmsi
1680  dm0k = -dm0i
1681  dm3k = -dm3i
1682  dm6k = -dm6i
1683  dmsk = -dmsi
1684  !inter sectional coagulation
1685  else
1686  call aero_inter(m0i ,m3i ,m6i , & !input unchanged
1687  m0j ,m3j ,m6j , & !input unchanged
1688  dm0i,dm3i,dm6i,dmsi, & !output
1689  dm0j,dm3j,dm6j,dmsj, & !output
1690  dm0k,dm3k,dm6k,dmsk, & !output
1691  dgi ,sgi, rhod_ae, & !input unchanged
1692  dgj ,sgj, rhod_ae, & !input unchanged
1693  rknc, rkfm, lambda, & !input unchanged
1694  deltt ) !input unchanged
1695  endif
1696 
1697  dt_m0(is_i(mc),ik_i(mc),ic_i(mc)) = dt_m0(is_i(mc),ik_i(mc),ic_i(mc)) + dm0i
1698  dt_m0(is_j(mc),ik_j(mc),ic_j(mc)) = dt_m0(is_j(mc),ik_j(mc),ic_j(mc)) + dm0j
1699  dt_m0(is_k(mc),ik_k(mc),ic_k(mc)) = dt_m0(is_k(mc),ik_k(mc),ic_k(mc)) + dm0k
1700  dt_m3(is_i(mc),ik_i(mc),ic_i(mc)) = dt_m3(is_i(mc),ik_i(mc),ic_i(mc)) + dm3i
1701  dt_m3(is_j(mc),ik_j(mc),ic_j(mc)) = dt_m3(is_j(mc),ik_j(mc),ic_j(mc)) + dm3j
1702  dt_m3(is_k(mc),ik_k(mc),ic_k(mc)) = dt_m3(is_k(mc),ik_k(mc),ic_k(mc)) + dm3k
1703  dt_m6(is_i(mc),ik_i(mc),ic_i(mc)) = dt_m6(is_i(mc),ik_i(mc),ic_i(mc)) + dm6i
1704  dt_m6(is_j(mc),ik_j(mc),ic_j(mc)) = dt_m6(is_j(mc),ik_j(mc),ic_j(mc)) + dm6j
1705  dt_m6(is_k(mc),ik_k(mc),ic_k(mc)) = dt_m6(is_k(mc),ik_k(mc),ic_k(mc)) + dm6k
1706  dt_ms(is_i(mc),ik_i(mc),ic_i(mc)) = dt_ms(is_i(mc),ik_i(mc),ic_i(mc)) + dmsi
1707  dt_ms(is_j(mc),ik_j(mc),ic_j(mc)) = dt_ms(is_j(mc),ik_j(mc),ic_j(mc)) + dmsj
1708  dt_ms(is_k(mc),ik_k(mc),ic_k(mc)) = dt_ms(is_k(mc),ik_k(mc),ic_k(mc)) + dmsk
1709 
1710  enddo !555 continue !mc(1:mcomb)
1711 
1712 
1713  !--- redistribution
1714  do ic = 1, n_ctg !aerosol category
1715  do ik = 1, n_kap(ic) !kappa bin
1716  do is0 = 1, n_siz(ic) !size bin
1717  if (aerosol_procs(ia_m0,is0,ik,ic) > 0._rp .AND. & !to avoid divided by zero
1718  dt_m0(is0,ik,ic)/aerosol_procs(ia_m0,is0,ik,ic) < 1._rp ) then !to avoid overflow in exp
1719  aerosol_procs(ia_m0,is0,ik,ic)=aerosol_procs(ia_m0,is0,ik,ic) &
1720  *exp( dt_m0(is0,ik,ic)/aerosol_procs(ia_m0,is0,ik,ic) )
1721  else !aerosol_procs = 0
1722  aerosol_procs(ia_m0,is0,ik,ic)=aerosol_procs(ia_m0,is0,ik,ic) + dt_m0(is0,ik,ic)
1723  endif
1724  if (aerosol_procs(ia_m3,is0,ik,ic) > 0._rp .AND. & !to avoid divided by zero
1725  dt_m3(is0,ik,ic)/aerosol_procs(ia_m3,is0,ik,ic) < 1._rp ) then !to avoid overflow in exp
1726  aerosol_procs(ia_m3,is0,ik,ic)=aerosol_procs(ia_m3,is0,ik,ic) &
1727  *exp( dt_m3(is0,ik,ic)/aerosol_procs(ia_m3,is0,ik,ic) )
1728  else !aerosol_procs = 0
1729  aerosol_procs(ia_m3,is0,ik,ic)=aerosol_procs(ia_m3,is0,ik,ic) + dt_m3(is0,ik,ic)
1730  endif
1731  if (sixth(is0,ik,ic) > 0._rp .AND. & !to avoid divided by zero
1732  dt_m6(is0,ik,ic)/sixth(is0,ik,ic) < 1._rp ) then !to avoid overflow in exp
1733  sixth(is0,ik,ic) =sixth(is0,ik,ic) &
1734  *exp( dt_m6(is0,ik,ic)/sixth(is0,ik,ic) )
1735  else !aerosol_procs = 0
1736  sixth(is0,ik,ic) =sixth(is0,ik,ic) + dt_m6(is0,ik,ic)
1737  endif
1738  if (aerosol_procs(ia_ms,is0,ik,ic) > 0._rp .AND. & !to avoid divided by zero
1739  dt_ms(is0,ik,ic)/aerosol_procs(ia_ms,is0,ik,ic) < 1._rp ) then !to avoid overflow in exp
1740  aerosol_procs(ia_ms,is0,ik,ic)=aerosol_procs(ia_ms,is0,ik,ic) &
1741  *exp( dt_ms(is0,ik,ic)/aerosol_procs(ia_ms,is0,ik,ic) )
1742  else !aerosol_procs = 0
1743  aerosol_procs(ia_ms,is0,ik,ic)=aerosol_procs(ia_ms,is0,ik,ic) + dt_ms(is0,ik,ic)
1744  endif
1745  m0t = aerosol_procs(ia_m0,is0,ik,ic)
1746  m3t = aerosol_procs(ia_m3,is0,ik,ic)
1747  m6t = sixth(is0,ik,ic)
1748  call diag_ds6(m0t,m3t,m6t, &
1749  m2t,dgt,sgt,dm2)
1750  if (dgt <= 0._rp) dgt = d_ct(is0,ic)
1751  if (sgt <= 0._rp) sgt = 1.3_rp
1752  if (m3t == 0._rp) aerosol_procs(ia_ms,is0,ik,ic)=0._rp
1753  aerosol_procs(ia_m2,is0,ik,ic)=m2t
1754  ! Y.Sato added
1755  aerosol_procs(ia_m0,is0,ik,ic)=m0t
1756  aerosol_procs(ia_m3,is0,ik,ic)=m3t
1757  enddo
1758  enddo
1759  enddo
1760 
1761  !== moving center
1762  do ic = 1, n_ctg !aerosol category
1763  do ik = 1, n_kap(ic) !kappa bin
1764 
1765  rm0_pls = 0._rp
1766  rm2_pls = 0._rp
1767  rm3_pls = 0._rp
1768  rms_pls = 0._rp
1769  rm0_mns = 0._rp
1770  rm2_mns = 0._rp
1771  rm3_mns = 0._rp
1772  rms_mns = 0._rp
1773 
1774  do is1 = 1, n_siz(ic) !size bin
1775 
1776  if (aerosol_procs(ia_m0,is1,ik,ic) &
1777  *aerosol_procs(ia_m2,is1,ik,ic) &
1778  *aerosol_procs(ia_m3,is1,ik,ic) > 0._rp) then
1779  m0t = aerosol_procs(ia_m0,is1,ik,ic)
1780  m2t = aerosol_procs(ia_m2,is1,ik,ic)
1781  m3t = aerosol_procs(ia_m3,is1,ik,ic)
1782  call diag_ds(m0t,m2t,m3t,dgt,sgt,dm2)
1783  if (dgt <= 0._rp) dgt = d_ct(is1,ic)
1784  if (sgt <= 0._rp) sgt = 1.3_rp
1785 
1786  do is2 = 1, n_siz(ic)
1787  if (dgt >= d_lw(is2,ic) .AND. dgt < d_up(is2,ic)) then !moving center
1788  rm0_pls(is2) = rm0_pls(is2) + aerosol_procs(ia_m0,is1,ik,ic) !is2 <=
1789  rm0_mns(is1) = aerosol_procs(ia_m0,is1,ik,ic) !is1 =>
1790  rm2_pls(is2) = rm2_pls(is2) + aerosol_procs(ia_m2,is1,ik,ic) !is2 <=
1791  rm2_mns(is1) = aerosol_procs(ia_m2,is1,ik,ic) !is1 =>
1792  rm3_pls(is2) = rm3_pls(is2) + aerosol_procs(ia_m3,is1,ik,ic) !is2 <=
1793  rm3_mns(is1) = aerosol_procs(ia_m3,is1,ik,ic) !is1 =>
1794  rms_pls(is2) = rms_pls(is2) + aerosol_procs(ia_ms,is1,ik,ic) !is2 <=
1795  rms_mns(is1) = aerosol_procs(ia_ms,is1,ik,ic) !is1 =>
1796  exit
1797  endif !d_lw(is2) < dgt < d_up(is2)
1798  enddo
1799 
1800  endif !aerosol number > 0
1801 
1802  enddo
1803 
1804  do is0 = 1, n_siz(ic)
1805  aerosol_procs(ia_m0,is0,ik,ic) = aerosol_procs(ia_m0,is0,ik,ic) &
1806  + rm0_pls(is0) - rm0_mns(is0)
1807  aerosol_procs(ia_m2,is0,ik,ic) = aerosol_procs(ia_m2,is0,ik,ic) &
1808  + rm2_pls(is0) - rm2_mns(is0)
1809  aerosol_procs(ia_m3,is0,ik,ic) = aerosol_procs(ia_m3,is0,ik,ic) &
1810  + rm3_pls(is0) - rm3_mns(is0)
1811  aerosol_procs(ia_ms,is0,ik,ic) = aerosol_procs(ia_ms,is0,ik,ic) &
1812  + rms_pls(is0) - rms_mns(is0)
1813  enddo !is(1:n_siz(ic))
1814 
1815  enddo
1816  enddo
1817 
1818 
1819  return
1820  end subroutine aerosol_coagulation
1821  !----------------------------------------------------------------------------------------
1822  ! subroutine 4. aerosol_activation
1823  ! Abdul-Razzak et al., JGR, 1998 [AR98]
1824  ! Abdul-Razzak and Ghan, JGR, 2000 [AR00]
1825  !----------------------------------------------------------------------------------------
1826  subroutine aerosol_activation(c_kappa, super, temp_k, ia_m0, ia_m2, ia_m3, &
1827  n_atr,n_siz_max,n_kap_max,n_ctg,n_siz,n_kap, &
1828  d_ct, aerosol_procs, aerosol_activ)
1829 
1830  implicit none
1831  !i/o variables
1832  real(rp),intent(in) :: super ! supersaturation [-]
1833  real(rp),intent(in) :: c_kappa ! hygroscopicity of condensable mass [-]
1834  real(rp),intent(in) :: temp_k ! temperature
1835  integer, intent(in) :: ia_m0, ia_m2, ia_m3
1836  integer, intent(in) :: n_atr
1837  integer, intent(in) :: n_siz_max
1838  integer, intent(in) :: n_kap_max
1839  integer, intent(in) :: n_ctg
1840  integer, intent(in) :: n_siz(n_ctg), n_kap(n_ctg)
1841  real(rp),intent(in) :: d_ct(n_siz_max,n_ctg)
1842  real(rp),intent(in) :: aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
1843 
1844  real(rp), intent(out) :: aerosol_activ(n_atr,n_siz_max,n_kap_max,n_ctg)
1845 
1846  !local variables
1847  real(rp),parameter :: two3 = 2._rp/3._rp
1848  real(rp),parameter :: rt2 = sqrt(2._rp)
1849  real(rp),parameter :: twort2 = rt2 ! 2/sqrt(2) = sqrt(2)
1850  real(rp),parameter :: thrrt2 = 3._rp/rt2 ! 3/sqrt(2)
1851  real(rp) :: smax_inv ! inverse
1852  real(rp) :: am,scrit_am,aa,tc,st,bb,ac
1853  real(rp) :: m0t,m2t,m3t,dgt,sgt,dm2
1854  real(rp) :: d_crit ! critical diameter
1855  real(rp) :: tmp1, tmp2, tmp3 !
1856  real(rp) :: ccn_frc,cca_frc,ccv_frc ! activated number,area,volume
1857  integer :: is0, ik, ic
1858 
1859  aerosol_activ(:,:,:,:) = 0._rp
1860 
1861  if (super<=0._rp) return
1862 
1863  smax_inv = 1._rp / super
1864 
1865  !--- surface tension of water
1866  tc = temp_k - stdtemp
1867  if (tc >= 0._rp ) then
1868  st = 75.94_rp-0.1365_rp*tc-0.3827e-3_rp*tc**2._rp !Gittens, JCIS, 69, Table 4.
1869  else !t[deg C]<0.
1870  st = 75.93_rp +0.115_rp*tc & !Eq.5-12, pp.130, PK97
1871  + 6.818e-2_rp*tc**2._rp+6.511e-3_rp*tc**3._rp &
1872  + 2.933e-4_rp*tc**4._rp+6.283e-6_rp*tc**5._rp &
1873  + 5.285e-8_rp*tc**6._rp
1874  endif
1875  st = st * 1.e-3_rp ![J/m2]
1876 
1877  !-- Kelvin effect
1878  ! [J m-2] [kg mol-1] [m3 kg-1] [mol K J-1] [K-1]
1879  aa = 2._rp * st * mwwat * 1.e-3_rp / (dnwat * rgas * temp_k ) ![m] Eq.5 in AR98
1880 
1881  do ic = 1, n_ctg
1882  do ik = 1, n_kap(ic)
1883  do is0 = 1, n_siz(ic)
1884  m0t = aerosol_procs(ia_m0,is0,ik,ic)
1885  m2t = aerosol_procs(ia_m2,is0,ik,ic)
1886  m3t = aerosol_procs(ia_m3,is0,ik,ic)
1887  call diag_ds(m0t,m2t,m3t,dgt,sgt,dm2)
1888  if (dgt <= 0._rp) dgt = d_ct(is0,ic)
1889  if (sgt <= 0._rp) sgt = 1.3_rp
1890  am = dgt * 0.5_rp !geometric dry mean radius [m]
1891  bb = c_kappa
1892  if (bb > 0._rp .AND. am > 0._rp ) then
1893  scrit_am = 2._rp/sqrt(bb)*(aa/(3._rp*am))**1.5_rp !AR00 Eq.9
1894  else
1895  scrit_am = 0._rp
1896  endif
1897  ac = am * (scrit_am * smax_inv)**two3 !AR00 Eq.12
1898  d_crit = ac * 2._rp
1899  tmp1 = log(d_crit) - log(dgt)
1900  tmp2 = 1._rp/(rt2*log(sgt))
1901  tmp3 = log(sgt)
1902  ccn_frc= 0.5_rp*(1._rp-erf(tmp1*tmp2))
1903  cca_frc= 0.5_rp*(1._rp-erf(tmp1*tmp2-twort2*tmp3))
1904  ccv_frc= 0.5_rp*(1._rp-erf(tmp1*tmp2-thrrt2*tmp3))
1905  ccn_frc=min(max(ccn_frc,0.0_rp),1.0_rp)
1906  cca_frc=min(max(cca_frc,0.0_rp),1.0_rp)
1907  ccv_frc=min(max(ccv_frc,0.0_rp),1.0_rp)
1908  aerosol_activ(ia_m0,is0,ik,ic) = ccn_frc * aerosol_procs(ia_m0,is0,ik,ic)
1909  aerosol_activ(ia_m2,is0,ik,ic) = cca_frc * aerosol_procs(ia_m2,is0,ik,ic)
1910  aerosol_activ(ia_m3,is0,ik,ic) = ccv_frc * aerosol_procs(ia_m3,is0,ik,ic)
1911  aerosol_activ(ia_ms,is0,ik,ic) = ccv_frc * aerosol_procs(ia_ms,is0,ik,ic)
1912  enddo !is(1:n_siz(ic))
1913  enddo !ik(1:n_kap(ic))
1914  enddo !ic(1:n_ctg)
1915 
1916  return
1917  end subroutine aerosol_activation
1918  !-----------------------------------------------------------------------------
1919  ! subroutine 5. diag_ds
1920  !----------------------------------------------------------------------------------------
1921  subroutine diag_ds(m0,m2,m3, & !i
1922  dg,sg,dm2) !o
1923  implicit none
1924  real(rp) :: m0,m2,m3,dg,sg,m3_bar,m2_bar
1925  real(rp) :: m2_new,m2_old,dm2
1926  real(rp), parameter :: sgmax=2.5_rp
1927  real(rp), parameter :: rk1=2._rp
1928  real(rp), parameter :: rk2=3._rp
1929  real(rp), parameter :: ratio =rk1/rk2
1930  real(rp), parameter :: rk1_hat=1._rp/(ratio*(rk2-rk1))
1931  real(rp), parameter :: rk2_hat=ratio/(rk1-rk2)
1932  real(dp), parameter :: tiny=1.e-50_dp
1933 
1934  dm2=0._rp
1935 
1936  if (m0 <= tiny .OR. m2 <= tiny .OR. m3 <= tiny) then
1937  m0=0._rp
1938  m2=0._rp
1939  m3=0._rp
1940  dg=-1._rp
1941  sg=-1._rp
1942  return
1943  endif
1944 
1945  m2_old = m2
1946  m3_bar = m3/m0
1947  m2_bar = m2/m0
1948  dg = m2_bar**rk1_hat*m3_bar**rk2_hat
1949 
1950  if (m2_bar/m3_bar**ratio < 1._rp) then !stdev > 1.
1951  sg = exp(sqrt(2._rp/(rk1*(rk1-rk2)) &
1952  * log(m2_bar/m3_bar**ratio) ))
1953  endif
1954 
1955  if (sg > sgmax) then
1956  ! print *,'sg=',sg
1957  sg = sgmax
1958  call diag_d2(m0,m3,sg, & !i
1959  m2,dg ) !o
1960  ! print *,'warning! modified as sg exceeded sgmax (diag_ds)'
1961  endif
1962 
1963  if (m2_bar/m3_bar**ratio >= 1._rp) then !stdev < 1.
1964  sg = 1._rp
1965  call diag_d2(m0,m3,sg, & !i
1966  m2,dg ) !o
1967  ! print *,'warning! modified as sg < 1. (diag_ds)'
1968  endif
1969 
1970  m2_new = m2
1971  dm2 = m2_old - m2_new !m2_pres - m2_diag
1972 
1973  return
1974  end subroutine diag_ds
1975  !----------------------------------------------------------------------------------------
1976  ! subroutine 6. diag_d2
1977  !----------------------------------------------------------------------------------------
1978  subroutine diag_d2(m0,m3,sg, & !i
1979  m2,dg ) !o
1980  implicit none
1981  real(rp) :: dg,sg,m0,m2,m3,aaa
1982  real(rp), parameter :: one3=1._rp/3._rp
1983  aaa = m0 * exp( 4.5_rp * (log(sg)**2._rp) )
1984  dg =(m3/aaa)**one3
1985  m2 = m0 * dg ** 2._rp * exp( 2.0_rp * (log(sg)**2._rp) )
1986 
1987  return
1988  end subroutine diag_d2
1989  !----------------------------------------------------------------------------------------
1990  ! subroutine 7. aero_intra :: intra_sectional Brownian coagulation
1991  !----------------------------------------------------------------------------------------
1992  subroutine aero_intra(m0, m3, & !input
1993  dm0,dm3,dm6,dmi, & !output
1994  dg, sg, lambda, & !input
1995  rknc,rkfm,dt) !input
1996 
1997  !---intra-modal coagulation for spheres
1998  implicit none
1999 
2000  !=== Input
2001  real(rp), intent(in) :: m0 ![/m3]
2002  real(rp), intent(in) :: m3 ![m2/m3]
2003  ! real(RP), intent(in) :: m6 ![m3/m3]
2004  real(rp), intent(out) :: dm0 ![/m3]
2005  real(rp), intent(out) :: dm3 ![m2/m3]
2006  real(rp), intent(out) :: dm6 ![m3/m3]
2007  real(rp), intent(out) :: dmi ! [ug/m3/s]
2008  real(rp), intent(in) :: dg ![m]
2009  real(rp), intent(in) :: sg ![-]
2010  real(rp), intent(in) :: lambda ! mean free path of air molecules [cm]
2011  real(rp), intent(in) :: rknc, rkfm ! constants for both regimes
2012  real(dp), intent(in) :: dt ! dt
2013 
2014  !=== Intermidiate quantities
2015  real(rp) :: mm2,mm1,m1,m4,m2,mm1p5,mm0p5,m0p5,m3p5,m1p5,m5,m2p5
2016  real(rp) :: dm0dt,dm6dt,dm0dt_nc,dm6dt_nc,dm0dt_fm,dm6dt_fm,bbb
2017 
2018  !=== Moment formulation for M0, M3 & M6 to extract Dg and Sg
2019  real(rp), parameter :: rk1=3._rp
2020  real(rp), parameter :: rk2=6._rp
2021  real(rp), parameter :: ratio=rk1/rk2
2022  real(rp), parameter :: rk1_hat=1._rp/(ratio*(rk2-rk1))
2023  real(rp), parameter :: rk2_hat=ratio/(rk1-rk2)
2024 
2025  !--- initialization
2026  dm0 = 0._rp
2027  dm3 = 0._rp
2028  dm6 = 0._rp
2029  dm0dt = 0._rp
2030  dm6dt = 0._rp
2031 
2032  !--- near continuum regime---------------------------------------------------------
2033  ! M(k) = M(0) * dmean ** k * exp( k**2.*0.5 *(dlog(stdev)**2.))
2034  mm2 =m0*dg**(-2._rp) *exp(2.0_rp *(log(sg)**2._rp)) !M-2 (dM0/dt)
2035  mm1 =m0*dg**(-1._rp) *exp(0.5_rp *(log(sg)**2._rp)) !M-1 (dM0/dt)
2036  m1 =m0*dg** 1._rp *exp(0.5_rp *(log(sg)**2._rp)) !M1 (dM0/dt dM6/dt)
2037  m2 =m0*dg** 2._rp *exp(2.0_rp *(log(sg)**2._rp)) !M2
2038  m4 =m0*dg** 4._rp *exp(8.0_rp *(log(sg)**2._rp)) !M4 (dM6/dt)
2039  ! m6 =m0*dg** 6._RP *exp(18._RP *(log(sg)**2._RP)) !M6
2040 
2041  dm0dt_nc=-1._rp*rknc*(m0*m0+m1*mm1+2.492e-2_rp*lambda*(m0*mm1+m1*mm2))
2042  dm6dt_nc= 2._rp*rknc*(m3*m3+m4*m2 +2.492e-2_rp*lambda*(m3*m2 +m4*m1 ))
2043 
2044  !--- free molecule regime----------------------------------------------------------
2045  ! M(k) = M(0 * dmean ** k * exp( k**2.*0.5 *(dlog(stdev)**2.) )
2046  mm1p5=m0*dg**(-1.5_rp)*exp(1.125_rp*(log(sg)**2._rp))!M-1.5(dM0/dt)
2047  mm0p5=m0*dg**(-0.5_rp)*exp(0.125_rp*(log(sg)**2._rp))!M-0.5(dM0/dt)
2048  m0p5 =m0*dg** 0.5_rp *exp(0.125_rp*(log(sg)**2._rp))!M0.5 (dM0/dt)
2049  m3p5 =m0*dg**( 3.5_rp)*exp(6.125_rp*(log(sg)**2._rp))!M3.5 (dM6/dt)
2050  m1p5 =m0*dg**( 1.5_rp)*exp(1.125_rp*(log(sg)**2._rp))!M1.5 (dM6/dt)
2051  m5 =m0*dg** 5._rp *exp(12.50_rp*(log(sg)**2._rp))!M5 (dM6/dt)
2052  m2p5 =m0*dg** 2.5_rp *exp(3.125_rp*(log(sg)**2._rp))!M2.5 (dM6/dt)
2053 
2054  bbb =1._rp+1.2_rp *exp(-2._rp*sg)-0.646_rp*exp(-0.35_rp*sg**2._rp) ! Eq.9 in Park et al. (1999) JAS
2055 
2056  dm0dt_fm= -1._rp*bbb*rkfm*(m0*m0p5+m2*mm1p5+2._rp*m1*mm0p5)
2057  dm6dt_fm= 2._rp*bbb*rkfm*(m3p5*m3+m1p5*m5 +2._rp*m2p5*m4 )
2058 
2059  !--- harmonic mean approach
2060  if (dm0dt_fm+dm0dt_nc /= 0._rp) dm0dt = dm0dt_fm*dm0dt_nc/(dm0dt_fm+dm0dt_nc)
2061  if (dm6dt_fm+dm6dt_nc /= 0._rp) dm6dt = dm6dt_fm*dm6dt_nc/(dm6dt_fm+dm6dt_nc)
2062 
2063  dm0 = dm0dt * dt
2064  dm3 = 0._rp
2065  dm6 = dm6dt * dt
2066 
2067  dmi = 0.0_rp !! should be check (S. Nishizawa 2018/3/2)
2068 
2069  return
2070  end subroutine aero_intra
2071  !----------------------------------------------------------------------------------------
2072  ! subroutine 8. aero_inter :: inter_sectional Brownian coagulation
2073  !----------------------------------------------------------------------------------------
2074  subroutine aero_inter(m0i ,m3i ,m6i , & !input unchanged
2075  m0j ,m3j ,m6j , & !input unchanged
2076  dm0i,dm3i,dm6i,dmsi, & !output
2077  dm0j,dm3j,dm6j,dmsj, & !output
2078  dm0k,dm3k,dm6k,dmsk, & !output
2079  dgi ,sgi, rhoi, & !input unchanged
2080  dgj ,sgj, rhoj, & !input unchanged
2081  rknc, rkfm, lambda, & !input unchanged
2082  dtrest ) !input unchanged
2083 
2084  !---inter-modal coagulation for spheres
2085  implicit none
2086 
2087  !=== input/output variables
2088  real(rp), intent(in) :: m0i,m0j ! [/m3]
2089  real(rp), intent(in) :: m3i,m3j ! [m3/m3]
2090  real(rp), intent(in) :: m6i,m6j ! [m6/m3]
2091  real(rp), intent(out) :: dm0i,dm0j,dm0k ! [/m3/s]
2092  real(rp), intent(out) :: dm3i,dm3j,dm3k ! [m3/m3/s]
2093  real(rp), intent(out) :: dm6i,dm6j,dm6k ! [m6/m3/s]
2094  real(rp), intent(out) :: dmsi,dmsj,dmsk ! [ug/m3/s]
2095  real(rp), intent(in) :: dgi,sgi,rhoi ! [m][-][g/cm3]
2096  real(rp), intent(in) :: dgj,sgj,rhoj ! [m][-][g/cm3]
2097  real(rp), intent(in) :: rknc, rkfm ! constants for both regimes
2098  real(rp), intent(in) :: lambda ! mean free path of air molecules [cm]
2099  real(dp) :: dtrest ! time[s]
2100 
2101  !=== local variables
2102  real(rp) :: dm0dt_i,dm3dt_i,dm6dt_i
2103  real(rp) :: dm0dt_j,dm3dt_j,dm6dt_j
2104  real(rp) :: dm0dt_k,dm3dt_k,dm6dt_k
2105  real(rp) :: mm2i,mm1i,m1i,m2i,m4i,m5i,m7i,m8i
2106  real(rp) :: mm2j,mm1j,m1j,m2j,m4j,m5j,m7j,m8j
2107  real(rp) :: mm1p5i,mm0p5i,m0p5i,m1p5i,m2p5i,m3p5i,m4p5i,m5p5i,m6p5i
2108  real(rp) :: mm1p5j,mm0p5j,m0p5j,m1p5j,m2p5j,m3p5j,m4p5j,m5p5j,m6p5j
2109  real(rp) :: dm6dt_nc_i,dm3dt_nc_i,dm0dt_nc_i
2110  real(rp) :: dm6dt_nc_j,dm3dt_nc_j,dm0dt_nc_j
2111  real(rp) :: dm6dt_nc_k,dm3dt_nc_k,dm0dt_nc_k
2112  real(rp) :: dm6dt_fm_i,dm3dt_fm_i,dm0dt_fm_i
2113  real(rp) :: dm6dt_fm_j,dm3dt_fm_j,dm0dt_fm_j
2114  real(rp) :: dm6dt_fm_k,dm3dt_fm_k,dm0dt_fm_k
2115  real(rp) :: bbb,gamma1,gamma2,alpha,beta
2116 
2117  !=== initialization
2118  dm0dt_i = 0._rp
2119  dm3dt_i = 0._rp
2120  dm6dt_i = 0._rp
2121  dm0dt_j = 0._rp
2122  dm3dt_j = 0._rp
2123  dm6dt_j = 0._rp
2124  dm0dt_k = 0._rp
2125  dm3dt_k = 0._rp
2126  dm6dt_k = 0._rp
2127  dm0i = 0._rp
2128  dm3i = 0._rp
2129  dm6i = 0._rp
2130  dmsi = 0._rp
2131  dm0j = 0._rp
2132  dm3j = 0._rp
2133  dm6j = 0._rp
2134  dmsj = 0._rp
2135  dm0k = 0._rp
2136  dm3k = 0._rp
2137  dm6k = 0._rp
2138  dmsk = 0._rp
2139 
2140  !---near continuum regime----------------------------------------------------------------------
2141  ! (mode i)
2142  mm2i=m0i*dgi**(-2._rp)*exp( 2.0_rp*(log(sgi)**2._rp))!M-2
2143  mm1i=m0i*dgi**(-1._rp)*exp( 0.5_rp*(log(sgi)**2._rp))!M-1
2144  m1i =m0i*dgi** 1._rp *exp( 0.5_rp*(log(sgi)**2._rp))!M1
2145  m2i =m0i*dgi** 2._rp *exp( 2.0_rp*(log(sgi)**2._rp))!M2
2146  ! m3i =m3i
2147  m4i =m0i*dgi** 4._rp *exp( 8.0_rp*(log(sgi)**2._rp))!M4
2148  m5i =m0i*dgi** 5._rp *exp(12.5_rp*(log(sgi)**2._rp))!M5
2149  ! m6i =m6i
2150  m7i =m0i*dgi** 7._rp *exp(24.5_rp*(log(sgi)**2._rp))!M7
2151  ! (mode j)
2152  mm2j=m0j*dgj**(-2._rp)*exp( 2.0_rp*(log(sgj)**2._rp))!M-2
2153  mm1j=m0j*dgj**(-1._rp)*exp( 0.5_rp*(log(sgj)**2._rp))!M-1
2154  m1j =m0j*dgj** 1._rp *exp( 0.5_rp*(log(sgj)**2._rp))!M1
2155  m2j =m0j*dgj** 2._rp *exp( 2.0_rp*(log(sgj)**2._rp))!M2
2156  ! m3j =m3j
2157  m4j =m0j*dgj** 4._rp *exp( 8.0_rp*(log(sgj)**2._rp))!M4
2158  m5j =m0j*dgj** 5._rp *exp(12.5_rp*(log(sgj)**2._rp))!M5
2159  ! m6j =m6j
2160  m7j =m0j*dgj** 7._rp *exp(24.5_rp*(log(sgj)**2._rp))!M7
2161 
2162  dm0dt_nc_i = -rknc*(2._rp*m0i*m0j+ mm1i*m1j &
2163  + mm1j*m1i &
2164  +2.492e-2_rp*lambda*(m0i *mm1j+m1i*mm2j &
2165  + m0j *mm1i+m1j*mm2i) )
2166  dm3dt_nc_i = -rknc*(2._rp*m3i*m0j+ m2i *m1j &
2167  + mm1j*m4i &
2168  +2.492e-2_rp*lambda*(m3i *mm1j+m4i*mm2j &
2169  + m0j *m2i +m1j*m1i ) )
2170  dm6dt_nc_i = -rknc*(2._rp*m6i*m0j+ m5i *m1j &
2171  + mm1j*m7i &
2172  +2.492e-2_rp*lambda*(m6i *mm1j+m7i*mm2j &
2173  + m0j *m5i +m1j*m4i ) )
2174  dm0dt_nc_j = -rknc*(2._rp*m0i*m0j+ mm1i*m1j & !=dm0dt_nc_i
2175  + mm1j*m1i &
2176  +2.492e-2_rp*lambda*(m0i *mm1j+m1i*mm2j &
2177  + m0j *mm1i+m1j*mm2i) )
2178  dm3dt_nc_j = -rknc*(2._rp*m0i*m3j+ mm1i*m4j &
2179  + m2j *m1i &
2180  +2.492e-2_rp*lambda*(m0i *m2j +m1i*m1j &
2181  + m3j *mm1i+m4j*mm2i) )
2182  dm6dt_nc_j = -rknc*(2._rp*m0i*m6j+ mm1i*m7j &
2183  + m5j *m1i &
2184  +2.492e-2_rp*lambda*(m0i *m5j +m1i*m4j &
2185  + m6j *mm1i+m7j*mm2i) )
2186  dm0dt_nc_k = -dm0dt_nc_i
2187  dm3dt_nc_k = -dm3dt_nc_i-dm3dt_nc_j
2188  dm6dt_nc_k = -dm6dt_nc_i-dm6dt_nc_j &
2189  +2._rp*rknc*(2._rp*m3i*m3j+ m2i *m4j &
2190  + m2j *m4i &
2191  +2.492e-2_rp*lambda*(m3i *m2j +m4i*m1j &
2192  + m3j *m2i +m4j*m1i) )
2193 
2194  !---free molecule regime-----------------------------------------------------------------------
2195  ! (mode i)
2196  mm1p5i=m0i*dgi**(-1.5_rp)*exp( 1.125_rp*(log(sgi)**2._rp) ) !M-1.5
2197  mm0p5i=m0i*dgi**(-0.5_rp)*exp( 0.125_rp*(log(sgi)**2._rp) ) !M-0.5
2198  m0p5i =m0i*dgi** 0.5_rp *exp( 0.125_rp*(log(sgi)**2._rp) ) !M0.5
2199  m1p5i =m0i*dgi** 1.5_rp *exp( 1.125_rp*(log(sgi)**2._rp) ) !M1.5
2200  m2p5i =m0i*dgi** 2.5_rp *exp( 3.125_rp*(log(sgi)**2._rp) ) !M2.5
2201  m3p5i =m0i*dgi** 3.5_rp *exp( 6.125_rp*(log(sgi)**2._rp) ) !M3.5
2202  m4p5i =m0i*dgi** 4.5_rp *exp(10.125_rp*(log(sgi)**2._rp) ) !M4.5
2203  m5p5i =m0i*dgi** 5.5_rp *exp(15.125_rp*(log(sgi)**2._rp) ) !M5.5
2204  m6p5i =m0i*dgi** 6.5_rp *exp(21.125_rp*(log(sgi)**2._rp) ) !M6.5
2205  m8i =m0i*dgi** 8._rp *exp(32.000_rp*(log(sgi)**2._rp) ) !M8
2206 
2207  ! (mode j)
2208  mm1p5j=m0j*dgj**(-1.5_rp)*exp( 1.125_rp*(log(sgj)**2._rp) ) !M-1.5
2209  mm0p5j=m0j*dgj**(-0.5_rp)*exp( 0.125_rp*(log(sgj)**2._rp) ) !M-0.5
2210  m0p5j =m0j*dgj** 0.5_rp *exp( 0.125_rp*(log(sgj)**2._rp) ) !M0.5
2211  m1p5j =m0j*dgj** 1.5_rp *exp( 1.125_rp*(log(sgj)**2._rp) ) !M1.5
2212  m2p5j =m0j*dgj** 2.5_rp *exp( 3.125_rp*(log(sgj)**2._rp) ) !M2.5
2213  m3p5j =m0j*dgj** 3.5_rp *exp( 6.125_rp*(log(sgj)**2._rp) ) !M3.5
2214  m4p5j =m0j*dgj** 4.5_rp *exp(10.125_rp*(log(sgj)**2._rp) ) !M4.5
2215  m5p5j =m0j*dgj** 5.5_rp *exp(15.125_rp*(log(sgj)**2._rp) ) !M5.5
2216  m6p5j =m0j*dgj** 6.5_rp *exp(21.125_rp*(log(sgj)**2._rp) ) !M6.5
2217  m8j =m0j*dgj** 8._rp *exp(32.000_rp*(log(sgj)**2._rp) ) !M8
2218 
2219  !---approximation function---------------------------------------------------------------------
2220  alpha = dgj/dgi
2221  beta =(1._rp-(sqrt(1._rp+alpha**3._rp)/(1._rp+sqrt(alpha**3._rp))))/(1._rp-1._rp/sqrt(2._rp))
2222  gamma1 =(sgi +alpha*sgj )/(1._rp+alpha)
2223  gamma2 =(sgi**2._rp +alpha*sgj**2._rp)/(1._rp+alpha)
2224  bbb = 1._rp+1.2_rp*beta*dexp(real(-2._rp*gamma1,kind=dp))-0.646_rp*beta*dexp(real(-0.35_rp*gamma2,kind=dp))! Kajino (2011) JAS
2225 
2226  dm0dt_fm_i=-bbb*rkfm*( &
2227  m0i *m0p5j +m2i *mm1p5j +2._rp*m1i *mm0p5j &
2228  + m0j *m0p5i +m2j *mm1p5i +2._rp*m1j *mm0p5i )
2229  dm3dt_fm_i=-bbb*rkfm*( &
2230  m3i *m0p5j +m5i *mm1p5j +2._rp*m4i *mm0p5j &
2231  + m0j *m3p5i +m2j *m1p5i +2._rp*m1j *m2p5i )
2232  dm6dt_fm_i=-bbb*rkfm*( &
2233  m6i *m0p5j +m8i *mm1p5j +2._rp*m7i *mm0p5j &
2234  + m0j *m6p5i +m2j *m4p5i +2._rp*m1j *m5p5i )
2235  dm0dt_fm_j= dm0dt_fm_i
2236  dm3dt_fm_j=-bbb*rkfm*( &
2237  m0i *m3p5j +m2i *m1p5j +2._rp*m1i *m2p5j &
2238  + m3j *m0p5i +m5j *mm1p5i +2._rp*m4j *mm0p5i )
2239  dm6dt_fm_j=-bbb*rkfm*( &
2240  m0i *m6p5j +m2i *m4p5j +2._rp*m1i *m5p5j &
2241  + m6j *m0p5i +m8j *mm1p5i +2._rp*m7j *mm0p5i )
2242  dm0dt_fm_k=-dm0dt_fm_i
2243  dm3dt_fm_k=-dm3dt_fm_i-dm3dt_fm_j
2244  dm6dt_fm_k=-dm6dt_fm_i-dm6dt_fm_j &
2245  +2._rp * bbb*rkfm*( &
2246  m3i *m3p5j +m5i *m1p5j +2._rp*m4i *m2p5j &
2247  + m3j *m3p5i +m5j *m1p5i +2._rp*m4j *m2p5i )
2248 
2249  !---harmonic mean approach---------------------------------------------------------------------
2250  if (dm0dt_fm_i+dm0dt_nc_i/=0._rp) &
2251  dm0dt_i = dm0dt_fm_i*dm0dt_nc_i/(dm0dt_fm_i+dm0dt_nc_i)
2252  if (dm3dt_fm_i+dm3dt_nc_i/=0._rp) &
2253  dm3dt_i = dm3dt_fm_i*dm3dt_nc_i/(dm3dt_fm_i+dm3dt_nc_i)
2254  if (dm6dt_fm_i+dm6dt_nc_i/=0._rp) &
2255  dm6dt_i = dm6dt_fm_i*dm6dt_nc_i/(dm6dt_fm_i+dm6dt_nc_i)
2256  if (dm0dt_fm_j+dm0dt_nc_j/=0._rp) &
2257  dm0dt_j = dm0dt_fm_j*dm0dt_nc_j/(dm0dt_fm_j+dm0dt_nc_j)
2258  if (dm3dt_fm_j+dm3dt_nc_j/=0._rp) &
2259  dm3dt_j = dm3dt_fm_j*dm3dt_nc_j/(dm3dt_fm_j+dm3dt_nc_j)
2260  if (dm6dt_fm_j+dm6dt_nc_j/=0._rp) &
2261  dm6dt_j = dm6dt_fm_j*dm6dt_nc_j/(dm6dt_fm_j+dm6dt_nc_j)
2262  if (dm0dt_fm_k+dm0dt_nc_k/=0._rp) &
2263  dm0dt_k = dm0dt_fm_k*dm0dt_nc_k/(dm0dt_fm_k+dm0dt_nc_k)
2264  ! if (dm3dt_fm_k+dm3dt_nc_k/=0._RP) &
2265  ! dm3dt_k = dm3dt_fm_k*dm3dt_nc_k/(dm3dt_fm_k+dm3dt_nc_k)
2266  if (dm6dt_fm_k+dm6dt_nc_k/=0._rp) &
2267  dm6dt_k = dm6dt_fm_k*dm6dt_nc_k/(dm6dt_fm_k+dm6dt_nc_k)
2268 
2269  dm0i = dm0dt_i * dtrest
2270  dm3i = dm3dt_i * dtrest
2271  dm6i = dm6dt_i * dtrest
2272  dmsi = dm3i * rhoi *pi6*1.e12_rp !3rd to mass !m3/m3=>ug/m3
2273  dm0j = dm0dt_j * dtrest
2274  dm3j = dm3dt_j * dtrest
2275  dm6j = dm6dt_j * dtrest
2276  dmsj = dm3j * rhoj *pi6*1.e12_rp !3rd to mass !m3/m3=>ug/m3
2277  dm0k = dm0dt_k * dtrest
2278  ! dm3k = dm3dt_k * dtrest
2279  dm6k = dm6dt_k * dtrest
2280  dm3k = -dm3i -dm3j
2281  dmsk = -dmsi -dmsj
2282 
2283  return
2284  end subroutine aero_inter
2285  !----------------------------------------------------------------------------------------
2286  ! subroutine 9. diag_ds6
2287  !----------------------------------------------------------------------------------------
2288  subroutine diag_ds6(m0,m3,m6,m2,dg,sg,dm2)
2289  implicit none
2290  real(rp) :: m0,m2,m3,m6,dg,sg,m3_bar,m6_bar
2291  real(rp) :: m2_new,m2_old,dm2
2292  real(rp), parameter :: sgmax=2.5_rp
2293  real(rp), parameter :: rk1=3._rp
2294  real(rp), parameter :: rk2=6._rp
2295  real(rp), parameter :: ratio =rk1/rk2
2296  real(rp), parameter :: rk1_hat=1._rp/(ratio*(rk2-rk1))
2297  real(rp), parameter :: rk2_hat=ratio/(rk1-rk2)
2298  real(dp), parameter :: tiny=1.e-50_dp
2299 
2300  dm2=0._rp
2301  if (m0 <= tiny .OR. m3 <= tiny .OR. m6 <= tiny) then
2302  m0=0._rp
2303  m2=0._rp
2304  m3=0._rp
2305  dg=-1._rp
2306  sg=-1._rp
2307  return
2308  endif
2309 
2310  m2_old = m2
2311  m3_bar = m3/m0
2312  m6_bar = m6/m0
2313  dg = m3_bar**rk1_hat*m6_bar**rk2_hat
2314 
2315  if (m3_bar/m6_bar**ratio < 1._rp) then !stdev > 1.
2316  sg = exp(dsqrt(real(2._rp/(rk1*(rk1-rk2)) &
2317  * log(m3_bar/m6_bar**ratio), kind=dp )))
2318  m2 = m0*dg**2._rp*exp(2._rp*(log(sg)**2._rp))
2319  m2_old = m2
2320  endif
2321 
2322  if (sg > sgmax) then
2323  ! print *,'sg=',sg
2324  sg = sgmax
2325  call diag_d2(m0,m3,sg, & !input
2326  m2,dg ) !output
2327  ! print *,'warning! modified as sg exceeded sgmax (diag_ds6)'
2328  endif
2329 
2330  if (m3_bar/m6_bar**ratio >= 1._rp) then !stdev < 1.
2331  sg = 1._rp
2332  call diag_d2(m0,m3,sg, & !input
2333  m2,dg ) !output
2334  ! print *,'warning! modified as sg < 1.(diag_ds6)'
2335  endif
2336 
2337  m2_new = m2
2338  dm2 = m2_old - m2_new !m2_pres - m2_diag
2339 
2340 
2341  return
2342  end subroutine diag_ds6
2343  !-----------------------------------------------------------------------------
2344  subroutine trans_ccn(aerosol_procs, aerosol_activ, t_ccn, t_cn, &
2345  n_ctg, n_kap_max, n_siz_max, n_atr, &
2346  ic_out, ia_m0, ia_m2, ia_m3, ik_out, n_siz, &
2347 ! rnum_out, nbins_out, dlog10d_out)
2348  rnum_out, nbins_out)
2349 
2350  implicit none
2351  !i/o variables
2352  integer, intent(in) :: n_ctg, n_kap_max, n_siz_max, n_atr
2353  integer, intent(in) :: ic_out, ik_out
2354  integer, intent(in) :: ia_m0, ia_m2, ia_m3
2355  integer, intent(in) :: n_siz(n_ctg)
2356  real(rp),intent(in) :: aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
2357  real(rp),intent(in) :: aerosol_activ(n_atr,n_siz_max,n_kap_max,n_ctg)
2358  integer, intent(in) :: nbins_out
2359  real(rp),intent(out):: rnum_out(nbins_out)
2360 ! real(RP),intent(out):: dlog10d_out
2361  real(rp),intent(out):: t_ccn !total ccn number concentration [#/m3]
2362  real(rp),intent(out):: t_cn !total cn number concentration [#/m3]
2363  !local variables
2364 ! real(RP) :: m0t,m2t,m3t,dgt,sgt,dm2
2365 ! real(RP) :: d_lw2,d_up2,dg2,sg2,fnum0,fnum1
2366 ! real(RP) :: dlogd
2367  integer :: is0!, is_out
2368 ! real(RP), parameter :: d_max = 1.E-5_RP
2369 ! real(RP), parameter :: d_min = 1.E-9_RP
2370 ! real(RP) :: d_lw(nbins_out), d_up(nbins_out)
2371 
2372  rnum_out(:) = 0._rp
2373 !
2374 ! dlogd = (log(d_max) - log(d_min))/float(nbins_out)
2375 ! dlog10d_out = (log10(d_max)-log10(d_min))/float(nbins_out)
2376 !
2377 ! do is_out = 1, nbins_out !size bin
2378 ! d_lw(is_out) = exp(log(d_min)+dlogd* float(is_out-1) )
2379 ! d_up(is_out) = exp(log(d_min)+dlogd* float(is_out) )
2380 ! enddo !is_out
2381 !
2382  t_ccn = 0._rp
2383  t_cn = 0._rp
2384  do is0 = 1, n_siz(ic_out)
2385 ! if (aerosol_procs(ia_m0,is0,ik_out,ic_out) > 0._RP) then
2386 ! m0t = aerosol_procs(ia_m0,is0,ik_out,ic_out)
2387 ! m2t = aerosol_procs(ia_m2,is0,ik_out,ic_out)
2388 ! m3t = aerosol_procs(ia_m3,is0,ik_out,ic_out)
2389 ! call diag_ds(m0t,m2t,m3t,dgt,sgt,dm2)
2390 ! if (dgt <= 0._RP) dgt = d_ct(is0,ic_out)
2391 ! if (sgt <= 0._RP) sgt = 1.3_RP
2392 !
2393 ! do is_out = 1, nbins_out
2394 ! sgt = max(sgt,1.0000001_RP) !to avoid floating divide by zero
2395 ! d_lw2 = log(d_lw(is_out))
2396 ! d_up2 = log(d_up(is_out))
2397 ! dg2 = log(dgt)
2398 ! sg2 = log(sgt)
2399 ! fnum0 = m0t*0.5_RP*(1._RP+erf((d_lw2-dg2)/(sqrt(2.0_RP)*sg2)))
2400 ! fnum1 = m0t*0.5_RP*(1._RP+erf((d_up2-dg2)/(sqrt(2.0_RP)*sg2)))
2401 ! rnum_out(is_out) = rnum_out(is_out) + fnum1 - fnum0
2402 ! enddo
2403 !
2404 ! endif !number>0
2405 
2406  t_ccn = t_ccn + aerosol_activ(ia_m0,is0,ik_out,ic_out)
2407  t_cn = t_cn + aerosol_procs(ia_m0,is0,ik_out,ic_out)
2408 
2409  enddo
2410 
2411  return
2412  end subroutine trans_ccn
2413 
2414 
2415 end module scale_atmos_phy_ae_kajino13
scale_atmos_phy_ae_kajino13
module atmosphere / physics / aerosol / Kajino13
Definition: scale_atmos_phy_ae_kajino13.F90:12
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_h2so4dt
real(rp), public atmos_phy_ae_kajino13_h2so4dt
Definition: scale_atmos_phy_ae_kajino13.F90:55
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_dg_reg
real(rp), public atmos_phy_ae_kajino13_dg_reg
Definition: scale_atmos_phy_ae_kajino13.F90:64
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_aerosol::n_ae
integer, parameter, public n_ae
Definition: scale_atmos_aerosol.F90:33
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_name
character(len=h_short), dimension(:), allocatable, public atmos_phy_ae_kajino13_name
Definition: scale_atmos_phy_ae_kajino13.F90:51
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:159
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_tracer_setup
subroutine, public atmos_phy_ae_kajino13_tracer_setup(QA_AE)
Tracer setup.
Definition: scale_atmos_phy_ae_kajino13.F90:183
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_flag_npf
logical, public atmos_phy_ae_kajino13_flag_npf
Definition: scale_atmos_phy_ae_kajino13.F90:59
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_effective_radius
subroutine, public atmos_phy_ae_kajino13_effective_radius(KA, IA, JA, QA_AE, QTRC, RH, Re)
Calculate Effective Radius.
Definition: scale_atmos_phy_ae_kajino13.F90:897
scale_file_history
module file_history
Definition: scale_file_history.F90:15
float
typedef float(real32_t)
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_tendency
subroutine, public atmos_phy_ae_kajino13_tendency(KA, KS, KE, IA, IS, IE, JA, JS, JE, QA_AE, TEMP, PRES, QDRY, NREG, DENS, QV, QTRC, EMIT, dt, RHOQ_t_AE, CN, CCN)
Aerosol Microphysics.
Definition: scale_atmos_phy_ae_kajino13.F90:619
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_flag_ccn_interactive
logical, public atmos_phy_ae_kajino13_flag_ccn_interactive
Definition: scale_atmos_phy_ae_kajino13.F90:62
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_setup
subroutine, public atmos_phy_ae_kajino13_setup
Setup.
Definition: scale_atmos_phy_ae_kajino13.F90:357
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_desc
character(len=h_mid), dimension(:), allocatable, public atmos_phy_ae_kajino13_desc
Definition: scale_atmos_phy_ae_kajino13.F90:52
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:44
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_const::const_mdry
real(rp), public const_mdry
mass weight (dry air) [g/mol]
Definition: scale_const.F90:54
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_unit
character(len=h_short), dimension(:), allocatable, public atmos_phy_ae_kajino13_unit
Definition: scale_atmos_phy_ae_kajino13.F90:53
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_ocgasdt
real(rp), public atmos_phy_ae_kajino13_ocgasdt
Definition: scale_atmos_phy_ae_kajino13.F90:56
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_negative_fixer
subroutine, public atmos_phy_ae_kajino13_negative_fixer(KA, KS, KE, IA, IS, IE, JA, JS, JE, QA_AE, QTRC)
Definition: scale_atmos_phy_ae_kajino13.F90:1143
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_mkinit
subroutine, public atmos_phy_ae_kajino13_mkinit(KA, KS, KE, IA, IS, IE, JA, JS, JE, QA_AE, DENS, TEMP, PRES, QDRY, QV, m0_init, dg_init, sg_init, d_min_inp, d_max_inp, k_min_inp, k_max_inp, n_kap_inp, QTRC, CCN)
Definition: scale_atmos_phy_ae_kajino13.F90:940
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_logk_aenucl
real(rp), public atmos_phy_ae_kajino13_logk_aenucl
Definition: scale_atmos_phy_ae_kajino13.F90:66
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_flag_regeneration
logical, public atmos_phy_ae_kajino13_flag_regeneration
Definition: scale_atmos_phy_ae_kajino13.F90:63
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_atmos_aerosol
module atmosphere / aerosol
Definition: scale_atmos_aerosol.F90:12
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_c_kappa
real(rp), public atmos_phy_ae_kajino13_c_kappa
Definition: scale_atmos_phy_ae_kajino13.F90:58
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_flag_cond
logical, public atmos_phy_ae_kajino13_flag_cond
Definition: scale_atmos_phy_ae_kajino13.F90:60
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_sg_reg
real(rp), public atmos_phy_ae_kajino13_sg_reg
Definition: scale_atmos_phy_ae_kajino13.F90:65
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_phy_ae_kajino13::atmos_phy_ae_kajino13_flag_coag
logical, public atmos_phy_ae_kajino13_flag_coag
Definition: scale_atmos_phy_ae_kajino13.F90:61
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
Generated by   doxygen 1.8.17