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