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