Line data Source code
1 : module microp_aero
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : ! CAM driver layer for aerosol activation processes.
6 : !
7 : ! ***N.B.*** This module is currently hardcoded to recognize only the aerosols/modes that
8 : ! affect the climate calculation. This is implemented by using list
9 : ! index 0 in all the calls to rad_constituent interfaces.
10 : !
11 : ! Author: Andrew Gettelman
12 : ! Based on code from: Hugh Morrison, Xiaohong Liu and Steve Ghan
13 : ! May 2010
14 : ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
15 : ! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
16 : ! for questions contact Andrew Gettelman (andrew@ucar.edu)
17 : ! Modifications: A. Gettelman Nov 2010 - changed to support separation of
18 : ! microphysics and macrophysics and concentrate aerosol information here
19 : ! B. Eaton, Sep 2014 - Refactored to move CAM interface code into the CAM
20 : ! interface modules and preserve just the driver layer functionality here.
21 : !
22 : !---------------------------------------------------------------------------------
23 :
24 : use shr_kind_mod, only: r8=>shr_kind_r8
25 : use spmd_utils, only: masterproc
26 : use ppgrid, only: pcols, pver, pverp, begchunk, endchunk
27 : use ref_pres, only: top_lev => trop_cloud_top_lev
28 : use physconst, only: rair
29 : use constituents, only: cnst_get_ind
30 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init, physics_ptend_sum, &
31 : physics_state_copy, physics_update
32 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field, &
33 : pbuf_get_chunk
34 : use phys_control, only: phys_getopts, use_hetfrz_classnuc
35 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, &
36 : rad_cnst_get_mode_num
37 :
38 : use nucleate_ice_cam, only: use_preexisting_ice, nucleate_ice_cam_readnl, nucleate_ice_cam_register, &
39 : nucleate_ice_cam_init, nucleate_ice_cam_calc
40 :
41 : use ndrop, only: ndrop_init, dropmixnuc
42 : use ndrop_bam, only: ndrop_bam_init, ndrop_bam_run, ndrop_bam_ccn
43 :
44 : use hetfrz_classnuc_cam, only: hetfrz_classnuc_cam_readnl, hetfrz_classnuc_cam_register, hetfrz_classnuc_cam_init, &
45 : hetfrz_classnuc_cam_calc
46 :
47 : use cam_history, only: addfld, add_default, outfld
48 : use cam_logfile, only: iulog
49 : use cam_abortutils, only: endrun
50 :
51 : use aerosol_properties_mod, only: aerosol_properties
52 : use modal_aerosol_properties_mod, only: modal_aerosol_properties
53 :
54 : use aerosol_state_mod, only: aerosol_state
55 : use modal_aerosol_state_mod, only: modal_aerosol_state
56 :
57 : implicit none
58 : private
59 : save
60 :
61 : public :: microp_aero_init, microp_aero_run, microp_aero_readnl, microp_aero_register
62 : public :: microp_aero_final
63 : public :: aerosol_state_object
64 : public :: aerosol_properties_object
65 :
66 : ! Private module data
67 : character(len=16) :: eddy_scheme
68 : real(r8), parameter :: unset_r8 = huge(1.0_r8)
69 :
70 : ! contact freezing due to dust
71 : ! dust number mean radius (m), Zender et al JGR 2003 assuming number mode radius of 0.6 micron, sigma=2
72 : real(r8), parameter :: rn_dst1 = 0.258e-6_r8
73 : real(r8), parameter :: rn_dst2 = 0.717e-6_r8
74 : real(r8), parameter :: rn_dst3 = 1.576e-6_r8
75 : real(r8), parameter :: rn_dst4 = 3.026e-6_r8
76 :
77 : ! Namelist parameters
78 : real(r8) :: bulk_scale ! prescribed aerosol bulk sulfur scale factor
79 : real(r8) :: npccn_scale ! scaling for activated number
80 : real(r8) :: wsub_scale ! scaling for sub-grid vertical velocity (liquid)
81 : real(r8) :: wsubi_scale ! scaling for sub-grid vertical velocity (ice)
82 : real(r8) :: wsub_min ! minimum sub-grid vertical velocity (liquid) before scale factor
83 : real(r8) :: wsub_min_asf ! minimum sub-grid vertical velocity (liquid) after scale factor
84 : real(r8) :: wsubi_min ! minimum sub-grid vertical velocity (ice)
85 :
86 : ! smallest mixing ratio considered in microphysics
87 : real(r8), parameter :: qsmall = 1.e-18_r8
88 :
89 : ! minimum allowed cloud fraction
90 : real(r8), parameter :: mincld = 0.0001_r8
91 :
92 : ! indices in state%q and pbuf structures
93 : integer :: cldliq_idx = -1
94 : integer :: cldice_idx = -1
95 : integer :: numliq_idx = -1
96 : integer :: numice_idx = -1
97 : integer :: kvh_idx = -1
98 : integer :: tke_idx = -1
99 : integer :: wp2_idx = -1
100 : integer :: ast_idx = -1
101 : integer :: cldo_idx = -1
102 : integer :: dgnumwet_idx = -1
103 :
104 : ! Bulk aerosols
105 : character(len=20), allocatable :: aername(:)
106 : real(r8), allocatable :: num_to_mass_aer(:)
107 :
108 : integer :: naer_all ! number of aerosols affecting climate
109 : integer :: idxsul = -1 ! index in aerosol list for sulfate
110 : integer :: idxdst2 = -1 ! index in aerosol list for dust2
111 : integer :: idxdst3 = -1 ! index in aerosol list for dust3
112 : integer :: idxdst4 = -1 ! index in aerosol list for dust4
113 :
114 : ! modal aerosols
115 : logical :: clim_modal_aero
116 :
117 : integer :: mode_accum_idx = -1 ! index of accumulation mode
118 : integer :: mode_aitken_idx = -1 ! index of aitken mode
119 : integer :: mode_coarse_idx = -1 ! index of coarse mode
120 : integer :: mode_coarse_dst_idx = -1 ! index of coarse dust mode
121 : integer :: mode_coarse_slt_idx = -1 ! index of coarse sea salt mode
122 : integer :: coarse_dust_idx = -1 ! index of dust in coarse mode
123 : integer :: coarse_nacl_idx = -1 ! index of nacl in coarse mode
124 : integer :: coarse_so4_idx = -1 ! index of sulfate in coarse mode
125 :
126 : integer :: npccn_idx, rndst_idx, nacon_idx
127 :
128 : logical :: separate_dust = .false.
129 :
130 : type aero_state_t
131 : class(aerosol_state), pointer :: obj=>null()
132 : end type aero_state_t
133 :
134 : class(aerosol_properties), pointer :: aero_props_obj=>null()
135 : type(aero_state_t), pointer :: aero_state(:) => null()
136 :
137 : !=========================================================================================
138 : contains
139 : !=========================================================================================
140 :
141 0 : subroutine microp_aero_register
142 : !-----------------------------------------------------------------------
143 : !
144 : ! Purpose:
145 : ! Register pbuf fields for aerosols needed by microphysics
146 : !
147 : ! Author: Cheryl Craig October 2012
148 : !
149 : !-----------------------------------------------------------------------
150 : use ppgrid, only: pcols
151 : use physics_buffer, only: pbuf_add_field, dtype_r8
152 :
153 0 : call pbuf_add_field('NPCCN', 'physpkg',dtype_r8,(/pcols,pver/), npccn_idx)
154 :
155 0 : call pbuf_add_field('RNDST', 'physpkg',dtype_r8,(/pcols,pver,4/), rndst_idx)
156 0 : call pbuf_add_field('NACON', 'physpkg',dtype_r8,(/pcols,pver,4/), nacon_idx)
157 :
158 0 : call nucleate_ice_cam_register()
159 0 : call hetfrz_classnuc_cam_register()
160 :
161 0 : end subroutine microp_aero_register
162 :
163 : !=========================================================================================
164 :
165 0 : subroutine microp_aero_init(phys_state,pbuf2d)
166 :
167 : !-----------------------------------------------------------------------
168 : !
169 : ! Purpose:
170 : ! Initialize constants for aerosols needed by microphysics
171 : !
172 : ! Author: Andrew Gettelman May 2010
173 : !
174 : !-----------------------------------------------------------------------
175 :
176 : type(physics_state), pointer :: phys_state(:)
177 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
178 :
179 : ! local variables
180 : integer :: iaer, ierr
181 : integer :: m, n, nmodes, nspec
182 :
183 : character(len=32) :: str32
184 : character(len=*), parameter :: routine = 'microp_aero_init'
185 : logical :: history_amwg
186 0 : type(physics_buffer_desc), pointer :: pbuf(:)
187 : integer :: c
188 :
189 : !-----------------------------------------------------------------------
190 :
191 : ! Query the PBL eddy scheme
192 : call phys_getopts(eddy_scheme_out = eddy_scheme, &
193 0 : history_amwg_out = history_amwg )
194 :
195 : ! Access the physical properties of the aerosols that are affecting the climate
196 : ! by using routines from the rad_constituents module.
197 :
198 : ! get indices into state and pbuf structures
199 0 : call cnst_get_ind('CLDLIQ', cldliq_idx)
200 0 : call cnst_get_ind('CLDICE', cldice_idx)
201 0 : call cnst_get_ind('NUMLIQ', numliq_idx)
202 0 : call cnst_get_ind('NUMICE', numice_idx)
203 :
204 0 : select case(trim(eddy_scheme))
205 : case ('diag_TKE')
206 0 : tke_idx = pbuf_get_index('tke')
207 : case ('CLUBB_SGS')
208 0 : wp2_idx = pbuf_get_index('WP2_nadv')
209 : case default
210 0 : kvh_idx = pbuf_get_index('kvh')
211 : end select
212 :
213 : ! clim_modal_aero determines whether modal aerosols are used in the climate calculation.
214 : ! The modal aerosols can be either prognostic or prescribed.
215 0 : call rad_cnst_get_info(0, nmodes=nmodes)
216 0 : clim_modal_aero = (nmodes > 0)
217 :
218 0 : ast_idx = pbuf_get_index('AST')
219 :
220 0 : if (clim_modal_aero) then
221 :
222 0 : cldo_idx = pbuf_get_index('CLDO')
223 0 : dgnumwet_idx = pbuf_get_index('DGNUMWET')
224 :
225 0 : aero_props_obj => modal_aerosol_properties()
226 0 : if (.not.associated(aero_props_obj)) then
227 0 : call endrun('ma_convproc_init: construction of modal_aerosol_properties object failed')
228 : end if
229 0 : call ndrop_init(aero_props_obj)
230 0 : call nucleate_ice_cam_init(mincld, bulk_scale, pbuf2d, aero_props=aero_props_obj)
231 :
232 0 : allocate(aero_state(begchunk:endchunk))
233 0 : do c = begchunk,endchunk
234 0 : pbuf => pbuf_get_chunk(pbuf2d, c)
235 0 : aero_state(c)%obj => modal_aerosol_state( phys_state(c), pbuf )
236 0 : if (.not.associated(aero_state(c)%obj)) then
237 0 : call endrun('microp_aero_init: construction of modal_aerosol_state object failed')
238 : end if
239 : end do
240 :
241 : ! Init indices for specific modes/species
242 :
243 : ! mode index for specified mode types
244 0 : do m = 1, nmodes
245 0 : call rad_cnst_get_info(0, m, mode_type=str32)
246 0 : select case (trim(str32))
247 : case ('accum')
248 0 : mode_accum_idx = m
249 : case ('aitken')
250 0 : mode_aitken_idx = m
251 : case ('coarse')
252 0 : mode_coarse_idx = m
253 : case ('coarse_dust')
254 0 : mode_coarse_dst_idx = m
255 : case ('coarse_seasalt')
256 0 : mode_coarse_slt_idx = m
257 : end select
258 : end do
259 :
260 : ! check if coarse dust is in separate mode
261 0 : separate_dust = mode_coarse_dst_idx > 0
262 :
263 : ! for 3-mode
264 0 : if ( mode_coarse_dst_idx<0 ) mode_coarse_dst_idx = mode_coarse_idx
265 0 : if ( mode_coarse_slt_idx<0 ) mode_coarse_slt_idx = mode_coarse_idx
266 :
267 : ! Check that required mode types were found
268 : if (mode_accum_idx == -1 .or. mode_aitken_idx == -1 .or. &
269 0 : mode_coarse_dst_idx == -1.or. mode_coarse_slt_idx == -1) then
270 0 : write(iulog,*) routine//': ERROR required mode type not found - mode idx:', &
271 0 : mode_accum_idx, mode_aitken_idx, mode_coarse_dst_idx, mode_coarse_slt_idx
272 0 : call endrun(routine//': ERROR required mode type not found')
273 : end if
274 :
275 : ! species indices for specified types
276 : ! find indices for the dust and seasalt species in the coarse mode
277 0 : call rad_cnst_get_info(0, mode_coarse_dst_idx, nspec=nspec)
278 0 : do n = 1, nspec
279 0 : call rad_cnst_get_info(0, mode_coarse_dst_idx, n, spec_type=str32)
280 0 : select case (trim(str32))
281 : case ('dust')
282 0 : coarse_dust_idx = n
283 : end select
284 : end do
285 0 : call rad_cnst_get_info(0, mode_coarse_slt_idx, nspec=nspec)
286 0 : do n = 1, nspec
287 0 : call rad_cnst_get_info(0, mode_coarse_slt_idx, n, spec_type=str32)
288 0 : select case (trim(str32))
289 : case ('seasalt')
290 0 : coarse_nacl_idx = n
291 : end select
292 : end do
293 0 : if (mode_coarse_idx>0) then
294 0 : call rad_cnst_get_info(0, mode_coarse_idx, nspec=nspec)
295 0 : do n = 1, nspec
296 0 : call rad_cnst_get_info(0, mode_coarse_idx, n, spec_type=str32)
297 0 : select case (trim(str32))
298 : case ('sulfate')
299 0 : coarse_so4_idx = n
300 : end select
301 : end do
302 : endif
303 :
304 : ! Check that required mode specie types were found
305 0 : if ( coarse_dust_idx == -1 .or. coarse_nacl_idx == -1 ) then
306 0 : write(iulog,*) routine//': ERROR required mode-species type not found - indicies:', &
307 0 : coarse_dust_idx, coarse_nacl_idx
308 0 : call endrun(routine//': ERROR required mode-species type not found')
309 : end if
310 :
311 : else
312 :
313 : ! Props needed for BAM number concentration calcs.
314 :
315 0 : call rad_cnst_get_info(0, naero=naer_all)
316 : allocate( &
317 0 : aername(naer_all), &
318 0 : num_to_mass_aer(naer_all) )
319 :
320 0 : do iaer = 1, naer_all
321 : call rad_cnst_get_aer_props(0, iaer, &
322 0 : aername = aername(iaer), &
323 0 : num_to_mass_aer = num_to_mass_aer(iaer) )
324 :
325 : ! Look for sulfate, dust, and soot in this list (Bulk aerosol only)
326 0 : if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
327 0 : if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer
328 0 : if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer
329 0 : if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer
330 : end do
331 :
332 0 : call ndrop_bam_init()
333 0 : call nucleate_ice_cam_init(mincld, bulk_scale, pbuf2d)
334 :
335 : end if
336 :
337 0 : call addfld('LCLOUD', (/ 'lev' /), 'A', ' ', 'Liquid cloud fraction used in stratus activation')
338 :
339 0 : call addfld('WSUB', (/ 'lev' /), 'A', 'm/s', 'Diagnostic sub-grid vertical velocity' )
340 0 : call addfld('WSUBI', (/ 'lev' /), 'A', 'm/s', 'Diagnostic sub-grid vertical velocity for ice' )
341 :
342 0 : if (history_amwg) then
343 0 : call add_default ('WSUB ', 1, ' ')
344 : end if
345 :
346 0 : if (use_hetfrz_classnuc) then
347 0 : if (associated(aero_props_obj)) then
348 0 : call hetfrz_classnuc_cam_init(mincld, aero_props_obj)
349 : else
350 0 : call endrun(routine//': cannot use hetfrz_classnuc without prognostic aerosols')
351 : endif
352 : endif
353 :
354 0 : end subroutine microp_aero_init
355 :
356 : !=========================================================================================
357 : ! returns a pointer to an aerosol state object for a given chunk index
358 0 : function aerosol_state_object(lchnk) result(obj)
359 :
360 : integer,intent(in) :: lchnk ! local chunk index
361 : class(aerosol_state), pointer :: obj ! aerosol state object pointer for local chunk
362 :
363 0 : obj => aero_state(lchnk)%obj
364 :
365 0 : end function aerosol_state_object
366 :
367 : !=========================================================================================
368 : ! returns a pointer to an aerosol properties object
369 1536 : function aerosol_properties_object() result(obj)
370 :
371 : class(aerosol_properties), pointer :: obj ! aerosol properties object pointer
372 :
373 1536 : obj => aero_props_obj
374 :
375 1536 : end function aerosol_properties_object
376 :
377 : !=========================================================================================
378 :
379 1536 : subroutine microp_aero_final
380 :
381 : integer :: c
382 :
383 1536 : if (associated(aero_props_obj)) then
384 0 : deallocate(aero_props_obj)
385 : end if
386 1536 : nullify(aero_props_obj)
387 :
388 1536 : if (associated(aero_state)) then
389 0 : do c = begchunk,endchunk
390 0 : deallocate(aero_state(c)%obj)
391 : end do
392 0 : deallocate(aero_state)
393 : nullify(aero_state)
394 : end if
395 :
396 1536 : end subroutine microp_aero_final
397 :
398 : !=========================================================================================
399 :
400 1536 : subroutine microp_aero_readnl(nlfile)
401 :
402 : use namelist_utils, only: find_group_name
403 : use units, only: getunit, freeunit
404 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8
405 :
406 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
407 :
408 : ! Namelist variables
409 : real(r8) :: microp_aero_bulk_scale = unset_r8 ! prescribed aerosol bulk sulfur scale factor
410 : real(r8) :: microp_aero_npccn_scale = unset_r8 ! prescribed aerosol bulk sulfur scale factor
411 : real(r8) :: microp_aero_wsub_scale = unset_r8 ! subgrid vertical velocity (liquid) scale factor
412 : real(r8) :: microp_aero_wsubi_scale = unset_r8 ! subgrid vertical velocity (ice) scale factor
413 : real(r8) :: microp_aero_wsub_min = unset_r8 ! subgrid vertical velocity (liquid) minimum (before scale factor)
414 : real(r8) :: microp_aero_wsub_min_asf = unset_r8 ! subgrid vertical velocity (liquid) minimum (after scale factor)
415 : real(r8) :: microp_aero_wsubi_min = unset_r8 ! subgrid vertical velocity (ice) minimum
416 :
417 : ! Local variables
418 : integer :: unitn, ierr
419 : character(len=*), parameter :: subname = 'microp_aero_readnl'
420 :
421 : namelist /microp_aero_nl/ microp_aero_bulk_scale, microp_aero_npccn_scale, microp_aero_wsub_min, &
422 : microp_aero_wsubi_min, microp_aero_wsub_scale, microp_aero_wsubi_scale, microp_aero_wsub_min_asf
423 : !-----------------------------------------------------------------------------
424 :
425 1536 : if (masterproc) then
426 2 : unitn = getunit()
427 2 : open( unitn, file=trim(nlfile), status='old' )
428 2 : call find_group_name(unitn, 'microp_aero_nl', status=ierr)
429 2 : if (ierr == 0) then
430 2 : read(unitn, microp_aero_nl, iostat=ierr)
431 2 : if (ierr /= 0) then
432 0 : call endrun(subname // ':: ERROR reading namelist')
433 : end if
434 : end if
435 2 : close(unitn)
436 2 : call freeunit(unitn)
437 : end if
438 :
439 : ! Broadcast namelist variables
440 1536 : call mpi_bcast(microp_aero_bulk_scale, 1, mpi_real8, mstrid, mpicom, ierr)
441 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_bulk_scale")
442 1536 : call mpi_bcast(microp_aero_npccn_scale, 1, mpi_real8, mstrid, mpicom, ierr)
443 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_npccn_scale")
444 1536 : call mpi_bcast(microp_aero_wsub_scale, 1, mpi_real8, mstrid, mpicom, ierr)
445 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsub_scale")
446 1536 : call mpi_bcast(microp_aero_wsubi_scale, 1, mpi_real8, mstrid, mpicom, ierr)
447 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsubi_scale")
448 1536 : call mpi_bcast(microp_aero_wsub_min, 1, mpi_real8, mstrid, mpicom, ierr)
449 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsub_min")
450 1536 : call mpi_bcast(microp_aero_wsub_min_asf, 1, mpi_real8, mstrid, mpicom, ierr)
451 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsub_min_asf")
452 1536 : call mpi_bcast(microp_aero_wsubi_min, 1, mpi_real8, mstrid, mpicom, ierr)
453 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsubi_min")
454 :
455 : ! set local variables
456 1536 : bulk_scale = microp_aero_bulk_scale
457 1536 : npccn_scale = microp_aero_npccn_scale
458 1536 : wsub_scale = microp_aero_wsub_scale
459 1536 : wsubi_scale = microp_aero_wsubi_scale
460 1536 : wsub_min = microp_aero_wsub_min
461 1536 : wsub_min_asf = microp_aero_wsub_min_asf
462 1536 : wsubi_min = microp_aero_wsubi_min
463 :
464 1536 : if(bulk_scale == unset_r8) call endrun(subname//": FATAL: bulk_scale is not set")
465 1536 : if(npccn_scale == unset_r8) call endrun(subname//": FATAL: npccn_scale is not set")
466 1536 : if(wsub_scale == unset_r8) call endrun(subname//": FATAL: wsub_scale is not set")
467 1536 : if(wsubi_scale == unset_r8) call endrun(subname//": FATAL: wsubi_scale is not set")
468 1536 : if(wsub_min == unset_r8) call endrun(subname//": FATAL: wsub_min is not set")
469 1536 : if(wsub_min_asf == unset_r8) call endrun(subname//": FATAL: wsub_min_asf is not set")
470 1536 : if(wsubi_min == unset_r8) call endrun(subname//": FATAL: wsubi_min is not set")
471 :
472 1536 : call nucleate_ice_cam_readnl(nlfile)
473 1536 : call hetfrz_classnuc_cam_readnl(nlfile)
474 :
475 1536 : end subroutine microp_aero_readnl
476 :
477 : !=========================================================================================
478 :
479 0 : subroutine microp_aero_run ( &
480 : state, ptend_all, deltatin, pbuf)
481 :
482 : ! input arguments
483 : type(physics_state), intent(in) :: state
484 : type(physics_ptend), intent(out) :: ptend_all
485 : real(r8), intent(in) :: deltatin ! time step (s)
486 : type(physics_buffer_desc), pointer :: pbuf(:)
487 :
488 : ! local workspace
489 : ! all units mks unless otherwise stated
490 :
491 : integer :: i, k, m
492 : integer :: itim_old
493 :
494 0 : type(physics_state), target :: state1 ! Local copy of state variable
495 0 : type(physics_ptend) :: ptend_loc
496 :
497 0 : real(r8), pointer :: ast(:,:)
498 :
499 0 : real(r8), pointer :: npccn(:,:) ! number of CCN (liquid activated)
500 :
501 0 : real(r8), pointer :: rndst(:,:,:) ! radius of 4 dust bins for contact freezing
502 0 : real(r8), pointer :: nacon(:,:,:) ! number in 4 dust bins for contact freezing
503 :
504 0 : real(r8), pointer :: num_coarse(:,:) ! number m.r. of coarse mode
505 0 : real(r8), pointer :: coarse_dust(:,:) ! mass m.r. of coarse dust
506 0 : real(r8), pointer :: coarse_nacl(:,:) ! mass m.r. of coarse nacl
507 0 : real(r8), pointer :: coarse_so4(:,:) ! mass m.r. of coarse sulfate
508 :
509 0 : real(r8), pointer :: kvh(:,:) ! vertical eddy diff coef (m2 s-1)
510 0 : real(r8), pointer :: tke(:,:) ! TKE from the UW PBL scheme (m2 s-2)
511 0 : real(r8), pointer :: wp2(:,:) ! CLUBB vertical velocity variance
512 :
513 0 : real(r8), pointer :: cldn(:,:) ! cloud fraction
514 0 : real(r8), pointer :: cldo(:,:) ! old cloud fraction
515 :
516 0 : real(r8), pointer :: dgnumwet(:,:,:) ! aerosol mode diameter
517 :
518 0 : real(r8), pointer :: aer_mmr(:,:) ! aerosol mass mixing ratio
519 :
520 : real(r8) :: rho(pcols,pver) ! air density (kg m-3)
521 :
522 : real(r8) :: lcldm(pcols,pver) ! liq cloud fraction
523 :
524 : real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud
525 : real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud
526 : real(r8) :: cldliqf(pcols,pver) ! fractional of total cloud that is liquid
527 : real(r8) :: qcld ! total cloud water
528 : real(r8) :: nctend_mixnuc(pcols,pver)
529 : real(r8) :: dum, dum2 ! temporary dummy variable
530 : real(r8) :: dmc, ssmc, so4mc ! variables for modal scheme.
531 :
532 : ! bulk aerosol variables
533 0 : real(r8), allocatable :: naer2(:,:,:) ! bulk aerosol number concentration (1/m3)
534 0 : real(r8), allocatable :: maerosol(:,:,:) ! bulk aerosol mass conc (kg/m3)
535 :
536 : real(r8) :: wsub(pcols,pver) ! diagnosed sub-grid vertical velocity st. dev. (m/s)
537 : real(r8) :: wsubi(pcols,pver) ! diagnosed sub-grid vertical velocity ice (m/s)
538 : real(r8) :: nucboas
539 :
540 : real(r8) :: wght
541 :
542 : integer :: lchnk, ncol, astat
543 :
544 0 : real(r8), allocatable :: factnum(:,:,:) ! activation fraction for aerosol number
545 :
546 : class(aerosol_state), pointer :: aero_state1_obj
547 :
548 : !-------------------------------------------------------------------------------
549 :
550 0 : nullify(aero_state1_obj)
551 :
552 0 : call physics_state_copy(state,state1)
553 :
554 0 : lchnk = state1%lchnk
555 0 : ncol = state1%ncol
556 :
557 0 : itim_old = pbuf_old_tim_idx()
558 0 : call pbuf_get_field(pbuf, ast_idx, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
559 :
560 0 : call pbuf_get_field(pbuf, npccn_idx, npccn)
561 :
562 0 : call pbuf_get_field(pbuf, nacon_idx, nacon)
563 0 : call pbuf_get_field(pbuf, rndst_idx, rndst)
564 :
565 0 : call physics_ptend_init(ptend_all, state%psetcols, 'microp_aero')
566 :
567 0 : if (clim_modal_aero) then
568 : ! create an aerosol state object specifically for cam state1
569 0 : aero_state1_obj => modal_aerosol_state( state1, pbuf )
570 0 : if (.not.associated(aero_state1_obj)) then
571 0 : call endrun('microp_aero_run: construction of aero_state1_obj modal_aerosol_state object failed')
572 : end if
573 :
574 0 : itim_old = pbuf_old_tim_idx()
575 :
576 0 : call pbuf_get_field(pbuf, ast_idx, cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
577 0 : call pbuf_get_field(pbuf, cldo_idx, cldo, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
578 0 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet)
579 : end if
580 :
581 : ! initialize output
582 0 : npccn(1:ncol,1:pver) = 0._r8
583 :
584 0 : nacon(1:ncol,1:pver,:) = 0._r8
585 :
586 : ! set default or fixed dust bins for contact freezing
587 0 : rndst(1:ncol,1:pver,1) = rn_dst1
588 0 : rndst(1:ncol,1:pver,2) = rn_dst2
589 0 : rndst(1:ncol,1:pver,3) = rn_dst3
590 0 : rndst(1:ncol,1:pver,4) = rn_dst4
591 :
592 : ! initialize time-varying parameters
593 0 : do k = top_lev, pver
594 0 : do i = 1, ncol
595 0 : rho(i,k) = state1%pmid(i,k)/(rair*state1%t(i,k))
596 : end do
597 : end do
598 :
599 0 : if (clim_modal_aero) then
600 : ! mode number mixing ratios
601 0 : call rad_cnst_get_mode_num(0, mode_coarse_dst_idx, 'a', state1, pbuf, num_coarse)
602 :
603 : ! mode specie mass m.r.
604 0 : call rad_cnst_get_aer_mmr(0, mode_coarse_dst_idx, coarse_dust_idx, 'a', state1, pbuf, coarse_dust)
605 0 : call rad_cnst_get_aer_mmr(0, mode_coarse_slt_idx, coarse_nacl_idx, 'a', state1, pbuf, coarse_nacl)
606 0 : if (mode_coarse_idx>0) then
607 0 : call rad_cnst_get_aer_mmr(0, mode_coarse_idx, coarse_so4_idx, 'a', state1, pbuf, coarse_so4)
608 : endif
609 :
610 : else
611 : ! init number/mass arrays for bulk aerosols
612 : allocate( &
613 0 : naer2(pcols,pver,naer_all), &
614 0 : maerosol(pcols,pver,naer_all))
615 :
616 0 : do m = 1, naer_all
617 0 : call rad_cnst_get_aer_mmr(0, m, state1, pbuf, aer_mmr)
618 0 : maerosol(:ncol,:,m) = aer_mmr(:ncol,:)*rho(:ncol,:)
619 :
620 0 : if (m .eq. idxsul) then
621 0 : naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)*bulk_scale
622 : else
623 0 : naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)
624 : end if
625 : end do
626 : end if
627 :
628 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
629 : ! More refined computation of sub-grid vertical velocity
630 : ! Set to be zero at the surface by initialization.
631 :
632 0 : select case (trim(eddy_scheme))
633 : case ('diag_TKE')
634 0 : call pbuf_get_field(pbuf, tke_idx, tke)
635 : case ('CLUBB_SGS')
636 0 : itim_old = pbuf_old_tim_idx()
637 0 : call pbuf_get_field(pbuf, wp2_idx, wp2, start=(/1,1,itim_old/),kount=(/pcols,pverp,1/))
638 0 : allocate(tke(pcols,pverp))
639 0 : tke(:ncol,:) = (3._r8/2._r8)*wp2(:ncol,:)
640 :
641 : case default
642 0 : call pbuf_get_field(pbuf, kvh_idx, kvh)
643 : end select
644 :
645 : ! Set minimum values above top_lev.
646 0 : wsub(:ncol,:top_lev-1) = wsub_min
647 0 : wsubi(:ncol,:top_lev-1) = wsubi_min
648 :
649 0 : do k = top_lev, pver
650 0 : do i = 1, ncol
651 :
652 0 : select case (trim(eddy_scheme))
653 : case ('diag_TKE', 'CLUBB_SGS')
654 0 : wsub(i,k) = sqrt(0.5_r8*(tke(i,k) + tke(i,k+1))*(2._r8/3._r8))
655 0 : wsub(i,k) = min(wsub(i,k),10._r8)
656 : case default
657 : ! get sub-grid vertical velocity from diff coef.
658 : ! following morrison et al. 2005, JAS
659 : ! assume mixing length of 30 m
660 0 : dum = (kvh(i,k) + kvh(i,k+1))/2._r8/30._r8
661 : ! use maximum sub-grid vertical vel of 10 m/s
662 0 : dum = min(dum, 10._r8)
663 : ! set wsub to value at current vertical level
664 0 : wsub(i,k) = dum
665 : end select
666 :
667 0 : wsubi(i,k) = max(wsubi_min, wsub(i,k)) * wsubi_scale
668 0 : if (.not. use_preexisting_ice) then
669 0 : wsubi(i,k) = min(wsubi(i,k), 0.2_r8)
670 : endif
671 :
672 0 : wsub(i,k) = max(wsub_min, wsub(i,k)) * wsub_scale
673 :
674 : end do
675 : end do
676 :
677 0 : call outfld('WSUB', wsub, pcols, lchnk)
678 0 : call outfld('WSUBI', wsubi, pcols, lchnk)
679 :
680 0 : if (trim(eddy_scheme) == 'CLUBB_SGS') deallocate(tke)
681 :
682 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
683 : !ICE Nucleation
684 :
685 0 : if (associated(aero_props_obj).and.associated(aero_state1_obj)) then
686 0 : call nucleate_ice_cam_calc(state1, wsubi, pbuf, deltatin, ptend_loc, aero_props_obj, aero_state1_obj)
687 : else
688 0 : call nucleate_ice_cam_calc(state1, wsubi, pbuf, deltatin, ptend_loc)
689 : end if
690 :
691 0 : call physics_ptend_sum(ptend_loc, ptend_all, ncol)
692 0 : call physics_update(state1, ptend_loc, deltatin)
693 :
694 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
695 : ! get liquid cloud fraction, check for minimum
696 :
697 0 : do k = top_lev, pver
698 0 : do i = 1, ncol
699 0 : lcldm(i,k) = max(ast(i,k), mincld)
700 : end do
701 : end do
702 :
703 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
704 : ! Droplet Activation
705 :
706 0 : if (clim_modal_aero) then
707 :
708 : ! for modal aerosol
709 :
710 : ! partition cloud fraction into liquid water part
711 0 : lcldn = 0._r8
712 0 : lcldo = 0._r8
713 0 : cldliqf = 0._r8
714 0 : do k = top_lev, pver
715 0 : do i = 1, ncol
716 0 : qcld = state1%q(i,k,cldliq_idx) + state1%q(i,k,cldice_idx)
717 0 : if (qcld > qsmall) then
718 0 : lcldn(i,k) = cldn(i,k)*state1%q(i,k,cldliq_idx)/qcld
719 0 : lcldo(i,k) = cldo(i,k)*state1%q(i,k,cldliq_idx)/qcld
720 0 : cldliqf(i,k) = state1%q(i,k,cldliq_idx)/qcld
721 : end if
722 : end do
723 : end do
724 :
725 0 : call outfld('LCLOUD', lcldn, pcols, lchnk)
726 :
727 0 : allocate(factnum(pcols,pver,aero_props_obj%nbins()),stat=astat)
728 0 : if (astat/=0) then
729 0 : call endrun('microp_aero_run: not able to allocate factnum')
730 : endif
731 :
732 : ! If not using preexsiting ice, then only use cloudbourne aerosol for the
733 : ! liquid clouds. This is the same behavior as CAM5.
734 0 : if (use_preexisting_ice) then
735 : call dropmixnuc( aero_props_obj, aero_state1_obj, &
736 : state1, ptend_loc, deltatin, pbuf, wsub, wsub_min_asf, &
737 0 : cldn, cldo, cldliqf, nctend_mixnuc, factnum)
738 : else
739 0 : cldliqf = 1._r8
740 : call dropmixnuc( aero_props_obj, aero_state1_obj, &
741 : state1, ptend_loc, deltatin, pbuf, wsub, wsub_min_asf, &
742 0 : lcldn, lcldo, cldliqf, nctend_mixnuc, factnum)
743 : end if
744 :
745 0 : npccn(:ncol,:) = nctend_mixnuc(:ncol,:)
746 :
747 0 : npccn(:ncol,:) = npccn(:ncol,:) * npccn_scale
748 :
749 : else
750 :
751 : ! for bulk aerosol
752 :
753 : ! no tendencies returned from ndrop_bam_run, so just init ptend here
754 0 : call physics_ptend_init(ptend_loc, state1%psetcols, 'none')
755 :
756 0 : do k = top_lev, pver
757 0 : do i = 1, ncol
758 :
759 0 : if (state1%q(i,k,cldliq_idx) >= qsmall) then
760 :
761 : ! get droplet activation rate
762 :
763 : call ndrop_bam_run( &
764 0 : wsub(i,k), state1%t(i,k), rho(i,k), naer2(i,k,:), naer_all, &
765 0 : naer_all, maerosol(i,k,:), &
766 0 : dum2)
767 0 : dum = dum2
768 : else
769 : dum = 0._r8
770 : end if
771 :
772 0 : npccn(i,k) = (dum*lcldm(i,k) - state1%q(i,k,numliq_idx))/deltatin
773 : end do
774 : end do
775 :
776 : end if
777 :
778 0 : call physics_ptend_sum(ptend_loc, ptend_all, ncol)
779 0 : call physics_update(state1, ptend_loc, deltatin)
780 :
781 :
782 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
783 : ! Contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
784 : ! estimate rndst and nanco for 4 dust bins here to pass to MG microphysics
785 :
786 0 : do k = top_lev, pver
787 0 : do i = 1, ncol
788 :
789 0 : if (state1%t(i,k) < 269.15_r8) then
790 :
791 0 : if (clim_modal_aero) then
792 :
793 : ! For modal aerosols:
794 : ! use size '3' for dust coarse mode...
795 : ! scale by dust fraction in coarse mode
796 :
797 0 : dmc = coarse_dust(i,k)
798 0 : ssmc = coarse_nacl(i,k)
799 :
800 0 : if ( separate_dust ) then
801 : ! 7-mode -- has separate dust and seasalt mode types and no need for weighting
802 : wght = 1._r8
803 : else
804 0 : so4mc = coarse_so4(i,k)
805 : ! 3-mode -- needs weighting for dust since dust, seasalt, and sulfate are combined in the "coarse" mode type
806 0 : wght = dmc/(ssmc + dmc + so4mc)
807 : endif
808 :
809 0 : if (dmc > 0.0_r8) then
810 0 : nacon(i,k,3) = wght*num_coarse(i,k)*rho(i,k)
811 : else
812 0 : nacon(i,k,3) = 0._r8
813 : end if
814 :
815 : !also redefine parameters based on size...
816 :
817 0 : rndst(i,k,3) = 0.5_r8*dgnumwet(i,k,mode_coarse_dst_idx)
818 0 : if (rndst(i,k,3) <= 0._r8) then
819 0 : rndst(i,k,3) = rn_dst3
820 : end if
821 :
822 : else
823 :
824 : !For Bulk Aerosols: set equal to aerosol number for dust for bins 2-4 (bin 1=0)
825 :
826 0 : if (idxdst2 > 0) then
827 0 : nacon(i,k,2) = naer2(i,k,idxdst2)
828 : end if
829 0 : if (idxdst3 > 0) then
830 0 : nacon(i,k,3) = naer2(i,k,idxdst3)
831 : end if
832 0 : if (idxdst4 > 0) then
833 0 : nacon(i,k,4) = naer2(i,k,idxdst4)
834 : end if
835 : end if
836 :
837 : end if
838 : end do
839 : end do
840 :
841 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
842 : !bulk aerosol ccn concentration (modal does it in ndrop, from dropmixnuc)
843 :
844 0 : if (.not. clim_modal_aero) then
845 :
846 : ! ccn concentration as diagnostic
847 0 : call ndrop_bam_ccn(lchnk, ncol, maerosol, naer2)
848 :
849 0 : deallocate( &
850 : naer2, &
851 0 : maerosol)
852 :
853 : end if
854 :
855 : ! heterogeneous freezing
856 0 : if (use_hetfrz_classnuc) then
857 :
858 0 : call hetfrz_classnuc_cam_calc(aero_props_obj, aero_state1_obj, state1, deltatin, factnum, pbuf)
859 :
860 : end if
861 :
862 0 : if (clim_modal_aero) then
863 0 : deallocate(factnum)
864 : end if
865 :
866 0 : if (associated(aero_state1_obj)) then
867 : ! destroy the aerosol state object
868 0 : deallocate(aero_state1_obj)
869 : nullify(aero_state1_obj)
870 : endif
871 :
872 0 : end subroutine microp_aero_run
873 :
874 : !=========================================================================================
875 :
876 0 : end module microp_aero
|
