PageRenderTime 80ms CodeModel.GetById 21ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/phys/module_mixactivate.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2743 lines | 1845 code | 334 blank | 564 comment | 153 complexity | aa80bcb1902394863ac73d83dbb3835e MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !************************************************************************
    2. 

  2. ! This computer software was prepared by Battelle Memorial Institute, 
    3. 

  3. ! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with
    4. 

  4. ! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE
    5. 

  5. ! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
    6. 

  6. ! LIABILITY FOR THE USE OF THIS SOFTWARE.
    7. 

  7. !
    8. 

  8. ! MOSAIC module: see chem/module_mosaic_driver.F for references and terms
    9. 

  9. ! of use
    10. 

  10. !************************************************************************
    11. 

  11. 

  12. MODULE module_mixactivate
    13. 

  13. PRIVATE
    14. 

  14. PUBLIC prescribe_aerosol_mixactivate, mixactivate
    15. 

  15. CONTAINS
    16. 

  16. 

  17. 

  18. !----------------------------------------------------------------------
    19. 

  19. !----------------------------------------------------------------------
    20. 

  20. ! 06-nov-2005 rce - grid_id & ktau added to arg list
    21. 

  21. ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
    22. 

  22.       subroutine prescribe_aerosol_mixactivate (                      &
    23. 

  23. 		grid_id, ktau, dtstep, naer,                          &
    24. 

  24. 		rho_phy, th_phy, pi_phy, w, cldfra, cldfra_old,       &
    25. 

  25. 		z, dz8w, p_at_w, t_at_w, exch_h,                      &
    26. 

  26.         qv, qc, qi, qndrop3d,                                 &
    27. 

  27.         nsource,                                              &
    28. 

  28. 		ids,ide, jds,jde, kds,kde,                            &
    29. 

  29. 		ims,ime, jms,jme, kms,kme,                            &
    30. 

  30. 		its,ite, jts,jte, kts,kte,                            &
    31. 

  31. 		f_qc, f_qi                                            )
    32. 

  32. 

  33. !        USE module_configure
    34. 

  34. 

  35. ! wrapper to call mixactivate for mosaic description of aerosol
    36. 

  36. 

  37. 	implicit none
    38. 

  38. 

  39. !   subr arguments
    40. 

  40. 	integer, intent(in) ::                  &
    41. 

  41. 		grid_id, ktau,                  &
    42. 

  42. 		ids, ide, jds, jde, kds, kde,   &
    43. 

  43. 		ims, ime, jms, jme, kms, kme,   &
    44. 

  44. 		its, ite, jts, jte, kts, kte
    45. 

  45. 

  46. 	real, intent(in) :: dtstep
    47. 

  47. 	real, intent(inout) :: naer ! aerosol number (/kg)
    48. 

  48. 

  49. 	real, intent(in),   &
    50. 

  50. 		dimension( ims:ime, kms:kme, jms:jme ) :: &
    51. 

  51. 		rho_phy, th_phy, pi_phy, w,  &
    52. 

  52. 		z, dz8w, p_at_w, t_at_w, exch_h
    53. 

  53. 

  54. 	real, intent(inout),   &
    55. 

  55. 		dimension( ims:ime, kms:kme, jms:jme ) :: cldfra, cldfra_old
    56. 

  56. 

  57. 	real, intent(in),   &
    58. 

  58. 		dimension( ims:ime, kms:kme, jms:jme ) :: &
    59. 

  59. 		qv, qc, qi
    60. 

  60. 

  61. 	real, intent(inout),   &
    62. 

  62. 		dimension( ims:ime, kms:kme, jms:jme ) :: &
    63. 

  63. 		qndrop3d
    64. 

  64. 

  65. 	real, intent(out),   &
    66. 

  66. 		dimension( ims:ime, kms:kme, jms:jme) :: nsource
    67. 

  67. 

  68.     LOGICAL, OPTIONAL :: f_qc, f_qi
    69. 

  69. 

  70. ! local vars
    71. 

  71. 	integer maxd_aphase, maxd_atype, maxd_asize, maxd_acomp, max_chem
    72. 

  72. 	parameter (maxd_aphase=2,maxd_atype=1,maxd_asize=1,maxd_acomp=1, max_chem=10)
    73. 

  73. 	real ddvel(its:ite, jts:jte, max_chem) ! dry deposition velosity
    74. 

  74. 	real qsrflx(ims:ime, jms:jme, max_chem) ! dry deposition flux of aerosol
    75. 

  75. 	real chem(ims:ime, kms:kme, jms:jme, max_chem) ! chem array
    76. 

  76. 	integer i,j,k,l,m,n,p
    77. 

  77. 	real hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk
    78. 

  78. 	integer ntype_aer, nsize_aer(maxd_atype),ncomp_aer(maxd_atype), nphase_aer
    79. 

  79.       	integer massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ),   &
    80. 

  80.       	  waterptr_aer( maxd_asize, maxd_atype ),   &
    81. 

  81.       	  numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
    82. 

  82. 	  ai_phase, cw_phase
    83. 

  83.         real dlo_sect( maxd_asize, maxd_atype ),   & ! minimum size of section (cm)
    84. 

  84.              dhi_sect( maxd_asize, maxd_atype ),   & ! maximum size of section (cm)
    85. 

  85. 	     sigmag_aer(maxd_asize, maxd_atype),   & ! geometric standard deviation of aerosol size dist
    86. 

  86. 	     dgnum_aer(maxd_asize, maxd_atype),    & ! median diameter (cm) of number distrib of mode
    87. 

  87. 	     dens_aer( maxd_acomp, maxd_atype),    & ! density (g/cm3) of material
    88. 

  88. 	     mw_aer( maxd_acomp, maxd_atype),      & ! molecular weight (g/mole)
    89. 

  89. 	     dpvolmean_aer(maxd_asize, maxd_atype)   ! mean-volume diameter (cm) of mode
    90. 

  90. ! terminology:  (pi/6) * (mean-volume diameter)**3 ==
    91. 

  91. !       (volume mixing ratio of section/mode)/(number mixing ratio)
    92. 

  92. 	real, dimension(ims:ime,kms:kme,jms:jme) :: &
    93. 

  93. 	     ccn1,ccn2,ccn3,ccn4,ccn5,ccn6  ! number conc of aerosols activated at supersat
    94. 

  94. 	integer idrydep_onoff
    95. 

  95. 	real, dimension(ims:ime,kms:kme,jms:jme) :: t_phy
    96. 

  96. 	integer msectional
    97. 

  97. 

  98. 

  99. 	  integer ptr
    100. 

  100. 	  real maer
    101. 

  101. 

  102.       if(naer.lt.1.)then
    103. 

  103. 	     naer=1000.e6 ! #/kg default value
    104. 

  104.       endif
    105. 

  105. 	  ai_phase=1
    106. 

  106. 	  cw_phase=2
    107. 

  107. 	  idrydep_onoff = 0
    108. 

  108. 	  msectional = 0
    109. 

  109. 

  110. 	  t_phy(its:ite,kts:kte,jts:jte)=th_phy(its:ite,kts:kte,jts:jte)*pi_phy(its:ite,kts:kte,jts:jte)
    111. 

  111. 

  112.       ntype_aer=maxd_atype
    113. 

  113.       do n=1,ntype_aer
    114. 

  114.          nsize_aer(n)=maxd_asize
    115. 

  115. 	 ncomp_aer(n)=maxd_acomp
    116. 

  116.       end do
    117. 

  117.       nphase_aer=maxd_aphase
    118. 

  118. 

  119. ! set properties for each type and size
    120. 

  120.        do n=1,ntype_aer
    121. 

  121.        do m=1,nsize_aer(n)
    122. 

  122.           dlo_sect( m,n )=0.01e-4    ! minimum size of section (cm)
    123. 

  123.           dhi_sect( m,n )=0.5e-4    ! maximum size of section (cm)
    124. 

  124. 	  sigmag_aer(m,n)=2.      ! geometric standard deviation of aerosol size dist
    125. 

  125. 	  dgnum_aer(m,n)=0.1e-4       ! median diameter (cm) of number distrib of mode
    126. 

  126. 	  dpvolmean_aer(m,n) = dgnum_aer(m,n) * exp( 1.5 * (log(sigmag_aer(m,n)))**2 )
    127. 

  127. 	  end do
    128. 

  128. 	  do l=1,ncomp_aer(n)
    129. 

  129. 	     dens_aer( l, n)=1.0   ! density (g/cm3) of material
    130. 

  130. 	     mw_aer( l, n)=132. ! molecular weight (g/mole)
    131. 

  131. 	  end do
    132. 

  132.       end do
    133. 

  133.        ptr=0
    134. 

  134.        do p=1,nphase_aer
    135. 

  135.        do n=1,ntype_aer
    136. 

  136.        do m=1,nsize_aer(n)
    137. 

  137.           ptr=ptr+1
    138. 

  138.           numptr_aer( m, n, p )=ptr
    139. 

  139. 	  if(p.eq.ai_phase)then
    140. 

  140. 	     chem(its:ite,kts:kte,jts:jte,ptr)=naer
    141. 

  141. 	  else
    142. 

  142. 	     chem(its:ite,kts:kte,jts:jte,ptr)=0.
    143. 

  143. 	  endif
    144. 

  144. 	end do ! size
    145. 

  145. 	end do ! type
    146. 

  146. 	end do ! phase
    147. 

  147.        do p=1,maxd_aphase
    148. 

  148.        do n=1,ntype_aer
    149. 

  149.        do m=1,nsize_aer(n)
    150. 

  150. 	  do l=1,ncomp_aer(n)
    151. 

  151.           ptr=ptr+1
    152. 

  152. 	     if(ptr.gt.max_chem)then
    153. 

  153. 	        write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
    154. 

  154. 	        call wrf_error_fatal("1")
    155. 

  155. 	     endif
    156. 

  156. 	     massptr_aer(l, m, n, p)=ptr
    157. 

  157. ! maer is ug/kg-air;  naer is #/kg-air;  dgnum is cm;  dens_aer is g/cm3
    158. 

  158. ! 1.e6 factor converts g to ug
    159. 

  159. 	     maer= 1.0e6 * naer * dens_aer(l,n) * ( (3.1416/6.) *   &
    160. 

  160.                  (dgnum_aer(m,n)**3) * exp( 4.5*((log(sigmag_aer(m,n)))**2) ) )
    161. 

  161. 	     if(p.eq.ai_phase)then
    162. 

  162. 	        chem(its:ite,kts:kte,jts:jte,ptr)=maer
    163. 

  163. 	     else
    164. 

  164. 	        chem(its:ite,kts:kte,jts:jte,ptr)=0.
    165. 

  165. 	     endif
    166. 

  166. 	  end do
    167. 

  167. 	end do ! size
    168. 

  168. 	end do ! type
    169. 

  169. 	end do ! phase
    170. 

  170.        do n=1,ntype_aer
    171. 

  171.        do m=1,nsize_aer(n)
    172. 

  172.           ptr=ptr+1
    173. 

  173. 	  if(ptr.gt.max_chem)then
    174. 

  174. 	     write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
    175. 

  175. 	     call wrf_error_fatal("1")
    176. 

  176. 	  endif
    177. 

  177. !wig	  waterptr_aer(m, n)=ptr
    178. 

  178. 	  waterptr_aer(m, n)=-1
    179. 

  179. 	end do ! size
    180. 

  180. 	end do ! type
    181. 

  181. 	ddvel(its:ite,jts:jte,:)=0.
    182. 

  182.     hygro(its:ite,kts:kte,jts:jte,:,:) = 0.5
    183. 

  183. 

  184. ! 06-nov-2005 rce - grid_id & ktau added to arg list
    185. 

  185.       call mixactivate(  msectional,     &
    186. 

  186.             chem,max_chem,qv,qc,qi,qndrop3d,        &
    187. 

  187.             t_phy, w, ddvel, idrydep_onoff,  &
    188. 

  188.             maxd_acomp, maxd_asize, maxd_atype, maxd_aphase,   &
    189. 

  189.             ncomp_aer, nsize_aer, ntype_aer, nphase_aer,  &
    190. 

  190.             numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dpvolmean_aer,  &
    191. 

  191.             dens_aer, mw_aer,           &
    192. 

  192.             waterptr_aer, hygro,  ai_phase, cw_phase,                &
    193. 

  193.             ids,ide, jds,jde, kds,kde,                            &
    194. 

  194.             ims,ime, jms,jme, kms,kme,                            &
    195. 

  195.             its,ite, jts,jte, kts,kte,                            &
    196. 

  196.             rho_phy, z, dz8w, p_at_w, t_at_w, exch_h,      &
    197. 

  197.             cldfra, cldfra_old, qsrflx,         &
    198. 

  198.             ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource,       &
    199. 

  199.             grid_id, ktau, dtstep, &
    200. 

  200.             F_QC=f_qc, F_QI=f_qi                              )
    201. 

  201. 

  202. 

  203.       end subroutine prescribe_aerosol_mixactivate
    204. 

  204. 

  205. !----------------------------------------------------------------------
    206. 

  206. !----------------------------------------------------------------------
    207. 

  207. !   nov-04 sg ! replaced amode with aer and expanded aerosol dimension to include type and phase
    208. 

  208. 

  209. ! 06-nov-2005 rce - grid_id & ktau added to arg list
    210. 

  210. ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
    211. 

  211. subroutine mixactivate(  msectional,            &
    212. 

  212.            chem, num_chem, qv, qc, qi, qndrop3d,         &
    213. 

  213.            temp, w, ddvel, idrydep_onoff,  &
    214. 

  214.            maxd_acomp, maxd_asize, maxd_atype, maxd_aphase,   &
    215. 

  215.            ncomp_aer, nsize_aer, ntype_aer, nphase_aer,  &
    216. 

  216.            numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dpvolmean_aer,  &
    217. 

  217.            dens_aer, mw_aer,               &
    218. 

  218.            waterptr_aer, hygro, ai_phase, cw_phase,              &
    219. 

  219.            ids,ide, jds,jde, kds,kde,                            &
    220. 

  220.            ims,ime, jms,jme, kms,kme,                            &
    221. 

  221.            its,ite, jts,jte, kts,kte,                            &
    222. 

  222.            rho, zm, dz8w, p_at_w, t_at_w, kvh,      &
    223. 

  223.            cldfra, cldfra_old, qsrflx,          &
    224. 

  224.            ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource,       &
    225. 

  225.            grid_id, ktau, dtstep, &
    226. 

  226.            f_qc, f_qi                       )
    227. 

  227. 

  228. 

  229. !     vertical diffusion and nucleation of cloud droplets
    230. 

  230. !     assume cloud presence controlled by cloud fraction
    231. 

  231. !     doesn't distinguish between warm, cold clouds
    232. 

    233. 

  233.   USE module_model_constants, only: g, rhowater, xlv, cp, rvovrd, r_d, r_v, mwdry, ep_2
    234. 

  234.   USE module_radiation_driver, only: cal_cldfra
    235. 

  235. 

  236.   implicit none
    237. 

  237. 

  238. !     input
    239. 

  239. 

  240.   INTEGER, intent(in) ::         grid_id, ktau
    241. 

  241.   INTEGER, intent(in) ::         num_chem
    242. 

  242.   integer, intent(in) ::         ids,ide, jds,jde, kds,kde,    &
    243. 

  243.                                  ims,ime, jms,jme, kms,kme,    &
    244. 

  244.                                  its,ite, jts,jte, kts,kte
    245. 

  245. 

  246.   integer maxd_aphase, nphase_aer, maxd_atype, ntype_aer
    247. 

  247.   integer maxd_asize, maxd_acomp, nsize_aer(maxd_atype)
    248. 

  248.   integer, intent(in) ::   &
    249. 

  249.        ncomp_aer( maxd_atype  ),   &
    250. 

  250.        massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ),   &
    251. 

  251.        waterptr_aer( maxd_asize, maxd_atype ),   &
    252. 

  252.        numptr_aer( maxd_asize, maxd_atype, maxd_aphase), &
    253. 

  253.        ai_phase, cw_phase
    254. 

  254.   integer, intent(in) :: msectional ! 1 for sectional, 0 for modal
    255. 

  255.   integer, intent(in) :: idrydep_onoff
    256. 

  256.   real, intent(in)  ::                       &
    257. 

  257.        dlo_sect( maxd_asize, maxd_atype ),   & ! minimum size of section (cm)
    258. 

  258.        dhi_sect( maxd_asize, maxd_atype ),   & ! maximum size of section (cm)
    259. 

  259.        sigmag_aer(maxd_asize, maxd_atype),   & ! geometric standard deviation of aerosol size dist
    260. 

  260.        dens_aer( maxd_acomp, maxd_atype),    & ! density (g/cm3) of material
    261. 

  261.        mw_aer( maxd_acomp, maxd_atype),      & ! molecular weight (g/mole)
    262. 

  262.        dpvolmean_aer(maxd_asize, maxd_atype)   ! mean-volume diameter (cm) of mode
    263. 

  263. ! terminology:  (pi/6) * (mean-volume diameter)**3 ==
    264. 

  264. !       (volume mixing ratio of section/mode)/(number mixing ratio)
    265. 

  265. 

  266. 

  267.   REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ) :: &
    268. 

  268.        chem ! aerosol molar mixing ratio (ug/kg or #/kg)
    269. 

  269. 

  270.   REAL, intent(in), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
    271. 

  271.        qv, qc, qi ! water species (vapor, cloud drops, cloud ice) mixing ratio (g/g)
    272. 

  272. 

  273.   LOGICAL, OPTIONAL :: f_qc, f_qi
    274. 

  274. 

  275.   REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
    276. 

  276.        qndrop3d    ! water species mixing ratio (g/g)
    277. 

  277. 

  278.   real, intent(in) :: dtstep             ! time step for microphysics (s)
    279. 

  279.   real, intent(in) :: temp(ims:ime, kms:kme, jms:jme)    ! temperature (K)
    280. 

  280.   real, intent(in) :: w(ims:ime, kms:kme, jms:jme)   ! vertical velocity (m/s)
    281. 

  281.   real, intent(in) :: rho(ims:ime, kms:kme, jms:jme)    ! density at mid-level  (kg/m3)
    282. 

  282.   REAL, intent(in) :: ddvel( its:ite, jts:jte, num_chem ) ! deposition velocity  (m/s)
    283. 

  283.   real, intent(in) :: zm(ims:ime, kms:kme, jms:jme)     ! geopotential height of level (m)
    284. 

  284.   real, intent(in) :: dz8w(ims:ime, kms:kme, jms:jme) ! layer thickness (m)
    285. 

  285.   real, intent(in) :: p_at_w(ims:ime, kms:kme, jms:jme) ! pressure at layer interface (Pa)
    286. 

  286.   real, intent(in) :: t_at_w(ims:ime, kms:kme, jms:jme) ! temperature at layer interface (K)
    287. 

  287.   real, intent(in) :: kvh(ims:ime, kms:kme, jms:jme)    ! vertical diffusivity (m2/s)
    288. 

  288.   real, intent(inout) :: cldfra_old(ims:ime, kms:kme, jms:jme)! cloud fraction on previous time step
    289. 

  289.   real, intent(inout) :: cldfra(ims:ime, kms:kme, jms:jme)    ! cloud fraction
    290. 

  290.   real, intent(in) :: hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk hygroscopicity   &
    291. 

  291. 

  292.   REAL, intent(out), DIMENSION( ims:ime, jms:jme, num_chem ) ::   qsrflx ! dry deposition rate for aerosol
    293. 

  293.   real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, &  ! droplet number source (#/kg/s)
    294. 

  294.        ccn1,ccn2,ccn3,ccn4,ccn5,ccn6  ! number conc of aerosols activated at supersat
    295. 

  295. 

  296. 

  297. !--------------------Local storage-------------------------------------
    298. 

  298. !
    299. 

  299.   real :: dgnum_aer(maxd_asize, maxd_atype) ! median diameter (cm) of number distrib of mode
    300. 

  300.   real :: qndrop(kms:kme)      ! cloud droplet number mixing ratio (#/kg)
    301. 

  301.   real :: lcldfra(kms:kme)     ! liquid cloud fraction
    302. 

  302.   real :: lcldfra_old(kms:kme) ! liquid cloud fraction for previous timestep
    303. 

  303.   real :: wtke(kms:kme)        ! turbulent vertical velocity at base of layer k (m2/s)
    304. 

  304.   real zn(kms:kme)             ! g/pdel (m2/g) for layer
    305. 

  305.   real zs(kms:kme)             ! inverse of distance between levels (m)
    306. 

  306.   real, parameter :: zkmin = 0.01
    307. 

  307.   real, parameter :: zkmax = 100.
    308. 

  308.   real cs(kms:kme)             ! air density (kg/m3) at layer center
    309. 

  309.   real csbot(kms:kme)          ! air density (kg/m3) at layer bottom
    310. 

  310.   real csbot_cscen(kms:kme)    ! csbot(k)/cs(k)
    311. 

  311.   real dz(kms:kme)             ! geometric thickness of layers (m)
    312. 

  312. 

  313.   real wdiab                   ! diabatic vertical velocity
    314. 

  314. !      real, parameter :: wmixmin = 0.1 ! minimum turbulence vertical velocity (m/s)
    315. 

  315.   real, parameter :: wmixmin = 0.2 ! minimum turbulence vertical velocity (m/s)
    316. 

  316. !      real, parameter :: wmixmin = 1.0 ! minimum turbulence vertical velocity (m/s)
    317. 

  317.   real :: qndrop_new(kms:kme)  ! droplet number nucleated on cloud boundaries
    318. 

  318.   real :: ekd(kms:kme)         ! diffusivity for droplets (m2/s)
    319. 

  319.   real :: ekk(kms:kme)         ! density*diffusivity for droplets (kg/m3 m2/s)
    320. 

  320.   real :: srcn(kms:kme)        ! droplet source rate (/s)
    321. 

  321.   real, parameter :: sq2pi = 2.5066282746
    322. 

  322.   real dtinv
    323. 

  323. 

  324.   integer km1,kp1
    325. 

  325.   real wbar,wmix,wmin,wmax
    326. 

  326.   real dum
    327. 

  327.   real tmpa, tmpb, tmpc, tmpc1, tmpc2, tmpd, tmpe, tmpf
    328. 

  328.   real tmpcourno
    329. 

  329.   real dact
    330. 

  330.   real fluxntot         ! (#/cm2/s)
    331. 

  331.   real fac_srflx
    332. 

  332.   real depvel_drop, depvel_tmp
    333. 

  333.   real, parameter :: depvel_uplimit = 1.0 ! upper limit for dep vels (m/s)
    334. 

  334.   real :: surfrate(num_chem) ! surface exchange rate (/s)
    335. 

  335.   real surfratemax      ! max surfrate for all species treated here
    336. 

  336.   real surfrate_drop    ! surfade exchange rate for droplelts
    337. 

  337.   real dtmin,tinv,dtt
    338. 

  338.   integer nsubmix,nsubmix_bnd
    339. 

  339.   integer i,j,k,m,n,nsub
    340. 

  340.   real dtmix
    341. 

  341.   real alogarg
    342. 

  342.   real qcld
    343. 

  343.   real pi
    344. 

  344.   integer nnew,nsav,ntemp
    345. 

  345.   real :: overlapp(kms:kme),overlapm(kms:kme) ! cloud overlap
    346. 

  346.   real ::  ekkp(kms:kme),ekkm(kms:kme) ! zn*zs*density*diffusivity
    347. 

  347. !  integer, save :: count_submix(100)=0 ! wig: Note that this is a no-no for tile threads with OMP
    348. 

  348. 

  349.   integer lnum,lnumcw,l,lmass,lmasscw,lsfc,lsfccw,ltype,lsig,lwater
    350. 

  350.   integer :: ntype(maxd_asize)
    351. 

  351. 

  352.   real ::  naerosol(maxd_asize, maxd_atype)    ! interstitial aerosol number conc (/m3)
    353. 

  353.   real ::  naerosolcw(maxd_asize, maxd_atype)  ! activated number conc (/m3)
    354. 

  354.   real ::   maerosol(maxd_acomp,maxd_asize, maxd_atype)   ! interstit mass conc (kg/m3)
    355. 

  355.   real ::   maerosolcw(maxd_acomp,maxd_asize, maxd_atype) ! activated mass conc (kg/m3)
    356. 

  356.   real ::   maerosol_tot(maxd_asize, maxd_atype)     ! species-total interstit mass conc (kg/m3)
    357. 

  357.   real ::   maerosol_totcw(maxd_asize, maxd_atype)   ! species-total activated mass conc (kg/m3)
    358. 

  358.   real ::   vaerosol(maxd_asize, maxd_atype) ! interstit+activated aerosol volume conc (m3/m3)
    359. 

  359.   real ::   vaerosolcw(maxd_asize, maxd_atype) ! activated aerosol volume conc (m3/m3)
    360. 

  360.   real ::   raercol(kms:kme,num_chem,2) ! aerosol mass, number mixing ratios
    361. 

  361.   real ::   source(kms:kme) !
    362. 

  362. 

  363.   real ::   fn(maxd_asize, maxd_atype)         ! activation fraction for aerosol number
    364. 

  364.   real ::   fs(maxd_asize, maxd_atype)         ! activation fraction for aerosol sfcarea
    365. 

  365.   real ::   fm(maxd_asize, maxd_atype)         ! activation fraction for aerosol mass
    366. 

  366.   integer ::   ncomp(maxd_atype)
    367. 

  367. 

  368.   real ::   fluxn(maxd_asize, maxd_atype)      ! number  activation fraction flux (m/s)
    369. 

  369.   real ::   fluxs(maxd_asize, maxd_atype)      ! sfcarea activation fraction flux (m/s)
    370. 

  370.   real ::   fluxm(maxd_asize, maxd_atype)      ! mass    activation fraction flux (m/s)
    371. 

  371.   real ::   flux_fullact(kms:kme)              ! 100%    activation fraction flux (m/s)
    372. 

  372. !     note:  activation fraction fluxes are defined as
    373. 

  373. !     fluxn = [flux of activated aero. number into cloud (#/m2/s)]
    374. 

  374. !           / [aero. number conc. in updraft, just below cloudbase (#/m3)]
    375. 

  375. 

  376.   real :: nact(kms:kme,maxd_asize, maxd_atype)  ! fractional aero. number  activation rate (/s)
    377. 

  377.   real :: mact(kms:kme,maxd_asize, maxd_atype)  ! fractional aero. mass    activation rate (/s)
    378. 

  378.   real :: npv(maxd_asize, maxd_atype) ! number per volume concentration (/m3)
    379. 

  379.   real scale
    380. 

  380. 

  381.   real :: hygro_aer(maxd_asize, maxd_atype)  ! hygroscopicity of aerosol mode
    382. 

  382.   real :: exp45logsig     ! exp(4.5*alogsig**2)
    383. 

  383.   real :: alogsig(maxd_asize, maxd_atype) ! natl log of geometric standard dev of aerosol
    384. 

  384.   integer, parameter :: psat=6  ! number of supersaturations to calc ccn concentration
    385. 

  385.   real ccn(kts:kte,psat)        ! number conc of aerosols activated at supersat
    386. 

  386.   real, parameter :: supersat(psat)= &! supersaturation (%) to determine ccn concentration
    387. 

  387.        (/0.02,0.05,0.1,0.2,0.5,1.0/)
    388. 

  388.   real super(psat) ! supersaturation
    389. 

  389.   real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m)
    390. 

  390.   real :: ccnfact(psat,maxd_asize, maxd_atype)
    391. 

  391.   real :: amcube(maxd_asize, maxd_atype) ! cube of dry mode radius (m)
    392. 

  392.   real :: argfactor(maxd_asize, maxd_atype)
    393. 

  393.   real aten ! surface tension parameter
    394. 

  394.   real t0 ! reference temperature
    395. 

  395.   real sm ! critical supersaturation
    396. 

  396.   real arg
    397. 

  397. 

  398.   integer,parameter :: icheck_colmass = 0
    399. 

  399.            ! icheck_colmass > 0 turns on mass/number conservation checking
    400. 

  400.            ! values of 1, 10, 100 produce less to more diagnostics
    401. 

  401.   integer :: colmass_worst_ij( 2, 0:maxd_acomp, maxd_asize, maxd_atype )
    402. 

  402.   integer :: colmass_maxworst_i(3)
    403. 

  403.   real :: colmass_bgn( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
    404. 

  404.   real :: colmass_end( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
    405. 

  405.   real :: colmass_sfc( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
    406. 

  406.   real :: colmass_worst( 0:maxd_acomp, maxd_asize, maxd_atype )
    407. 

  407.   real :: colmass_maxworst_r
    408. 

  408.   real :: rhodz( kts:kte ), rhodzsum
    409. 

  409. 

  410. !!$#if (defined AIX)
    411. 

  411. !!$#define ERF erf
    412. 

  412. !!$#define ERFC erfc
    413. 

  413. !!$#else
    414. 

  414. !!$#define ERF erf
    415. 

  415. !!$    real erf
    416. 

  416. !!$#define ERFC erfc
    417. 

  417. !!$    real erfc
    418. 

  418. !!$#endif
    419. 

  419. 

  420.   character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
    421. 

  421. 

  422. 

  423.   colmass_worst(:,:,:) = 0.0
    424. 

  424.   colmass_worst_ij(:,:,:,:) = -1
    425. 

  425. 

  426. 

  427.   arg = 1.0
    428. 

  428.   if (abs(0.8427-ERF_ALT(arg))/0.8427>0.001) then
    429. 

  429.      write (6,*) 'erf_alt(1.0) = ',ERF_ALT(arg)
    430. 

  430.      call wrf_error_fatal('dropmixnuc: Error function error')
    431. 

  431.   endif
    432. 

  432.   arg = 0.0
    433. 

  433.   if (ERF_ALT(arg) /= 0.0) then
    434. 

  434.      write (6,*) 'erf_alt(0.0) = ',ERF_ALT(arg)
    435. 

  435.      call wrf_error_fatal('dropmixnuc: Error function error')
    436. 

  436.   endif
    437. 

  437. 

  438.   pi = 4.*atan(1.0)
    439. 

  439.   dtinv=1./dtstep
    440. 

  440. 

  441.   depvel_drop =  0.1 ! prescribed here rather than getting it from dry_dep_driver
    442. 

  442.   if (idrydep_onoff .le. 0) depvel_drop =  0.0
    443. 

  443.   depvel_drop =  min(depvel_drop,depvel_uplimit)
    444. 

  444. 

  445.   do n=1,ntype_aer
    446. 

  446.      do m=1,nsize_aer(n)
    447. 

  447.         ncomp(n)=ncomp_aer(n)
    448. 

  448.         alogsig(m,n)=alog(sigmag_aer(m,n))
    449. 

  449.         dgnum_aer(m,n) = dpvolmean_aer(m,n) * exp( -1.5*alogsig(m,n)*alogsig(m,n) )
    450. 

  450. !    print *,'sigmag_aer,dgnum_aer=',sigmag_aer(m,n),dgnum_aer(m,n)
    451. 

  451.         ! npv is used only if number is diagnosed from volume
    452. 

  452.         npv(m,n)=6./(pi*(0.01*dgnum_aer(m,n))**3*exp(4.5*alogsig(m,n)*alogsig(m,n)))
    453. 

  453.      end do
    454. 

  454.   end do
    455. 

  455.   t0=273.15   !wig, 1-Mar-2009: Added .15
    456. 

  456.   aten=2.*surften/(r_v*t0*rhowater)
    457. 

  457.   super(:)=0.01*supersat(:)
    458. 

  458.   do n=1,ntype_aer
    459. 

  459.      do m=1,nsize_aer(n)
    460. 

  460.         exp45logsig=exp(4.5*alogsig(m,n)*alogsig(m,n))
    461. 

  461.         argfactor(m,n)=2./(3.*sqrt(2.)*alogsig(m,n))
    462. 

  462.         amcube(m,n)=3./(4.*pi*exp45logsig*npv(m,n))
    463. 

  463.      enddo
    464. 

  464.   enddo
    465. 

  465. 

  466.   IF( PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
    467. 

  467.      CALL cal_cldfra(CLDFRA,qc,qi,f_qc,f_qi,      &
    468. 

  468.           ids,ide, jds,jde, kds,kde,              &
    469. 

  469.           ims,ime, jms,jme, kms,kme,              &
    470. 

  470.           its,ite, jts,jte, kts,kte               )
    471. 

  471.   END IF
    472. 

  472. 

  473.   qsrflx(its:ite,jts:jte,:) = 0.
    474. 

  474. 

  475. !     start loop over columns
    476. 

  476. 

  477. OVERALL_MAIN_J_LOOP: do j=jts,jte
    478. 

  478. OVERALL_MAIN_I_LOOP: do i=its,ite
    479. 

  479. 

  480. !        load number nucleated into qndrop on cloud boundaries
    481. 

  481. 

  482. ! initialization for current i .........................................
    483. 

  483. 

  484.      do k=kts+1,kte
    485. 

  485. 	    zs(k)=1./(zm(i,k,j)-zm(i,k-1,j))
    486. 

  486. 	 enddo
    487. 

  487. 	 zs(kts)=zs(kts+1)
    488. 

  488.      zs(kte+1)=0.
    489. 

  489. 

  490.      do k=kts,kte
    491. 

  491. !!$	    if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
    492. 

  492. !!$!	       call wrf_error_fatal("1")
    493. 

  493. !!$	    endif
    494. 

  494.         if(f_qi)then
    495. 

  495.            qcld=qc(i,k,j)+qi(i,k,j)
    496. 

  496.         else
    497. 

  497.            qcld=qc(i,k,j)
    498. 

  498.         endif
    499. 

  499.         if(qcld.lt.-1..or.qcld.gt.1.)then
    500. 

  500.            write(6,'(a,g12.2,a,3i5)')'qcld=',qcld,' for i,k,j=',i,k,j
    501. 

  501.            call wrf_error_fatal("1")
    502. 

  502.         endif
    503. 

  503.         if(qcld.gt.1.e-20)then
    504. 

  504.            lcldfra(k)=cldfra(i,k,j)*qc(i,k,j)/qcld
    505. 

  505.            lcldfra_old(k)=cldfra_old(i,k,j)*qc(i,k,j)/qcld
    506. 

  506.         else
    507. 

  507.            lcldfra(k)=0.
    508. 

  508.            lcldfra_old(k)=0.
    509. 

  509.         endif
    510. 

  510.         qndrop(k)=qndrop3d(i,k,j)
    511. 

  511. !	    qndrop(k)=1.e5
    512. 

  512.         cs(k)=rho(i,k,j) ! air density (kg/m3)
    513. 

  513.         dz(k)=dz8w(i,k,j)
    514. 

  514.         do n=1,ntype_aer
    515. 

  515.            do m=1,nsize_aer(n)
    516. 

  516.               nact(k,m,n)=0.
    517. 

  517.               mact(k,m,n)=0.
    518. 

  518.            enddo
    519. 

  519.         enddo
    520. 

  520.         zn(k)=1./(cs(k)*dz(k))
    521. 

  521.         if(k>kts)then
    522. 

  522.            ekd(k)=kvh(i,k,j)
    523. 

  523.            ekd(k)=max(ekd(k),zkmin)
    524. 

  524.            ekd(k)=min(ekd(k),zkmax)
    525. 

  525.         else
    526. 

  526.            ekd(k)=0
    527. 

  527.         endif
    528. 

  528. !           diagnose subgrid vertical velocity from diffusivity
    529. 

  529.         if(k.eq.kts)then
    530. 

  530.            wtke(k)=sq2pi*depvel_drop
    531. 

  531. !               wtke(k)=sq2pi*kvh(i,k,j)
    532. 

  532. !               wtke(k)=max(wtke(k),wmixmin)
    533. 

  533.         else
    534. 

  534.            wtke(k)=sq2pi*ekd(k)/dz(k)
    535. 

  535.         endif
    536. 

  536.         wtke(k)=max(wtke(k),wmixmin)
    537. 

  537.         nsource(i,k,j)=0.
    538. 

  538.      enddo
    539. 

  539.      nsource(i,kte+1,j) = 0.
    540. 

  540.      qndrop(kte+1)      = 0.
    541. 

  541.      zn(kte+1)          = 0.
    542. 

  542. 

  543.      do k = kts+1, kte
    544. 

  544.         tmpa = dz(k-1) ; tmpb = dz(k)
    545. 

  545.         tmpc = tmpa/(tmpa + tmpb)
    546. 

  546.         csbot(k) = cs(k-1)*(1.0-tmpc) + cs(k)*tmpc
    547. 

  547.         csbot_cscen(k) = csbot(k)/cs(k)
    548. 

  548.      end do
    549. 

  549.      csbot(kts) = cs(kts)
    550. 

  550.      csbot_cscen(kts) = 1.0
    551. 

  551.      csbot(kte+1) = cs(kte)
    552. 

  552.      csbot_cscen(kte+1) = 1.0
    553. 

  553. 

  554.     !  calculate surface rate and mass mixing ratio for aerosol
    555. 

  555. 

  556.      surfratemax = 0.0
    557. 

  557.      nsav=1
    558. 

  558.      nnew=2
    559. 

  559.      surfrate_drop=depvel_drop/dz(kts)
    560. 

  560.      surfratemax = max( surfratemax, surfrate_drop )
    561. 

  561.      do n=1,ntype_aer
    562. 

  562.         do m=1,nsize_aer(n)
    563. 

  563.            lnum=numptr_aer(m,n,ai_phase)
    564. 

  564.            lnumcw=numptr_aer(m,n,cw_phase)
    565. 

  565.            if(lnum>0)then
    566. 

  566.               depvel_tmp = max( 0.0, min( ddvel(i,j,lnum), depvel_uplimit ) )
    567. 

  567.               surfrate(lnum)=depvel_tmp/dz(kts)
    568. 

  568.               surfrate(lnumcw)=surfrate_drop
    569. 

  569.               surfratemax = max( surfratemax, surfrate(lnum) )
    570. 

  570. !             scale = 1000./mwdry ! moles/kg
    571. 

  571.               scale = 1.
    572. 

  572.               raercol(kts:kte,lnumcw,nsav)=chem(i,kts:kte,j,lnumcw)*scale ! #/kg
    573. 

  573.               raercol(kts:kte,lnum,nsav)=chem(i,kts:kte,j,lnum)*scale
    574. 

  574.            endif
    575. 

  575.            do l=1,ncomp(n)
    576. 

  576.               lmass=massptr_aer(l,m,n,ai_phase)
    577. 

  577.               lmasscw=massptr_aer(l,m,n,cw_phase)
    578. 

  578. !             scale = mw_aer(l,n)/mwdry
    579. 

  579.               scale = 1.e-9 ! kg/ug
    580. 

  580.               depvel_tmp = max( 0.0, min( ddvel(i,j,lmass), depvel_uplimit ) )
    581. 

  581.               surfrate(lmass)=depvel_tmp/dz(kts)
    582. 

  582.               surfrate(lmasscw)=surfrate_drop
    583. 

  583.               surfratemax = max( surfratemax, surfrate(lmass) )
    584. 

  584.               raercol(kts:kte,lmasscw,nsav)=chem(i,kts:kte,j,lmasscw)*scale ! kg/kg
    585. 

  585.               raercol(kts:kte,lmass,nsav)=chem(i,kts:kte,j,lmass)*scale ! kg/kg
    586. 

  586.            enddo
    587. 

  587.            lwater=waterptr_aer(m,n)
    588. 

  588.            if(lwater>0)then
    589. 

  589.               depvel_tmp = max( 0.0, min( ddvel(i,j,lwater), depvel_uplimit ) )
    590. 

  590.               surfrate(lwater)=depvel_tmp/dz(kts)
    591. 

  591.               surfratemax = max( surfratemax, surfrate(lwater) )
    592. 

  592.               raercol(kts:kte,lwater,nsav)=chem(i,kts:kte,j,lwater) ! don't bother to convert units,
    593. 

  593.              ! because it doesn't contribute to aerosol mass
    594. 

  594.            endif
    595. 

  595.         enddo ! size
    596. 

  596.      enddo ! type
    597. 

  597. 

  598. 

  599. ! mass conservation checking
    600. 

  600.      if (icheck_colmass > 0) then
    601. 

  601. 

  602. ! calc initial column burdens
    603. 

  603.         colmass_bgn(:,:,:,:) = 0.0
    604. 

  604.         colmass_end(:,:,:,:) = 0.0
    605. 

  605.         colmass_sfc(:,:,:,:) = 0.0
    606. 

  606.         rhodz(kts:kte) = 1.0/zn(kts:kte)
    607. 

  607.         rhodzsum = sum( rhodz(kts:kte) )
    608. 

  608.         do n=1,ntype_aer
    609. 

  609.            do m=1,nsize_aer(n)
    610. 

  610.               lnum=numptr_aer(m,n,ai_phase)
    611. 

  611.               lnumcw=numptr_aer(m,n,cw_phase)
    612. 

  612.               if(lnum>0)then
    613. 

  613.                  colmass_bgn(0,m,n,1) = sum( chem(i,kts:kte,j,lnum  )*rhodz(kts:kte) )
    614. 

  614.                  colmass_bgn(0,m,n,2) = sum( chem(i,kts:kte,j,lnumcw)*rhodz(kts:kte) )
    615. 

  615.               endif
    616. 

  616.               do l=1,ncomp(n)
    617. 

  617.                  lmass=massptr_aer(l,m,n,ai_phase)
    618. 

  618.                  lmasscw=massptr_aer(l,m,n,cw_phase)
    619. 

  619.                  colmass_bgn(l,m,n,1) = sum( chem(i,kts:kte,j,lmass  )*rhodz(kts:kte) )
    620. 

  620.                  colmass_bgn(l,m,n,2) = sum( chem(i,kts:kte,j,lmasscw)*rhodz(kts:kte) )
    621. 

  621.               enddo
    622. 

  622.            enddo ! size
    623. 

  623.         enddo ! type
    624. 

  624.      endif ! (icheck_colmass > 0)
    625. 

  625. 

  626. 

  627. !        droplet nucleation/aerosol activation
    628. 

  628. 

  629. ! k-loop for growing/shrinking cloud calcs .............................
    630. 

  630. GROW_SHRINK_MAIN_K_LOOP: do k=kts,kte
    631. 

  631.         km1=max0(k-1,1)
    632. 

  632.         kp1=min0(k+1,kde-1)
    633. 

  633. 

  634. 

  635. !       if(lcldfra(k)-lcldfra_old(k).gt.0.01)then   ! this line is the "old" criterion
    636. 

  636. !       go to 10
    637. 

  637. 

  638. !       growing cloud PLUS
    639. 

  639. !       upwards vertical advection when lcldfra(k-1) < lcldfra(k)
    640. 

  640. !
    641. 

  641. ! tmpc1 = cloud fraction increase from previous time step
    642. 

  642.         tmpc1 = max( (lcldfra(k)-lcldfra_old(k)), 0.0 )
    643. 

  643.         if (k > kts) then
    644. 

  644. ! tmpc2 = fraction of layer for which vertical advection from below
    645. 

  645. !             (over dtstep) displaces cloudy air with clear air
    646. 

  646. !       = (courant number using upwards w at layer bottom)*(difference in cloud fraction)
    647. 

  647.            tmpcourno = dtstep*max(w(i,k,j),0.0)/dz(k)
    648. 

  648.            tmpc2 = max( (lcldfra(k)-lcldfra(km1)), 0.0 ) * tmpcourno
    649. 

  649.            tmpc2 = min( tmpc2, 1.0 )
    650. 

  650. !          tmpc2 = 0.0   ! this turns off the vertical advect part
    651. 

  651.         else
    652. 

  652.            tmpc2 = 0.0
    653. 

  653.         endif
    654. 

  654. 

  655.         if ((tmpc1 > 0.001) .or. (tmpc2 > 0.001)) then
    656. 

  656. 

  657. !                wmix=wtke(k)
    658. 

  658.            wbar=w(i,k,j)+wtke(k)
    659. 

  659.            wmix=0.
    660. 

  660.            wmin=0.
    661. 

  661. ! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
    662. 

  662.            wmax=50.
    663. 

  663.            wdiab=0
    664. 

  664. 

  665. !                load aerosol properties, assuming external mixtures
    666. 

  666.            do n=1,ntype_aer
    667. 

  667.               do m=1,nsize_aer(n)
    668. 

  668.                  call loadaer(raercol(1,1,nsav),k,kms,kme,num_chem,    &
    669. 

  669.                       cs(k), npv(m,n), dlo_sect(m,n),dhi_sect(m,n),             &
    670. 

  670.                       maxd_acomp, ncomp(n), &
    671. 

  671.                       grid_id, ktau, i, j, m, n,   &
    672. 

  672.                       numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase),  &
    673. 

  673.                       dens_aer(1,n),    &
    674. 

  674.                       massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase),  &
    675. 

  675.                       maerosol(1,m,n), maerosolcw(1,m,n),          &
    676. 

  676.                       maerosol_tot(m,n), maerosol_totcw(m,n),      &
    677. 

  677.                       naerosol(m,n), naerosolcw(m,n),                  &
    678. 

  678.                       vaerosol(m,n), vaerosolcw(m,n) )
    679. 

  679. 

  680.                  hygro_aer(m,n)=hygro(i,k,j,m,n)
    681. 

  681.               enddo
    682. 

  682.            enddo
    683. 

  683. 

  684. ! 06-nov-2005 rce - grid_id & ktau added to arg list
    685. 

  685.            call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
    686. 

  686.                 msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
    687. 

  687.                 naerosol, vaerosol,  &
    688. 

  688.                 dlo_sect,dhi_sect,sigmag_aer,hygro_aer,              &
    689. 

  689.                 fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k )
    690. 

  690. 

  691.            do n = 1,ntype_aer
    692. 

  692.               do m = 1,nsize_aer(n)
    693. 

  693.                  lnum   = numptr_aer(m,n,ai_phase)
    694. 

  694.                  lnumcw = numptr_aer(m,n,cw_phase)
    695. 

  695.                  if (tmpc1 > 0.0) then
    696. 

  696.                     dact = tmpc1*fn(m,n)*raercol(k,lnum,nsav) ! interstitial only
    697. 

  697.                  else
    698. 

  698.                     dact = 0.0
    699. 

  699.                  endif
    700. 

  700.                  if (tmpc2 > 0.0) then
    701. 

  701.                     dact = dact + tmpc2*fn(m,n)*raercol(km1,lnum,nsav) ! interstitial only
    702. 

  702.                  endif
    703. 

  703.                  dact = min( dact, 0.99*raercol(k,lnum,nsav) )
    704. 

  704.                  raercol(k,lnumcw,nsav) = raercol(k,lnumcw,nsav)+dact
    705. 

  705.                  raercol(k,lnum,  nsav) = raercol(k,lnum,  nsav)-dact
    706. 

  706.                  qndrop(k) = qndrop(k)+dact
    707. 

  707.                  nsource(i,k,j) = nsource(i,k,j)+dact*dtinv
    708. 

  708.                  do l = 1,ncomp(n)
    709. 

  709.                     lmass   = massptr_aer(l,m,n,ai_phase)
    710. 

  710.                     lmasscw = massptr_aer(l,m,n,cw_phase)
    711. 

  711.                     if (tmpc1 > 0.0) then
    712. 

  712.                        dact = tmpc1*fm(m,n)*raercol(k,lmass,nsav) ! interstitial only
    713. 

  713.                     else
    714. 

  714.                        dact = 0.0
    715. 

  715.                     endif
    716. 

  716.                     if (tmpc2 > 0.0) then
    717. 

  717.                        dact = dact + tmpc2*fm(m,n)*raercol(km1,lmass,nsav) ! interstitial only
    718. 

  718.                     endif
    719. 

  719.                     dact = min( dact, 0.99*raercol(k,lmass,nsav) )
    720. 

  720.                     raercol(k,lmasscw,nsav)  =  raercol(k,lmasscw,nsav)+dact
    721. 

  721.                     raercol(k,lmass,  nsav)  =  raercol(k,lmass,  nsav)-dact
    722. 

  722.                  enddo
    723. 

  723.               enddo
    724. 

  724.            enddo
    725. 

  725. !   10 continue
    726. 

  726.         endif   ! ((tmpc1 > 0.001) .or. (tmpc2 > 0.001))
    727. 

  727. 

  728. 

  729.         if(lcldfra(k) < lcldfra_old(k) .and. lcldfra_old(k) > 1.e-20)then   ! this line is the "old" criterion
    730. 

  730. !         go to 20
    731. 

  731. 

  732. !       shrinking cloud ......................................................
    733. 

  733. 

  734. !                droplet loss in decaying cloud
    735. 

  735.            nsource(i,k,j)=nsource(i,k,j)+qndrop(k)*(lcldfra(k)-lcldfra_old(k))*dtinv
    736. 

  736.            qndrop(k)=qndrop(k)*(1.+lcldfra(k)-lcldfra_old(k))
    737. 

  737. !                 convert activated aerosol to interstitial in decaying cloud
    738. 

  738. 

  739.            tmpc = (lcldfra(k)-lcldfra_old(k))/lcldfra_old(k)
    740. 

  740.            do n=1,ntype_aer
    741. 

  741.               do m=1,nsize_aer(n)
    742. 

  742.                  lnum=numptr_aer(m,n,ai_phase)
    743. 

  743.                  lnumcw=numptr_aer(m,n,cw_phase)
    744. 

  744.                  if(lnum.gt.0)then
    745. 

  745.                     dact=raercol(k,lnumcw,nsav)*tmpc
    746. 

  746.                     raercol(k,lnumcw,nsav)=raercol(k,lnumcw,nsav)+dact
    747. 

  747.                     raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
    748. 

  748.                  endif
    749. 

  749.                  do l=1,ncomp(n)
    750. 

  750.                     lmass=massptr_aer(l,m,n,ai_phase)
    751. 

  751.                     lmasscw=massptr_aer(l,m,n,cw_phase)
    752. 

  752.                     dact=raercol(k,lmasscw,nsav)*tmpc
    753. 

  753.                     raercol(k,lmasscw,nsav)=raercol(k,lmasscw,nsav)+dact
    754. 

  754.                     raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
    755. 

  755.                  enddo
    756. 

  756.               enddo
    757. 

  757.            enddo
    758. 

  758. !             20 continue
    759. 

  759.         endif
    760. 

  760. 

  761.      enddo GROW_SHRINK_MAIN_K_LOOP
    762. 

  762. ! end of k-loop for growing/shrinking cloud calcs ......................
    763. 

  763. 

  764. 

  765. 

  766. ! ......................................................................
    767. 

  767. ! start of main k-loop for calc of old cloud activation tendencies ..........
    768. 

  768. ! this loop does "set up" for the nsubmix loop
    769. 

  769. !
    770. 

  770. ! rce-comment
    771. 

  771. !    changed this part of code to use current cloud fraction (lcldfra) exclusively
    772. 

  772. 

  773. OLD_CLOUD_MAIN_K_LOOP: do k=kts,kte
    774. 

  774.         km1=max0(k-1,kts)
    775. 

  775.         kp1=min0(k+1,kde-1)
    776. 

  776.         flux_fullact(k) = 0.0
    777. 

  777.         if(lcldfra(k).gt.0.01)then
    778. 

  778. !          go to 30
    779. 

  779. 

  780. !               old cloud
    781. 

  781.               if(lcldfra(k)-lcldfra(km1).gt.0.01.or.k.eq.kts)then
    782. 

  782. 

  783. !                   interior cloud
    784. 

  784. !                   cloud base
    785. 

  785. 

  786.                  wdiab=0
    787. 

  787.                  wmix=wtke(k) ! spectrum of updrafts
    788. 

  788.                  wbar=w(i,k,j) ! spectrum of updrafts
    789. 

  789. !                    wmix=0. ! single updraft
    790. 

  790. !               wbar=wtke(k) ! single updraft
    791. 

  791. ! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
    792. 

  792.                  wmax=50.
    793. 

  793.                  ekd(k)=wtke(k)*dz(k)/sq2pi
    794. 

  794.                  alogarg=max(1.e-20,1/lcldfra(k)-1.)
    795. 

  795.                  wmin=wbar+wmix*0.25*sq2pi*alog(alogarg)
    796. 

  796. 

  797.                  do n=1,ntype_aer
    798. 

  798.                     do m=1,nsize_aer(n)
    799. 

  799.                        call loadaer(raercol(1,1,nsav),km1,kms,kme,num_chem,    &
    800. 

  800.                             cs(k), npv(m,n),dlo_sect(m,n),dhi_sect(m,n),               &
    801. 

  801.                             maxd_acomp, ncomp(n), &
    802. 

  802.                             grid_id, ktau, i, j, m, n,   &
    803. 

  803.                             numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase),  &
    804. 

  804.                             dens_aer(1,n),   &
    805. 

  805.                             massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase),  &
    806. 

  806.                             maerosol(1,m,n), maerosolcw(1,m,n),          &
    807. 

  807.                             maerosol_tot(m,n), maerosol_totcw(m,n),      &
    808. 

  808.                             naerosol(m,n), naerosolcw(m,n),                  &
    809. 

  809.                             vaerosol(m,n), vaerosolcw(m,n) )
    810. 

  810. 

  811.                        hygro_aer(m,n)=hygro(i,k,j,m,n)
    812. 

  812. 

  813.                     enddo
    814. 

  814.                  enddo
    815. 

  815. !          print *,'old cloud wbar,wmix=',wbar,wmix
    816. 

  816. 

  817.                  call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
    818. 

  818.                       msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
    819. 

  819.                       naerosol, vaerosol,  &
    820. 

  820.                       dlo_sect,dhi_sect, sigmag_aer,hygro_aer,                    &
    821. 

  821.                       fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k )
    822. 

  822.                  
    823. 

  823. ! rce-comment
    824. 

  824. !    the activation-fraction fluxes (fluxn, fluxm) from subr activate assume that
    825. 

  825. !       wbar << wmix, which is valid for global-model scale but not mesoscale
    826. 

  826. !    for wrf-chem application, divide these by flux_fullact to get a 
    827. 

  827. !       "flux-weighted-average" activation fraction, then multiply by (ekd(k)*zs(k)) 
    828. 

  828. !       which is the local "turbulent vertical-mixing velocity"
    829. 

  829.                  if (k > kts) then
    830. 

  830.                     if (flux_fullact(k) > 1.0e-20) then
    831. 

  831.                        tmpa = ekd(k)*zs(k)
    832. 

  832.                        tmpf = flux_fullact(k)
    833. 

  833.                        do n=1,ntype_aer
    834. 

  834.                        do m=1,nsize_aer(n)
    835. 

  835.                           tmpb = max( fluxn(m,n), 0.0 ) / max( fluxn(m,n), tmpf )
    836. 

  836.                           fluxn(m,n) = tmpa*tmpb
    837. 

  837.                           tmpb = max( fluxm(m,n), 0.0 ) / max( fluxm(m,n), tmpf )
    838. 

  838.                           fluxm(m,n) = tmpa*tmpb
    839. 

  839.                        enddo
    840. 

  840.                        enddo
    841. 

  841.                     else
    842. 

  842.                        fluxn(:,:) = 0.0
    843. 

  843.                        fluxm(:,:) = 0.0
    844. 

  844.                     endif
    845. 

  845.                  endif
    846. 

  846. 

  847.                  if(k.gt.kts)then
    848. 

  848.                     tmpc = lcldfra(k)-lcldfra(km1)
    849. 

  849.                  else
    850. 

  850.                     tmpc=lcldfra(k)
    851. 

  851.                  endif
    852. 

  852. ! rce-comment
    853. 

  853. !    flux of activated mass into layer k (in kg/m2/s)
    854. 

  854. !       = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
    855. 

  855. !    source of activated mass (in kg/kg/s) = flux divergence
    856. 

  856. !       = actmassflux/(cs(i,k)*dz(i,k))
    857. 

  857. !    so need factor of csbot_cscen = csbot(k)/cs(i,k)
    858. 

  858. !                tmpe=1./(dz(k))
    859. 

  859.                  tmpe = csbot_cscen(k)/(dz(k))
    860. 

  860.                  fluxntot=0.
    861. 

  861.                  do n=1,ntype_aer
    862. 

  862.                  do m=1,nsize_aer(n)
    863. 

  863.                     fluxn(m,n)=fluxn(m,n)*tmpc
    864. 

  864. !                   fluxs(m,n)=fluxs(m,n)*tmpc
    865. 

  865.                     fluxm(m,n)=fluxm(m,n)*tmpc
    866. 

  866.                     lnum=numptr_aer(m,n,ai_phase)
    867. 

  867.                     fluxntot=fluxntot+fluxn(m,n)*raercol(km1,lnum,nsav)
    868. 

  868. !             print *,'fn=',fn(m,n),' for m,n=',m,n
    869. 

  869. !             print *,'old cloud tmpc=',tmpc,' fn=',fn(m,n),' for m,n=',m,n
    870. 

  870.                     nact(k,m,n)=nact(k,m,n)+fluxn(m,n)*tmpe
    871. 

  871.                     mact(k,m,n)=mact(k,m,n)+fluxm(m,n)*tmpe
    872. 

  872.                  enddo
    873. 

  873.                  enddo
    874. 

  874.                  flux_fullact(k) = flux_fullact(k)*tmpe
    875. 

  875.                  nsource(i,k,j)=nsource(i,k,j)+fluxntot*zs(k)
    876. 

  876.                  fluxntot=fluxntot*cs(k)
    877. 

  877.               endif
    878. 

  878. !       30 continue
    879. 

  879. 

  880.         else
    881. 

  881. !       go to 40
    882. 

  882. !              no cloud
    883. 

  883.            if(qndrop(k).gt.10000.e6)then
    884. 

  884.               print *,'i,k,j,lcldfra,qndrop=',i,k,j,lcldfra(k),qndrop(k)
    885. 

  885.               print *,'cldfra,ql,qi',cldfra(i,k,j),qc(i,k,j),qi(i,k,j)
    886. 

  886.            endif
    887. 

  887.            nsource(i,k,j)=nsource(i,k,j)-qndrop(k)*dtinv
    888. 

  888.            qndrop(k)=0.
    889. 

  889. !              convert activated aerosol to interstitial in decaying cloud
    890. 

  890.            do n=1,ntype_aer
    891. 

  891.               do m=1,nsize_aer(n)
    892. 

  892.                  lnum=numptr_aer(m,n,ai_phase)
    893. 

  893.                  lnumcw=numptr_aer(m,n,cw_phase)
    894. 

  894.                  if(lnum.gt.0)then
    895. 

  895.                     raercol(k,lnum,nsav)=raercol(k,lnum,nsav)+raercol(k,lnumcw,nsav)
    896. 

  896.                     raercol(k,lnumcw,nsav)=0.
    897. 

  897.                  endif
    898. 

  898.                  do l=1,ncomp(n)
    899. 

  899.                     lmass=massptr_aer(l,m,n,ai_phase)
    900. 

  900.                     lmasscw=massptr_aer(l,m,n,cw_phase)
    901. 

  901.                     raercol(k,lmass,nsav)=raercol(k,lmass,nsav)+raercol(k,lmasscw,nsav)
    902. 

  902.                     raercol(k,lmasscw,nsav)=0.
    903. 

  903.                  enddo
    904. 

  904.               enddo
    905. 

  905.            enddo
    906. 

  906. !      40 continue
    907. 

  907.         endif
    908. 

  908. 

  909.      enddo OLD_CLOUD_MAIN_K_LOOP
    910. 

  910. 

  911. !    cycle OVERALL_MAIN_I_LOOP
    912. 

  912. 

  913. 

  914. !    switch nsav, nnew so that nnew is the updated aerosol
    915. 

  915. 

  916.      ntemp=nsav
    917. 

  917.      nsav=nnew
    918. 

  918.      nnew=ntemp
    919. 

  919. 

  920. !    load new droplets in layers above, below clouds
    921. 

  921. 

  922.      dtmin=dtstep
    923. 

  923.      ekk(kts)=0.0
    924. 

  924. ! rce-comment -- ekd(k) is eddy-diffusivity at k/k-1 interface
    925. 

  925. !   want ekk(k) = ekd(k) * (density at k/k-1 interface)
    926. 

  926.      do k=kts+1,kte
    927. 

  927.         ekk(k)=ekd(k)*csbot(k)
    928. 

  928.      enddo
    929. 

  929.      ekk(kte+1)=0.0
    930. 

  930.      do k=kts,kte
    931. 

  931.         ekkp(k)=zn(k)*ekk(k+1)*zs(k+1)
    932. 

  932.         ekkm(k)=zn(k)*ekk(k)*zs(k)
    933. 

  933.         tinv=ekkp(k)+ekkm(k)
    934. 

  934.         if(k.eq.kts)tinv=tinv+surfratemax
    935. 

  935.         if(tinv.gt.1.e-6)then
    936. 

  936.            dtt=1./tinv
    937. 

  937.            dtmin=min(dtmin,dtt)
    938. 

  938.         endif
    939. 

  939.      enddo
    940. 

  940.      dtmix=0.9*dtmin
    941. 

  941.      nsubmix=dtstep/dtmix+1
    942. 

  942.      if(nsubmix>100)then
    943. 

  943.         nsubmix_bnd=100
    944. 

  944.      else
    945. 

  945.         nsubmix_bnd=nsubmix
    946. 

  946.      endif
    947. 

  947. !     count_submix(nsubmix_bnd)=count_submix(nsubmix_bnd)+1
    948. 

  948.      dtmix=dtstep/nsubmix
    949. 

  949.      fac_srflx = -1.0/(zn(1)*nsubmix)
    950. 

  950.      
    951. 

  951.      do k=kts,kte
    952. 

  952.         kp1=min(k+1,kde-1)
    953. 

  953.         km1=max(k-1,1)
    954. 

  954.         if(lcldfra(kp1).gt.0)then
    955. 

  955.            overlapp(k)=min(lcldfra(k)/lcldfra(kp1),1.)
    956. 

  956.         else
    957. 

  957.            overlapp(k)=1.
    958. 

  958.         endif
    959. 

  959.         if(lcldfra(km1).gt.0)then
    960. 

  960.            overlapm(k)=min(lcldfra(k)/lcldfra(km1),1.)
    961. 

  961.         else
    962. 

  962.            overlapm(k)=1.
    963. 

  963.         endif
    964. 

  964.      enddo
    965. 

  965. 

  966. 

  967. 

  968. ! ......................................................................
    969. 

  969. ! start of nsubmix-loop for calc of old cloud activation tendencies ....
    970. 

  970. OLD_CLOUD_NSUBMIX_LOOP: do nsub=1,nsubmix
    971. 

  971.         qndrop_new(kts:kte)=qndrop(kts:kte)
    972. 

  972. !           switch nsav, nnew so that nsav is the updated aerosol
    973. 

  973.         ntemp=nsav
    974. 

  974.         nsav=nnew
    975. 

  975.         nnew=ntemp
    976. 

  976.         srcn(:)=0.0
    977. 

  977.         do n=1,ntype_aer
    978. 

  978.            do m=1,nsize_aer(n)
    979. 

  979.               lnum=numptr_aer(m,n,ai_phase)
    980. 

  980. !              update droplet source
    981. 

  981. ! rce-comment - activation source in layer k involves particles from k-1
    982. 

  982. !             srcn(kts  :kte)=srcn(kts  :kte)+nact(kts  :kte,m,n)*(raercol(kts:kte  ,lnum,nsav))
    983. 

  983.               srcn(kts+1:kte)=srcn(kts+1:kte)+nact(kts+1:kte,m,n)*(raercol(kts:kte-1,lnum,nsav))
    984. 

  984. ! rce-comment - new formulation for k=kts should be implemented
    985. 

  985.               srcn(kts      )=srcn(kts      )+nact(kts      ,m,n)*(raercol(kts      ,lnum,nsav))
    986. 

  986.            enddo
    987. 

  987.         enddo
    988. 

  988.         call explmix(qndrop,srcn,ekkp,ekkm,overlapp,overlapm,   &
    989. 

  989.              qndrop_new,surfrate_drop,kms,kme,kts,kte,dtmix,.false.)
    990. 

  990.         do n=1,ntype_aer
    991. 

  991.            do m=1,nsize_aer(n)
    992. 

  992.               lnum=numptr_aer(m,n,ai_phase)
    993. 

  993.               lnumcw=numptr_aer(m,n,cw_phase)
    994. 

  994.               if(lnum>0)then
    995. 

  995. ! rce-comment - activation source in layer k involves particles from k-1
    996. 

  996. !                source(kts  :kte)= nact(kts  :kte,m,n)*(raercol(kts:kte  ,lnum,nsav))
    997. 

  997.                  source(kts+1:kte)= nact(kts+1:kte,m,n)*(raercol(kts:kte-1,lnum,nsav))
    998. 

  998. ! rce-comment - new formulation for k=kts should be implemented
    999. 

  999.                  source(kts      )= nact(kts      ,m,n)*(raercol(kts      ,lnum,nsav))
    1000. 

  1000.                  call explmix(raercol(1,lnumcw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
    1001. 

  1001.                       raercol(1,lnumcw,nsav),surfrate(lnumcw),kms,kme,kts,kte,dtmix,&
    1002. 

  1002.                       .false.)
    1003. 

  1003.                  call explmix(raercol(1,lnum,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
    1004. 

  1004.                       raercol(1,lnum,nsav),surfrate(lnum),kms,kme,kts,kte,dtmix, &
    1005. 

  1005.                       .true.,raercol(1,lnumcw,nsav))
    1006. 

  1006.                  qsrflx(i,j,lnum) = qsrflx(i,j,lnum) + fac_srflx*            &
    1007. 

  1007.                       raercol(kts,lnum,nsav)*surfrate(lnum)
    1008. 

  1008.                  qsrflx(i,j,lnumcw) = qsrflx(i,j,lnumcw) + fac_srflx*        &
    1009. 

  1009.                       raercol(kts,lnumcw,nsav)*surfrate(lnumcw)
    1010. 

  1010.                  if (icheck_colmass > 0) then
    1011. 

  1011.                     tmpf = dtmix*rhodz(kts)
    1012. 

  1012.                     colmass_sfc(0,m,n,1) = colmass_sfc(0,m,n,1) &
    1013. 

  1013.                           + raercol(kts,lnum  ,nsav)*surfrate(lnum  )*tmpf
    1014. 

  1014.                     colmass_sfc(0,m,n,2) = colmass_sfc(0,m,n,2) &
    1015. 

  1015.                           + raercol(kts,lnumcw,nsav)*surfrate(lnumcw)*tmpf
    1016. 

  1016.                  endif
    1017. 

  1017.               endif
    1018. 

  1018.               do l=1,ncomp(n)
    1019. 

  1019.                  lmass=massptr_aer(l,m,n,ai_phase)
    1020. 

  1020.                  lmasscw=massptr_aer(l,m,n,cw_phase)
    1021. 

  1021. ! rce-comment - activation source in layer k involves particles from k-1
    1022. 

  1022. !                source(kts  :kte)= mact(kts  :kte,m,n)*(raercol(kts:kte  ,lmass,nsav))
    1023. 

  1023.                  source(kts+1:kte)= mact(kts+1:kte,m,n)*(raercol(kts:kte-1,lmass,nsav))
    1024. 

  1024. ! rce-comment - new formulation for k=kts should be implemented
    1025. 

  1025.                  source(kts      )= mact(kts      ,m,n)*(raercol(kts      ,lmass,nsav))
    1026. 

  1026.                  call explmix(raercol(1,lmasscw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
    1027. 

  1027.                       raercol(1,lmasscw,nsav),surfrate(lmasscw),kms,kme,kts,kte,dtmix,  &
    1028. 

  1028.                       .false.)
    1029. 

  1029.                  call explmix(raercol(1,lmass,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
    1030. 

  1030.                       raercol(1,lmass,nsav),surfrate(lmass),kms,kme,kts,kte,dtmix,  &
    1031. 

  1031.                       .true.,raercol(1,lmasscw,nsav))
    1032. 

  1032.                  qsrflx(i,j,lmass) = qsrflx(i,j,lmass) + fac_srflx*          &
    1033. 

  1033.                       raercol(kts,lmass,nsav)*surfrate(lmass)
    1034. 

  1034.                  qsrflx(i,j,lmasscw) = qsrflx(i,j,lmasscw) + fac_srflx*      &
    1035. 

  1035.                       raercol(kts,lmasscw,nsav)*surfrate(lmasscw)
    1036. 

  1036.                  if (icheck_colmass > 0) then
    1037. 

  1037.                     ! colmass_sfc calculation
    1038. 

  1038.                     !    colmass_bgn/end = bgn/end column burden = sum.over.k.of{ rho(k)*dz(k)*chem(k,l) }
    1039. 

  1039.                     !    colmass_sfc = surface loss over dtstep
    1040. 

  1040.                     !       = sum.over.nsubmix.substeps{ depvel(l)*rho(kts)*chem(kts,l)*dtmix }
    1041. 

  1041.                     !    surfrate(l) = depvel(l)/dz(kts) so need to multiply by dz(kts)
    1042. 

  1042.                     !    for mass, raercol(k,l) = chem(k,l)*1.0e-9, so need to multiply by 1.0e9
    1043. 

  1043.                     tmpf = dtmix*rhodz(kts)*1.0e9
    1044. 

  1044.                     colmass_sfc(l,m,n,1) = colmass_sfc(l,m,n,1) &
    1045. 

  1045.                           + raercol(kts,lmass  ,nsav)*surfrate(lmass  )*tmpf
    1046. 

  1046.                     colmass_sfc(l,m,n,2) = colmass_sfc(l,m,n,2) &
    1047. 

  1047.                           + raercol(kts,lmasscw,nsav)*surfrate(lmasscw)*tmpf
    1048. 

  1048.                  endif
    1049. 

  1049.               enddo
    1050. 

  1050.               lwater=waterptr_aer(m,n)  ! aerosol water
    1051. 

  1051.               if(lwater>0)then
    1052. 

  1052.                  source(:)=0.
    1053. 

  1053.                  call explmix(   raercol(1,lwater,nnew),source,ekkp,ekkm,overlapp,overlapm,   &
    1054. 

  1054.                       raercol(1,lwater,nsav),surfrate(lwater),kms,kme,kts,kte,dtmix,  &
    1055. 

  1055.                       .true.,source)
    1056. 

  1056.               endif
    1057. 

  1057.            enddo ! size
    1058. 

  1058.         enddo ! type
    1059. 

  1059. 

  1060.      enddo OLD_CLOUD_NSUBMIX_LOOP
    1061. 

  1061. 

  1062. !    cycle OVERALL_MAIN_I_LOOP
    1063. 

  1063. 

  1064. !        evaporate particles again if no cloud
    1065. 

  1065. 

  1066.      do k=kts,kte
    1067. 

  1067.         if(lcldfra(k).eq.0.)then
    1068. 

  1068. 

  1069. !              no cloud
    1070. 

  1070. 

  1071.            qndrop(k)=0.
    1072. 

  1072. !              convert activated aerosol to interstitial in decaying cloud
    1073. 

  1073.            do n=1,ntype_aer
    1074. 

  1074.               do m=1,nsize_aer(n)
    1075. 

  1075.                  lnum=numptr_aer(m,n,ai_phase)
    1076. 

  1076.                  lnumcw=numptr_aer(m,n,cw_phase)
    1077. 

  1077.                  if(lnum.gt.0)then
    1078. 

  1078.                     raercol(k,lnum,nnew)=raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew)
    1079. 

  1079.                     raercol(k,lnumcw,nnew)=0.
    1080. 

  1080.                  endif
    1081. 

  1081.                  do l=1,ncomp(n)
    1082. 

  1082.                     lmass=massptr_aer(l,m,n,ai_phase)
    1083. 

  1083.                     lmasscw=massptr_aer(l,m,n,cw_phase)
    1084. 

  1084.                     raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol(k,lmasscw,nnew)
    1085. 

  1085.                     raercol(k,lmasscw,nnew)=0.
    1086. 

  1086.                  enddo
    1087. 

  1087.               enddo
    1088. 

  1088.            enddo
    1089. 

  1089.         endif
    1090. 

  1090.      enddo
    1091. 

  1091. 

  1092. !         cycle OVERALL_MAIN_I_LOOP
    1093. 

  1093. 

  1094. !        droplet number
    1095. 

  1095. 

  1096.      do k=kts,kte
    1097. 

  1097. !       if(lcldfra(k).gt.0.1)then
    1098. 

  1098. !           write(6,'(a,3i5,f12.1)')'i,j,k,qndrop=',i,j,k,qndrop(k)
    1099. 

  1099. !       endif
    1100. 

  1100.         if(qndrop(k).lt.-10.e6.or.qndrop(k).gt.1.e12)then
    1101. 

  1101.            write(6,'(a,g12.2,a,3i5)')'after qndrop=',qndrop(k),' for i,k,j=',i,k,j
    1102. 

  1102.         endif
    1103. 

  1103. 

  1104.         qndrop3d(i,k,j) = max(qndrop(k),1.e-6)
    1105. 

  1105. 

  1106.         if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
    1107. 

  1107.            write(6,'(a,g12.2,a,3i5)')'after qndrop3d=',qndrop3d(i,k,j),' for i,k,j=',i,k,j
    1108. 

  1108.         endif
    1109. 

  1109.         if(qc(i,k,j).lt.-1..or.qc(i,k,j).gt.1.)then
    1110. 

  1110.            write(6,'(a,g12.2,a,3i5)')'qc=',qc(i,k,j),' for i,k,j=',i,k,j
    1111. 

  1111.            call wrf_error_fatal("1")
    1112. 

  1112.         endif
    1113. 

  1113.         if(qi(i,k,j).lt.-1..or.qi(i,k,j).gt.1.)then
    1114. 

  1114.            write(6,'(a,g12.2,a,3i5)')'qi=',qi(i,k,j),' for i,k,j=',i,k,j
    1115. 

  1115.            call wrf_error_fatal("1")
    1116. 

  1116.         endif
    1117. 

  1117.         if(qv(i,k,j).lt.-1..or.qv(i,k,j).gt.1.)then
    1118. 

  1118.            write(6,'(a,g12.2,a,3i5)')'qv=',qv(i,k,j),' for i,k,j=',i,k,j
    1119. 

  1119.            call wrf_error_fatal("1")
    1120. 

  1120.         endif
    1121. 

  1121.         cldfra_old(i,k,j) = cldfra(i,k,j)
    1122. 

  1122. !       if(k.gt.6.and.k.lt.11)cldfra_old(i,k,j)=1.
    1123. 

  1123.      enddo
    1124. 

  1124. 

  1125. !    cycle OVERALL_MAIN_I_LOOP
    1126. 

  1126. 

  1127. !        update chem and convert back to mole/mole
    1128. 

  1128. 

  1129.      ccn(:,:) = 0.
    1130. 

  1130.      do n=1,ntype_aer
    1131. 

  1131.         do m=1,nsize_aer(n)
    1132. 

  1132.            lnum=numptr_aer(m,n,ai_phase)
    1133. 

  1133.            lnumcw=numptr_aer(m,n,cw_phase)
    1134. 

  1134.            if(lnum.gt.0)then
    1135. 

  1135.               !          scale=mwdry*0.001
    1136. 

  1136.               scale = 1.
    1137. 

  1137.               chem(i,kts:kte,j,lnumcw)= raercol(kts:kte,lnumcw,nnew)*scale
    1138. 

  1138.               chem(i,kts:kte,j,lnum)= raercol(kts:kte,lnum,nnew)*scale
    1139. 

  1139.            endif
    1140. 

  1140.            do l=1,ncomp(n)
    1141. 

  1141.               lmass=massptr_aer(l,m,n,ai_phase)
    1142. 

  1142.               lmasscw=massptr_aer(l,m,n,cw_phase)
    1143. 

  1143. !          scale = mwdry/mw_aer(l,n)
    1144. 

  1144.               scale = 1.e9
    1145. 

  1145.               chem(i,kts:kte,j,lmasscw)=raercol(kts:kte,lmasscw,nnew)*scale ! ug/kg
    1146. 

  1146.               chem(i,kts:kte,j,lmass)=raercol(kts:kte,lmass,nnew)*scale ! ug/kg
    1147. 

  1147.            enddo
    1148. 

  1148.            lwater=waterptr_aer(m,n)
    1149. 

  1149.            if(lwater>0)chem(i,kts:kte,j,lwater)=raercol(kts:kte,lwater,nnew) ! don't convert units
    1150. 

  1150.            do k=kts,kte
    1151. 

  1151.               sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n)))
    1152. 

  1152.               do l=1,psat
    1153. 

  1153.                  arg=argfactor(m,n)*log(sm/super(l))
    1154. 

  1154.                  if(arg<2)then
    1155. 

  1155.                     if(arg<-2)then
    1156. 

  1156.                        ccnfact(l,m,n)=1.e-6 ! convert from #/m3 to #/cm3
    1157. 

  1157.                     else
    1158. 

  1158.                        ccnfact(l,m,n)=1.e-6*0.5*ERFC_NUM_RECIPES(arg)
    1159. 

  1159.                     endif
    1160. 

  1160.                  else
    1161. 

  1161.                     ccnfact(l,m,n) = 0.
    1162. 

  1162.                  endif
    1163. 

  1163. !                 ccn concentration as diagnostic
    1164. 

  1164. !                 assume same hygroscopicity and ccnfact for cloud-phase and aerosol phase particles
    1165. 

  1165.                  ccn(k,l)=ccn(k,l)+(raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew))*cs(k)*ccnfact(l,m,n)
    1166. 

  1166.               enddo
    1167. 

  1167.            enddo
    1168. 

  1168.         enddo
    1169. 

  1169.      enddo
    1170. 

  1170.      do l=1,psat
    1171. 

  1171.         !wig, 22-Nov-2006: added vertical bounds to prevent out-of-bounds at top
    1172. 

  1172.         if(l.eq.1)ccn1(i,kts:kte,j)=ccn(:,l)
    1173. 

  1173.         if(l.eq.2)ccn2(i,kts:kte,j)=ccn(:,l)
    1174. 

  1174.         if(l.eq.3)ccn3(i,kts:kte,j)=ccn(:,l)
    1175. 

  1175.         if(l.eq.4)ccn4(i,kts:kte,j)=ccn(:,l)
    1176. 

  1176.         if(l.eq.5)ccn5(i,kts:kte,j)=ccn(:,l)
    1177. 

  1177.         if(l.eq.6)ccn6(i,kts:kte,j)=ccn(:,l)
    1178. 

  1178.      end do
    1179. 

  1179. 

  1180. ! mass conservation checking
    1181. 

  1181.      if (icheck_colmass > 0) then
    1182. 

  1182. ! calc final column burdens
    1183. 

  1183.         do n=1,ntype_aer
    1184. 

  1184.         do m=1,nsize_aer(n)
    1185. 

  1185.            lnum=numptr_aer(m,n,ai_phase)
    1186. 

  1186.            lnumcw=numptr_aer(m,n,cw_phase)
    1187. 

  1187.            if(lnum>0)then
    1188. 

  1188.               colmass_end(0,m,n,1) = sum( chem(i,kts:kte,j,lnum  )*rhodz(kts:kte) )
    1189. 

  1189.               colmass_end(0,m,n,2) = sum( chem(i,kts:kte,j,lnumcw)*rhodz(kts:kte) )
    1190. 

  1190.            endif
    1191. 

  1191.            do l=1,ncomp(n)
    1192. 

  1192.               lmass=massptr_aer(l,m,n,ai_phase)
    1193. 

  1193.               lmasscw=massptr_aer(l,m,n,cw_phase)
    1194. 

  1194.               colmass_end(l,m,n,1) = sum( chem(i,kts:kte,j,lmass  )*rhodz(kts:kte) )
    1195. 

  1195.               colmass_end(l,m,n,2) = sum( chem(i,kts:kte,j,lmasscw)*rhodz(kts:kte) )
    1196. 

  1196.            enddo
    1197. 

  1197.         enddo ! size
    1198. 

  1198.         enddo ! type
    1199. 

  1199. ! calc burden change errors for each interstitial/activated pair
    1200. 

  1200.         do n=1,ntype_aer
    1201. 

  1201.         do m=1,nsize_aer(n)
    1202. 

  1202.            do l=0,ncomp(n)
    1203. 

  1203.               ! tmpa & tmpb = beginning & ending column burden divided by rhodzsum,
    1204. 

  1204.               !             = beginning & ending column-mean mixing ratios
    1205. 

  1205.               ! tmpc = loss to surface divided by rhodzsum,
    1206. 

  1206.               tmpa = ( colmass_bgn(l,m,n,1) + colmass_bgn(l,m,n,2) )/rhodzsum
    1207. 

  1207.               tmpb = ( colmass_end(l,m,n,1) + colmass_end(l,m,n,2) )/rhodzsum
    1208. 

  1208.               tmpc = ( colmass_sfc(l,m,n,1) + colmass_sfc(l,m,n,2) )/rhodzsum
    1209. 

  1209. 

  1210.               ! tmpd = ((final burden) + (sfc loss)) - (initial burden)
    1211. 

  1211.               !      = burden change error
    1212. 

  1212.               tmpd = (tmpb + tmpc) - tmpa
    1213. 

  1213.               tmpe = max( tmpa, 1.0e-20 )
    1214. 

  1214. 

  1215.               ! tmpf = (burden change error) / (initial burden)
    1216. 

  1216.               if (abs(tmpd) < 1.0e5*tmpe) then
    1217. 

  1217.                  tmpf = tmpd/tmpe
    1218. 

  1218.               else if (tmpf < 0.0) then
    1219. 

  1219.                  tmpf = -1.0e5
    1220. 

  1220.               else
    1221. 

  1221.                  tmpf = 1.0e5
    1222. 

  1222.               end if
    1223. 

  1223.               if (abs(tmpf) > abs(colmass_worst(l,m,n))) then
    1224. 

  1224.                  colmass_worst(l,m,n) = tmpf
    1225. 

  1225.                  colmass_worst_ij(1,l,m,n) = i
    1226. 

  1226.                  colmass_worst_ij(2,l,m,n) = j
    1227. 

  1227.               endif
    1228. 

  1228.            enddo
    1229. 

  1229.         enddo ! size
    1230. 

  1230.         enddo ! type
    1231. 

  1231.      endif ! (icheck_colmass > 0)
    1232. 

  1232. 

  1233. 

  1234.      enddo OVERALL_MAIN_I_LOOP ! end of main loop over i
    1235. 

  1235.      enddo OVERALL_MAIN_J_LOOP ! end of main loop over j
    1236. 

  1236. 

  1237. 

  1238. ! mass conservation checking
    1239. 

  1239.      if (icheck_colmass > 0) then
    1240. 

  1240.         if (icheck_colmass >= 100) write(*,'(a)') &
    1241. 

  1241.              'mixactivate colmass worst errors bgn - type, size, comp, err, i, j'
    1242. 

  1242.         colmass_maxworst_r = 0.0
    1243. 

  1243.         colmass_maxworst_i(:) = -1
    1244. 

  1244.         do n=1,ntype_aer
    1245. 

  1245.         do m=1,nsize_aer(n)
    1246. 

  1246.            do l=0,ncomp(n)
    1247. 

  1247.               if (icheck_colmass >= 100) &
    1248. 

  1248.                  write(*,'(3i3,1p,e10.2,2i4)') n, m, l, &
    1249. 

  1249.                  colmass_worst(l,m,n), colmass_worst_ij(1:2,l,m,n) 
    1250. 

  1250.               if (abs(colmass_worst(l,m,n)) > abs(colmass_maxworst_r)) then
    1251. 

  1251.                  colmass_maxworst_r = colmass_worst(l,m,n) 
    1252. 

  1252.                  colmass_maxworst_i(1) = n
    1253. 

  1253.                  colmass_maxworst_i(2) = m
    1254. 

  1254.                  colmass_maxworst_i(3) = l
    1255. 

  1255.               end if
    1256. 

  1256.            enddo
    1257. 

  1257.         enddo ! size
    1258. 

  1258.         enddo ! type
    1259. 

  1259.         if ((icheck_colmass >= 10) .or. (abs(colmass_maxworst_r) >= 1.0e-6)) &
    1260. 

  1260.              write(*,'(a,3i3,1p,e10.2)') 'mixactivate colmass maxworst', &
    1261. 

  1261.              colmass_maxworst_i(1:3), colmass_maxworst_r
    1262. 

  1262.      endif ! (icheck_colmass > 0)
    1263. 

  1263. 

  1264. 

  1265.      return
    1266. 

  1266.    end subroutine mixactivate
    1267. 

  1267. 

  1268. 

  1269. !----------------------------------------------------------------------
    1270. 

  1270. !----------------------------------------------------------------------
    1271. 

  1271.    subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, &
    1272. 

  1272.                        qold, surfrate, kms, kme, kts, kte, dt, &
    1273. 

  1273.                        is_unact, qactold )
    1274. 

  1274. 

  1275. !  explicit integration of droplet/aerosol mixing
    1276. 

  1276. !     with source due to activation/nucleation
    1277. 

  1277. 

  1278. 

  1279.    implicit none
    1280. 

  1280.    integer, intent(in) :: kms,kme ! number of levels for array definition
    1281. 

  1281.    integer, intent(in) :: kts,kte ! number of levels for looping
    1282. 

  1282.    real, intent(inout) :: q(kms:kme) ! mixing ratio to be updated
    1283. 

  1283.    real, intent(in) :: qold(kms:kme) ! mixing ratio from previous time step
    1284. 

  1284.    real, intent(in) :: src(kms:kme) ! source due to activation/nucleation (/s)
    1285. 

  1285.    real, intent(in) :: ekkp(kms:kme) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
    1286. 

  1286.                       ! below layer k  (k,k+1 interface)
    1287. 

  1287.    real, intent(in) :: ekkm(kms:kme) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
    1288. 

  1288.                       ! above layer k  (k,k+1 interface)
    1289. 

  1289.    real, intent(in) :: overlapp(kms:kme) ! cloud overlap below
    1290. 

  1290.    real, intent(in) :: overlapm(kms:kme) ! cloud overlap above
    1291. 

  1291.    real, intent(in) :: surfrate ! surface exchange rate (/s)
    1292. 

  1292.    real, intent(in) :: dt ! time step (s)
    1293. 

  1293.    logical, intent(in) :: is_unact ! true if this is an unactivated species
    1294. 

  1294.    real, intent(in),optional :: qactold(kms:kme)
    1295. 

  1295.           ! mixing ratio of ACTIVATED species from previous step
    1296. 

  1296.           ! *** this should only be present
    1297. 

  1297.           !     if the current species is unactivated number/sfc/mass
    1298. 

  1298. 

  1299.    integer k,kp1,km1
    1300. 

  1300. 

  1301.    if ( is_unact ) then
    1302. 

  1302. !     the qactold*(1-overlap) terms are resuspension of activated material
    1303. 

  1303.       do k=kts,kte
    1304. 

  1304.          kp1=min(k+1,kte)
    1305. 

  1305.          km1=max(k-1,kts)
    1306. 

  1306.          q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) +  &
    1307. 

  1307.                            qactold(kp1)*(1.0-overlapp(k)))               &
    1308. 

  1308.                                   + ekkm(k)*(qold(km1) - qold(k) +     &
    1309. 

  1309.                            qactold(km1)*(1.0-overlapm(k))) )
    1310. 

  1310. !         if(q(k)<-1.e-30)then ! force to non-negative
    1311. 

  1311. !            print *,'q=',q(k),' in explmix'
    1312. 

  1312.              q(k)=max(q(k),0.)
    1313. 

  1313. !         endif
    1314. 

  1314.       end do
    1315. 

  1315. 

  1316.    else
    1317. 

  1317.       do k=kts,kte
    1318. 

  1318.          kp1=min(k+1,kte)
    1319. 

  1319.          km1=max(k-1,kts)
    1320. 

  1320.          q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) +  &
    1321. 

  1321.                                     ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
    1322. 

  1322. !        if(q(k)<-1.e-30)then ! force to non-negative
    1323. 

  1323. !           print *,'q=',q(k),' in explmix'
    1324. 

  1324.             q(k)=max(q(k),0.)
    1325. 

  1325. !        endif
    1326. 

  1326.       end do
    1327. 

  1327.    end if
    1328. 

  1328. 

  1329. !  dry deposition loss at base of lowest layer
    1330. 

  1330.    q(kts)=q(kts)-surfrate*qold(kts)*dt
    1331. 

  1331. !  if(q(kts)<-1.e-30)then ! force to non-negative
    1332. 

  1332. !     print *,'q=',q(kts),' in explmix'
    1333. 

  1333.       q(kts)=max(q(kts),0.)
    1334. 

  1334. !  endif
    1335. 

  1335. 

  1336.    return
    1337. 

  1337.    end subroutine explmix
    1338. 

  1338. 

  1339. !----------------------------------------------------------------------
    1340. 

  1340. !----------------------------------------------------------------------
    1341. 

  1341. ! 06-nov-2005 rce - grid_id & ktau added to arg list
    1342. 

  1342.       subroutine activate(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,  &
    1343. 

  1343.                       msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
    1344. 

  1344.                       na, volc, dlo_sect,dhi_sect,sigman, hygro, &
    1345. 

  1345.                       fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, &
    1346. 

  1346.                       grid_id, ktau, ii, jj, kk )
    1347. 

  1347. 

  1348. !      calculates number, surface, and mass fraction of aerosols activated as CCN
    1349. 

  1349. !      calculates flux of cloud droplets, surface area, and aerosol mass into cloud
    1350. 

  1350. !      assumes an internal mixture within each of aerosol mode.
    1351. 

  1351. !      A sectional treatment within each type is assumed if ntype_aer >7.
    1352. 

  1352. !      A gaussiam spectrum of updrafts can be treated.
    1353. 

  1353. 

  1354. !      mks units
    1355. 

  1355. 

  1356. !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
    1357. 

  1357. !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
    1358. 

  1358. 

  1359.       USE module_model_constants, only: g,rhowater, xlv, cp, rvovrd, r_d, r_v, &
    1360. 

  1360.               mwdry,svp1,svp2,svp3,ep_2
    1361. 

  1361. 

  1362.       implicit none
    1363. 

  1363. 

  1364. 

  1365. !      input
    1366. 

  1366. 

  1367.       integer,intent(in) :: maxd_atype      ! dimension of types
    1368. 

  1368.       integer,intent(in) :: maxd_asize      ! dimension of sizes
    1369. 

  1369.       integer,intent(in) :: ntype_aer       ! number of types
    1370. 

  1370.       integer,intent(in) :: nsize_aer(maxd_atype) ! number of sizes for type
    1371. 

  1371.       integer,intent(in) :: msectional      ! 1 for sectional, 0 for modal
    1372. 

  1372.       integer,intent(in) :: grid_id         ! WRF grid%id
    1373. 

  1373.       integer,intent(in) :: ktau            ! WRF time step count
    1374. 

  1374.       integer,intent(in) :: ii, jj, kk      ! i,j,k of current grid cell
    1375. 

  1375.       real,intent(in) :: wbar          ! grid cell mean vertical velocity (m/s)
    1376. 

  1376.       real,intent(in) :: sigw          ! subgrid standard deviation of vertical vel (m/s)
    1377. 

  1377.       real,intent(in) :: wdiab         ! diabatic vertical velocity (0 if adiabatic)
    1378. 

  1378.       real,intent(in) :: wminf         ! minimum updraft velocity for integration (m/s)
    1379. 

  1379.       real,intent(in) :: wmaxf         ! maximum updraft velocity for integration (m/s)
    1380. 

  1380.       real,intent(in) :: tair          ! air temperature (K)
    1381. 

  1381.       real,intent(in) :: rhoair        ! air density (kg/m3)
    1382. 

  1382.       real,intent(in) :: na(maxd_asize,maxd_atype)     ! aerosol number concentration (/m3)
    1383. 

  1383.       real,intent(in) :: sigman(maxd_asize,maxd_atype) ! geometric standard deviation of aerosol size distribution
    1384. 

  1384.       real,intent(in) :: hygro(maxd_asize,maxd_atype)  ! bulk hygroscopicity of aerosol mode
    1385. 

  1385.       real,intent(in) :: volc(maxd_asize,maxd_atype)   ! total aerosol volume  concentration (m3/m3)
    1386. 

  1386.       real,intent(in) :: dlo_sect( maxd_asize, maxd_atype ), &  ! minimum size of section (cm)
    1387. 

  1387.            dhi_sect( maxd_asize, maxd_atype )     ! maximum size of section (cm)
    1388. 

  1388. 

  1389. !      output
    1390. 

  1390. 

  1391.       real,intent(inout) :: fn(maxd_asize,maxd_atype)    ! number fraction of aerosols activated
    1392. 

  1392.       real,intent(inout) :: fs(maxd_asize,maxd_atype)    ! surface fraction of aerosols activated
    1393. 

  1393.       real,intent(inout) :: fm(maxd_asize,maxd_atype)    ! mass fraction of aerosols activated
    1394. 

  1394.       real,intent(inout) :: fluxn(maxd_asize,maxd_atype) ! flux of activated aerosol number fraction into cloud (m/s)
    1395. 

  1395.       real,intent(inout) :: fluxs(maxd_asize,maxd_atype) ! flux of activated aerosol surface fraction (m/s)
    1396. 

  1396.       real,intent(inout) :: fluxm(maxd_asize,maxd_atype) ! flux of activated aerosol mass fraction into cloud (m/s)
    1397. 

  1397.       real,intent(inout) :: flux_fullact                 ! flux when activation fraction = 100% (m/s)
    1398. 

  1398. 

  1399. !      local
    1400. 

  1400. 

  1401. !!$      external erf,erfc
    1402. 

  1402. !!$      real erf,erfc
    1403. 

  1403. !      external qsat_water
    1404. 

  1404.       integer, parameter:: nx=200
    1405. 

  1405.       integer iquasisect_option, isectional
    1406. 

  1406.       real integ,integf
    1407. 

  1407.       real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m)
    1408. 

  1408.       real, parameter :: p0 = 1013.25e2  ! reference pressure (Pa)
    1409. 

  1409.       real, parameter :: t0 = 273.15     ! reference temperature (K)
    1410. 

  1410.       real ylo(maxd_asize,maxd_atype),yhi(maxd_asize,maxd_atype) ! 1-particle volume at section interfaces
    1411. 

  1411.       real ymean(maxd_asize,maxd_atype) ! 1-particle volume at r=rmean
    1412. 

  1412.       real ycut, lnycut, betayy, betayy2, gammayy, phiyy
    1413. 

  1413.       real surfc(maxd_asize,maxd_atype) ! surface concentration (m2/m3)
    1414. 

  1414.       real sign(maxd_asize,maxd_atype)    ! geometric standard deviation of size distribution
    1415. 

  1415.       real alnsign(maxd_asize,maxd_atype) ! natl log of geometric standard dev of aerosol
    1416. 

  1416.       real am(maxd_asize,maxd_atype) ! number mode radius of dry aerosol (m)
    1417. 

  1417.       real lnhygro(maxd_asize,maxd_atype) ! ln(b)
    1418. 

  1418.       real f1(maxd_asize,maxd_atype) ! array to hold parameter for maxsat
    1419. 

  1419.       real pres ! pressure (Pa)
    1420. 

  1420.       real path ! mean free path (m)
    1421. 

  1421.       real diff ! diffusivity (m2/s)
    1422. 

  1422.       real conduct ! thermal conductivity (Joule/m/sec/deg)
    1423. 

  1423.       real diff0,conduct0
    1424. 

  1424.       real es ! saturation vapor pressure
    1425. 

  1425.       real qs ! water vapor saturation mixing ratio
    1426. 

  1426.       real dqsdt ! change in qs with temperature
    1427. 

  1427.       real dqsdp ! change in qs with pressure
    1428. 

  1428.       real gg ! thermodynamic function (m2/s)
    1429. 

  1429.       real sqrtg ! sqrt(gg)
    1430. 

  1430.       real sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius
    1431. 

  1431.       real lnsm(maxd_asize,maxd_atype) ! ln( sm )
    1432. 

  1432.       real zeta, eta(maxd_asize,maxd_atype)
    1433. 

  1433.       real lnsmax ! ln(smax)
    1434. 

  1434.       real alpha
    1435. 

  1435.       real gamma
    1436. 

  1436.       real beta
    1437. 

  1437.       real gaus
    1438. 

  1438.       logical :: top        ! true if cloud top, false if cloud base or new cloud
    1439. 

  1439.       real asub(maxd_asize,maxd_atype),bsub(maxd_asize,maxd_atype) ! coefficients of submode size distribution N=a+bx
    1440. 

  1440.       real totn(maxd_atype) ! total aerosol number concentration
    1441. 

  1441.       real aten ! surface tension parameter
    1442. 

  1442.       real gmrad(maxd_atype) ! geometric mean radius
    1443. 

  1443.       real gmradsq(maxd_atype) ! geometric mean of radius squared
    1444. 

  1444.       real gmlnsig(maxd_atype) ! geometric standard deviation
    1445. 

  1445.       real gmsm(maxd_atype) ! critical supersaturation at radius gmrad
    1446. 

  1446.       real sumflxn(maxd_asize,maxd_atype)
    1447. 

  1447.       real sumflxs(maxd_asize,maxd_atype)
    1448. 

  1448.       real sumflxm(maxd_asize,maxd_atype)
    1449. 

  1449.       real sumflx_fullact
    1450. 

  1450.       real sumfn(maxd_asize,maxd_atype)
    1451. 

  1451.       real sumfs(maxd_asize,maxd_atype)
    1452. 

  1452.       real sumfm(maxd_asize,maxd_atype)
    1453. 

  1453.       real sumns(maxd_atype)
    1454. 

  1454.       real fnold(maxd_asize,maxd_atype)   ! number fraction activated
    1455. 

  1455.       real fsold(maxd_asize,maxd_atype)   ! surface fraction activated
    1456. 

  1456.       real fmold(maxd_asize,maxd_atype)   ! mass fraction activated
    1457. 

  1457.       real wold,gold
    1458. 

  1458.       real alogten,alog2,alog3,alogaten
    1459. 

  1459.       real alogam
    1460. 

  1460.       real rlo(maxd_asize,maxd_atype), rhi(maxd_asize,maxd_atype)
    1461. 

  1461.       real rmean(maxd_asize,maxd_atype)
    1462. 

  1462.                   ! mean radius (m) for the section (not used with modal)
    1463. 

  1463.                   ! calculated from current volume & number
    1464. 

  1464.       real ccc
    1465. 

  1465.       real dumaa,dumbb
    1466. 

  1466.       real wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
    1467. 

  1467.       real dfmin,dfmax,fnew,fold,fnmin,fnbar,fsbar,fmbar
    1468. 

  1468.       real alw,sqrtalw
    1469. 

  1469.       real smax
    1470. 

  1470.       real x,arg
    1471. 

  1471.       real xmincoeff,xcut
    1472. 

  1472.       real z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
    1473. 

  1473.       real etafactor1,etafactor2(maxd_asize,maxd_atype),etafactor2max
    1474. 

  1474.       integer m,n,nw,nwmax
    1475. 

  1475. 

  1476. !      numerical integration parameters
    1477. 

  1477.       real, parameter :: eps = 0.3
    1478. 

  1478.       real, parameter :: fmax = 0.99
    1479. 

  1479.       real, parameter :: sds = 3.
    1480. 

  1480. 

  1481. !      mathematical constants
    1482. 

  1482.       real third, twothird, sixth, zero, one, two, three
    1483. 

  1483. 

  1484.       real, parameter :: sq2  = 1.4142135624
    1485. 

  1485.       real, parameter :: sqpi = 1.7724538509
    1486. 

  1486.       real, parameter :: pi   = 3.1415926536
    1487. 

  1487. 

  1488. !      integer, save :: ndist(nx)  ! accumulates frequency distribution of integration bins required
    1489. 

  1489. !      data ndist/nx*0/
    1490. 

  1490. 

  1491. !     for nsize_aer>7, a sectional approach is used and isectional = iquasisect_option
    1492. 

  1492. !     activation fractions (fn,fs,fm) are computed as follows
    1493. 

  1493. !     iquasisect_option = 1,3 - each section treated as a narrow lognormal
    1494. 

  1494. !     iquasisect_option = 2,4 - within-section dn/dx = a + b*x,  x = ln(r)
    1495. 

  1495. !     smax is computed as follows (when explicit activation is OFF)
    1496. 

  1496. !     iquasisect_option = 1,2 - razzak-ghan modal parameterization with
    1497. 

  1497. !     single mode having same ntot, dgnum, sigmag as the combined sections
    1498. 

  1498. !     iquasisect_option = 3,4 - razzak-ghan sectional parameterization
    1499. 

  1499. !     for nsize_aer=<9, a modal approach is used and isectional = 0
    1500. 

  1500. 

  1501. ! rce 08-jul-2005
    1502. 

  1502. ! if either (na(n,m) < nsmall) or (volc(n,m) < vsmall)
    1503. 

  1503. ! then treat bin/mode (n,m) as being empty, and set its fn/fs/fm=0.0
    1504. 

  1504. !     (for single precision, gradual underflow starts around 1.0e-38,
    1505. 

  1505. !      and strange things can happen when in that region)
    1506. 

  1506.       real, parameter :: nsmall = 1.0e-20    ! aer number conc in #/m3
    1507. 

  1507.       real, parameter :: vsmall = 1.0e-37    ! aer volume conc in m3/m3
    1508. 

  1508.       logical bin_is_empty(maxd_asize,maxd_atype), all_bins_empty
    1509. 

  1509.       logical bin_is_narrow(maxd_asize,maxd_atype)
    1510. 

  1510. 

  1511.       integer idiagaa, ipass_nwloop
    1512. 

  1512.       integer idiag_dndy_neg, idiag_fnsm_prob
    1513. 

  1513. 

  1514. ! The flag for cloud top is no longer used so set it to false. This is an
    1515. 

  1515. ! antiquated feature related to radiation enhancing mass fluxes at cloud
    1516. 

  1516. ! top. It is currently, as of Feb. 2009, set to false in the CAM version
    1517. 

  1517. ! as well.
    1518. 

  1518.       top = .false.
    1519. 

  1519. 

  1520. !.......................................................................
    1521. 

  1521. !
    1522. 

  1522. !   start calc. of modal or sectional activation properties (start of section 1)
    1523. 

  1523. !
    1524. 

  1524. !.......................................................................
    1525. 

  1525.       idiag_dndy_neg = 1      ! set this to 0 to turn off 
    1526. 

  1526.                               !     warnings about dn/dy < 0
    1527. 

  1527.       idiag_fnsm_prob = 1     ! set this to 0 to turn off 
    1528. 

  1528.                               !     warnings about fn/fs/fm misbehavior
    1529. 

  1529. 

  1530.       iquasisect_option = 2
    1531. 

  1531.       if(msectional.gt.0)then
    1532. 

  1532.          isectional = iquasisect_option
    1533. 

  1533.       else
    1534. 

  1534.          isectional = 0
    1535. 

  1535.       endif
    1536. 

  1536. 

  1537.       do n=1,ntype_aer
    1538. 

  1538. !         print *,'ntype_aer,n,nsize_aer(n)=',ntype_aer,n,nsize_aer(n)
    1539. 

  1539. 

  1540.         if(ntype_aer.eq.1.and.nsize_aer(n).eq.1.and.na(1,1).lt.1.e-20)then
    1541. 

  1541.          fn(1,1)=0.
    1542. 

  1542.          fs(1,1)=0.
    1543. 

  1543.          fm(1,1)=0.
    1544. 

  1544.          fluxn(1,1)=0.
    1545. 

  1545.          fluxs(1,1)=0.
    1546. 

  1546.          fluxm(1,1)=0.
    1547. 

  1547.          flux_fullact=0.
    1548. 

  1548.          return
    1549. 

  1549.         endif
    1550. 

  1550.       enddo
    1551. 

  1551. 

  1552.       zero = 0.0
    1553. 

  1553.       one = 1.0
    1554. 

  1554.       two = 2.0
    1555. 

  1555.       three = 3.0
    1556. 

  1556.       third = 1.0/3.0
    1557. 

  1557.       twothird = 2.0/3.0 !wig, 1-Mar-2009: Corrected value from 2/6
    1558. 

  1558.       sixth = 1.0/6.0
    1559. 

  1559. 

  1560.       pres=r_d*rhoair*tair
    1561. 

  1561.       diff0=0.211e-4*(p0/pres)*(tair/t0)**1.94
    1562. 

  1562.       conduct0=(5.69+0.017*(tair-t0))*4.186e2*1.e-5 ! convert to J/m/s/deg
    1563. 

  1563.       es=1000.*svp1*exp( svp2*(tair-t0)/(tair-svp3) )
    1564. 

  1564.       qs=ep_2*es/(pres-es)
    1565. 

  1565.       dqsdt=xlv/(r_v*tair*tair)*qs
    1566. 

  1566.       alpha=g*(xlv/(cp*r_v*tair*tair)-1./(r_d*tair))
    1567. 

  1567.       gamma=(1+xlv/cp*dqsdt)/(rhoair*qs)
    1568. 

  1568.       gg=1./(rhowater/(diff0*rhoair*qs)+xlv*rhowater/(conduct0*tair)*(xlv/(r_v*tair)-1.))
    1569. 

  1569.       sqrtg=sqrt(gg)
    1570. 

  1570.       beta=4.*pi*rhowater*gg*gamma
    1571. 

  1571.       aten=2.*surften/(r_v*tair*rhowater)
    1572. 

  1572.       alogaten=log(aten)
    1573. 

  1573.       alog2=log(two)
    1574. 

  1574.       alog3=log(three)
    1575. 

  1575.       ccc=4.*pi*third
    1576. 

  1576.       etafactor2max=1.e10/(alpha*wmaxf)**1.5 ! this should make eta big if na is very small.
    1577. 

  1577. 

  1578.       all_bins_empty = .true.
    1579. 

  1579.       do n=1,ntype_aer
    1580. 

  1580.       totn(n)=0.
    1581. 

  1581.       gmrad(n)=0.
    1582. 

  1582.       gmradsq(n)=0.
    1583. 

  1583.       sumns(n)=0.
    1584. 

  1584.       do m=1,nsize_aer(n)
    1585. 

  1585.          alnsign(m,n)=log(sigman(m,n))
    1586. 

  1586. !         internal mixture of aerosols
    1587. 

  1587. 

  1588.          bin_is_empty(m,n) = .true.
    1589. 

  1589.          if (volc(m,n).gt.vsmall .and. na(m,n).gt.nsmall) then
    1590. 

  1590.             bin_is_empty(m,n) = .false.
    1591. 

  1591.             all_bins_empty = .false.
    1592. 

  1592.             lnhygro(m,n)=log(hygro(m,n))
    1593. 

  1593. !            number mode radius (m,n)
    1594. 

  1594. !           write(6,*)'alnsign,volc,na=',alnsign(m,n),volc(m,n),na(m,n)
    1595. 

  1595.             am(m,n)=exp(-1.5*alnsign(m,n)*alnsign(m,n))*              &
    1596. 

  1596.               (3.*volc(m,n)/(4.*pi*na(m,n)))**third
    1597. 

  1597. 

  1598.             if (isectional .gt. 0) then
    1599. 

  1599. !               sectional model.
    1600. 

  1600. !               need to use bulk properties because parameterization doesn't
    1601. 

  1601. !               work well for narrow bins.
    1602. 

  1602.                totn(n)=totn(n)+na(m,n)
    1603. 

  1603.                alogam=log(am(m,n))
    1604. 

  1604.                gmrad(n)=gmrad(n)+na(m,n)*alogam
    1605. 

  1605.                gmradsq(n)=gmradsq(n)+na(m,n)*alogam*alogam
    1606. 

  1606.             endif
    1607. 

  1607.             etafactor2(m,n)=1./(na(m,n)*beta*sqrtg)
    1608. 

  1608. 

  1609.             if(hygro(m,n).gt.1.e-10)then
    1610. 

  1610.                sm(m,n)=2.*aten/(3.*am(m,n))*sqrt(aten/(3.*hygro(m,n)*am(m,n)))
    1611. 

  1611.             else
    1612. 

  1612.                sm(m,n)=100.
    1613. 

  1613.             endif
    1614. 

  1614. !           write(6,*)'sm,hygro,am=',sm(m,n),hygro(m,n),am(m,n)
    1615. 

  1615.          else
    1616. 

  1616.             sm(m,n)=1.
    1617. 

  1617.             etafactor2(m,n)=etafactor2max ! this should make eta big if na is very small.
    1618. 

  1618. 

  1619.          endif
    1620. 

  1620.          lnsm(m,n)=log(sm(m,n))
    1621. 

  1621.          if ((isectional .eq. 3) .or. (isectional .eq. 4)) then
    1622. 

  1622.             sumns(n)=sumns(n)+na(m,n)/sm(m,n)**twothird
    1623. 

  1623.          endif
    1624. 

  1624. !        write(6,'(a,i4,6g12.2)')'m,na,am,hygro,lnhygro,sm,lnsm=',m,na(m,n),am(m,n),hygro(m,n),lnhygro(m,n),sm(m,n),lnsm(m,n)
    1625. 

  1625.       end do ! size
    1626. 

  1626.       end do ! type
    1627. 

  1627. 

  1628. !  if all bins are empty, set all activation fractions to zero and exit
    1629. 

  1629.          if ( all_bins_empty ) then
    1630. 

  1630.             do n=1,ntype_aer
    1631. 

  1631.             do m=1,nsize_aer(n)
    1632. 

  1632.                fluxn(m,n)=0.
    1633. 

  1633.                fn(m,n)=0.
    1634. 

  1634.                fluxs(m,n)=0.
    1635. 

  1635.                fs(m,n)=0.
    1636. 

  1636.                fluxm(m,n)=0.
    1637. 

  1637.                fm(m,n)=0.
    1638. 

  1638.             end do
    1639. 

  1639.             end do
    1640. 

  1640.             flux_fullact=0.
    1641. 

  1641.             return
    1642. 

  1642.          endif
    1643. 

  1643. 

  1644. 

  1645. 

  1646.          if (isectional .le. 0) then
    1647. 

  1647.             ! Initialize maxsat at this cell and timestep for the
    1648. 

  1648.             ! modal setup (the sectional case is handled below).
    1649. 

  1649.             call maxsat_init(maxd_atype, ntype_aer, &
    1650. 

  1650.                  maxd_asize, nsize_aer, alnsign, f1)
    1651. 

  1651. 

  1652.             goto 30000
    1653. 

  1653.          end if
    1654. 

  1654. 

  1655.          do n=1,ntype_aer
    1656. 

  1656.             !wig 19-Oct-2006: Add zero trap based May 2006 e-mail from
    1657. 

  1657.             !Ghan. Transport can clear out a cell leading to
    1658. 

  1658.             !inconsistencies with the mass.
    1659. 

  1659.             gmrad(n)=gmrad(n)/max(totn(n),1e-20)
    1660. 

  1660.             gmlnsig=gmradsq(n)/totn(n)-gmrad(n)*gmrad(n)    ! [ln(sigmag)]**2
    1661. 

  1661.             gmlnsig(n)=sqrt( max( 1.e-4, gmlnsig(n) ) )
    1662. 

  1662.             gmrad(n)=exp(gmrad(n))
    1663. 

  1663.             if ((isectional .eq. 3) .or. (isectional .eq. 4)) then
    1664. 

  1664.                gmsm(n)=totn(n)/sumns(n)
    1665. 

  1665.                gmsm(n)=gmsm(n)*gmsm(n)*gmsm(n)
    1666. 

  1666.                gmsm(n)=sqrt(gmsm(n))
    1667. 

  1667.             else
    1668. 

  1668. !              gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(1,n)*gmrad(n)))
    1669. 

  1669.                gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(nsize_aer(n),n)*gmrad(n)))
    1670. 

  1670.             endif
    1671. 

  1671.          enddo
    1672. 

  1672.          
    1673. 

  1673.          ! Initialize maxsat at this cell and timestep for the
    1674. 

  1674.          ! sectional setup (the modal case is handled above)...
    1675. 

  1675.          call maxsat_init(maxd_atype, ntype_aer, &
    1676. 

  1676.               maxd_asize, (/1/), gmlnsig, f1)
    1677. 

  1677. 

  1678. !.......................................................................
    1679. 

  1679. !   calculate sectional "sub-bin" size distribution
    1680. 

  1680. !
    1681. 

  1681. !   dn/dy = nt*( a + b*y )   for  ylo < y < yhi
    1682. 

  1682. !
    1683. 

  1683. !   nt = na(m,n) = number mixing ratio of the bin
    1684. 

  1684. !   y = v/vhi
    1685. 

  1685. !       v = (4pi/3)*r**3 = particle volume
    1686. 

  1686. !       vhi = v at r=rhi (upper bin boundary)
    1687. 

  1687. !   ylo = y at lower bin boundary = vlo/vhi = (rlo/rhi)**3
    1688. 

  1688. !   yhi = y at upper bin boundary = 1.0
    1689. 

  1689. !
    1690. 

  1690. !   dv/dy = v * dn/dy = nt*vhi*( a*y + b*y*y )
    1691. 

  1691. !
    1692. 

  1692. !.......................................................................
    1693. 

  1693. ! 02-may-2006 - this dn/dy replaces the previous
    1694. 

  1694. !       dn/dx = a + b*x   where l = ln(r)
    1695. 

  1695. !    the old dn/dx was overly complicated for cases of rmean near rlo or rhi
    1696. 

  1696. !    the new dn/dy is consistent with that used in the movesect routine,
    1697. 

  1697. !       which does continuous growth by condensation and aqueous chemistry
    1698. 

  1698. !.......................................................................
    1699. 

  1699.              do 25002 n = 1,ntype_aer
    1700. 

  1700.              do 25000 m = 1,nsize_aer(n)
    1701. 

  1701. 

  1702. ! convert from diameter in cm to radius in m
    1703. 

  1703.                 rlo(m,n) = 0.5*0.01*dlo_sect(m,n)
    1704. 

  1704.                 rhi(m,n) = 0.5*0.01*dhi_sect(m,n)
    1705. 

  1705.                 ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
    1706. 

  1706.                 yhi(m,n) = 1.0
    1707. 

  1707. 

  1708. ! 04-nov-2005 - extremely narrow bins will be treated using 0/1 activation
    1709. 

  1709. !    this is to avoid potential numerical problems
    1710. 

  1710.                 bin_is_narrow(m,n) = .false.
    1711. 

  1711.                 if ((rhi(m,n)/rlo(m,n)) .le. 1.01) bin_is_narrow(m,n) = .true.
    1712. 

  1712. 

  1713. ! rmean is mass mean radius for the bin; xmean = log(rmean)
    1714. 

  1714. ! just use section midpoint if bin is empty
    1715. 

  1715.                 if ( bin_is_empty(m,n) ) then
    1716. 

  1716.                    rmean(m,n) = sqrt(rlo(m,n)*rhi(m,n)) 
    1717. 

  1717.                    ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
    1718. 

  1718.                    goto 25000
    1719. 

  1719.                 end if
    1720. 

  1720. 

  1721.                 rmean(m,n) = (volc(m,n)/(ccc*na(m,n)))**third
    1722. 

  1722.                 rmean(m,n) = max( rlo(m,n), min( rhi(m,n), rmean(m,n) ) )
    1723. 

  1723.                 ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
    1724. 

  1724.                 if ( bin_is_narrow(m,n) ) goto 25000
    1725. 

  1725. 

  1726. ! if rmean is extremely close to either rlo or rhi, 
    1727. 

  1727. ! treat the bin as extremely narrow
    1728. 

  1728.                 if ((rhi(m,n)/rmean(m,n)) .le. 1.01) then
    1729. 

  1729.                    bin_is_narrow(m,n) = .true.
    1730. 

  1730.                    rlo(m,n) = min( rmean(m,n), (rhi(m,n)/1.01) )
    1731. 

  1731.                    ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
    1732. 

  1732.                    goto 25000
    1733. 

  1733.                 else if ((rmean(m,n)/rlo(m,n)) .le. 1.01) then
    1734. 

  1734.                    bin_is_narrow(m,n) = .true.
    1735. 

  1735.                    rhi(m,n) = max( rmean(m,n), (rlo(m,n)*1.01) )
    1736. 

  1736.                    ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
    1737. 

  1737.                    ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
    1738. 

  1738.                    goto 25000
    1739. 

  1739.                 endif
    1740. 

  1740. 

  1741. ! if rmean is somewhat close to either rlo or rhi, then dn/dy will be 
    1742. 

  1742. !    negative near the upper or lower bin boundary
    1743. 

  1743. ! in these cases, assume that all the particles are in a subset of the full bin,
    1744. 

  1744. !    and adjust rlo or rhi so that rmean will be near the center of this subset
    1745. 

  1745. ! note that the bin is made narrower LOCALLY/TEMPORARILY, 
    1746. 

  1746. !    just for the purposes of the activation calculation
    1747. 

  1747.                 gammayy = (ymean(m,n)-ylo(m,n)) / (yhi(m,n)-ylo(m,n))
    1748. 

  1748.                 if (gammayy .lt. 0.34) then
    1749. 

  1749.                    dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*(gammayy/0.34)
    1750. 

  1750.                    rhi(m,n) = rhi(m,n)*(dumaa**third)
    1751. 

  1751.                    ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
    1752. 

  1752.                    ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
    1753. 

  1753.                 else if (gammayy .ge. 0.66) then
    1754. 

  1754.                    dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*((gammayy-0.66)/0.34)
    1755. 

  1755.                    ylo(m,n) = dumaa
    1756. 

  1756.                    rlo(m,n) = rhi(m,n)*(dumaa**third)
    1757. 

  1757.                 end if
    1758. 

  1758.                 if ((rhi(m,n)/rlo(m,n)) .le. 1.01) then
    1759. 

  1759.                    bin_is_narrow(m,n) = .true.
    1760. 

  1760.                    goto 25000
    1761. 

  1761.                 end if
    1762. 

  1762. 

  1763.                 betayy = ylo(m,n)/yhi(m,n)
    1764. 

  1764.                 betayy2 = betayy*betayy
    1765. 

  1765.                 bsub(m,n) = (12.0*ymean(m,n) - 6.0*(1.0+betayy)) /   &
    1766. 

  1766.                    (4.0*(1.0-betayy2*betayy) - 3.0*(1.0-betayy2)*(1.0+betayy))
    1767. 

  1767.                 asub(m,n) = (1.0 - bsub(m,n)*(1.0-betayy2)*0.5) / (1.0-betayy)
    1768. 

  1768. 

  1769.                 if ( asub(m,n)+bsub(m,n)*ylo(m,n) .lt. 0. ) then
    1770. 

  1770.                   if (idiag_dndy_neg .gt. 0) then
    1771. 

  1771.                     print *,'dndy<0 at lower boundary'
    1772. 

  1772.                     print *,'n,m=',n,m
    1773. 

  1773.                     print *,'na=',na(m,n),' volc=',volc(m,n)
    1774. 

  1774.                     print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc))
    1775. 

  1775.                     print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n)
    1776. 

  1776.                     print *,'dlo_sect/2,dhi_sect/2=',   &
    1777. 

  1777.                              (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n))
    1778. 

  1778.                     print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n)
    1779. 

  1779.                     print *,'asub+bsub*ylo=',   &
    1780. 

  1780.                              (asub(m,n)+bsub(m,n)*ylo(m,n))
    1781. 

  1781.                     print *,'subr activate error 11 - i,j,k =', ii, jj, kk
    1782. 

  1782.                   endif
    1783. 

  1783.                 endif
    1784. 

  1784.                 if ( asub(m,n)+bsub(m,n)*yhi(m,n) .lt. 0. ) then
    1785. 

  1785.                   if (idiag_dndy_neg .gt. 0) then
    1786. 

  1786.                     print *,'dndy<0 at upper boundary'
    1787. 

  1787.                     print *,'n,m=',n,m
    1788. 

  1788.                     print *,'na=',na(m,n),' volc=',volc(m,n)
    1789. 

  1789.                     print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc))
    1790. 

  1790.                     print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n)
    1791. 

  1791.                     print *,'dlo_sect/2,dhi_sect/2=',   &
    1792. 

  1792.                              (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n))
    1793. 

  1793.                     print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n)
    1794. 

  1794.                     print *,'asub+bsub*yhi=',   &
    1795. 

  1795.                              (asub(m,n)+bsub(m,n)*yhi(m,n))
    1796. 

  1796.                     print *,'subr activate error 12 - i,j,k =', ii, jj, kk
    1797. 

  1797.                   endif
    1798. 

  1798.                 endif
    1799. 

  1799. 

  1800. 25000        continue      ! m=1,nsize_aer(n)
    1801. 

  1801. 25002        continue      ! n=1,ntype_aer
    1802. 

  1802. 

  1803. 

  1804. 30000    continue
    1805. 

  1805. !.......................................................................
    1806. 

  1806. !
    1807. 

  1807. !   end calc. of modal or sectional activation properties (end of section 1)
    1808. 

  1808. !
    1809. 

  1809. !.......................................................................
    1810. 

  1810. 

  1811. 

  1812. 

  1813. !      sjg 7-16-98  upward
    1814. 

  1814. !      print *,'wbar,sigw=',wbar,sigw
    1815. 

  1815. 

  1816.       if(sigw.le.1.e-5) goto 50000
    1817. 

  1817. 

  1818. !.......................................................................
    1819. 

  1819. !
    1820. 

  1820. !   start calc. of activation fractions/fluxes
    1821. 

  1821. !   for spectrum of updrafts (start of section 2)
    1822. 

  1822. !
    1823. 

  1823. !.......................................................................
    1824. 

  1824.          ipass_nwloop = 1
    1825. 

  1825.          idiagaa = 0
    1826. 

  1826. ! 06-nov-2005 rce - set idiagaa=1 for testing/debugging
    1827. 

  1827. !        if ((grid_id.eq.1) .and. (ktau.eq.167) .and.   &
    1828. 

  1828. !            (ii.eq.24) .and. (jj.eq. 1) .and. (kk.eq.14)) idiagaa = 1
    1829. 

  1829. 

  1830. 40000    continue
    1831. 

  1831.          if(top)then
    1832. 

  1832.            wmax=0.
    1833. 

  1833.            wmin=min(zero,-wdiab)
    1834. 

  1834.          else
    1835. 

  1835.            wmax=min(wmaxf,wbar+sds*sigw)
    1836. 

  1836.            wmin=max(wminf,-wdiab)
    1837. 

  1837.          endif
    1838. 

  1838.          wmin=max(wmin,wbar-sds*sigw)
    1839. 

  1839.          w=wmin
    1840. 

  1840.          dwmax=eps*sigw
    1841. 

  1841.          dw=dwmax
    1842. 

  1842.          dfmax=0.2
    1843. 

  1843.          dfmin=0.1
    1844. 

  1844.          if(wmax.le.w)then
    1845. 

  1845.             do n=1,ntype_aer
    1846. 

  1846.             do m=1,nsize_aer(n)
    1847. 

  1847.                fluxn(m,n)=0.
    1848. 

  1848.                fn(m,n)=0.
    1849. 

  1849.                fluxs(m,n)=0.
    1850. 

  1850.                fs(m,n)=0.
    1851. 

  1851.                fluxm(m,n)=0.
    1852. 

  1852.                fm(m,n)=0.
    1853. 

  1853.             end do
    1854. 

  1854.             end do
    1855. 

  1855.             flux_fullact=0.
    1856. 

  1856.             return
    1857. 

  1857.          endif
    1858. 

  1858.          do n=1,ntype_aer
    1859. 

  1859.          do m=1,nsize_aer(n)
    1860. 

  1860.             sumflxn(m,n)=0.
    1861. 

  1861.             sumfn(m,n)=0.
    1862. 

  1862.             fnold(m,n)=0.
    1863. 

  1863.             sumflxs(m,n)=0.
    1864. 

  1864.             sumfs(m,n)=0.
    1865. 

  1865.             fsold(m,n)=0.
    1866. 

  1866.             sumflxm(m,n)=0.
    1867. 

  1867.             sumfm(m,n)=0.
    1868. 

  1868.             fmold(m,n)=0.
    1869. 

  1869.          enddo
    1870. 

  1870.          enddo
    1871. 

  1871.          sumflx_fullact=0.
    1872. 

  1872. 

  1873.          fold=0
    1874. 

  1874.          gold=0
    1875. 

  1875. ! 06-nov-2005 rce - set wold=w here
    1876. 

  1876. !        wold=0
    1877. 

  1877.          wold=w
    1878. 

  1878. 

  1879. 

  1880. ! 06-nov-2005 rce - define nwmax; calc dwmin from nwmax
    1881. 

  1881.          nwmax = 200
    1882. 

  1882. !        dwmin = min( dwmax, 0.01 )
    1883. 

  1883.          dwmin = (wmax - wmin)/(nwmax-1)
    1884. 

  1884.          dwmin = min( dwmax, dwmin )
    1885. 

  1885.          dwmin = max( 0.01,  dwmin )
    1886. 

  1886. 

  1887. !
    1888. 

  1888. ! loop over updrafts, incrementing sums as you go
    1889. 

  1889. ! the "200" is (arbitrary) upper limit for number of updrafts
    1890. 

  1890. ! if integration finishes before this, OK; otherwise, ERROR
    1891. 

  1891. !
    1892. 

  1892.          if (idiagaa.gt.0) then
    1893. 

  1893.              write(*,94700) ktau, grid_id, ii, jj, kk, nwmax
    1894. 

  1894.              write(*,94710) 'wbar,sigw,wdiab=', wbar, sigw, wdiab
    1895. 

  1895.              write(*,94710) 'wmin,wmax,dwmin,dwmax=', wmin, wmax, dwmin, dwmax
    1896. 

  1896.              write(*,94720) -1, w, wold, dw
    1897. 

  1897.          end if
    1898. 

  1898. 94700    format( / 'activate 47000 - ktau,id,ii,jj,kk,nwmax=', 6i5 )
    1899. 

  1899. 94710    format( 'activate 47000 - ', a, 6(1x,f11.5) )
    1900. 

  1900. 94720    format( 'activate 47000 - nw,w,wold,dw=', i5, 3(1x,f11.5) )
    1901. 

  1901. 

  1902.          do 47000 nw = 1, nwmax
    1903. 

  1903. 41000       wnuc=w+wdiab
    1904. 

  1904. 

  1905.             if (idiagaa.gt.0) write(*,94720) nw, w, wold, dw
    1906. 

  1906. 

  1907. !           write(6,*)'wnuc=',wnuc
    1908. 

  1908.             alw=alpha*wnuc
    1909. 

  1909.             sqrtalw=sqrt(alw)
    1910. 

  1910.             zeta=2.*sqrtalw*aten/(3.*sqrtg)
    1911. 

  1911.             etafactor1=2.*alw*sqrtalw
    1912. 

  1912.             if (isectional .gt. 0) then
    1913. 

  1913. !              sectional model.
    1914. 

  1914. !              use bulk properties
    1915. 

  1915. 

  1916.               do n=1,ntype_aer
    1917. 

  1917.                  if(totn(n).gt.1.e-10)then
    1918. 

  1918.                     eta(1,n)=etafactor1/(totn(n)*beta*sqrtg)
    1919. 

  1919.                  else
    1920. 

  1920.                     eta(1,n)=1.e10
    1921. 

  1921.                  endif
    1922. 

  1922.               enddo
    1923. 

  1923.               call maxsat(zeta,eta,maxd_atype,ntype_aer, &
    1924. 

  1924.                    maxd_asize,(/1/),gmsm,gmlnsig,f1,smax)
    1925. 

  1925.               lnsmax=log(smax)
    1926. 

  1926.               x=2*(log(gmsm(1))-lnsmax)/(3*sq2*gmlnsig(1))
    1927. 

  1927.               fnew=0.5*(1.-ERF_ALT(x))
    1928. 

  1928. 

  1929.             else
    1930. 

  1930. 

  1931.               do n=1,ntype_aer
    1932. 

  1932.               do m=1,nsize_aer(n)
    1933. 

  1933.                  eta(m,n)=etafactor1*etafactor2(m,n)
    1934. 

  1934.               enddo
    1935. 

  1935.               enddo
    1936. 

  1936. 

  1937.               call maxsat(zeta,eta,maxd_atype,ntype_aer, &
    1938. 

  1938.                    maxd_asize,nsize_aer,sm,alnsign,f1,smax)
    1939. 

  1939. !             write(6,*)'w,smax=',w,smax
    1940. 

  1940. 

  1941.               lnsmax=log(smax)
    1942. 

  1942. 

  1943.               x=2*(lnsm(nsize_aer(1),1)-lnsmax)/(3*sq2*alnsign(nsize_aer(1),1))
    1944. 

  1944.               fnew=0.5*(1.-ERF_ALT(x))
    1945. 

  1945. 

  1946.             endif
    1947. 

  1947. 

  1948.             dwnew = dw
    1949. 

  1949. ! 06-nov-2005 rce - "n" here should be "nw" (?) 
    1950. 

  1950. !           if(fnew-fold.gt.dfmax.and.n.gt.1)then
    1951. 

  1951.             if(fnew-fold.gt.dfmax.and.nw.gt.1)then
    1952. 

  1952. !              reduce updraft increment for greater accuracy in integration
    1953. 

  1953.                if (dw .gt. 1.01*dwmin) then
    1954. 

  1954.                   dw=0.7*dw
    1955. 

  1955.                   dw=max(dw,dwmin)
    1956. 

  1956.                   w=wold+dw
    1957. 

  1957.                   go to 41000
    1958. 

  1958.                else
    1959. 

  1959.                   dwnew = dwmin
    1960. 

  1960.                endif
    1961. 

  1961.             endif
    1962. 

  1962. 

  1963.             if(fnew-fold.lt.dfmin)then
    1964. 

  1964. !              increase updraft increment to accelerate integration
    1965. 

  1965.                dwnew=min(1.5*dw,dwmax)
    1966. 

  1966.             endif
    1967. 

  1967.             fold=fnew
    1968. 

  1968. 

  1969.             z=(w-wbar)/(sigw*sq2)
    1970. 

  1970.             gaus=exp(-z*z)
    1971. 

  1971.             fnmin=1.
    1972. 

  1972.             xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
    1973. 

  1973. !           write(6,*)'xmincoeff=',xmincoeff
    1974. 

  1974. 

  1975. 

  1976.             do 44002 n=1,ntype_aer
    1977. 

  1977.             do 44000 m=1,nsize_aer(n)
    1978. 

  1978.                if ( bin_is_empty(m,n) ) then
    1979. 

  1979.                    fn(m,n)=0.
    1980. 

  1980.                    fs(m,n)=0.
    1981. 

  1981.                    fm(m,n)=0.
    1982. 

  1982.                else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then
    1983. 

  1983. !                 sectional
    1984. 

  1984. !                  within-section dn/dx = a + b*x
    1985. 

  1985.                   xcut=xmincoeff-third*lnhygro(m,n)
    1986. 

  1986. !                 ycut=(exp(xcut)/rhi(m,n))**3
    1987. 

  1987. ! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20
    1988. 

  1988. ! if (ycut > yhi), then actual value of ycut is unimportant, 
    1989. 

  1989. ! so do the following to avoid overflow
    1990. 

  1990.                   lnycut = 3.0 * ( xcut - log(rhi(m,n)) )
    1991. 

  1991.                   lnycut = min( lnycut, log(yhi(m,n)*1.0e5) )
    1992. 

  1992.                   ycut=exp(lnycut)
    1993. 

  1993. !                 write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n)
    1994. 

  1994. !                   if(lnsmax.lt.lnsmn(m,n))then
    1995. 

  1995.                   if(ycut.gt.yhi(m,n))then
    1996. 

  1996.                      fn(m,n)=0.
    1997. 

  1997.                      fs(m,n)=0.
    1998. 

  1998.                      fm(m,n)=0.
    1999. 

  1999.                   elseif(ycut.lt.ylo(m,n))then
    2000. 

  2000.                      fn(m,n)=1.
    2001. 

  2001.                      fs(m,n)=1.
    2002. 

  2002.                      fm(m,n)=1.
    2003. 

  2003.                   elseif ( bin_is_narrow(m,n) ) then
    2004. 

  2004. ! 04-nov-2005 rce - for extremely narrow bins, 
    2005. 

  2005. ! do zero activation if xcut>xmean, 100% activation otherwise
    2006. 

  2006.                      if (ycut.gt.ymean(m,n)) then
    2007. 

  2007.                         fn(m,n)=0.
    2008. 

  2008.                         fs(m,n)=0.
    2009. 

  2009.                         fm(m,n)=0.
    2010. 

  2010.                      else
    2011. 

  2011.                         fn(m,n)=1.
    2012. 

  2012.                         fs(m,n)=1.
    2013. 

  2013.                         fm(m,n)=1.
    2014. 

  2014.                      endif
    2015. 

  2015.                   else
    2016. 

  2016.                      phiyy=ycut/yhi(m,n)
    2017. 

  2017.                      fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy)
    2018. 

  2018.                      if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then
    2019. 

  2019.                       if (idiag_fnsm_prob .gt. 0) then
    2020. 

  2020.                         print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21'
    2021. 

  2021.                         print *,'na,volc       =', na(m,n), volc(m,n)
    2022. 

  2022.                         print *,'asub,bsub     =', asub(m,n), bsub(m,n)
    2023. 

  2023.                         print *,'yhi,ycut      =', yhi(m,n), ycut
    2024. 

  2024.                       endif
    2025. 

  2025.                      endif
    2026. 

  2026. 

  2027.                      if (fn(m,n) .le. zero) then
    2028. 

  2028. ! 10-nov-2005 rce - if fn=0, then fs & fm must be 0
    2029. 

  2029.                         fn(m,n)=zero
    2030. 

  2030.                         fs(m,n)=zero
    2031. 

  2031.                         fm(m,n)=zero
    2032. 

  2032.                      else if (fn(m,n) .ge. one) then
    2033. 

  2033. ! 10-nov-2005 rce - if fn=1, then fs & fm must be 1
    2034. 

  2034.                         fn(m,n)=one
    2035. 

  2035.                         fs(m,n)=one
    2036. 

  2036.                         fm(m,n)=one
    2037. 

  2037.                      else
    2038. 

  2038. ! 10-nov-2005 rce - otherwise, calc fm and check it
    2039. 

  2039.                         fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) +   &
    2040. 

  2040.                                   third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy))
    2041. 

  2041.                         if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then
    2042. 

  2042.                          if (idiag_fnsm_prob .gt. 0) then
    2043. 

  2043.                            print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22'
    2044. 

  2044.                            print *,'na,volc,fn    =', na(m,n), volc(m,n), fn(m,n)
    2045. 

  2045.                            print *,'asub,bsub     =', asub(m,n), bsub(m,n)
    2046. 

  2046.                            print *,'yhi,ycut     =', yhi(m,n), ycut
    2047. 

  2047.                          endif
    2048. 

  2048.                         endif
    2049. 

  2049.                         if (fm(m,n) .le. fn(m,n)) then
    2050. 

  2050. ! 10-nov-2005 rce - if fm=fn, then fs must =fn
    2051. 

  2051.                            fm(m,n)=fn(m,n)
    2052. 

  2052.                            fs(m,n)=fn(m,n)
    2053. 

  2053.                         else if (fm(m,n) .ge. one) then
    2054. 

  2054. ! 10-nov-2005 rce - if fm=1, then fs & fn must be 1
    2055. 

  2055.                            fm(m,n)=one
    2056. 

  2056.                            fs(m,n)=one
    2057. 

  2057.                            fn(m,n)=one
    2058. 

  2058.                         else
    2059. 

  2059. ! 10-nov-2005 rce - these two checks assure that the mean size
    2060. 

  2060. ! of the activated & interstitial particles will be between rlo & rhi
    2061. 

  2061.                            dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n)) 
    2062. 

  2062.                            fm(m,n) = min( fm(m,n), dumaa )
    2063. 

  2063.                            dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n)) 
    2064. 

  2064.                            fm(m,n) = min( fm(m,n), dumaa )
    2065. 

  2065. ! 10-nov-2005 rce - now calculate fs and bound it by fn, fm
    2066. 

  2066.                            betayy = ylo(m,n)/yhi(m,n)
    2067. 

  2067.                            dumaa = phiyy**twothird
    2068. 

  2068.                            dumbb = betayy**twothird
    2069. 

  2069.                            fs(m,n) =   &
    2070. 

  2070.                               (asub(m,n)*(1.0-phiyy*dumaa) +   &
    2071. 

  2071.                                   0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) /   &
    2072. 

  2072.                               (asub(m,n)*(1.0-betayy*dumbb) +   &
    2073. 

  2073.                                   0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb))
    2074. 

  2074.                            fs(m,n)=max(fs(m,n),fn(m,n))
    2075. 

  2075.                            fs(m,n)=min(fs(m,n),fm(m,n))
    2076. 

  2076.                         endif
    2077. 

  2077.                      endif
    2078. 

  2078.                   endif
    2079. 

  2079. 

  2080.                else
    2081. 

  2081. !                 modal
    2082. 

  2082.                   x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n))
    2083. 

  2083.                   fn(m,n)=0.5*(1.-ERF_ALT(x))
    2084. 

  2084.                   arg=x-sq2*alnsign(m,n)
    2085. 

  2085.                   fs(m,n)=0.5*(1.-ERF_ALT(arg))
    2086. 

  2086.                   arg=x-1.5*sq2*alnsign(m,n)
    2087. 

  2087.                   fm(m,n)=0.5*(1.-ERF_ALT(arg))
    2088. 

  2088. !                 print *,'w,x,fn,fs,fm=',w,x,fn(m,n),fs(m,n),fm(m,n)
    2089. 

  2089.                endif
    2090. 

  2090. 

  2091. !                     fn(m,n)=1.  !test
    2092. 

  2092. !                     fs(m,n)=1.
    2093. 

  2093. !                     fm(m,n)=1.
    2094. 

  2094.                fnmin=min(fn(m,n),fnmin)
    2095. 

  2095. !               integration is second order accurate
    2096. 

  2096. !               assumes linear variation of f*gaus with w
    2097. 

  2097.                wb=(w+wold)
    2098. 

  2098.                fnbar=(fn(m,n)*gaus+fnold(m,n)*gold)
    2099. 

  2099.                fsbar=(fs(m,n)*gaus+fsold(m,n)*gold)
    2100. 

  2100.                fmbar=(fm(m,n)*gaus+fmold(m,n)*gold)
    2101. 

  2101.                if((top.and.w.lt.0.).or.(.not.top.and.w.gt.0.))then
    2102. 

  2102.                   sumflxn(m,n)=sumflxn(m,n)+sixth*(wb*fnbar           &
    2103. 

  2103.                       +(fn(m,n)*gaus*w+fnold(m,n)*gold*wold))*dw
    2104. 

  2104.                   sumflxs(m,n)=sumflxs(m,n)+sixth*(wb*fsbar           &
    2105. 

  2105.                       +(fs(m,n)*gaus*w+fsold(m,n)*gold*wold))*dw
    2106. 

  2106.                   sumflxm(m,n)=sumflxm(m,n)+sixth*(wb*fmbar           &
    2107. 

  2107.                       +(fm(m,n)*gaus*w+fmold(m,n)*gold*wold))*dw
    2108. 

  2108.                endif
    2109. 

  2109.                sumfn(m,n)=sumfn(m,n)+0.5*fnbar*dw
    2110. 

  2110. !              write(6,'(a,9g10.2)')'lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw=', &
    2111. 

  2111. !                lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw
    2112. 

  2112.                fnold(m,n)=fn(m,n)
    2113. 

  2113.                sumfs(m,n)=sumfs(m,n)+0.5*fsbar*dw
    2114. 

  2114.                fsold(m,n)=fs(m,n)
    2115. 

  2115.                sumfm(m,n)=sumfm(m,n)+0.5*fmbar*dw
    2116. 

  2116.                fmold(m,n)=fm(m,n)
    2117. 

  2117. 

  2118. 44000       continue      ! m=1,nsize_aer(n)
    2119. 

  2119. 44002       continue      ! n=1,ntype_aer
    2120. 

  2120. 

  2121. !           same form as sumflxm(m,n) but replace the fm/fmold(m,n) with 1.0
    2122. 

  2122.             sumflx_fullact = sumflx_fullact &
    2123. 

  2123.                            + sixth*(wb*(gaus+gold) + (gaus*w + gold*wold))*dw
    2124. 

  2124. !            sumg=sumg+0.5*(gaus+gold)*dw
    2125. 

  2125.             gold=gaus
    2126. 

  2126.             wold=w
    2127. 

  2127.             dw=dwnew
    2128. 

  2128. 

  2129.             if(nw.gt.1.and.(w.gt.wmax.or.fnmin.gt.fmax))go to 48000
    2130. 

  2130.             w=w+dw
    2131. 

  2131. 

  2132. 47000    continue      ! nw = 1, nwmax
    2133. 

  2133. 

  2134. 

  2135.          print *,'do loop is too short in activate'
    2136. 

  2136.          print *,'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
    2137. 

  2137.          print *,'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
    2138. 

  2138.          print *,'wnuc=',wnuc
    2139. 

  2139.          do n=1,ntype_aer
    2140. 

  2140.             print *,'ntype=',n
    2141. 

  2141.             print *,'na=',(na(m,n),m=1,nsize_aer(n))
    2142. 

  2142.             print *,'fn=',(fn(m,n),m=1,nsize_aer(n))
    2143. 

  2143.          end do
    2144. 

  2144. !   dump all subr parameters to allow testing with standalone code
    2145. 

  2145. !   (build a driver that will read input and call activate)
    2146. 

  2146.          print *,'top,wbar,sigw,wdiab,tair,rhoair,ntype_aer='
    2147. 

  2147.          print *, top,wbar,sigw,wdiab,tair,rhoair,ntype_aer
    2148. 

  2148.          print *,'na='
    2149. 

  2149.          print *, na
    2150. 

  2150.          print *,'volc='
    2151. 

  2151.          print *, volc
    2152. 

  2152.          print *,'sigman='
    2153. 

  2153.          print *, sigman
    2154. 

  2154.          print *,'hygro='
    2155. 

  2155.          print *, hygro
    2156. 

  2156. 

  2157.          print *,'subr activate error 31 - i,j,k =', ii, jj, kk
    2158. 

  2158. ! 06-nov-2005 rce - if integration fails, repeat it once with additional diagnostics
    2159. 

  2159.          if (ipass_nwloop .eq. 1) then
    2160. 

  2160.              ipass_nwloop = 2
    2161. 

  2161.              idiagaa = 2
    2162. 

  2162.              goto 40000
    2163. 

  2163.          end if
    2164. 

  2164.          call wrf_error_fatal("STOP: activate before 48000")
    2165. 

  2165. 

  2166. 48000    continue
    2167. 

  2167. 

  2168. 

  2169. !         ndist(n)=ndist(n)+1
    2170. 

  2170.          if(.not.top.and.w.lt.wmaxf)then
    2171. 

  2171. 

  2172. !            contribution from all updrafts stronger than wmax
    2173. 

  2173. !            assuming constant f (close to fmax)
    2174. 

  2174.             wnuc=w+wdiab
    2175. 

  2175. 

  2176.             z1=(w-wbar)/(sigw*sq2)
    2177. 

  2177.             z2=(wmaxf-wbar)/(sigw*sq2)
    2178. 

  2178.             integ=sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(z1)-ERFC_NUM_RECIPES(z2))
    2179. 

  2179. !            consider only upward flow into cloud base when estimating flux
    2180. 

  2180.             wf1=max(w,zero)
    2181. 

  2181.             zf1=(wf1-wbar)/(sigw*sq2)
    2182. 

  2182.             gf1=exp(-zf1*zf1)
    2183. 

  2183.             wf2=max(wmaxf,zero)
    2184. 

  2184.             zf2=(wf2-wbar)/(sigw*sq2)
    2185. 

  2185.             gf2=exp(-zf2*zf2)
    2186. 

  2186.             gf=(gf1-gf2)
    2187. 

  2187.             integf=wbar*sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(zf1)-ERFC_NUM_RECIPES(zf2))+sigw*sigw*gf
    2188. 

  2188. 

  2189.             do n=1,ntype_aer
    2190. 

  2190.             do m=1,nsize_aer(n)
    2191. 

  2191.                sumflxn(m,n)=sumflxn(m,n)+integf*fn(m,n)
    2192. 

  2192.                sumfn(m,n)=sumfn(m,n)+fn(m,n)*integ
    2193. 

  2193.                sumflxs(m,n)=sumflxs(m,n)+integf*fs(m,n)
    2194. 

  2194.                sumfs(m,n)=sumfs(m,n)+fs(m,n)*integ
    2195. 

  2195.                sumflxm(m,n)=sumflxm(m,n)+integf*fm(m,n)
    2196. 

  2196.                sumfm(m,n)=sumfm(m,n)+fm(m,n)*integ
    2197. 

  2197.             end do
    2198. 

  2198.             end do
    2199. 

  2199. !           same form as sumflxm(m,n) but replace the fm(m,n) with 1.0
    2200. 

  2200.             sumflx_fullact = sumflx_fullact + integf
    2201. 

  2201. !            sumg=sumg+integ
    2202. 

  2202.          endif
    2203. 

  2203. 

  2204. 

  2205.          do n=1,ntype_aer
    2206. 

  2206.          do m=1,nsize_aer(n)
    2207. 

  2207. 

  2208. !           fn(m,n)=sumfn(m,n)/(sumg)
    2209. 

  2209.             fn(m,n)=sumfn(m,n)/(sq2*sqpi*sigw)
    2210. 

  2210.             fluxn(m,n)=sumflxn(m,n)/(sq2*sqpi*sigw)
    2211. 

  2211.             if(fn(m,n).gt.1.01)then
    2212. 

  2212.              if (idiag_fnsm_prob .gt. 0) then
    2213. 

  2213.                print *,'fn=',fn(m,n),' > 1 - activate err41'
    2214. 

  2214.                print *,'w,m,n,na,am=',w,m,n,na(m,n),am(m,n)
    2215. 

  2215.                print *,'integ,sumfn,sigw=',integ,sumfn(m,n),sigw
    2216. 

  2216.                print *,'subr activate error - i,j,k =', ii, jj, kk
    2217. 

  2217.              endif
    2218. 

  2218.              fluxn(m,n) = fluxn(m,n)/fn(m,n)
    2219. 

  2219.             endif
    2220. 

  2220. 

  2221.             fs(m,n)=sumfs(m,n)/(sq2*sqpi*sigw)
    2222. 

  2222.             fluxs(m,n)=sumflxs(m,n)/(sq2*sqpi*sigw)
    2223. 

  2223.             if(fs(m,n).gt.1.01)then
    2224. 

  2224.              if (idiag_fnsm_prob .gt. 0) then
    2225. 

  2225.                print *,'fs=',fs(m,n),' > 1 - activate err42'
    2226. 

  2226.                print *,'m,n,isectional=',m,n,isectional
    2227. 

  2227.                print *,'alnsign(m,n)=',alnsign(m,n)
    2228. 

  2228.                print *,'rcut,rlo(m,n),rhi(m,n)',exp(xcut),rlo(m,n),rhi(m,n)
    2229. 

  2229.                print *,'w,m,na,am=',w,m,na(m,n),am(m,n)
    2230. 

  2230.                print *,'integ,sumfs,sigw=',integ,sumfs(m,n),sigw
    2231. 

  2231.              endif
    2232. 

  2232.              fluxs(m,n) = fluxs(m,n)/fs(m,n)
    2233. 

  2233.             endif
    2234. 

  2234. 

  2235. !           fm(m,n)=sumfm(m,n)/(sumg)
    2236. 

  2236.             fm(m,n)=sumfm(m,n)/(sq2*sqpi*sigw)
    2237. 

  2237.             fluxm(m,n)=sumflxm(m,n)/(sq2*sqpi*sigw)
    2238. 

  2238.             if(fm(m,n).gt.1.01)then
    2239. 

  2239.              if (idiag_fnsm_prob .gt. 0) then
    2240. 

  2240.                print *,'fm(',m,n,')=',fm(m,n),' > 1 - activate err43'
    2241. 

  2241.              endif
    2242. 

  2242.              fluxm(m,n) = fluxm(m,n)/fm(m,n)
    2243. 

  2243.             endif
    2244. 

  2244. 

  2245.          end do
    2246. 

  2246.          end do
    2247. 

  2247. !        same form as fluxm(m,n)
    2248. 

  2248.          flux_fullact = sumflx_fullact/(sq2*sqpi*sigw)
    2249. 

  2249. 

  2250.       goto 60000
    2251. 

  2251. !.......................................................................
    2252. 

  2252. !
    2253. 

  2253. !   end calc. of activation fractions/fluxes
    2254. 

  2254. !   for spectrum of updrafts (end of section 2)
    2255. 

  2255. !
    2256. 

  2256. !.......................................................................
    2257. 

  2257. 

  2258. !.......................................................................
    2259. 

  2259. !
    2260. 

  2260. !   start calc. of activation fractions/fluxes
    2261. 

  2261. !   for (single) uniform updraft (start of section 3)
    2262. 

  2262. !
    2263. 

  2263. !.......................................................................
    2264. 

  2264. 50000 continue
    2265. 

  2265. 

  2266.          wnuc=wbar+wdiab
    2267. 

  2267. !         write(6,*)'uniform updraft =',wnuc
    2268. 

  2268. 

  2269. ! 04-nov-2005 rce - moved the code for "wnuc.le.0" code to here
    2270. 

  2270.          if(wnuc.le.0.)then
    2271. 

  2271.             do n=1,ntype_aer
    2272. 

  2272.             do m=1,nsize_aer(n)
    2273. 

  2273.                fn(m,n)=0
    2274. 

  2274.                fluxn(m,n)=0
    2275. 

  2275.                fs(m,n)=0
    2276. 

  2276.                fluxs(m,n)=0
    2277. 

  2277.                fm(m,n)=0
    2278. 

  2278.                fluxm(m,n)=0
    2279. 

  2279.             end do
    2280. 

  2280.             end do
    2281. 

  2281.             flux_fullact=0.
    2282. 

  2282.             return
    2283. 

  2283.          endif
    2284. 

  2284. 

  2285.             w=wbar
    2286. 

  2286.             alw=alpha*wnuc
    2287. 

  2287.             sqrtalw=sqrt(alw)
    2288. 

  2288.             zeta=2.*sqrtalw*aten/(3.*sqrtg)
    2289. 

  2289. 

  2290.             if (isectional .gt. 0) then
    2291. 

  2291. !              sectional model.
    2292. 

  2292. !              use bulk properties
    2293. 

  2293.               do n=1,ntype_aer
    2294. 

  2294.               if(totn(n).gt.1.e-10)then
    2295. 

  2295.                  eta(1,n)=2*alw*sqrtalw/(totn(n)*beta*sqrtg)
    2296. 

  2296.               else
    2297. 

  2297.                  eta(1,n)=1.e10
    2298. 

  2298.               endif
    2299. 

  2299.               end do
    2300. 

  2300.                call maxsat(zeta,eta,maxd_atype,ntype_aer, &
    2301. 

  2301.                     maxd_asize,(/1/),gmsm,gmlnsig,f1,smax)
    2302. 

  2302. 

  2303.             else
    2304. 

  2304. 

  2305.               do n=1,ntype_aer
    2306. 

  2306.               do m=1,nsize_aer(n)
    2307. 

  2307.                  if(na(m,n).gt.1.e-10)then
    2308. 

  2308.                     eta(m,n)=2*alw*sqrtalw/(na(m,n)*beta*sqrtg)
    2309. 

  2309.                  else
    2310. 

  2310.                     eta(m,n)=1.e10
    2311. 

  2311.                  endif
    2312. 

  2312.               end do
    2313. 

  2313.               end do
    2314. 

  2314. 

  2315.               call maxsat(zeta,eta,maxd_atype,ntype_aer, &
    2316. 

  2316.                    maxd_asize,nsize_aer,sm,alnsign,f1,smax)
    2317. 

  2317. 

  2318.             endif
    2319. 

  2319. 

  2320.             lnsmax=log(smax)
    2321. 

  2321.             xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
    2322. 

  2322. 

  2323.             do 55002 n=1,ntype_aer
    2324. 

  2324.             do 55000 m=1,nsize_aer(n)
    2325. 

  2325. 

  2326. ! 04-nov-2005 rce - check for bin_is_empty here too, just like earlier
    2327. 

  2327.                if ( bin_is_empty(m,n) ) then
    2328. 

  2328.                    fn(m,n)=0.
    2329. 

  2329.                    fs(m,n)=0.
    2330. 

  2330.                    fm(m,n)=0.
    2331. 

  2331. 

  2332.                else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then
    2333. 

  2333. !                 sectional
    2334. 

  2334. !                  within-section dn/dx = a + b*x
    2335. 

  2335.                   xcut=xmincoeff-third*lnhygro(m,n)
    2336. 

  2336. !                 ycut=(exp(xcut)/rhi(m,n))**3
    2337. 

  2337. ! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20
    2338. 

  2338. ! if (ycut > yhi), then actual value of ycut is unimportant, 
    2339. 

  2339. ! so do the following to avoid overflow
    2340. 

  2340.                   lnycut = 3.0 * ( xcut - log(rhi(m,n)) )
    2341. 

  2341.                   lnycut = min( lnycut, log(yhi(m,n)*1.0e5) )
    2342. 

  2342.                   ycut=exp(lnycut)
    2343. 

  2343. !                 write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n)
    2344. 

  2344. !                   if(lnsmax.lt.lnsmn(m,n))then
    2345. 

  2345.                   if(ycut.gt.yhi(m,n))then
    2346. 

  2346.                      fn(m,n)=0.
    2347. 

  2347.                      fs(m,n)=0.
    2348. 

  2348.                      fm(m,n)=0.
    2349. 

  2349. !                   elseif(lnsmax.gt.lnsmx(m,n))then
    2350. 

  2350.                   elseif(ycut.lt.ylo(m,n))then
    2351. 

  2351.                      fn(m,n)=1.
    2352. 

  2352.                      fs(m,n)=1.
    2353. 

  2353.                      fm(m,n)=1.
    2354. 

  2354.                   elseif ( bin_is_narrow(m,n) ) then
    2355. 

  2355. ! 04-nov-2005 rce - for extremely narrow bins, 
    2356. 

  2356. ! do zero activation if xcut>xmean, 100% activation otherwise
    2357. 

  2357.                      if (ycut.gt.ymean(m,n)) then
    2358. 

  2358.                         fn(m,n)=0.
    2359. 

  2359.                         fs(m,n)=0.
    2360. 

  2360.                         fm(m,n)=0.
    2361. 

  2361.                      else
    2362. 

  2362.                         fn(m,n)=1.
    2363. 

  2363.                         fs(m,n)=1.
    2364. 

  2364.                         fm(m,n)=1.
    2365. 

  2365.                      endif
    2366. 

  2366.                   else
    2367. 

  2367.                      phiyy=ycut/yhi(m,n)
    2368. 

  2368.                      fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy)
    2369. 

  2369.                      if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then
    2370. 

  2370.                       if (idiag_fnsm_prob .gt. 0) then
    2371. 

  2371.                         print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21'
    2372. 

  2372.                         print *,'na,volc       =', na(m,n), volc(m,n)
    2373. 

  2373.                         print *,'asub,bsub     =', asub(m,n), bsub(m,n)
    2374. 

  2374.                         print *,'yhi,ycut      =', yhi(m,n), ycut
    2375. 

  2375.                       endif
    2376. 

  2376.                      endif
    2377. 

  2377. 

  2378.                      if (fn(m,n) .le. zero) then
    2379. 

  2379. ! 10-nov-2005 rce - if fn=0, then fs & fm must be 0
    2380. 

  2380.                         fn(m,n)=zero
    2381. 

  2381.                         fs(m,n)=zero
    2382. 

  2382.                         fm(m,n)=zero
    2383. 

  2383.                      else if (fn(m,n) .ge. one) then
    2384. 

  2384. ! 10-nov-2005 rce - if fn=1, then fs & fm must be 1
    2385. 

  2385.                         fn(m,n)=one
    2386. 

  2386.                         fs(m,n)=one
    2387. 

  2387.                         fm(m,n)=one
    2388. 

  2388.                      else
    2389. 

  2389. ! 10-nov-2005 rce - otherwise, calc fm and check it
    2390. 

  2390.                         fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) +   &
    2391. 

  2391.                                   third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy))
    2392. 

  2392.                         if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then
    2393. 

  2393.                          if (idiag_fnsm_prob .gt. 0) then
    2394. 

  2394.                            print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22'
    2395. 

  2395.                            print *,'na,volc,fn    =', na(m,n), volc(m,n), fn(m,n)
    2396. 

  2396.                            print *,'asub,bsub     =', asub(m,n), bsub(m,n)
    2397. 

  2397.                            print *,'yhi,ycut      =', yhi(m,n), ycut
    2398. 

  2398.                          endif
    2399. 

  2399.                         endif
    2400. 

  2400.                         if (fm(m,n) .le. fn(m,n)) then
    2401. 

  2401. ! 10-nov-2005 rce - if fm=fn, then fs must =fn
    2402. 

  2402.                            fm(m,n)=fn(m,n)
    2403. 

  2403.                            fs(m,n)=fn(m,n)
    2404. 

  2404.                         else if (fm(m,n) .ge. one) then
    2405. 

  2405. ! 10-nov-2005 rce - if fm=1, then fs & fn must be 1
    2406. 

  2406.                            fm(m,n)=one
    2407. 

  2407.                            fs(m,n)=one
    2408. 

  2408.                            fn(m,n)=one
    2409. 

  2409.                         else
    2410. 

  2410. ! 10-nov-2005 rce - these two checks assure that the mean size
    2411. 

  2411. ! of the activated & interstitial particles will be between rlo & rhi
    2412. 

  2412.                            dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n)) 
    2413. 

  2413.                            fm(m,n) = min( fm(m,n), dumaa )
    2414. 

  2414.                            dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n))
    2415. 

  2415.                            fm(m,n) = min( fm(m,n), dumaa )
    2416. 

  2416. ! 10-nov-2005 rce - now calculate fs and bound it by fn, fm
    2417. 

  2417.                            betayy = ylo(m,n)/yhi(m,n)
    2418. 

  2418.                            dumaa = phiyy**twothird
    2419. 

  2419.                            dumbb = betayy**twothird
    2420. 

  2420.                            fs(m,n) =   &
    2421. 

  2421.                               (asub(m,n)*(1.0-phiyy*dumaa) +   &
    2422. 

  2422.                                   0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) /   &
    2423. 

  2423.                               (asub(m,n)*(1.0-betayy*dumbb) +   &
    2424. 

  2424.                                   0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb))
    2425. 

  2425.                            fs(m,n)=max(fs(m,n),fn(m,n))
    2426. 

  2426.                            fs(m,n)=min(fs(m,n),fm(m,n))
    2427. 

  2427.                         endif
    2428. 

  2428.                      endif
    2429. 

  2429. 

  2430.                   endif
    2431. 

  2431. 

  2432.                else
    2433. 

  2433. !                 modal
    2434. 

  2434.                   x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n))
    2435. 

  2435.                   fn(m,n)=0.5*(1.-ERF_ALT(x))
    2436. 

  2436.                   arg=x-sq2*alnsign(m,n)
    2437. 

  2437.                   fs(m,n)=0.5*(1.-ERF_ALT(arg))
    2438. 

  2438.                   arg=x-1.5*sq2*alnsign(m,n)
    2439. 

  2439.                   fm(m,n)=0.5*(1.-ERF_ALT(arg))
    2440. 

  2440.                endif
    2441. 

  2441. 

  2442. !                     fn(m,n)=1. ! test
    2443. 

  2443. !                     fs(m,n)=1.
    2444. 

  2444. !                     fm(m,n)=1.
    2445. 

  2445.                 if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then
    2446. 

  2446.                    fluxn(m,n)=fn(m,n)*w
    2447. 

  2447.                    fluxs(m,n)=fs(m,n)*w
    2448. 

  2448.                    fluxm(m,n)=fm(m,n)*w
    2449. 

  2449.                 else
    2450. 

  2450.                    fluxn(m,n)=0
    2451. 

  2451.                    fluxs(m,n)=0
    2452. 

  2452.                    fluxm(m,n)=0
    2453. 

  2453.                endif
    2454. 

  2454. 

  2455. 55000       continue      ! m=1,nsize_aer(n)
    2456. 

  2456. 55002       continue      ! n=1,ntype_aer
    2457. 

  2457. 

  2458.             if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then
    2459. 

  2459.                flux_fullact = w
    2460. 

  2460.             else
    2461. 

  2461.                flux_fullact = 0.0
    2462. 

  2462.             endif
    2463. 

  2463. 

  2464. ! 04-nov-2005 rce - moved the code for "wnuc.le.0" from here 
    2465. 

  2465. ! to near the start the uniform undraft section
    2466. 

  2466. 

  2467. !.......................................................................
    2468. 

  2468. !
    2469. 

  2469. !   end calc. of activation fractions/fluxes 
    2470. 

  2470. !   for (single) uniform updraft (end of section 3)
    2471. 

  2471. !
    2472. 

  2472. !.......................................................................
    2473. 

  2473. 

  2474. 

  2475. 

  2476. 60000 continue
    2477. 

  2477. 

  2478. 

  2479. !            do n=1,ntype_aer
    2480. 

  2480. !            do m=1,nsize_aer(n)
    2481. 

  2481. !                write(6,'(a,2i3,5e10.1)')'n,m,na,wbar,sigw,fn,fm=',n,m,na(m,n),wbar,sigw,fn(m,n),fm(m,n)
    2482. 

  2482. !            end do
    2483. 

  2483. !            end do
    2484. 

  2484. 

  2485. 

  2486.       return
    2487. 

  2487.       end subroutine activate
    2488. 

  2488. 

  2489. 

  2490. 

  2491. !----------------------------------------------------------------------
    2492. 

  2492. !----------------------------------------------------------------------
    2493. 

  2493.       subroutine maxsat(zeta,eta, &
    2494. 

  2494.                         maxd_atype,ntype_aer,maxd_asize,nsize_aer, &
    2495. 

  2495.                         sm,alnsign,f1,smax)
    2496. 

  2496. 

  2497. !      Calculates maximum supersaturation for multiple competing aerosol
    2498. 

  2498. !      modes. Note that maxsat_init must be called before calling this
    2499. 

  2499. !      subroutine.
    2500. 

  2500. 

  2501. !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
    2502. 

  2502. !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
    2503. 

  2503. 

  2504.       implicit none
    2505. 

  2505. 

  2506.       integer, intent(in) :: maxd_atype
    2507. 

  2507.       integer, intent(in) :: ntype_aer
    2508. 

  2508.       integer, intent(in) :: maxd_asize
    2509. 

  2509.       integer, intent(in) :: nsize_aer(maxd_atype) ! number of size bins
    2510. 

  2510.       real, intent(in) :: sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius
    2511. 

  2511.       real, intent(in) :: zeta, eta(maxd_asize,maxd_atype)
    2512. 

  2512.       real, intent(in) :: alnsign(maxd_asize,maxd_atype) ! ln(sigma)
    2513. 

  2513.       real, intent(in) :: f1(maxd_asize,maxd_atype)
    2514. 

  2514.       real, intent(out) :: smax ! maximum supersaturation
    2515. 

  2515. 

  2516.       real :: g1, g2
    2517. 

  2517.       real thesum
    2518. 

  2518.       integer m ! size index
    2519. 

  2519.       integer n ! type index
    2520. 

  2520. 

  2521.       do n=1,ntype_aer
    2522. 

  2522.       do m=1,nsize_aer(n)
    2523. 

  2523.          if(zeta.gt.1.e5*eta(m,n) .or. &
    2524. 

  2524.               sm(m,n)*sm(m,n).gt.1.e5*eta(m,n))then
    2525. 

  2525. !           weak forcing. essentially none activated
    2526. 

  2526.             smax=1.e-20
    2527. 

  2527.          else
    2528. 

  2528. !           significant activation of this mode. calc activation all modes.
    2529. 

  2529.             go to 1
    2530. 

  2530.          endif
    2531. 

  2531.       end do
    2532. 

  2532.       end do
    2533. 

  2533. 

  2534.       return
    2535. 

  2535. 

  2536.   1   continue
    2537. 

  2537. 

  2538.       thesum=0
    2539. 

  2539.       do n=1,ntype_aer
    2540. 

  2540.       do m=1,nsize_aer(n)
    2541. 

  2541.          if(eta(m,n).gt.1.e-20)then
    2542. 

  2542.             g1=sqrt(zeta/eta(m,n))
    2543. 

  2543.             g1=g1*g1*g1
    2544. 

  2544.             g2=sm(m,n)/sqrt(eta(m,n)+3*zeta)
    2545. 

  2545.             g2=sqrt(g2)
    2546. 

  2546.             g2=g2*g2*g2
    2547. 

  2547.             thesum=thesum + &
    2548. 

  2548.                  (f1(m,n)*g1+(1.+0.25*alnsign(m,n))*g2)/(sm(m,n)*sm(m,n))
    2549. 

  2549.          else
    2550. 

  2550.             thesum=1.e20
    2551. 

  2551.          endif
    2552. 

  2552.       end do
    2553. 

  2553.       end do
    2554. 

  2554. 

  2555.       smax=1./sqrt(thesum)
    2556. 

  2556. 

  2557.       return
    2558. 

  2558.       end subroutine maxsat
    2559. 

  2559. 

  2560. 

  2561. 

  2562. !----------------------------------------------------------------------
    2563. 

  2563. !----------------------------------------------------------------------
    2564. 

  2564.       subroutine maxsat_init(maxd_atype, ntype_aer, &
    2565. 

  2565.            maxd_asize, nsize_aer, alnsign, f1)
    2566. 

  2566. 

  2567. !     Calculates the f1 paramter needed by maxsat.
    2568. 

  2568. 

  2569. !     Abdul-Razzak and Ghan, A parameterization of aerosol activation.
    2570. 

  2570. !     2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
    2571. 

  2571. 

  2572.       implicit none
    2573. 

  2573. 

  2574.       integer, intent(in)  :: maxd_atype
    2575. 

  2575.       integer, intent(in)  :: ntype_aer ! number of aerosol types
    2576. 

  2576.       integer, intent(in)  :: maxd_asize
    2577. 

  2577.       integer, intent(in)  :: nsize_aer(maxd_atype) ! number of size bins
    2578. 

  2578.       real,    intent(in)  :: alnsign(maxd_asize,maxd_atype) ! ln(sigma)
    2579. 

  2579.       real,    intent(out) :: f1(maxd_asize,maxd_atype)
    2580. 

  2580. 

  2581.       integer m ! size index
    2582. 

  2582.       integer n ! type index
    2583. 

  2583. 

  2584. !  calculate and save f1(sigma), assumes sigma is invariant
    2585. 

  2585. !  between calls to this init routine
    2586. 

  2586. 

  2587.       do n=1,ntype_aer
    2588. 

  2588.          do m=1,nsize_aer(n)
    2589. 

  2589.             f1(m,n)=0.5*exp(2.5*alnsign(m,n)*alnsign(m,n))
    2590. 

  2590.          end do
    2591. 

  2591.       end do
    2592. 

  2592. 

  2593.       end subroutine maxsat_init
    2594. 

  2594. 

  2595. 

  2596. 

  2597. !----------------------------------------------------------------------
    2598. 

  2598. !----------------------------------------------------------------------
    2599. 

  2599. ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3);
    2600. 

  2600. !     grid_id, ktau, i, j, isize, itype added to arg list to assist debugging
    2601. 

  2601.        subroutine loadaer(chem,k,kmn,kmx,num_chem,cs,npv, &
    2602. 

  2602.                           dlo_sect,dhi_sect,maxd_acomp, ncomp,                &
    2603. 

  2603.                           grid_id, ktau, i, j, isize, itype,   &
    2604. 

  2604.                           numptr_aer, numptrcw_aer, dens_aer,   &
    2605. 

  2605.                           massptr_aer, massptrcw_aer,   &
    2606. 

  2606.                           maerosol, maerosolcw,                 &
    2607. 

  2607.                           maerosol_tot, maerosol_totcw,         &
    2608. 

  2608.                           naerosol, naerosolcw,                 &
    2609. 

  2609.                           vaerosol, vaerosolcw)
    2610. 

  2610. 

  2611.       implicit none
    2612. 

  2612. 

  2613. !      load aerosol number, surface, mass concentrations
    2614. 

  2614. 

  2615. !      input
    2616. 

  2616. 

  2617.        integer, intent(in) ::  num_chem ! maximum number of consituents
    2618. 

  2618.        integer, intent(in) ::  k,kmn,kmx
    2619. 

  2619.        real,    intent(in) ::  chem(kmn:kmx,num_chem) ! aerosol mass, number mixing ratios
    2620. 

  2620.        real,    intent(in) ::  cs  ! air density (kg/m3)
    2621. 

  2621.        real,    intent(in) ::  npv ! number per volume concentration (/m3)
    2622. 

  2622.        integer, intent(in) ::  maxd_acomp,ncomp
    2623. 

  2623.        integer, intent(in) ::  numptr_aer,numptrcw_aer
    2624. 

  2624.        integer, intent(in) ::  massptr_aer(maxd_acomp), massptrcw_aer(maxd_acomp)
    2625. 

  2625.        real,    intent(in) ::  dens_aer(maxd_acomp) ! aerosol material density (g/cm3)
    2626. 

  2626.        real,    intent(in) ::  dlo_sect,dhi_sect ! minimum, maximum diameter of section (cm)
    2627. 

  2627.        integer, intent(in) ::  grid_id, ktau, i, j, isize, itype
    2628. 

  2628. 

  2629. !      output
    2630. 

  2630. 

  2631.        real, intent(out) ::  naerosol                ! interstitial number conc (/m3)
    2632. 

  2632.        real, intent(out) ::  naerosolcw              ! activated    number conc (/m3)
    2633. 

  2633.        real, intent(out) ::  maerosol(maxd_acomp)   ! interstitial mass conc (kg/m3)
    2634. 

  2634.        real, intent(out) ::  maerosolcw(maxd_acomp) ! activated    mass conc (kg/m3)
    2635. 

  2635.        real, intent(out) ::  maerosol_tot   ! total-over-species interstitial mass conc (kg/m3)
    2636. 

  2636.        real, intent(out) ::  maerosol_totcw ! total-over-species activated    mass conc (kg/m3)
    2637. 

  2637.        real, intent(out) ::  vaerosol       ! interstitial volume conc (m3/m3)
    2638. 

  2638.        real, intent(out) ::  vaerosolcw     ! activated volume conc (m3/m3)
    2639. 

  2639. 

  2640. !      internal
    2641. 

  2641. 

  2642.        integer lnum,lnumcw,l,ltype,lmass,lmasscw,lsfc,lsfccw
    2643. 

  2643.        real num_at_dhi, num_at_dlo
    2644. 

  2644.        real npv_at_dhi, npv_at_dlo
    2645. 

  2645.        real, parameter :: pi = 3.1415926526
    2646. 

  2646.        real specvol ! inverse aerosol material density (m3/kg)
    2647. 

  2647. 

  2648.           lnum=numptr_aer
    2649. 

  2649.           lnumcw=numptrcw_aer
    2650. 

  2650.           maerosol_tot=0.
    2651. 

  2651.           maerosol_totcw=0.
    2652. 

  2652.           vaerosol=0.
    2653. 

  2653.           vaerosolcw=0.
    2654. 

  2654.           do l=1,ncomp
    2655. 

  2655.              lmass=massptr_aer(l)
    2656. 

  2656.              lmasscw=massptrcw_aer(l)
    2657. 

  2657.              maerosol(l)=chem(k,lmass)*cs
    2658. 

  2658.              maerosol(l)=max(maerosol(l),0.)
    2659. 

  2659.              maerosolcw(l)=chem(k,lmasscw)*cs
    2660. 

  2660.              maerosolcw(l)=max(maerosolcw(l),0.)
    2661. 

  2661.              maerosol_tot=maerosol_tot+maerosol(l)
    2662. 

  2662.              maerosol_totcw=maerosol_totcw+maerosolcw(l)
    2663. 

  2663. ! [ 1.e-3 factor because dens_aer is (g/cm3), specvol is (m3/kg) ]
    2664. 

  2664.              specvol=1.0e-3/dens_aer(l)
    2665. 

  2665.              vaerosol=vaerosol+maerosol(l)*specvol
    2666. 

  2666.              vaerosolcw=vaerosolcw+maerosolcw(l)*specvol
    2667. 

  2667. !            write(6,'(a,3e12.2)')'maerosol,dens_aer,vaerosol=',maerosol(l),dens_aer(l),vaerosol
    2668. 

  2668.           enddo
    2669. 

  2669. 

  2670.           if(lnum.gt.0)then
    2671. 

  2671. !            aerosol number predicted
    2672. 

  2672. ! [ 1.0e6 factor because because dhi_ & dlo_sect are (cm), vaerosol is (m3) ]
    2673. 

  2673.              npv_at_dhi = 6.0e6/(pi*dhi_sect*dhi_sect*dhi_sect)
    2674. 

  2674.              npv_at_dlo = 6.0e6/(pi*dlo_sect*dlo_sect*dlo_sect)
    2675. 

  2675. 

  2676.              naerosol=chem(k,lnum)*cs
    2677. 

  2677.              naerosolcw=chem(k,lnumcw)*cs
    2678. 

  2678.              num_at_dhi = vaerosol*npv_at_dhi
    2679. 

  2679.              num_at_dlo = vaerosol*npv_at_dlo
    2680. 

  2680.              naerosol = max( num_at_dhi, min( num_at_dlo, naerosol ) )
    2681. 

  2681. 

  2682. !            write(6,'(a,5e10.1)')'naerosol,num_at_dhi,num_at_dlo,dhi_sect,dlo_sect', &
    2683. 

  2683. !                          naerosol,num_at_dhi,num_at_dlo,dhi_sect,dlo_sect
    2684. 

  2684.              num_at_dhi = vaerosolcw*npv_at_dhi
    2685. 

  2685.              num_at_dlo = vaerosolcw*npv_at_dlo
    2686. 

  2686.              naerosolcw = max( num_at_dhi, min( num_at_dlo, naerosolcw ) )
    2687. 

  2687.           else
    2688. 

  2688. !            aerosol number diagnosed from mass and prescribed size
    2689. 

  2689.              naerosol=vaerosol*npv
    2690. 

  2690.              naerosol=max(naerosol,0.)
    2691. 

  2691.              naerosolcw=vaerosolcw*npv
    2692. 

  2692.              naerosolcw=max(naerosolcw,0.)
    2693. 

  2693.           endif
    2694. 

  2694. 

  2695. 

  2696.        return
    2697. 

  2697.        end subroutine loadaer
    2698. 

  2698. 

  2699. 

  2700. 

  2701. !-----------------------------------------------------------------------
    2702. 

  2702.         real function erfc_num_recipes( x )
    2703. 

  2703. !
    2704. 

  2704. !   from press et al, numerical recipes, 1990, page 164
    2705. 

  2705. !
    2706. 

  2706.         implicit none
    2707. 

  2707.         real x
    2708. 

  2708.         double precision erfc_dbl, dum, t, zz
    2709. 

  2709. 

  2710.         zz = abs(x)
    2711. 

  2711.         t = 1.0/(1.0 + 0.5*zz)
    2712. 

  2712. 

  2713. !       erfc_num_recipes =
    2714. 

  2714. !     &   t*exp( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
    2715. 

  2715. !     &   t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
    2716. 

  2716. !     &                                    t*(-1.13520398 +
    2717. 

  2717. !     &   t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
    2718. 

  2718. 

  2719.         dum =  ( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +   &
    2720. 

  2720.           t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +   &
    2721. 

  2721.                                            t*(-1.13520398 +   &
    2722. 

  2722.           t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
    2723. 

  2723. 

  2724.         erfc_dbl = t * exp(dum)
    2725. 

  2725.         if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl
    2726. 

  2726. 

  2727.         erfc_num_recipes = erfc_dbl
    2728. 

  2728. 

  2729.         return
    2730. 

  2730.         end function erfc_num_recipes     
    2731. 

  2731. 

  2732. !-----------------------------------------------------------------------
    2733. 

  2733.     real function erf_alt( x )
    2734. 

  2734. 

  2735.     implicit none
    2736. 

  2736. 

  2737.     real,intent(in) :: x
    2738. 

  2738. 

  2739.     erf_alt = 1. - erfc_num_recipes(x)
    2740. 

  2740. 

  2741.     end function erf_alt
    2742. 

  2742. 

  2743. END MODULE module_mixactivate
    2744. 

  2744.