Line data Source code
1 : module nucleate_ice_cam
2 :
3 : !---------------------------------------------------------------------------------
4 : !
5 : ! CAM Interfaces for nucleate_ice module.
6 : !
7 : ! B. Eaton - Sept 2014
8 : !---------------------------------------------------------------------------------
9 :
10 : use shr_kind_mod, only: r8=>shr_kind_r8
11 : use spmd_utils, only: masterproc
12 : use ppgrid, only: pcols, pver
13 : use physconst, only: pi, rair, tmelt
14 : use constituents, only: pcnst, cnst_get_ind
15 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
16 : use physics_buffer, only: physics_buffer_desc
17 : use phys_control, only: use_hetfrz_classnuc
18 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props
19 :
20 : use physics_buffer, only: pbuf_add_field, dtype_r8, pbuf_old_tim_idx, &
21 : pbuf_get_index, pbuf_get_field, &
22 : pbuf_set_field
23 : use cam_history, only: addfld, add_default, outfld
24 :
25 : use ref_pres, only: top_lev => trop_cloud_top_lev
26 : use wv_saturation, only: qsat_water
27 :
28 : use cam_logfile, only: iulog
29 : use cam_abortutils, only: endrun
30 :
31 : use nucleate_ice, only: nucleati_init, nucleati
32 :
33 : use aerosol_properties_mod, only: aerosol_properties
34 : use aerosol_state_mod, only: aerosol_state
35 :
36 : use phys_control, only: cam_physpkg_is
37 :
38 : implicit none
39 : private
40 : save
41 :
42 : public :: &
43 : nucleate_ice_cam_readnl, &
44 : nucleate_ice_cam_register, &
45 : nucleate_ice_cam_init, &
46 : nucleate_ice_cam_calc
47 :
48 : ! Namelist variables
49 : logical, public, protected :: use_preexisting_ice = .false.
50 : logical :: hist_preexisting_ice = .false.
51 : logical :: nucleate_ice_incloud = .false.
52 : logical :: nucleate_ice_use_troplev = .false.
53 : real(r8) :: nucleate_ice_subgrid = -1._r8
54 : real(r8) :: nucleate_ice_subgrid_strat = -1._r8
55 : real(r8) :: nucleate_ice_strat = 0.0_r8
56 :
57 : ! Vars set via init method.
58 : real(r8) :: mincld ! minimum allowed cloud fraction
59 : real(r8) :: bulk_scale ! prescribed aerosol bulk sulfur scale factor
60 :
61 : ! constituent indices
62 : integer :: &
63 : cldliq_idx = -1, &
64 : cldice_idx = -1, &
65 : numice_idx = -1
66 :
67 : integer :: &
68 : naai_idx = -1, &
69 : naai_hom_idx = -1
70 :
71 : integer :: &
72 : ast_idx = -1
73 :
74 : integer :: &
75 : qsatfac_idx = -1
76 :
77 : ! Bulk aerosols
78 : character(len=20), allocatable :: aername(:)
79 : real(r8), allocatable :: num_to_mass_aer(:)
80 :
81 : integer :: naer_all = -1 ! number of aerosols affecting climate
82 : integer :: idxsul = -1 ! index in aerosol list for sulfate
83 : integer :: idxdst1 = -1 ! index in aerosol list for dust1
84 : integer :: idxdst2 = -1 ! index in aerosol list for dust2
85 : integer :: idxdst3 = -1 ! index in aerosol list for dust3
86 : integer :: idxdst4 = -1 ! index in aerosol list for dust4
87 : integer :: idxbcphi = -1 ! index in aerosol list for Soot (BCPHIL)
88 :
89 : ! modal aerosols
90 : logical :: clim_modal_aero = .false.
91 : logical :: prog_modal_aero = .false.
92 :
93 : logical :: lq(pcnst) = .false. ! set flags true for constituents with non-zero tendencies
94 :
95 : integer, allocatable :: aer_cnst_idx(:,:)
96 :
97 : !===============================================================================
98 : contains
99 : !===============================================================================
100 :
101 178008 : subroutine nucleate_ice_cam_readnl(nlfile)
102 :
103 : use namelist_utils, only: find_group_name
104 : use units, only: getunit, freeunit
105 : use spmd_utils, only: mpicom, masterprocid, mpi_logical, mpi_real8
106 :
107 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
108 :
109 : ! Local variables
110 : integer :: unitn, ierr
111 : character(len=*), parameter :: subname = 'nucleate_ice_cam_readnl'
112 :
113 : namelist /nucleate_ice_nl/ use_preexisting_ice, hist_preexisting_ice, &
114 : nucleate_ice_subgrid, nucleate_ice_subgrid_strat, nucleate_ice_strat, &
115 : nucleate_ice_incloud, nucleate_ice_use_troplev
116 :
117 : !-----------------------------------------------------------------------------
118 :
119 1536 : if (masterproc) then
120 2 : unitn = getunit()
121 2 : open( unitn, file=trim(nlfile), status='old' )
122 2 : call find_group_name(unitn, 'nucleate_ice_nl', status=ierr)
123 2 : if (ierr == 0) then
124 2 : read(unitn, nucleate_ice_nl, iostat=ierr)
125 2 : if (ierr /= 0) then
126 0 : call endrun(subname // ':: ERROR reading namelist')
127 : end if
128 : end if
129 2 : close(unitn)
130 2 : call freeunit(unitn)
131 :
132 : end if
133 :
134 : ! Broadcast namelist variables
135 1536 : call mpi_bcast(use_preexisting_ice, 1, mpi_logical,masterprocid, mpicom, ierr)
136 1536 : call mpi_bcast(hist_preexisting_ice, 1, mpi_logical,masterprocid, mpicom, ierr)
137 1536 : call mpi_bcast(nucleate_ice_subgrid, 1, mpi_real8, masterprocid, mpicom, ierr)
138 1536 : call mpi_bcast(nucleate_ice_subgrid_strat, 1, mpi_real8, masterprocid, mpicom, ierr)
139 1536 : call mpi_bcast(nucleate_ice_strat, 1, mpi_real8, masterprocid, mpicom, ierr)
140 1536 : call mpi_bcast(nucleate_ice_incloud, 1, mpi_logical,masterprocid, mpicom, ierr)
141 1536 : call mpi_bcast(nucleate_ice_use_troplev, 1, mpi_logical,masterprocid, mpicom, ierr)
142 :
143 1536 : end subroutine nucleate_ice_cam_readnl
144 :
145 : !================================================================================================
146 :
147 1536 : subroutine nucleate_ice_cam_register()
148 :
149 : ! global scope for NAAI needed when clubb_do_icesuper=.true.
150 1536 : call pbuf_add_field('NAAI', 'global', dtype_r8, (/pcols,pver/), naai_idx)
151 1536 : call pbuf_add_field('NAAI_HOM', 'physpkg', dtype_r8, (/pcols,pver/), naai_hom_idx)
152 :
153 1536 : end subroutine nucleate_ice_cam_register
154 :
155 : !================================================================================================
156 :
157 1536 : subroutine nucleate_ice_cam_init(mincld_in, bulk_scale_in, pbuf2d, aero_props)
158 : use phys_control, only: phys_getopts
159 : use time_manager, only: is_first_step
160 :
161 : real(r8), intent(in) :: mincld_in
162 : real(r8), intent(in) :: bulk_scale_in
163 : class(aerosol_properties), optional, intent(in) :: aero_props
164 :
165 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
166 :
167 : ! local variables
168 : integer :: iaer
169 : integer :: ierr
170 : integer :: ispc, ibin
171 : integer :: idxtmp
172 : integer :: nmodes
173 :
174 : character(len=*), parameter :: routine = 'nucleate_ice_cam_init'
175 : logical :: history_cesm_forcing
176 :
177 : character(len=32) :: tmpname
178 :
179 : !--------------------------------------------------------------------------------------------
180 1536 : call phys_getopts(prog_modal_aero_out = prog_modal_aero, history_cesm_forcing_out = history_cesm_forcing)
181 :
182 1536 : mincld = mincld_in
183 1536 : bulk_scale = bulk_scale_in
184 :
185 1536 : lq(:) = .false.
186 :
187 1536 : if (prog_modal_aero.and.use_preexisting_ice) then
188 :
189 1536 : if (.not. present(aero_props)) then
190 0 : call endrun(routine//' : aero_props must be present')
191 : end if
192 :
193 : ! constituent tendencies are calculated only if use_preexisting_ice is TRUE
194 : ! set lq for constituent tendencies --
195 :
196 13824 : allocate(aer_cnst_idx(aero_props%nbins(),0:maxval(aero_props%nspecies())), stat=ierr)
197 1536 : if( ierr /= 0 ) then
198 0 : call endrun(routine//': aer_cnst_idx allocation failed')
199 : end if
200 55296 : aer_cnst_idx = -1
201 :
202 9216 : do ibin = 1, aero_props%nbins()
203 7680 : if (aero_props%icenuc_updates_num(ibin)) then
204 :
205 : ! constituents of this bin will need to be updated
206 :
207 1536 : if (aero_props%icenuc_updates_mmr(ibin,0)) then ! species 0 indicates bin MMR
208 0 : call aero_props%amb_mmr_name( ibin, 0, tmpname)
209 : else
210 1536 : call aero_props%amb_num_name( ibin, tmpname)
211 : end if
212 :
213 1536 : call cnst_get_ind(tmpname, idxtmp, abort=.false.)
214 1536 : aer_cnst_idx(ibin,0) = idxtmp
215 1536 : if (idxtmp>0) then
216 1536 : lq(idxtmp) = .true.
217 : end if
218 :
219 : ! iterate over the species within the bin
220 6144 : do ispc = 1, aero_props%nspecies(ibin)
221 6144 : if (aero_props%icenuc_updates_mmr(ibin,ispc)) then
222 : ! this aerosol constituent will be updated
223 1536 : call aero_props%amb_mmr_name( ibin, ispc, tmpname)
224 1536 : call cnst_get_ind(tmpname, idxtmp, abort=.false.)
225 1536 : aer_cnst_idx(ibin,ispc) = idxtmp
226 1536 : if (idxtmp>0) then
227 1536 : lq(idxtmp) = .true.
228 : end if
229 : end if
230 : end do
231 :
232 : end if
233 : end do
234 :
235 : end if
236 :
237 : ! Initialize naai.
238 1536 : if (is_first_step()) then
239 768 : call pbuf_set_field(pbuf2d, naai_idx, 0.0_r8)
240 : end if
241 :
242 1536 : if( masterproc ) then
243 2 : write(iulog,*) 'nucleate_ice parameters:'
244 2 : write(iulog,*) ' mincld = ', mincld_in
245 2 : write(iulog,*) ' bulk_scale = ', bulk_scale_in
246 2 : write(iulog,*) ' use_preexisiting_ice = ', use_preexisting_ice
247 2 : write(iulog,*) ' hist_preexisiting_ice = ', hist_preexisting_ice
248 2 : write(iulog,*) ' nucleate_ice_subgrid = ', nucleate_ice_subgrid
249 2 : write(iulog,*) ' nucleate_ice_subgrid_strat = ', nucleate_ice_subgrid_strat
250 2 : write(iulog,*) ' nucleate_ice_strat = ', nucleate_ice_strat
251 2 : write(iulog,*) ' nucleate_ice_incloud = ', nucleate_ice_incloud
252 2 : write(iulog,*) ' nucleate_ice_use_troplev = ', nucleate_ice_use_troplev
253 : end if
254 :
255 1536 : call cnst_get_ind('CLDLIQ', cldliq_idx)
256 1536 : call cnst_get_ind('CLDICE', cldice_idx)
257 1536 : call cnst_get_ind('NUMICE', numice_idx)
258 1536 : qsatfac_idx = pbuf_get_index('QSATFAC', ierr)
259 :
260 1536 : if (((nucleate_ice_subgrid .eq. -1._r8) .or. (nucleate_ice_subgrid_strat .eq. -1._r8)) .and. (qsatfac_idx .eq. -1)) then
261 0 : call endrun(routine//': ERROR qsatfac is required when subgrid = -1 or subgrid_strat = -1')
262 : end if
263 :
264 1536 : if (cam_physpkg_is("cam7")) then
265 : ! Updates for PUMAS v1.21+
266 3072 : call addfld('NIHFTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to homogenous freezing')
267 3072 : call addfld('NIDEPTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to deposition nucleation')
268 3072 : call addfld('NIIMMTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to immersion freezing')
269 3072 : call addfld('NIMEYTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to meyers deposition')
270 : else
271 0 : call addfld('NIHF', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to homogenous freezing')
272 0 : call addfld('NIDEP', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to deposition nucleation')
273 0 : call addfld('NIIMM', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to immersion freezing')
274 0 : call addfld('NIMEY', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to meyers deposition')
275 : endif
276 :
277 3072 : call addfld('NIREGM',(/ 'lev' /), 'A', 'C', 'Ice Nucleation Temperature Threshold for Regime')
278 3072 : call addfld('NISUBGRID',(/ 'lev' /), 'A', '', 'Ice Nucleation subgrid saturation factor')
279 3072 : call addfld('NITROP_PD',(/ 'lev' /), 'A', '', 'Chemical Tropopause probability')
280 1536 : if ( history_cesm_forcing ) then
281 0 : call add_default('NITROP_PD',8,' ')
282 : endif
283 :
284 1536 : if (use_preexisting_ice) then
285 3072 : call addfld('fhom', (/ 'lev' /), 'A','fraction', 'Fraction of cirrus where homogeneous freezing occur' )
286 3072 : call addfld ('WICE', (/ 'lev' /), 'A','m/s','Vertical velocity Reduction caused by preexisting ice' )
287 3072 : call addfld ('WEFF', (/ 'lev' /), 'A','m/s','Effective Vertical velocity for ice nucleation' )
288 :
289 1536 : if (cam_physpkg_is("cam7")) then
290 : ! Updates for PUMAS v1.21+
291 3072 : call addfld ('INnso4TEN', (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency so4 (in) to ice_nucleation')
292 3072 : call addfld ('INnbcTEN', (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency bc (in) to ice_nucleation')
293 3072 : call addfld ('INndustTEN', (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency dust (in) ice_nucleation')
294 3072 : call addfld ('INondustTEN', (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency dust (out) from ice_nucleation')
295 : call addfld ('INhetTEN', (/ 'lev' /), 'A','1/m3/s', &
296 3072 : 'Tendency for contribution for in-cloud ice number density increase by het nucleation in ice cloud')
297 : call addfld ('INhomTEN', (/ 'lev' /), 'A','1/m3/s', &
298 3072 : 'Tendency for contribution for in-cloud ice number density increase by hom nucleation in ice cloud')
299 : else
300 0 : call addfld ('INnso4', (/ 'lev' /), 'A','1/m3','Number Concentration so4 (in) to ice_nucleation')
301 0 : call addfld ('INnbc', (/ 'lev' /), 'A','1/m3','Number Concentration bc (in) to ice_nucleation')
302 0 : call addfld ('INndust', (/ 'lev' /), 'A','1/m3','Number Concentration dust (in) ice_nucleation')
303 0 : call addfld ('INondust', (/ 'lev' /), 'A','1/m3','Number Concentration dust (out) from ice_nucleation')
304 : call addfld ('INhet', (/ 'lev' /), 'A','1/m3', &
305 0 : 'contribution for in-cloud ice number density increase by het nucleation in ice cloud')
306 : call addfld ('INhom', (/ 'lev' /), 'A','1/m3', &
307 0 : 'contribution for in-cloud ice number density increase by hom nucleation in ice cloud')
308 : endif
309 :
310 3072 : call addfld ('INFrehom', (/ 'lev' /), 'A','frequency','hom IN frequency ice cloud')
311 3072 : call addfld ('INFreIN', (/ 'lev' /), 'A','frequency','frequency of ice nucleation occur')
312 :
313 1536 : if (hist_preexisting_ice) then
314 0 : call add_default ('WSUBI ', 1, ' ') ! addfld/outfld calls are in microp_aero
315 :
316 0 : call add_default ('fhom ', 1, ' ')
317 0 : call add_default ('WICE ', 1, ' ')
318 0 : call add_default ('WEFF ', 1, ' ')
319 0 : call add_default ('INnso4 ', 1, ' ')
320 0 : call add_default ('INnbc ', 1, ' ')
321 0 : call add_default ('INndust ', 1, ' ')
322 0 : call add_default ('INhet ', 1, ' ')
323 0 : call add_default ('INhom ', 1, ' ')
324 0 : call add_default ('INFrehom', 1, ' ')
325 0 : call add_default ('INFreIN ', 1, ' ')
326 : end if
327 : end if
328 :
329 : ! clim_modal_aero determines whether modal aerosols are used in the climate calculation.
330 : ! The modal aerosols can be either prognostic or prescribed.
331 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
332 :
333 1536 : clim_modal_aero = (nmodes > 0)
334 :
335 1536 : if (.not. clim_modal_aero) then
336 :
337 : ! Props needed for BAM number concentration calcs.
338 :
339 0 : call rad_cnst_get_info(0, naero=naer_all)
340 : allocate( &
341 0 : aername(naer_all), &
342 0 : num_to_mass_aer(naer_all) )
343 :
344 0 : do iaer = 1, naer_all
345 : call rad_cnst_get_aer_props(0, iaer, &
346 0 : aername = aername(iaer), &
347 0 : num_to_mass_aer = num_to_mass_aer(iaer))
348 : ! Look for sulfate, dust, and soot in this list (Bulk aerosol only)
349 0 : if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
350 0 : if (trim(aername(iaer)) == 'DUST1') idxdst1 = iaer
351 0 : if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer
352 0 : if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer
353 0 : if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer
354 0 : if (trim(aername(iaer)) == 'BCPHI') idxbcphi = iaer
355 : end do
356 : end if
357 :
358 :
359 : call nucleati_init(use_preexisting_ice, use_hetfrz_classnuc, nucleate_ice_incloud, iulog, pi, &
360 1536 : mincld)
361 :
362 : ! get indices for fields in the physics buffer
363 1536 : ast_idx = pbuf_get_index('AST')
364 :
365 1536 : end subroutine nucleate_ice_cam_init
366 :
367 : !================================================================================================
368 :
369 7588296 : subroutine nucleate_ice_cam_calc( &
370 176472 : state, wsubi, pbuf, dtime, ptend, aero_props, aero_state )
371 :
372 1536 : use tropopause, only: tropopause_findChemTrop
373 :
374 : ! arguments
375 : type(physics_state), target, intent(in) :: state
376 : real(r8), intent(in) :: wsubi(:,:)
377 : type(physics_buffer_desc), pointer :: pbuf(:)
378 : real(r8), intent(in) :: dtime
379 : type(physics_ptend), intent(out) :: ptend
380 : class(aerosol_properties),optional, intent(in) :: aero_props
381 : class(aerosol_state),optional, intent(in) :: aero_state
382 :
383 : ! local workspace
384 :
385 : ! naai and naai_hom are the outputs shared with the microphysics
386 176472 : real(r8), pointer :: naai(:,:) ! number of activated aerosol for ice nucleation
387 176472 : real(r8), pointer :: naai_hom(:,:) ! number of activated aerosol for ice nucleation (homogeneous freezing only)
388 :
389 : integer :: lchnk, ncol
390 : integer :: itim_old
391 : integer :: i, k, l, m
392 :
393 : character(len=32) :: spectype
394 :
395 176472 : real(r8), pointer :: t(:,:) ! input temperature (K)
396 176472 : real(r8), pointer :: qn(:,:) ! input water vapor mixing ratio (kg/kg)
397 176472 : real(r8), pointer :: qc(:,:) ! cloud water mixing ratio (kg/kg)
398 176472 : real(r8), pointer :: qi(:,:) ! cloud ice mixing ratio (kg/kg)
399 176472 : real(r8), pointer :: ni(:,:) ! cloud ice number conc (1/kg)
400 176472 : real(r8), pointer :: pmid(:,:) ! pressure at layer midpoints (pa)
401 :
402 176472 : real(r8), pointer :: aer_mmr(:,:) ! aerosol mass mixing ratio
403 :
404 176472 : real(r8), pointer :: ast(:,:)
405 : real(r8) :: icecldf(pcols,pver) ! ice cloud fraction
406 176472 : real(r8), pointer :: qsatfac(:,:) ! Subgrid cloud water saturation scaling factor.
407 :
408 : real(r8) :: rho(pcols,pver) ! air density (kg m-3)
409 :
410 176472 : real(r8), allocatable :: naer2(:,:,:) ! bulk aerosol number concentration (1/m3)
411 176472 : real(r8), allocatable :: maerosol(:,:,:) ! bulk aerosol mass conc (kg/m3)
412 :
413 : real(r8) :: qs(pcols) ! liquid-ice weighted sat mixing rat (kg/kg)
414 : real(r8) :: es(pcols) ! liquid-ice weighted sat vapor press (pa)
415 : real(r8) :: gammas(pcols) ! parameter for cond/evap of cloud water
416 : integer :: troplev(pcols) ! tropopause level
417 :
418 : real(r8) :: relhum(pcols,pver) ! relative humidity
419 : real(r8) :: icldm(pcols,pver) ! ice cloud fraction
420 :
421 : real(r8) :: dst_num ! total dust aerosol number (#/cm^3)
422 : real(r8) :: dso4_num ! so4 aerosol number (#/cm^3)
423 : real(r8) :: so4_num ! so4 aerosol number (#/cm^3)
424 : real(r8) :: soot_num ! soot (hydrophilic) aerosol number (#/cm^3)
425 : real(r8) :: wght
426 : real(r8) :: oso4_num
427 : real(r8) :: odst_num
428 : real(r8) :: osoot_num
429 : real(r8) :: so4_num_st_cr_tot
430 : real(r8) :: ramp
431 :
432 : real(r8) :: subgrid(pcols,pver)
433 : real(r8) :: trop_pd(pcols,pver)
434 :
435 : ! For pre-existing ice
436 : real(r8) :: fhom(pcols,pver) ! how much fraction of cloud can reach Shom
437 : real(r8) :: wice(pcols,pver) ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at Shom
438 : real(r8) :: weff(pcols,pver) ! effective Vertical velocity for ice nucleation (m/s); weff=wsubi-wice
439 : real(r8) :: INnso4(pcols,pver) ! #/m3, so4 aerosol number used for ice nucleation
440 : real(r8) :: INnbc(pcols,pver) ! #/m3, bc aerosol number used for ice nucleation
441 : real(r8) :: INndust(pcols,pver) ! #/m3, dust aerosol number used for ice nucleation
442 : real(r8) :: INondust(pcols,pver) ! #/m3, dust aerosol number used for ice nucleation
443 : real(r8) :: INhet(pcols,pver) ! #/m3, ice number from het freezing
444 : real(r8) :: INhom(pcols,pver) ! #/m3, ice number from hom freezing
445 : real(r8) :: INFrehom(pcols,pver) ! hom freezing occurence frequency. 1 occur, 0 not occur.
446 : real(r8) :: INFreIN(pcols,pver) ! ice nucleation occerence frequency. 1 occur, 0 not occur.
447 :
448 : ! history output for ice nucleation
449 : real(r8) :: nihf(pcols,pver) !output number conc of ice nuclei due to heterogenous freezing (1/m3)
450 : real(r8) :: niimm(pcols,pver) !output number conc of ice nuclei due to immersion freezing (hetero nuc) (1/m3)
451 : real(r8) :: nidep(pcols,pver) !output number conc of ice nuclei due to deoposion nucleation (hetero nuc) (1/m3)
452 : real(r8) :: nimey(pcols,pver) !output number conc of ice nuclei due to meyers deposition (1/m3)
453 : real(r8) :: regm(pcols,pver) !output temperature thershold for nucleation regime
454 :
455 : real(r8) :: size_wghts(pcols,pver)
456 : real(r8) :: type_wghts(pcols,pver)
457 : real(r8), pointer :: num_col(:,:)
458 : real(r8) :: dust_num_col(pcols,pver)
459 : real(r8) :: sulf_num_col(pcols,pver)
460 : real(r8) :: soot_num_col(pcols,pver)
461 : real(r8) :: sulf_num_tot_col(pcols,pver)
462 :
463 : integer :: idxtmp
464 176472 : real(r8), pointer :: amb_num(:,:)
465 176472 : real(r8), pointer :: amb_mmr(:,:)
466 176472 : real(r8), pointer :: cld_num(:,:)
467 176472 : real(r8), pointer :: cld_mmr(:,:)
468 :
469 : real(r8) :: delmmr, delmmr_sum
470 : real(r8) :: delnum, delnum_sum
471 :
472 : real(r8), parameter :: per_cm3 = 1.e-6_r8 ! factor for m-3 to cm-3 conversions
473 :
474 : !-------------------------------------------------------------------------------
475 :
476 176472 : lchnk = state%lchnk
477 176472 : ncol = state%ncol
478 176472 : t => state%t
479 176472 : qn => state%q(:,:,1)
480 176472 : qc => state%q(:,:,cldliq_idx)
481 176472 : qi => state%q(:,:,cldice_idx)
482 176472 : ni => state%q(:,:,numice_idx)
483 176472 : pmid => state%pmid
484 :
485 274216968 : rho(:ncol,:) = pmid(:ncol,:)/(rair*t(:ncol,:))
486 :
487 176472 : if (clim_modal_aero) then
488 :
489 176472 : call physics_ptend_init(ptend, state%psetcols, 'nucleatei', lq=lq)
490 :
491 : else
492 : ! init number/mass arrays for bulk aerosols
493 : allocate( &
494 : naer2(pcols,pver,naer_all), &
495 0 : maerosol(pcols,pver,naer_all))
496 :
497 0 : do m = 1, naer_all
498 0 : call rad_cnst_get_aer_mmr(0, m, state, pbuf, aer_mmr)
499 0 : maerosol(:ncol,:,m) = aer_mmr(:ncol,:)*rho(:ncol,:)
500 :
501 0 : if (m .eq. idxsul) then
502 0 : naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)*bulk_scale
503 : else
504 0 : naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)
505 : end if
506 : end do
507 :
508 0 : call physics_ptend_init(ptend, state%psetcols, 'nucleatei')
509 : end if
510 :
511 176472 : itim_old = pbuf_old_tim_idx()
512 705888 : call pbuf_get_field(pbuf, ast_idx, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
513 :
514 274216968 : icecldf(:ncol,:pver) = ast(:ncol,:pver)
515 :
516 : ! naai and naai_hom are the outputs from this parameterization
517 176472 : call pbuf_get_field(pbuf, naai_idx, naai)
518 176472 : call pbuf_get_field(pbuf, naai_hom_idx, naai_hom)
519 274216968 : naai(1:ncol,1:pver) = 0._r8
520 274216968 : naai_hom(1:ncol,1:pver) = 0._r8
521 :
522 : ! Use the same criteria that is used in chemistry and in CLUBB (for cloud fraction)
523 : ! to determine whether to use tropospheric or stratospheric settings. Include the
524 : ! tropopause level so that the cold point tropopause will use the stratospheric values.
525 176472 : call tropopause_findChemTrop(state, troplev)
526 :
527 176472 : if ((nucleate_ice_subgrid .eq. -1._r8) .or. (nucleate_ice_subgrid_strat .eq. -1._r8)) then
528 0 : call pbuf_get_field(pbuf, qsatfac_idx, qsatfac)
529 : end if
530 :
531 176472 : trop_pd(:,:) = 0._r8
532 :
533 15000120 : do k = top_lev, pver
534 247696920 : do i = 1, ncol
535 232696800 : trop_pd(i, troplev(i)) = 1._r8
536 :
537 247520448 : if (k <= troplev(i)) then
538 119318698 : if (nucleate_ice_subgrid_strat .eq. -1._r8) then
539 0 : subgrid(i, k) = 1._r8 / qsatfac(i, k)
540 : else
541 119318698 : subgrid(i, k) = nucleate_ice_subgrid_strat
542 : end if
543 : else
544 113378102 : if (nucleate_ice_subgrid .eq. -1._r8) then
545 0 : subgrid(i, k) = 1._r8 / qsatfac(i, k)
546 : else
547 113378102 : subgrid(i, k) = nucleate_ice_subgrid
548 : end if
549 : end if
550 : end do
551 : end do
552 :
553 :
554 : ! initialize history output fields for ice nucleation
555 274216968 : nihf(1:ncol,1:pver) = 0._r8
556 274216968 : niimm(1:ncol,1:pver) = 0._r8
557 274216968 : nidep(1:ncol,1:pver) = 0._r8
558 274216968 : nimey(1:ncol,1:pver) = 0._r8
559 :
560 274216968 : regm(1:ncol,1:pver) = 0._r8
561 :
562 176472 : if (use_preexisting_ice) then
563 176472 : fhom(:,:) = 0.0_r8
564 176472 : wice(:,:) = 0.0_r8
565 176472 : weff(:,:) = 0.0_r8
566 176472 : INnso4(:,:) = 0.0_r8
567 176472 : INnbc(:,:) = 0.0_r8
568 176472 : INndust(:,:) = 0.0_r8
569 176472 : INondust(:,:) = 0.0_r8
570 176472 : INhet(:,:) = 0.0_r8
571 176472 : INhom(:,:) = 0.0_r8
572 176472 : INFrehom(:,:) = 0.0_r8
573 176472 : INFreIN(:,:) = 0.0_r8
574 : endif
575 :
576 15000120 : do k = top_lev, pver
577 : ! Get humidity and saturation vapor pressures
578 14823648 : call qsat_water(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol, gam=gammas(1:ncol))
579 :
580 247696920 : do i = 1, ncol
581 :
582 232696800 : relhum(i,k) = qn(i,k)/qs(i)
583 :
584 : ! get cloud fraction, check for minimum
585 247520448 : icldm(i,k) = max(icecldf(i,k), mincld)
586 :
587 : end do
588 : end do
589 :
590 176472 : dust_num_col = 0._r8
591 176472 : sulf_num_col = 0._r8
592 176472 : sulf_num_tot_col = 0._r8
593 176472 : soot_num_col = 0._r8
594 :
595 176472 : if (clim_modal_aero) then
596 :
597 529416 : if (.not.(present(aero_props).and.present(aero_state))) then
598 0 : call endrun('nucleate_ice_cam_calc: aero_props and aero_state must be present')
599 : end if
600 :
601 : ! collect number densities (#/cm^3) for dust, sulfate, and soot
602 : call aero_state%nuclice_get_numdens( aero_props, use_preexisting_ice, ncol, pver, rho, &
603 176472 : dust_num_col, sulf_num_col, soot_num_col, sulf_num_tot_col )
604 :
605 : else
606 : ! for bulk model
607 0 : dust_num_col(:ncol,:) = naer2(:ncol,:,idxdst1)/25._r8 * per_cm3 & ! #/cm3
608 0 : + naer2(:ncol,:,idxdst2)/25._r8 * per_cm3 &
609 0 : + naer2(:ncol,:,idxdst3)/25._r8 * per_cm3 &
610 0 : + naer2(:ncol,:,idxdst4)/25._r8 * per_cm3
611 0 : sulf_num_col(:ncol,:) = naer2(:ncol,:,idxsul)/25._r8 * per_cm3
612 0 : soot_num_col(:ncol,:) = naer2(:ncol,:,idxbcphi)/25._r8 * per_cm3
613 : endif
614 :
615 15000120 : kloop: do k = top_lev, pver
616 247696920 : iloop: do i = 1, ncol
617 :
618 232696800 : so4_num_st_cr_tot = 0._r8
619 :
620 247520448 : freezing: if (t(i,k) < tmelt - 5._r8) then
621 :
622 : ! set aerosol number for so4, soot, and dust with units #/cm^3
623 180451908 : so4_num = sulf_num_col(i,k)
624 180451908 : dst_num = dust_num_col(i,k)
625 180451908 : so4_num_st_cr_tot=sulf_num_tot_col(i,k)
626 :
627 : ! *** Turn off soot nucleation ***
628 180451908 : soot_num = 0.0_r8
629 :
630 180451908 : if (cam_physpkg_is("cam7")) then
631 :
632 : call nucleati( &
633 360903816 : wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k), &
634 180451908 : qc(i,k), qi(i,k), ni(i,k), rho(i,k), &
635 : so4_num, dst_num, soot_num, subgrid(i,k), &
636 0 : naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), &
637 : wice(i,k), weff(i,k), fhom(i,k), regm(i,k), &
638 : oso4_num, odst_num, osoot_num, &
639 541355724 : call_frm_zm_in = .false., add_preexisting_ice_in = .false.)
640 :
641 : else
642 :
643 : call nucleati( &
644 0 : wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k), &
645 0 : qc(i,k), qi(i,k), ni(i,k), rho(i,k), &
646 : so4_num, dst_num, soot_num, subgrid(i,k), &
647 0 : naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), &
648 : wice(i,k), weff(i,k), fhom(i,k), regm(i,k), &
649 0 : oso4_num, odst_num, osoot_num)
650 :
651 : end if
652 :
653 : ! Move aerosol used for nucleation from interstial to cloudborne,
654 : ! otherwise the same coarse mode aerosols will be available again
655 : ! in the next timestep and will supress homogeneous freezing.
656 :
657 :
658 180451908 : if (prog_modal_aero .and. use_preexisting_ice) then
659 :
660 : ! compute tendencies for transported aerosol constituents
661 : ! and update not-transported constituents
662 :
663 902259540 : do m = 1, aero_props%nbins()
664 :
665 902259540 : if (aero_props%icenuc_updates_num(m)) then
666 :
667 : ! constituents of this bin will need to be updated
668 :
669 180451908 : call aero_state%get_ambient_num(m, amb_num)
670 180451908 : call aero_state%get_cldbrne_num(m, cld_num)
671 :
672 180451908 : if (amb_num(i,k)>0._r8) then
673 180451908 : delmmr_sum = 0._r8
674 180451908 : delnum_sum = 0._r8
675 :
676 : ! iterate over the species within the bin
677 721807632 : do l = 1, aero_props%nspecies(m)
678 721807632 : if (aero_props%icenuc_updates_mmr(m,l)) then
679 :
680 180451908 : call aero_props%species_type(m, l, spectype)
681 180451908 : call aero_state%icenuc_size_wght( m, i,k, spectype, use_preexisting_ice, wght)
682 :
683 180451908 : if (wght>0._r8) then
684 :
685 : ! this aerosol constituent will be updated
686 :
687 180451908 : idxtmp = aer_cnst_idx(m,l)
688 :
689 180451908 : call aero_state%get_ambient_mmr(l,m,amb_mmr)
690 180451908 : call aero_state%get_cldbrne_mmr(l,m,cld_mmr)
691 :
692 : ! determine change in aerosol mass
693 180451908 : delmmr = 0._r8
694 180451908 : delnum = 0._r8
695 180451908 : if (trim(spectype)=='dust') then
696 180451908 : if (dst_num>0._r8) then
697 180451908 : delmmr = (odst_num / dst_num) * icldm(i,k) * amb_mmr(i,k) * wght
698 180451908 : delnum = (odst_num * icldm(i,k)) /rho(i,k)/per_cm3
699 : endif
700 0 : elseif (trim(spectype)=='sulfate') then
701 0 : if (so4_num>0._r8) then
702 0 : delmmr = (oso4_num / so4_num) * icldm(i,k) * amb_mmr(i,k) * wght
703 0 : delnum = (oso4_num * icldm(i,k)) /rho(i,k)/per_cm3
704 : endif
705 : endif
706 :
707 180451908 : if (idxtmp>0) then
708 : ! constituent tendency (for transported species)
709 180451908 : ptend%q(i,k,idxtmp) = -delmmr/dtime
710 : else
711 : ! apply change of mass to not-transported species
712 0 : amb_mmr(i,k) = amb_mmr(i,k) - delmmr
713 : endif
714 180451908 : cld_mmr(i,k) = cld_mmr(i,k) + delmmr
715 :
716 180451908 : delmmr_sum = delmmr_sum + delmmr
717 180451908 : delnum_sum = delnum_sum + delnum
718 : end if
719 : end if
720 : end do
721 :
722 180451908 : idxtmp = aer_cnst_idx(m,0)
723 :
724 : ! update aerosol state bin and tendency for grid box i,k
725 180451908 : call aero_state%update_bin( m,i,k, delmmr_sum, delnum_sum, idxtmp, dtime, ptend%q )
726 :
727 : end if
728 :
729 : end if
730 : end do
731 :
732 : end if
733 :
734 :
735 : ! Liu&Penner does not generate enough nucleation in the polar winter
736 : ! stratosphere, which affects surface area density, dehydration and
737 : ! ozone chemistry. Part of this is that there are a larger number of
738 : ! particles in the accumulation mode than in the Aitken mode. In volcanic
739 : ! periods, the coarse mode may also be important. As a short
740 : ! term work around, include the accumulation and coarse mode particles
741 : ! and assume a larger fraction of the sulfates nucleate in the polar
742 : ! stratosphere.
743 : !
744 : ! Do not include the tropopause level, as stratospheric aerosols
745 : ! only exist above the tropopause level.
746 : !
747 : ! NOTE: This may still not represent the proper particles that
748 : ! participate in nucleation, because it doesn't include STS and NAT
749 : ! particles. It may not represent the proper saturation threshold for
750 : ! nucleation, and wsubi from CLUBB is probably not representative of
751 : ! wave driven varaibility in the polar stratosphere.
752 180451908 : if (nucleate_ice_use_troplev .and. clim_modal_aero) then
753 180451908 : if ((k < troplev(i)) .and. (nucleate_ice_strat > 0._r8) .and. (oso4_num > 0._r8)) then
754 100 : dso4_num = max(0._r8, (nucleate_ice_strat*so4_num_st_cr_tot - oso4_num) * 1e6_r8 / rho(i,k))
755 100 : naai(i,k) = naai(i,k) + dso4_num
756 100 : nihf(i,k) = nihf(i,k) + dso4_num
757 : endif
758 : else
759 : ! This maintains backwards compatibility with the previous version.
760 0 : if (pmid(i,k) <= 12500._r8 .and. pmid(i,k) > 100._r8 .and. abs(state%lat(i)) >= 60._r8 * pi / 180._r8) then
761 0 : ramp = 1._r8 - min(1._r8, max(0._r8, (pmid(i,k) - 10000._r8) / 2500._r8))
762 :
763 0 : if (oso4_num > 0._r8) then
764 0 : dso4_num = (max(oso4_num, ramp * nucleate_ice_strat * so4_num) - oso4_num) * 1e6_r8 / rho(i,k)
765 0 : naai(i,k) = naai(i,k) + dso4_num
766 0 : nihf(i,k) = nihf(i,k) + dso4_num
767 : end if
768 : end if
769 : end if
770 :
771 180451908 : if (cam_physpkg_is("cam7")) then
772 : !Updates for pumas v1.21+
773 :
774 180451908 : naai_hom(i,k) = nihf(i,k)/dtime
775 180451908 : naai(i,k)= naai(i,k)/dtime
776 :
777 : ! output activated ice (convert from #/kg -> #/m3/s)
778 180451908 : nihf(i,k) = nihf(i,k) *rho(i,k)/dtime
779 180451908 : niimm(i,k) = niimm(i,k)*rho(i,k)/dtime
780 180451908 : nidep(i,k) = nidep(i,k)*rho(i,k)/dtime
781 180451908 : nimey(i,k) = nimey(i,k)*rho(i,k)/dtime
782 :
783 180451908 : if (use_preexisting_ice) then
784 180451908 : INnso4(i,k) =so4_num*1e6_r8/dtime ! (convert from #/cm3 -> #/m3/s)
785 180451908 : INnbc(i,k) =soot_num*1e6_r8/dtime
786 180451908 : INndust(i,k)=dst_num*1e6_r8/dtime
787 180451908 : INondust(i,k)=odst_num*1e6_r8/dtime
788 180451908 : INFreIN(i,k)=1.0_r8 ! 1,ice nucleation occur
789 180451908 : INhet(i,k) = (niimm(i,k) + nidep(i,k)) ! #/m3/s, nimey not in cirrus
790 180451908 : INhom(i,k) = nihf(i,k) ! #/m3/s
791 180451908 : if (INhom(i,k).gt.1e3_r8) then ! > 1/L
792 1495 : INFrehom(i,k)=1.0_r8 ! 1, hom freezing occur
793 : endif
794 :
795 : ! exclude no ice nucleaton
796 180451908 : if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8)) then
797 180323503 : INnso4(i,k) =0.0_r8
798 180323503 : INnbc(i,k) =0.0_r8
799 180323503 : INndust(i,k)=0.0_r8
800 180323503 : INondust(i,k)=0.0_r8
801 180323503 : INFreIN(i,k)=0.0_r8
802 180323503 : INhet(i,k) = 0.0_r8
803 180323503 : INhom(i,k) = 0.0_r8
804 180323503 : INFrehom(i,k)=0.0_r8
805 180323503 : wice(i,k) = 0.0_r8
806 180323503 : weff(i,k) = 0.0_r8
807 180323503 : fhom(i,k) = 0.0_r8
808 : endif
809 : endif
810 :
811 : else ! Not cam7
812 :
813 0 : naai_hom(i,k) = nihf(i,k)
814 :
815 : ! output activated ice (convert from #/kg -> #/m3/s)
816 0 : nihf(i,k) = nihf(i,k) *rho(i,k)
817 0 : niimm(i,k) = niimm(i,k)*rho(i,k)
818 0 : nidep(i,k) = nidep(i,k)*rho(i,k)
819 0 : nimey(i,k) = nimey(i,k)*rho(i,k)
820 :
821 0 : if (use_preexisting_ice) then
822 0 : INnso4(i,k) =so4_num*1e6_r8 ! (convert from #/cm3 -> #/m3/s)
823 0 : INnbc(i,k) =soot_num*1e6_r8
824 0 : INndust(i,k)=dst_num*1e6_r8
825 0 : INondust(i,k)=odst_num*1e6_r8
826 0 : INFreIN(i,k)=1.0_r8 ! 1,ice nucleation occur
827 0 : INhet(i,k) = (niimm(i,k) + nidep(i,k)) ! #/m3, nimey not in cirrus
828 0 : INhom(i,k) = nihf(i,k) ! #/m3
829 0 : if (INhom(i,k).gt.1e3_r8) then ! > 1/L
830 0 : INFrehom(i,k)=1.0_r8 ! 1, hom freezing occur
831 : endif
832 :
833 : ! exclude no ice nucleaton
834 0 : if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8)) then
835 0 : INnso4(i,k) =0.0_r8
836 0 : INnbc(i,k) =0.0_r8
837 0 : INndust(i,k)=0.0_r8
838 0 : INondust(i,k)=0.0_r8
839 0 : INFreIN(i,k)=0.0_r8
840 0 : INhet(i,k) = 0.0_r8
841 0 : INhom(i,k) = 0.0_r8
842 0 : INFrehom(i,k)=0.0_r8
843 0 : wice(i,k) = 0.0_r8
844 0 : weff(i,k) = 0.0_r8
845 0 : fhom(i,k) = 0.0_r8
846 : endif
847 : end if
848 :
849 : end if ! cam7
850 : end if freezing
851 : end do iloop
852 : end do kloop
853 :
854 176472 : if (.not. clim_modal_aero) then
855 0 : deallocate( &
856 : naer2, &
857 0 : maerosol)
858 : end if
859 :
860 176472 : if (cam_physpkg_is("cam7")) then
861 : ! Updates for PUMAS v1.21+
862 176472 : call outfld('NIHFTEN', nihf, pcols, lchnk)
863 176472 : call outfld('NIIMMTEN', niimm, pcols, lchnk)
864 176472 : call outfld('NIDEPTEN', nidep, pcols, lchnk)
865 176472 : call outfld('NIMEYTEN', nimey, pcols, lchnk)
866 : else
867 0 : call outfld('NIHF', nihf, pcols, lchnk)
868 0 : call outfld('NIIMM', niimm, pcols, lchnk)
869 0 : call outfld('NIDEP', nidep, pcols, lchnk)
870 0 : call outfld('NIMEY', nimey, pcols, lchnk)
871 : end if
872 176472 : call outfld('NIREGM', regm, pcols, lchnk)
873 176472 : call outfld('NISUBGRID', subgrid, pcols, lchnk)
874 176472 : call outfld('NITROP_PD', trop_pd, pcols, lchnk)
875 :
876 176472 : if (use_preexisting_ice) then
877 176472 : call outfld( 'fhom' , fhom, pcols, lchnk)
878 176472 : call outfld( 'WICE' , wice, pcols, lchnk)
879 176472 : call outfld( 'WEFF' , weff, pcols, lchnk)
880 176472 : if (cam_physpkg_is("cam7")) then
881 : ! Updates for PUMAS v1.21+
882 176472 : call outfld('INnso4TEN',INnso4 , pcols,lchnk)
883 176472 : call outfld('INnbcTEN',INnbc , pcols,lchnk)
884 176472 : call outfld('INndustTEN',INndust, pcols,lchnk)
885 176472 : call outfld('INondustTEN',INondust, pcols,lchnk)
886 176472 : call outfld('INhetTEN',INhet , pcols,lchnk)
887 176472 : call outfld('INhomTEN',INhom , pcols,lchnk)
888 : else
889 0 : call outfld('INnso4 ',INnso4 , pcols,lchnk)
890 0 : call outfld('INnbc ',INnbc , pcols,lchnk)
891 0 : call outfld('INndust ',INndust, pcols,lchnk)
892 0 : call outfld('INondust ',INondust, pcols,lchnk)
893 0 : call outfld('INhet ',INhet , pcols,lchnk)
894 0 : call outfld('INhom ',INhom , pcols,lchnk)
895 : end if
896 176472 : call outfld('INFrehom',INFrehom,pcols,lchnk)
897 176472 : call outfld('INFreIN ',INFreIN, pcols,lchnk)
898 : end if
899 :
900 352944 : end subroutine nucleate_ice_cam_calc
901 :
902 : !================================================================================================
903 :
904 : end module nucleate_ice_cam
|
