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 4469064 : 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 192103704 : subroutine nucleate_ice_cam_calc( &
370 4467528 : 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 4467528 : real(r8), pointer :: naai(:,:) ! number of activated aerosol for ice nucleation
387 4467528 : 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 4467528 : real(r8), pointer :: t(:,:) ! input temperature (K)
396 4467528 : real(r8), pointer :: qn(:,:) ! input water vapor mixing ratio (kg/kg)
397 4467528 : real(r8), pointer :: qc(:,:) ! cloud water mixing ratio (kg/kg)
398 4467528 : real(r8), pointer :: qi(:,:) ! cloud ice mixing ratio (kg/kg)
399 4467528 : real(r8), pointer :: ni(:,:) ! cloud ice number conc (1/kg)
400 4467528 : real(r8), pointer :: pmid(:,:) ! pressure at layer midpoints (pa)
401 :
402 4467528 : real(r8), pointer :: aer_mmr(:,:) ! aerosol mass mixing ratio
403 :
404 4467528 : real(r8), pointer :: ast(:,:)
405 : real(r8) :: icecldf(pcols,pver) ! ice cloud fraction
406 4467528 : real(r8), pointer :: qsatfac(:,:) ! Subgrid cloud water saturation scaling factor.
407 :
408 : real(r8) :: rho(pcols,pver) ! air density (kg m-3)
409 :
410 4467528 : real(r8), allocatable :: naer2(:,:,:) ! bulk aerosol number concentration (1/m3)
411 4467528 : 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 4467528 : real(r8), pointer :: amb_num(:,:)
465 4467528 : real(r8), pointer :: amb_mmr(:,:)
466 4467528 : real(r8), pointer :: cld_num(:,:)
467 4467528 : 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 4467528 : lchnk = state%lchnk
477 4467528 : ncol = state%ncol
478 4467528 : t => state%t
479 4467528 : qn => state%q(:,:,1)
480 4467528 : qc => state%q(:,:,cldliq_idx)
481 4467528 : qi => state%q(:,:,cldice_idx)
482 4467528 : ni => state%q(:,:,numice_idx)
483 4467528 : pmid => state%pmid
484 :
485 6942019032 : rho(:ncol,:) = pmid(:ncol,:)/(rair*t(:ncol,:))
486 :
487 4467528 : if (clim_modal_aero) then
488 :
489 4467528 : 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 4467528 : itim_old = pbuf_old_tim_idx()
512 17870112 : call pbuf_get_field(pbuf, ast_idx, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
513 :
514 6942019032 : icecldf(:ncol,:pver) = ast(:ncol,:pver)
515 :
516 : ! naai and naai_hom are the outputs from this parameterization
517 4467528 : call pbuf_get_field(pbuf, naai_idx, naai)
518 4467528 : call pbuf_get_field(pbuf, naai_hom_idx, naai_hom)
519 6942019032 : naai(1:ncol,1:pver) = 0._r8
520 6942019032 : 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 : !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
526 4467528 : troplev(:) = 0
527 : !REMOVECAM_END
528 4467528 : call tropopause_findChemTrop(state, troplev)
529 :
530 4467528 : if ((nucleate_ice_subgrid .eq. -1._r8) .or. (nucleate_ice_subgrid_strat .eq. -1._r8)) then
531 0 : call pbuf_get_field(pbuf, qsatfac_idx, qsatfac)
532 : end if
533 :
534 4467528 : trop_pd(:,:) = 0._r8
535 :
536 379739880 : do k = top_lev, pver
537 6270643080 : do i = 1, ncol
538 5890903200 : trop_pd(i, troplev(i)) = 1._r8
539 :
540 6266175552 : if (k <= troplev(i)) then
541 3033110807 : if (nucleate_ice_subgrid_strat .eq. -1._r8) then
542 0 : subgrid(i, k) = 1._r8 / qsatfac(i, k)
543 : else
544 3033110807 : subgrid(i, k) = nucleate_ice_subgrid_strat
545 : end if
546 : else
547 2857792393 : if (nucleate_ice_subgrid .eq. -1._r8) then
548 0 : subgrid(i, k) = 1._r8 / qsatfac(i, k)
549 : else
550 2857792393 : subgrid(i, k) = nucleate_ice_subgrid
551 : end if
552 : end if
553 : end do
554 : end do
555 :
556 :
557 : ! initialize history output fields for ice nucleation
558 6942019032 : nihf(1:ncol,1:pver) = 0._r8
559 6942019032 : niimm(1:ncol,1:pver) = 0._r8
560 6942019032 : nidep(1:ncol,1:pver) = 0._r8
561 6942019032 : nimey(1:ncol,1:pver) = 0._r8
562 :
563 6942019032 : regm(1:ncol,1:pver) = 0._r8
564 :
565 4467528 : if (use_preexisting_ice) then
566 4467528 : fhom(:,:) = 0.0_r8
567 4467528 : wice(:,:) = 0.0_r8
568 4467528 : weff(:,:) = 0.0_r8
569 4467528 : INnso4(:,:) = 0.0_r8
570 4467528 : INnbc(:,:) = 0.0_r8
571 4467528 : INndust(:,:) = 0.0_r8
572 4467528 : INondust(:,:) = 0.0_r8
573 4467528 : INhet(:,:) = 0.0_r8
574 4467528 : INhom(:,:) = 0.0_r8
575 4467528 : INFrehom(:,:) = 0.0_r8
576 4467528 : INFreIN(:,:) = 0.0_r8
577 : endif
578 :
579 379739880 : do k = top_lev, pver
580 : ! Get humidity and saturation vapor pressures
581 375272352 : call qsat_water(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol, gam=gammas(1:ncol))
582 :
583 6270643080 : do i = 1, ncol
584 :
585 5890903200 : relhum(i,k) = qn(i,k)/qs(i)
586 :
587 : ! get cloud fraction, check for minimum
588 6266175552 : icldm(i,k) = max(icecldf(i,k), mincld)
589 :
590 : end do
591 : end do
592 :
593 4467528 : dust_num_col = 0._r8
594 4467528 : sulf_num_col = 0._r8
595 4467528 : sulf_num_tot_col = 0._r8
596 4467528 : soot_num_col = 0._r8
597 :
598 4467528 : if (clim_modal_aero) then
599 :
600 13402584 : if (.not.(present(aero_props).and.present(aero_state))) then
601 0 : call endrun('nucleate_ice_cam_calc: aero_props and aero_state must be present')
602 : end if
603 :
604 : ! collect number densities (#/cm^3) for dust, sulfate, and soot
605 : call aero_state%nuclice_get_numdens( aero_props, use_preexisting_ice, ncol, pver, rho, &
606 4467528 : dust_num_col, sulf_num_col, soot_num_col, sulf_num_tot_col )
607 :
608 : else
609 : ! for bulk model
610 0 : dust_num_col(:ncol,:) = naer2(:ncol,:,idxdst1)/25._r8 * per_cm3 & ! #/cm3
611 0 : + naer2(:ncol,:,idxdst2)/25._r8 * per_cm3 &
612 0 : + naer2(:ncol,:,idxdst3)/25._r8 * per_cm3 &
613 0 : + naer2(:ncol,:,idxdst4)/25._r8 * per_cm3
614 0 : sulf_num_col(:ncol,:) = naer2(:ncol,:,idxsul)/25._r8 * per_cm3
615 0 : soot_num_col(:ncol,:) = naer2(:ncol,:,idxbcphi)/25._r8 * per_cm3
616 : endif
617 :
618 379739880 : kloop: do k = top_lev, pver
619 6270643080 : iloop: do i = 1, ncol
620 :
621 5890903200 : so4_num_st_cr_tot = 0._r8
622 :
623 6266175552 : freezing: if (t(i,k) < tmelt - 5._r8) then
624 :
625 : ! set aerosol number for so4, soot, and dust with units #/cm^3
626 4575144061 : so4_num = sulf_num_col(i,k)
627 4575144061 : dst_num = dust_num_col(i,k)
628 4575144061 : so4_num_st_cr_tot=sulf_num_tot_col(i,k)
629 :
630 : ! *** Turn off soot nucleation ***
631 4575144061 : soot_num = 0.0_r8
632 :
633 4575144061 : if (cam_physpkg_is("cam7")) then
634 :
635 : call nucleati( &
636 9150288122 : wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k), &
637 4575144061 : qc(i,k), qi(i,k), ni(i,k), rho(i,k), &
638 : so4_num, dst_num, soot_num, subgrid(i,k), &
639 0 : naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), &
640 : wice(i,k), weff(i,k), fhom(i,k), regm(i,k), &
641 : oso4_num, odst_num, osoot_num, &
642 13725432183 : call_frm_zm_in = .false., add_preexisting_ice_in = .false.)
643 :
644 : else
645 :
646 : call nucleati( &
647 0 : wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k), &
648 0 : qc(i,k), qi(i,k), ni(i,k), rho(i,k), &
649 : so4_num, dst_num, soot_num, subgrid(i,k), &
650 0 : naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), &
651 : wice(i,k), weff(i,k), fhom(i,k), regm(i,k), &
652 0 : oso4_num, odst_num, osoot_num)
653 :
654 : end if
655 :
656 : ! Move aerosol used for nucleation from interstial to cloudborne,
657 : ! otherwise the same coarse mode aerosols will be available again
658 : ! in the next timestep and will supress homogeneous freezing.
659 :
660 :
661 4575144061 : if (prog_modal_aero .and. use_preexisting_ice) then
662 :
663 : ! compute tendencies for transported aerosol constituents
664 : ! and update not-transported constituents
665 :
666 22875720305 : do m = 1, aero_props%nbins()
667 :
668 22875720305 : if (aero_props%icenuc_updates_num(m)) then
669 :
670 : ! constituents of this bin will need to be updated
671 :
672 4575144061 : call aero_state%get_ambient_num(m, amb_num)
673 4575144061 : call aero_state%get_cldbrne_num(m, cld_num)
674 :
675 4575144061 : if (amb_num(i,k)>0._r8) then
676 4575144061 : delmmr_sum = 0._r8
677 4575144061 : delnum_sum = 0._r8
678 :
679 : ! iterate over the species within the bin
680 18300576244 : do l = 1, aero_props%nspecies(m)
681 18300576244 : if (aero_props%icenuc_updates_mmr(m,l)) then
682 :
683 4575144061 : call aero_props%species_type(m, l, spectype)
684 4575144061 : call aero_state%icenuc_size_wght( m, i,k, spectype, use_preexisting_ice, wght)
685 :
686 4575144061 : if (wght>0._r8) then
687 :
688 : ! this aerosol constituent will be updated
689 :
690 4575144061 : idxtmp = aer_cnst_idx(m,l)
691 :
692 4575144061 : call aero_state%get_ambient_mmr(l,m,amb_mmr)
693 4575144061 : call aero_state%get_cldbrne_mmr(l,m,cld_mmr)
694 :
695 : ! determine change in aerosol mass
696 4575144061 : delmmr = 0._r8
697 4575144061 : delnum = 0._r8
698 4575144061 : if (trim(spectype)=='dust') then
699 4575144061 : if (dst_num>0._r8) then
700 4575144061 : delmmr = (odst_num / dst_num) * icldm(i,k) * amb_mmr(i,k) * wght
701 4575144061 : delnum = (odst_num * icldm(i,k)) /rho(i,k)/per_cm3
702 : endif
703 0 : elseif (trim(spectype)=='sulfate') then
704 0 : if (so4_num>0._r8) then
705 0 : delmmr = (oso4_num / so4_num) * icldm(i,k) * amb_mmr(i,k) * wght
706 0 : delnum = (oso4_num * icldm(i,k)) /rho(i,k)/per_cm3
707 : endif
708 : endif
709 :
710 4575144061 : if (idxtmp>0) then
711 : ! constituent tendency (for transported species)
712 4575144061 : ptend%q(i,k,idxtmp) = -delmmr/dtime
713 : else
714 : ! apply change of mass to not-transported species
715 0 : amb_mmr(i,k) = amb_mmr(i,k) - delmmr
716 : endif
717 4575144061 : cld_mmr(i,k) = cld_mmr(i,k) + delmmr
718 :
719 4575144061 : delmmr_sum = delmmr_sum + delmmr
720 4575144061 : delnum_sum = delnum_sum + delnum
721 : end if
722 : end if
723 : end do
724 :
725 4575144061 : idxtmp = aer_cnst_idx(m,0)
726 :
727 : ! update aerosol state bin and tendency for grid box i,k
728 4575144061 : call aero_state%update_bin( m,i,k, delmmr_sum, delnum_sum, idxtmp, dtime, ptend%q )
729 :
730 : end if
731 :
732 : end if
733 : end do
734 :
735 : end if
736 :
737 :
738 : ! Liu&Penner does not generate enough nucleation in the polar winter
739 : ! stratosphere, which affects surface area density, dehydration and
740 : ! ozone chemistry. Part of this is that there are a larger number of
741 : ! particles in the accumulation mode than in the Aitken mode. In volcanic
742 : ! periods, the coarse mode may also be important. As a short
743 : ! term work around, include the accumulation and coarse mode particles
744 : ! and assume a larger fraction of the sulfates nucleate in the polar
745 : ! stratosphere.
746 : !
747 : ! Do not include the tropopause level, as stratospheric aerosols
748 : ! only exist above the tropopause level.
749 : !
750 : ! NOTE: This may still not represent the proper particles that
751 : ! participate in nucleation, because it doesn't include STS and NAT
752 : ! particles. It may not represent the proper saturation threshold for
753 : ! nucleation, and wsubi from CLUBB is probably not representative of
754 : ! wave driven varaibility in the polar stratosphere.
755 4575144061 : if (nucleate_ice_use_troplev .and. clim_modal_aero) then
756 4575144061 : if ((k < troplev(i)) .and. (nucleate_ice_strat > 0._r8) .and. (oso4_num > 0._r8)) then
757 2560 : dso4_num = max(0._r8, (nucleate_ice_strat*so4_num_st_cr_tot - oso4_num) * 1e6_r8 / rho(i,k))
758 2560 : naai(i,k) = naai(i,k) + dso4_num
759 2560 : nihf(i,k) = nihf(i,k) + dso4_num
760 : endif
761 : else
762 : ! This maintains backwards compatibility with the previous version.
763 0 : if (pmid(i,k) <= 12500._r8 .and. pmid(i,k) > 100._r8 .and. abs(state%lat(i)) >= 60._r8 * pi / 180._r8) then
764 0 : ramp = 1._r8 - min(1._r8, max(0._r8, (pmid(i,k) - 10000._r8) / 2500._r8))
765 :
766 0 : if (oso4_num > 0._r8) then
767 0 : dso4_num = (max(oso4_num, ramp * nucleate_ice_strat * so4_num) - oso4_num) * 1e6_r8 / rho(i,k)
768 0 : naai(i,k) = naai(i,k) + dso4_num
769 0 : nihf(i,k) = nihf(i,k) + dso4_num
770 : end if
771 : end if
772 : end if
773 :
774 4575144061 : if (cam_physpkg_is("cam7")) then
775 : !Updates for pumas v1.21+
776 :
777 4575144061 : naai_hom(i,k) = nihf(i,k)/dtime
778 4575144061 : naai(i,k)= naai(i,k)/dtime
779 :
780 : ! output activated ice (convert from #/kg -> #/m3/s)
781 4575144061 : nihf(i,k) = nihf(i,k) *rho(i,k)/dtime
782 4575144061 : niimm(i,k) = niimm(i,k)*rho(i,k)/dtime
783 4575144061 : nidep(i,k) = nidep(i,k)*rho(i,k)/dtime
784 4575144061 : nimey(i,k) = nimey(i,k)*rho(i,k)/dtime
785 :
786 4575144061 : if (use_preexisting_ice) then
787 4575144061 : INnso4(i,k) =so4_num*1e6_r8/dtime ! (convert from #/cm3 -> #/m3/s)
788 4575144061 : INnbc(i,k) =soot_num*1e6_r8/dtime
789 4575144061 : INndust(i,k)=dst_num*1e6_r8/dtime
790 4575144061 : INondust(i,k)=odst_num*1e6_r8/dtime
791 4575144061 : INFreIN(i,k)=1.0_r8 ! 1,ice nucleation occur
792 4575144061 : INhet(i,k) = (niimm(i,k) + nidep(i,k)) ! #/m3/s, nimey not in cirrus
793 4575144061 : INhom(i,k) = nihf(i,k) ! #/m3/s
794 4575144061 : if (INhom(i,k).gt.1e3_r8) then ! > 1/L
795 39701 : INFrehom(i,k)=1.0_r8 ! 1, hom freezing occur
796 : endif
797 :
798 : ! exclude no ice nucleaton
799 4575144061 : if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8)) then
800 4565036680 : INnso4(i,k) =0.0_r8
801 4565036680 : INnbc(i,k) =0.0_r8
802 4565036680 : INndust(i,k)=0.0_r8
803 4565036680 : INondust(i,k)=0.0_r8
804 4565036680 : INFreIN(i,k)=0.0_r8
805 4565036680 : INhet(i,k) = 0.0_r8
806 4565036680 : INhom(i,k) = 0.0_r8
807 4565036680 : INFrehom(i,k)=0.0_r8
808 4565036680 : wice(i,k) = 0.0_r8
809 4565036680 : weff(i,k) = 0.0_r8
810 4565036680 : fhom(i,k) = 0.0_r8
811 : endif
812 : endif
813 :
814 : else ! Not cam7
815 :
816 0 : naai_hom(i,k) = nihf(i,k)
817 :
818 : ! output activated ice (convert from #/kg -> #/m3/s)
819 0 : nihf(i,k) = nihf(i,k) *rho(i,k)
820 0 : niimm(i,k) = niimm(i,k)*rho(i,k)
821 0 : nidep(i,k) = nidep(i,k)*rho(i,k)
822 0 : nimey(i,k) = nimey(i,k)*rho(i,k)
823 :
824 0 : if (use_preexisting_ice) then
825 0 : INnso4(i,k) =so4_num*1e6_r8 ! (convert from #/cm3 -> #/m3/s)
826 0 : INnbc(i,k) =soot_num*1e6_r8
827 0 : INndust(i,k)=dst_num*1e6_r8
828 0 : INondust(i,k)=odst_num*1e6_r8
829 0 : INFreIN(i,k)=1.0_r8 ! 1,ice nucleation occur
830 0 : INhet(i,k) = (niimm(i,k) + nidep(i,k)) ! #/m3, nimey not in cirrus
831 0 : INhom(i,k) = nihf(i,k) ! #/m3
832 0 : if (INhom(i,k).gt.1e3_r8) then ! > 1/L
833 0 : INFrehom(i,k)=1.0_r8 ! 1, hom freezing occur
834 : endif
835 :
836 : ! exclude no ice nucleaton
837 0 : if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8)) then
838 0 : INnso4(i,k) =0.0_r8
839 0 : INnbc(i,k) =0.0_r8
840 0 : INndust(i,k)=0.0_r8
841 0 : INondust(i,k)=0.0_r8
842 0 : INFreIN(i,k)=0.0_r8
843 0 : INhet(i,k) = 0.0_r8
844 0 : INhom(i,k) = 0.0_r8
845 0 : INFrehom(i,k)=0.0_r8
846 0 : wice(i,k) = 0.0_r8
847 0 : weff(i,k) = 0.0_r8
848 0 : fhom(i,k) = 0.0_r8
849 : endif
850 : end if
851 :
852 : end if ! cam7
853 : end if freezing
854 : end do iloop
855 : end do kloop
856 :
857 4467528 : if (.not. clim_modal_aero) then
858 0 : deallocate( &
859 : naer2, &
860 0 : maerosol)
861 : end if
862 :
863 4467528 : if (cam_physpkg_is("cam7")) then
864 : ! Updates for PUMAS v1.21+
865 4467528 : call outfld('NIHFTEN', nihf, pcols, lchnk)
866 4467528 : call outfld('NIIMMTEN', niimm, pcols, lchnk)
867 4467528 : call outfld('NIDEPTEN', nidep, pcols, lchnk)
868 4467528 : call outfld('NIMEYTEN', nimey, pcols, lchnk)
869 : else
870 0 : call outfld('NIHF', nihf, pcols, lchnk)
871 0 : call outfld('NIIMM', niimm, pcols, lchnk)
872 0 : call outfld('NIDEP', nidep, pcols, lchnk)
873 0 : call outfld('NIMEY', nimey, pcols, lchnk)
874 : end if
875 4467528 : call outfld('NIREGM', regm, pcols, lchnk)
876 4467528 : call outfld('NISUBGRID', subgrid, pcols, lchnk)
877 4467528 : call outfld('NITROP_PD', trop_pd, pcols, lchnk)
878 :
879 4467528 : if (use_preexisting_ice) then
880 4467528 : call outfld( 'fhom' , fhom, pcols, lchnk)
881 4467528 : call outfld( 'WICE' , wice, pcols, lchnk)
882 4467528 : call outfld( 'WEFF' , weff, pcols, lchnk)
883 4467528 : if (cam_physpkg_is("cam7")) then
884 : ! Updates for PUMAS v1.21+
885 4467528 : call outfld('INnso4TEN',INnso4 , pcols,lchnk)
886 4467528 : call outfld('INnbcTEN',INnbc , pcols,lchnk)
887 4467528 : call outfld('INndustTEN',INndust, pcols,lchnk)
888 4467528 : call outfld('INondustTEN',INondust, pcols,lchnk)
889 4467528 : call outfld('INhetTEN',INhet , pcols,lchnk)
890 4467528 : call outfld('INhomTEN',INhom , pcols,lchnk)
891 : else
892 0 : call outfld('INnso4 ',INnso4 , pcols,lchnk)
893 0 : call outfld('INnbc ',INnbc , pcols,lchnk)
894 0 : call outfld('INndust ',INndust, pcols,lchnk)
895 0 : call outfld('INondust ',INondust, pcols,lchnk)
896 0 : call outfld('INhet ',INhet , pcols,lchnk)
897 0 : call outfld('INhom ',INhom , pcols,lchnk)
898 : end if
899 4467528 : call outfld('INFrehom',INFrehom,pcols,lchnk)
900 4467528 : call outfld('INFreIN ',INFreIN, pcols,lchnk)
901 : end if
902 :
903 8935056 : end subroutine nucleate_ice_cam_calc
904 :
905 : !================================================================================================
906 :
907 : end module nucleate_ice_cam
|
