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