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