mam  v1.0
A Modal Aerosol Model
mode.F90
Go to the documentation of this file.
1! Copyright (C) 2021 National Center for Atmospheric Research
2! SPDX-License-Identifier: Apache-2.0
3!
4!> \file
5!> The mam_mode module
6
7!> The mode_t type and related functions
8module mam_mode
9
10 use ai_aerosol, only : aerosol_t
11 use ai_aerosol_state, only : aerosol_state_t
12 use mam_optics_lookup, only : optics_lookup_t
13 use mam_species, only : species_t
14 use musica_constants, only : musica_dk
15
16 implicit none
17 private
18
19 public :: mode_t, mode_state_t
20
21 !> An aerosol mode
22 type, extends(aerosol_t) :: mode_t
23 private
24 !> Geometric mean diameter [m] of number distribution
25 real(kind=musica_dk) :: geometric_mean_diameter__m_
26 !> Geometric standard deviation of number distribution
27 real(kind=musica_dk) :: geometric_standard_deviation_
28 !> Chemical species that can be present in the mode
29 type(species_t), allocatable :: species_(:)
30 !> Shortwave optical property lookup table
31 type(optics_lookup_t) :: shortwave_lookup_
32 !> Longwave optical property lookup table
33 type(optics_lookup_t) :: longwave_lookup_
34 contains
35 procedure :: new_state
36 procedure :: new_optics
37 procedure, private :: shortwave_optics_scalar
38 procedure, private :: shortwave_optics_array
39 procedure, private :: longwave_optics_scalar
40 procedure, private :: longwave_optics_array
41 procedure :: shortwave_grid
42 procedure :: longwave_grid
43 procedure :: print_state
44 procedure :: add_shortwave_optics
45 procedure :: add_longwave_optics
46 procedure :: geometric_mean_diameter_of_number_distribution__m
47 procedure :: geometric_standard_deviation_of_number_distribution
48 procedure :: wet_volume_to_mass_mixing_ratio__m3_kg
49 procedure :: wet_number_mode_radius__m
50 procedure :: wet_number_mode_diameter__m
51 procedure :: wet_surface_mode_radius__m
52 procedure :: wet_surface_mode_diameter__m
53 procedure :: number_mixing_ratio__num_mol
54 procedure :: net_shortwave_refractive_index
55 procedure :: net_longwave_refractive_index
56 procedure :: specific_absorption__m2_kg
57 procedure :: specific_extinction__m2_kg
58 procedure :: asymmetry_factor
59 end type mode_t
60
61 interface mode_t
62 module procedure :: constructor
63 end interface mode_t
64
65 !> Aerosol mode state
66 type, extends(aerosol_state_t) :: mode_state_t
67 private
68 !> Mode diameter of the wet aerosol number distribution [m]
69 real(kind=musica_dk) :: wet_number_mode_diameter__m_
70 !> Particle number mixing ratio per mol_air [# mol-1]
71 !! \todo Is this dry air or wet air?
72 real(kind=musica_dk) :: number_mixing_ratio__num_mol_
73 !> Mass mixing ratio of each mode species relative to dry air [kg kg-1]
74 real(kind=musica_dk), allocatable :: mass_mixing_ratio__kg_kg_(:)
75 contains
76 procedure :: raw_size
77 procedure :: load_state
78 procedure :: dump_state
79 procedure :: randomize
80 end type mode_state_t
81
82contains
83
84!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
85
86 !> Constructor of mode_t objects
87 function constructor( config ) result( new_obj )
88
89 use musica_config, only : config_t
90 use musica_iterator, only : iterator_t
91
92 type(mode_t) :: new_obj
93 class(config_t), intent(inout) :: config
94
95 character(len=*), parameter :: my_name = "MAM mode_t constructor"
96 type(config_t) :: optics_config, species_set, species
97 class(iterator_t), pointer :: iter
98 integer :: i_species
99
100 call config%get( "geometric mean diameter of number distribution", "m", &
101 new_obj%geometric_mean_diameter__m_, my_name )
102 call config%get( "geometric standard deviation of number distribution", &
103 new_obj%geometric_standard_deviation_, my_name )
104 call config%get( "species", species_set, my_name )
105 allocate( new_obj%species_( species_set%number_of_children( ) ) )
106 iter => species_set%get_iterator( )
107 i_species = 1
108 do while( iter%next( ) )
109 call species_set%get( iter, species, my_name )
110 call species%add( "name", species_set%key( iter ), my_name )
111 new_obj%species_( i_species ) = species_t( species )
112 i_species = i_species + 1
113 end do
114 call config%get( "shortwave lookup tables", optics_config, my_name )
115 new_obj%shortwave_lookup_ = optics_lookup_t( optics_config )
116 call config%get( "longwave lookup tables", optics_config, my_name )
117 new_obj%longwave_lookup_ = optics_lookup_t( optics_config )
118 deallocate( iter )
119
120 end function constructor
121
122!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
123
124 !> Creates a new state object for the mode
125 function new_state( this )
126
127 use ai_aerosol_state, only : aerosol_state_t
128
129 class(aerosol_state_t), pointer :: new_state
130 class(mode_t), intent(in) :: this
131
132 allocate( mode_state_t :: new_state )
133 select type( new_state )
134 type is( mode_state_t )
135 allocate( new_state%mass_mixing_ratio__kg_kg_( size( this%species_ ) ) )
136 end select
137
138 end function new_state
139
140!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141
142 !> Creates an optics_t object for a given property
143 function new_optics( this, property, output_grid, interpolation_strategy )
144
145 use ai_optics, only : optics_t
146 use ai_wavelength_grid, only : wavelength_grid_t
147 use mam_optics_util, only : create_optics
148 use musica_interpolator, only : interpolation_strategy_i
149 use musica_property, only : property_t
150
151 class(optics_t), pointer :: new_optics
152 class(mode_t), intent(in) :: this
153 class(property_t), intent(in) :: property
154 type(wavelength_grid_t), intent(in) :: output_grid
155 procedure(interpolation_strategy_i), optional :: interpolation_strategy
156
157 new_optics => &
158 create_optics( property, this%shortwave_lookup_%grid( ), &
159 this%longwave_lookup_%grid( ), output_grid, &
160 interpolation_strategy )
161
162 end function new_optics
163
164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165
166 !> Returns shortwave optical properties
167 subroutine shortwave_optics_scalar( this, environmental_state, &
168 aerosol_state, optics )
169
170 use ai_environmental_state, only : environmental_state_t
171 use ai_optics, only : optics_t, optics_ptr
172
173 class(mode_t), intent(in) :: this
174 class(environmental_state_t), intent(in) :: environmental_state
175 class(aerosol_state_t), intent(in) :: aerosol_state
176 class(optics_t), target, intent(inout) :: optics
177
178 type(optics_ptr) :: optics_set(1)
179
180 optics_set(1)%ptr_ => optics
181 call this%shortwave_optics( environmental_state, aerosol_state, &
182 optics_set )
183 nullify( optics_set(1)%ptr_ )
184
185 end subroutine shortwave_optics_scalar
186
187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188
189 !> Returns shortwave optical properties
190 subroutine shortwave_optics_array( this, environmental_state, &
191 aerosol_state, optics )
192
193 use ai_environmental_state, only : environmental_state_t
194 use ai_optics, only : optics_ptr
195 use musica_assert, only : die_msg
196
197 class(mode_t), intent(in) :: this
198 class(environmental_state_t), intent(in) :: environmental_state
199 class(aerosol_state_t), intent(in) :: aerosol_state
200 type(optics_ptr), intent(inout) :: optics(:)
201
202 integer :: i_prop
203
204 select type( aerosol_state )
205 class is( mode_state_t )
206 do i_prop = 1, size( optics )
207 call optics( i_prop )%ptr_%reset_values( )
208 end do
209 call this%add_shortwave_optics( environmental_state, aerosol_state, &
210 optics )
211 class default
212 call die_msg( 713383720, "Invalid state passed to MAM mode shortwave "//&
213 "optics calculator" )
214 end select
215
216 end subroutine shortwave_optics_array
217
218!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
219
220 !> Returns longwave optical properties
221 subroutine longwave_optics_scalar( this, environmental_state, &
222 aerosol_state, optics )
223
224 use ai_environmental_state, only : environmental_state_t
225 use ai_optics, only : optics_t, optics_ptr
226
227 class(mode_t), intent(in) :: this
228 class(environmental_state_t), intent(in) :: environmental_state
229 class(aerosol_state_t), intent(in) :: aerosol_state
230 class(optics_t), target, intent(inout) :: optics
231
232 type(optics_ptr) :: optics_set(1)
233
234 optics_set(1)%ptr_ => optics
235 call this%longwave_optics( environmental_state, aerosol_state, &
236 optics_set )
237 nullify( optics_set(1)%ptr_ )
238
239 end subroutine longwave_optics_scalar
240
241!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
242
243 !> Returns longwave optical properties
244 subroutine longwave_optics_array( this, environmental_state, &
245 aerosol_state, optics )
246
247 use ai_environmental_state, only : environmental_state_t
248 use ai_optics, only : optics_ptr
249 use musica_assert, only : die_msg
250
251 class(mode_t), intent(in) :: this
252 class(environmental_state_t), intent(in) :: environmental_state
253 class(aerosol_state_t), intent(in) :: aerosol_state
254 type(optics_ptr), intent(inout) :: optics(:)
255
256 integer :: i_prop
257
258 select type( aerosol_state )
259 class is( mode_state_t )
260 do i_prop = 1, size( optics )
261 call optics( i_prop )%ptr_%reset_values( )
262 end do
263 call this%add_longwave_optics( environmental_state, aerosol_state, &
264 optics )
265 class default
266 call die_msg( 713383720, "Invalid state passed to MAM mode longwave "// &
267 "optics calculator" )
268 end select
269
270 end subroutine longwave_optics_array
271
272!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273
274 !> Returns the native shortwave optics wavelength grid
275 function shortwave_grid( this )
276
277 use ai_wavelength_grid, only : wavelength_grid_t
278
279 type(wavelength_grid_t) :: shortwave_grid
280 class(mode_t), intent(in) :: this
281
282 shortwave_grid = this%shortwave_lookup_%grid( )
283
284 end function shortwave_grid
285
286!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
287
288 !> Returns the native longwave optics wavelength grid
289 function longwave_grid( this )
290
291 use ai_wavelength_grid, only : wavelength_grid_t
292
293 type(wavelength_grid_t) :: longwave_grid
294 class(mode_t), intent(in) :: this
295
296 longwave_grid = this%longwave_lookup_%grid( )
297
298 end function longwave_grid
299
300!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
301
302 !> Outputs the current mode state
303 subroutine print_state( this, aerosol_state, io_unit )
304
305 use musica_assert, only : die_msg
306 use musica_string, only : to_char
307
308 !> MAM mode
309 class(mode_t), intent(in) :: this
310 !> MAM mode state
311 class(aerosol_state_t), intent(in) :: aerosol_state
312 !> Optional output unit (defaults to 6)
313 integer, optional, intent(in) :: io_unit
314
315 integer :: lunit, i_species
316
317 lunit = 6
318 if( present( io_unit ) ) lunit = io_unit
319 select type( aerosol_state )
320 class is( mode_state_t )
321 write(lunit,*) "** MAM Mode State **"
322 if( .not. allocated( this%species_ ) ) then
323 write(lunit,*) "--- Uninitialized MAM Mode ---"
324 write(lunit,*) "** End MAM Mode State **"
325 return
326 end if
327 if( .not. allocated( aerosol_state%mass_mixing_ratio__kg_kg_ ) ) then
328 write(lunit,*) "--- Uninitialized MAM Mode State ---"
329 write(lunit,*) "** End MAM Mode State **"
330 return
331 end if
332 write(lunit,*) "wet mode diameter of number distribution [m]: "// &
333 trim( to_char( &
334 aerosol_state%wet_number_mode_diameter__m_ ) )
335 write(lunit,*) "number mixing ratio [# mol]: "//trim( to_char( &
336 aerosol_state%number_mixing_ratio__num_mol_ ) )
337 do i_species = 1, size( this%species_ )
338 write(lunit,*) "species '"//this%species_( i_species )%name( )// &
339 "' mass mixing ratio [kg kg-1]: "// &
340 trim( to_char( &
341 aerosol_state%mass_mixing_ratio__kg_kg_( i_species ) ) )
342 end do
343 class default
344 call die_msg( 285756314, "Invalid state passed to MAM mode" )
345 end select
346
347 end subroutine print_state
348
349!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
350
351 !> Calculates shortwave optical properties for a given state and adds them
352 !! to the set of optical properties
353 subroutine add_shortwave_optics( this, environmental_state, mode_state, &
354 optics )
355
356 use ai_environmental_state, only : environmental_state_t
357 use ai_optics, only : optics_ptr
358 use mam_constants, only : kaccellerationbygravity
359 use mam_optics_util, only : add_shortwave_property
360 use musica_math, only : chebyshev_function
361
362 class(mode_t), intent(in) :: this
363 class(environmental_state_t), intent(in) :: environmental_state
364 class(mode_state_t), intent(in) :: mode_state
365 class(optics_ptr), intent(inout) :: optics(:)
366
367 integer :: i_prop, n_cheby, n_bands
368 real(kind=musica_dk), allocatable :: absorption_coefficients(:,:)
369 real(kind=musica_dk), allocatable :: extinction_coefficients(:,:)
370 real(kind=musica_dk), allocatable :: asymmetry_factor_coefficients(:,:)
371 real(kind=musica_dk), allocatable :: size_function(:)
372 complex(kind=musica_dk), allocatable :: net_refractive_index(:)
373 real(kind=musica_dk), allocatable :: absorption(:)
374 real(kind=musica_dk), allocatable :: extinction(:)
375 real(kind=musica_dk), allocatable :: asymmetry_factor(:)
376 real(kind=musica_dk), allocatable :: ssa(:) ! single scattering albedo
377 reaL(kind=musica_dk), allocatable :: layer_aod(:) ! layer aerosol optical depth
378 real(kind=musica_dk) :: normalized_radius
379
380 n_cheby = this%shortwave_lookup_%number_of_chebyshev_coefficients( )
381 n_bands = this%shortwave_lookup_%number_of_wavelength_bands( )
382 allocate( absorption_coefficients( n_cheby, n_bands ) )
383 allocate( extinction_coefficients( n_cheby, n_bands ) )
384 allocate( asymmetry_factor_coefficients( n_cheby, n_bands ) )
385 allocate( size_function( n_cheby ) )
386 allocate( net_refractive_index( n_bands ) )
387 allocate( absorption( n_bands ) )
388 allocate( extinction( n_bands ) )
389 allocate( asymmetry_factor( n_bands ) )
390 allocate( ssa( n_bands ) )
391 allocate( layer_aod( n_bands ) )
392 net_refractive_index = &
393 this%net_shortwave_refractive_index( mode_state, n_bands )
394 ! lookup Chebyshev coefficients for optical properties
395 call this%shortwave_lookup_%get_optics( net_refractive_index, &
396 absorption = absorption_coefficients, &
397 extinction = extinction_coefficients, &
398 asymmetry_factor = asymmetry_factor_coefficients )
399 ! get Chebyshev function for normalized wet number mode diameter
400 normalized_radius = this%shortwave_lookup_%normalize_radius( &
401 this%wet_surface_mode_radius__m( mode_state ) )
402 call chebyshev_function( normalized_radius, size_function )
403 ! calculate optical properties
404 extinction = this%specific_extinction__m2_kg( mode_state, n_bands, &
405 n_cheby, &
406 extinction_coefficients, &
407 size_function, &
408 this%shortwave_lookup_ )
409 absorption = this%specific_absorption__m2_kg( mode_state, n_bands, &
410 n_cheby, &
411 absorption_coefficients, &
412 size_function, &
413 max_absorption = extinction )
414 absorption(:) = min( absorption(:), extinction(:) )
415 asymmetry_factor = this%asymmetry_factor( mode_state, n_bands, n_cheby, &
416 asymmetry_factor_coefficients, &
417 size_function )
418 !> \todo Are single-scatter albedo and layer_aod the correct terms to use?
419 ssa(:) = 1.0_musica_dk - absorption(:) &
420 / ( max( extinction(:), 1.0e-40_musica_dk ) )
421 layer_aod = extinction(:) * environmental_state%layer_thickness__Pa( ) &
422 / kaccellerationbygravity
423 do i_prop = 1, size( optics )
424 call add_shortwave_property( layer_aod, ssa, asymmetry_factor, &
425 optics( i_prop )%ptr_ )
426 end do
427
428 end subroutine add_shortwave_optics
429
430!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
431! to the set of optical properties
432 !> Calculates longwave optical properties for a given state and adds them
433 !!
434 subroutine add_longwave_optics( this, environmental_state, mode_state, &
435 optics )
436
437 use ai_environmental_state, only : environmental_state_t
438 use ai_optics, only : optics_ptr
439 use mam_constants, only : kaccellerationbygravity
440 use mam_optics_util, only : add_longwave_property
441 use musica_math, only : chebyshev_function
442
443 class(mode_t), intent(in) :: this
444 class(environmental_state_t), intent(in) :: environmental_state
445 class(mode_state_t), intent(in) :: mode_state
446 class(optics_ptr), intent(inout) :: optics(:)
447
448 integer :: i_prop, n_cheby, n_bands
449 real(kind=musica_dk), allocatable :: absorption_coefficients(:,:)
450 real(kind=musica_dk), allocatable :: chebyshev_coefficients(:)
451 real(kind=musica_dk), allocatable :: size_function(:)
452 complex(kind=musica_dk), allocatable :: net_refractive_index(:)
453 real(kind=musica_dk), allocatable :: absorption(:)
454 real(kind=musica_dk), allocatable :: layer_aod(:)
455 real(kind=musica_dk) :: normalized_radius
456
457 n_cheby = this%longwave_lookup_%number_of_chebyshev_coefficients( )
458 n_bands = this%longwave_lookup_%number_of_wavelength_bands( )
459 allocate( absorption_coefficients( n_cheby, n_bands ) )
460 allocate( chebyshev_coefficients( n_cheby ) )
461 allocate( size_function( n_cheby ) )
462 allocate( net_refractive_index( n_bands ) )
463 allocate( absorption( n_bands ) )
464 allocate( layer_aod( n_bands ) )
465 net_refractive_index = &
466 this%net_longwave_refractive_index( mode_state, n_bands )
467 ! lookup Chebyshev coefficients for optical properties
468 call this%longwave_lookup_%get_optics( net_refractive_index, &
469 absorption = absorption_coefficients )
470 ! get Chebyshev function for normalized wet number mode diameter
471 normalized_radius = this%longwave_lookup_%normalize_radius( &
472 this%wet_surface_mode_radius__m( mode_state ) )
473 call chebyshev_function( normalized_radius, size_function )
474 ! calculate optical properties
475 absorption = this%specific_absorption__m2_kg( mode_state, n_bands, &
476 n_cheby, &
477 absorption_coefficients, &
478 size_function )
479 layer_aod(:) = absorption(:) * environmental_state%layer_thickness__Pa( ) &
480 / kaccellerationbygravity
481 do i_prop = 1, size( optics )
482 call add_longwave_property( layer_aod, optics( i_prop )%ptr_ )
483 end do
484
485 end subroutine add_longwave_optics
486
487!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
488
489 !> Returns the geometric mean diameter [m] of the number distribution of the
490 !! mode
491 elemental function geometric_mean_diameter_of_number_distribution__m( this, &
492 mode_state ) result( gmd )
493
494 real(kind=musica_dk) :: gmd
495 class(mode_t), intent(in) :: this
496 class(mode_state_t), intent(in) :: mode_state
497
498 gmd = this%geometric_mean_diameter__m_
499
500 end function geometric_mean_diameter_of_number_distribution__m
501
502!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
503
504 !> Returns the geometric standard deviation of the number distribution of
505 !! the mode
506 elemental function geometric_standard_deviation_of_number_distribution( &
507 this, mode_state ) result( gsd )
508
509 real(kind=musica_dk) :: gsd
510 class(mode_t), intent(in) :: this
511 class(mode_state_t), intent(in) :: mode_state
512
513 gsd = this%geometric_standard_deviation_
514
515 end function geometric_standard_deviation_of_number_distribution
516
517!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
518
519 !> Returns the wet volume to mass mixing ratio for the mode [m3 kg_dryair-1]
520 elemental function wet_volume_to_mass_mixing_ratio__m3_kg( this, &
521 mode_state ) result( vmmr )
522
523 real(kind=musica_dk) :: vmmr
524 class(mode_t), intent(in) :: this
525 class(mode_state_t), intent(in) :: mode_state
526
527 integer :: i_species
528
529 vmmr = 0.0
530 do i_species = 1, size( this%species_ )
531 associate( species => this%species_( i_species ), &
532 species_mmr => &
533 mode_state%mass_mixing_ratio__kg_kg_( i_species ) )
534 vmmr = vmmr + species%volume__m3( species_mmr )
535 end associate
536 end do
537
538 end function wet_volume_to_mass_mixing_ratio__m3_kg
539
540!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
541
542 !> Returns the mode radius [m] of the number distribution for the mode
543 real(kind=musica_dk) elemental function wet_number_mode_radius__m( this, &
544 mode_state ) result( radius )
545
546 class(mode_t), intent(in) :: this
547 class(mode_state_t), intent(in) :: mode_state
548
549 radius = mode_state%wet_number_mode_diameter__m_ * 0.5_musica_dk
550
551 end function wet_number_mode_radius__m
552
553!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
554
555 !> Returns the mode diameter [m] of the number distribution for the mode
556 real(kind=musica_dk) elemental function wet_number_mode_diameter__m( this, &
557 mode_state ) result( diameter )
558
559 class(mode_t), intent(in) :: this
560 class(mode_state_t), intent(in) :: mode_state
561
562 diameter = mode_state%wet_number_mode_diameter__m_
563
564 end function wet_number_mode_diameter__m
565
566!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
567
568 !> Returns the mode diameter [m] of the surface area distribution for the
569 !! mode
570 real(kind=musica_dk) elemental function wet_surface_mode_diameter__m( this, &
571 mode_state ) result( radius )
572
573 class(mode_t), intent(in) :: this
574 class(mode_state_t), intent(in) :: mode_state
575
576 radius = 2.0_musica_dk * this%wet_surface_mode_radius__m( mode_state )
577
578 end function wet_surface_mode_diameter__m
579
580!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581
582 !> Returns the mode radius [m] of the surface area distribution for the mode
583 real(kind=musica_dk) elemental function wet_surface_mode_radius__m( this, &
584 mode_state ) result( radius )
585
586 class(mode_t), intent(in) :: this
587 class(mode_state_t), intent(in) :: mode_state
588
589 associate( sigma => &
590 this%geometric_standard_deviation_of_number_distribution( &
591 mode_state ) )
592 radius = this%wet_number_mode_radius__m( mode_state ) * &
593 exp( 2.0_musica_dk * ( log( sigma )**2 ) )
594 end associate
595
596 end function wet_surface_mode_radius__m
597
598!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
599
600 !> Returns the particle number mixing ratio per mole air [# mol-1]
601 !! \todo is this per mole dry or wet air?
602 real(kind=musica_dk) elemental function number_mixing_ratio__num_mol( this, &
603 mode_state ) result( mixing_ratio )
604
605 class(mode_t), intent(in) :: this
606 class(mode_state_t), intent(in) :: mode_state
607
608 mixing_ratio = mode_state%number_mixing_ratio__num_mol_
609
610 end function number_mixing_ratio__num_mol
611
612!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
613
614 !> Calculates the net refractive index for the mode based on weighted
615 !! contributions from constituent species in the shortwave region
616 pure function net_shortwave_refractive_index( this, mode_state, &
617 number_of_bands ) result( refractive_index )
618
619 complex(kind=musica_dk) :: refractive_index( number_of_bands )
620 class(mode_t), intent(in) :: this
621 class(mode_state_t), intent(in) :: mode_state
622 integer, intent(in) :: number_of_bands
623
624 integer :: i_species, i_band
625 real(kind=musica_dk) :: total_vmmr, species_vmmr
626
627 refractive_index(:) = ( 0.0_musica_dk, 0.0_musica_dk )
628 total_vmmr = 0.0_musica_dk
629 do i_species = 1, size( this%species_ )
630 associate( species => this%species_( i_species ), &
631 species_mmr => &
632 mode_state%mass_mixing_ratio__kg_kg_( i_species ) )
633 species_vmmr = species%volume__m3( species_mmr )
634 do i_band = 1, number_of_bands
635 refractive_index( i_band ) = refractive_index( i_band ) &
636 + species%shortwave_refractive_index( i_band ) * species_vmmr
637 end do
638 total_vmmr = total_vmmr + species_vmmr
639 end associate
640 end do
641 refractive_index(:) = refractive_index(:) &
642 / max( total_vmmr, 1.0e-60_musica_dk )
643
644 end function net_shortwave_refractive_index
645
646!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
647
648 !> Calculates the net refractive index for the mode based on weighted
649 !! contributions from constituent species in the longwave region
650 pure function net_longwave_refractive_index( this, mode_state, &
651 number_of_bands ) result( refractive_index )
652
653 complex(kind=musica_dk) :: refractive_index( number_of_bands )
654 class(mode_t), intent(in) :: this
655 class(mode_state_t), intent(in) :: mode_state
656 integer, intent(in) :: number_of_bands
657
658 integer :: i_species, i_band
659 real(kind=musica_dk) :: total_vmmr, species_vmmr
660
661 refractive_index(:) = ( 0.0_musica_dk, 0.0_musica_dk )
662 total_vmmr = 0.0_musica_dk
663 do i_species = 1, size( this%species_ )
664 associate( species => this%species_( i_species ), &
665 species_mmr => &
666 mode_state%mass_mixing_ratio__kg_kg_( i_species ) )
667 species_vmmr = species%volume__m3( species_mmr )
668 do i_band = 1, number_of_bands
669 refractive_index( i_band ) = refractive_index( i_band ) &
670 + species%longwave_refractive_index( i_band ) * species_vmmr
671 end do
672 total_vmmr = total_vmmr + species_vmmr
673 end associate
674 end do
675 refractive_index(:) = refractive_index(:) &
676 / max( total_vmmr, 1.0e-60_musica_dk )
677
678 end function net_longwave_refractive_index
679
680!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
681
682 !> Calculates the specific absorption [m2 kg-1] for a given Chebyshev
683 !! function
684 !!
685 !! \todo is "specific absorption" (m2 kg-1) the correct name/units for what
686 !! is returned from this function?
687 pure function specific_absorption__m2_kg( this, mode_state, number_of_bands,&
688 number_of_coefficients, coefficients, size_function, max_absorption )
689
690 use mam_constants, only : kwaterdensitystp
691 use musica_math, only : weighted_chebyshev
692
693 real(kind=musica_dk) :: specific_absorption__m2_kg( number_of_bands )
694 class(mode_t), intent(in) :: this
695 class(mode_state_t), intent(in) :: mode_state
696 integer, intent(in) :: number_of_bands
697 integer, intent(in) :: number_of_coefficients
698 !> Chebyshev coefficients for absorption calculation
699 real(kind=musica_dk), intent(in) :: coefficients( number_of_coefficients, &
700 number_of_bands )
701 !> Chebyshev function for the current surface mode radius
702 real(kind=musica_dk), intent(in) :: size_function( number_of_coefficients )
703 !> Maximum values for specific absorption for each band [m2 kg-1]
704 real(kind=musica_dk), intent(in), optional :: max_absorption( &
705 number_of_bands )
706
707 integer :: i_band
708 do i_band = 1, number_of_bands
709 associate( absorp => specific_absorption__m2_kg( i_band ) )
710 absorp = weighted_chebyshev( number_of_coefficients, &
711 coefficients( :, i_band ), &
712 size_function )
713 absorp = absorp * &
714 this%wet_volume_to_mass_mixing_ratio__m3_kg( mode_state ) * &
715 kwaterdensitystp
716 absorp = max( 0.0_musica_dk, absorp )
717 if( present( max_absorption ) ) then
718 absorp = min( max_absorption( i_band ), absorp )
719 end if
720 end associate
721 end do
722
723 end function specific_absorption__m2_kg
724
725!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
726
727 !> Calculates the specific extinction [m2 kg-1] for a given Chebyshev
728 !! function
729 !!
730 !! \todo is "specific extinction" (m2 kg-1) the correct name/units for what
731 !! is returned from this function?
732 pure function specific_extinction__m2_kg( this, mode_state, number_of_bands,&
733 number_of_coefficients, coefficients, size_function, optics_lookup )
734
735 use mam_constants, only : kwaterdensitystp
736 use musica_math, only : weighted_chebyshev
737
738 !> Specific extinction by wavelength
739 real(kind=musica_dk) :: specific_extinction__m2_kg( number_of_bands )
740 !> MAM mode
741 class(mode_t), intent(in) :: this
742 !> MAM mode state
743 class(mode_state_t), intent(in) :: mode_state
744 !> Number of wavelength bands
745 integer, intent(in) :: number_of_bands
746 !> Number of Chebyshev coefficients
747 integer, intent(in) :: number_of_coefficients
748 !> Chebyshev coefficients for extinction calculation
749 real(kind=musica_dk), intent(in) :: coefficients( number_of_coefficients, &
750 number_of_bands )
751 !> Chebyshev function for the current surface mode radius
752 real(kind=musica_dk), intent(in) :: size_function( number_of_coefficients )
753 !> Optics lookup table
754 class(optics_lookup_t), intent(in) :: optics_lookup
755
756 integer :: i_band
757 real(kind=musica_dk) :: surf_rad
758
759 surf_rad = this%wet_surface_mode_radius__m( mode_state )
760 if( surf_rad .le. optics_lookup%maximum_radius__m( ) ) then
761 do i_band = 1, number_of_bands
762 associate( ext => specific_extinction__m2_kg( i_band ) )
763 ext = weighted_chebyshev( number_of_coefficients, &
764 coefficients( :, i_band ), &
765 size_function )
766 !> \todo is the data returned from the lookup table in ln( m2/kg )?
767 ext = exp( ext )
768 end associate
769 end do
770 else
771 specific_extinction__m2_kg(:) = 1.5_musica_dk &
772 / ( surf_rad * kwaterdensitystp )
773 end if
774 !> \todo in the original code, this is labelled as converting "m2/kg
775 !! water to m2/kg aerosol", but it appears to convert instead to
776 !! m2/kg dry air - is this right???
777 specific_extinction__m2_kg(:) = specific_extinction__m2_kg(:) &
778 * this%wet_volume_to_mass_mixing_ratio__m3_kg( mode_state ) &
779 * kwaterdensitystp
780
781 end function specific_extinction__m2_kg
782
783!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
784
785 !> Calculates the asymmetry parameter [unitless ratio] for a given Chebyshev
786 !! function
787 pure function asymmetry_factor( this, mode_state, number_of_bands, &
788 number_of_coefficients, coefficients, size_function )
789
790 use mam_constants, only : kwaterdensitystp
791 use musica_math, only : weighted_chebyshev
792
793 !> Asymmetry factor by wavelength
794 real(kind=musica_dk) :: asymmetry_factor( number_of_bands )
795 !> MAM mode
796 class(mode_t), intent(in) :: this
797 !> MAM mode state
798 class(mode_state_t), intent(in) :: mode_state
799 !> Number of wavelength bands
800 integer, intent(in) :: number_of_bands
801 !> Number of Chebyshev coefficients
802 integer, intent(in) :: number_of_coefficients
803 !> Chebyshev coefficients for asymmetry parameter calculation
804 real(kind=musica_dk), intent(in) :: coefficients( number_of_coefficients, &
805 number_of_bands )
806 !> Chebyshev function for the current surface mode radius
807 real(kind=musica_dk), intent(in) :: size_function( number_of_coefficients )
808
809 integer :: i_band
810
811 do i_band = 1, number_of_bands
812 associate( asym => asymmetry_factor( i_band ) )
813 asym = weighted_chebyshev( number_of_coefficients, &
814 coefficients( :, i_band ), &
815 size_function )
816 end associate
817 end do
818
819 end function asymmetry_factor
820
821!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
822!!
823!! mode_state_t functions
824!!
825!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
826
827 !> Returns the number of doubles needed to hold the raw mode state
828 integer function raw_size( this )
829
830 class(mode_state_t), intent(in) :: this
831
832 raw_size = 2 + size( this%mass_mixing_ratio__kg_kg_ )
833
834 end function raw_size
835
836!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
837
838 !> Loads raw mode state data to the mode_state_t object
839 subroutine load_state( this, raw_state, index )
840
841 class(mode_state_t), intent(inout) :: this
842 real(kind=musica_dk), intent(in) :: raw_state(:)
843 integer, optional, intent(inout) :: index
844
845 integer :: id, last_id
846
847 id = 1
848 if( present( index ) ) id = index
849 this%wet_number_mode_diameter__m_ = raw_state( id )
850 this%number_mixing_ratio__num_mol_ = raw_state( id + 1 )
851 id = id + 2
852 last_id = id + size( this%mass_mixing_ratio__kg_kg_ ) - 1
853 this%mass_mixing_ratio__kg_kg_(:) = raw_state( id : last_id )
854 if( present( index ) ) index = last_id + 1
855
856 end subroutine load_state
857
858!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859
860 !> Dumps the aerosol state into the raw state array
861 subroutine dump_state( this, raw_state, index )
862
863 class(mode_state_t), intent(in) :: this
864 real(kind=musica_dk), intent(inout) :: raw_state(:)
865 integer, optional, intent(inout) :: index
866
867 integer :: id, last_id
868
869 id = 1
870 if( present( index ) ) id = index
871 raw_state( id ) = this%wet_number_mode_diameter__m_
872 raw_state( id + 1 ) = this%number_mixing_ratio__num_mol_
873 id = id + 2
874 last_id = id + size( this%mass_mixing_ratio__kg_kg_ ) - 1
875 raw_state( id : last_id ) = this%mass_mixing_ratio__kg_kg_(:)
876 if( present( index ) ) index = last_id + 1
877
878 end subroutine dump_state
879
880!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
881
882 !> Set the mode state to a random, but reasonable, state. For testing only.
883 subroutine randomize( this )
884
885 class(mode_state_t), intent(inout) :: this
886
887 real(kind=musica_dk) :: rand_val
888 integer :: i_species
889
890 call random_number( rand_val )
891 this%wet_number_mode_diameter__m_ = 10.0**( rand_val * 3 - 5 )
892 this%number_mixing_ratio__num_mol_ = 10.0**( rand_val * 6 )
893 do i_species = 1, size( this%mass_mixing_ratio__kg_kg_ )
894 call random_number( rand_val )
895 this%mass_mixing_ratio__kg_kg_( i_species ) = rand_val * 2.0e-9
896 end do
897
898 end subroutine randomize
899
900!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
901
902end module mam_mode
mam_constants
Constants used in MAM calculations.
Definition: constants.F90:8
mam_constants::kaccellerationbygravity
real(kind=musica_dk), parameter kaccellerationbygravity
AccellerationByGravity [m s-2].
Definition: constants.F90:17
mam_constants::kwaterdensitystp
real(kind=musica_dk), parameter kwaterdensitystp
Water density at STP [kg m-3].
Definition: constants.F90:15
mam_mode
The mode_t type and related functions.
Definition: mode.F90:8
mam_mode::longwave_grid
type(wavelength_grid_t) function longwave_grid(this)
Returns the native longwave optics wavelength grid.
Definition: mode.F90:290
mam_mode::net_longwave_refractive_index
pure complex(kind=musica_dk) function, dimension(number_of_bands) net_longwave_refractive_index(this, mode_state, number_of_bands)
Calculates the net refractive index for the mode based on weighted contributions from constituent spe...
Definition: mode.F90:652
mam_mode::wet_number_mode_radius__m
real(kind=musica_dk) elemental function wet_number_mode_radius__m(this, mode_state)
Returns the mode radius [m] of the number distribution for the mode.
Definition: mode.F90:545
mam_mode::print_state
subroutine print_state(this, aerosol_state, io_unit)
Outputs the current mode state.
Definition: mode.F90:304
mam_mode::wet_surface_mode_diameter__m
real(kind=musica_dk) elemental function wet_surface_mode_diameter__m(this, mode_state)
Returns the mode diameter [m] of the surface area distribution for the mode.
Definition: mode.F90:572
mam_mode::raw_size
integer function raw_size(this)
Returns the number of doubles needed to hold the raw mode state.
Definition: mode.F90:829
mam_mode::constructor
type(mode_t) function constructor(config)
Constructor of mode_t objects.
Definition: mode.F90:88
mam_mode::wet_number_mode_diameter__m
real(kind=musica_dk) elemental function wet_number_mode_diameter__m(this, mode_state)
Returns the mode diameter [m] of the number distribution for the mode.
Definition: mode.F90:558
mam_mode::geometric_standard_deviation_of_number_distribution
elemental real(kind=musica_dk) function geometric_standard_deviation_of_number_distribution(this, mode_state)
Returns the geometric standard deviation of the number distribution of the mode.
Definition: mode.F90:508
mam_mode::specific_absorption__m2_kg
pure real(kind=musica_dk) function, dimension(number_of_bands) specific_absorption__m2_kg(this, mode_state, number_of_bands, number_of_coefficients, coefficients, size_function, max_absorption)
Calculates the specific absorption [m2 kg-1] for a given Chebyshev function.
Definition: mode.F90:689
mam_mode::randomize
subroutine randomize(this)
Set the mode state to a random, but reasonable, state. For testing only.
Definition: mode.F90:884
mam_mode::wet_surface_mode_radius__m
real(kind=musica_dk) elemental function wet_surface_mode_radius__m(this, mode_state)
Returns the mode radius [m] of the surface area distribution for the mode.
Definition: mode.F90:585
mam_mode::dump_state
subroutine dump_state(this, raw_state, index)
Dumps the aerosol state into the raw state array.
Definition: mode.F90:862
mam_mode::number_mixing_ratio__num_mol
real(kind=musica_dk) elemental function number_mixing_ratio__num_mol(this, mode_state)
Returns the particle number mixing ratio per mole air [# mol-1].
Definition: mode.F90:604
mam_mode::specific_extinction__m2_kg
pure real(kind=musica_dk) function, dimension(number_of_bands) specific_extinction__m2_kg(this, mode_state, number_of_bands, number_of_coefficients, coefficients, size_function, optics_lookup)
Calculates the specific extinction [m2 kg-1] for a given Chebyshev function.
Definition: mode.F90:734
mam_mode::shortwave_optics_array
subroutine shortwave_optics_array(this, environmental_state, aerosol_state, optics)
Returns shortwave optical properties.
Definition: mode.F90:192
mam_mode::new_state
class(aerosol_state_t) function, pointer new_state(this)
Creates a new state object for the mode.
Definition: mode.F90:126
mam_mode::net_shortwave_refractive_index
pure complex(kind=musica_dk) function, dimension(number_of_bands) net_shortwave_refractive_index(this, mode_state, number_of_bands)
Calculates the net refractive index for the mode based on weighted contributions from constituent spe...
Definition: mode.F90:618
mam_mode::shortwave_optics_scalar
subroutine shortwave_optics_scalar(this, environmental_state, aerosol_state, optics)
Returns shortwave optical properties.
Definition: mode.F90:169
mam_mode::add_shortwave_optics
subroutine add_shortwave_optics(this, environmental_state, mode_state, optics)
Calculates shortwave optical properties for a given state and adds them to the set of optical propert...
Definition: mode.F90:355
mam_mode::shortwave_grid
type(wavelength_grid_t) function shortwave_grid(this)
Returns the native shortwave optics wavelength grid.
Definition: mode.F90:276
mam_mode::longwave_optics_scalar
subroutine longwave_optics_scalar(this, environmental_state, aerosol_state, optics)
Returns longwave optical properties.
Definition: mode.F90:223
mam_mode::new_optics
class(optics_t) function, pointer new_optics(this, property, output_grid, interpolation_strategy)
Creates an optics_t object for a given property.
Definition: mode.F90:144
mam_mode::add_longwave_optics
subroutine add_longwave_optics(this, environmental_state, mode_state, optics)
Calculates longwave optical properties for a given state and adds them to the set of optical properti...
Definition: mode.F90:436
mam_mode::load_state
subroutine load_state(this, raw_state, index)
Loads raw mode state data to the mode_state_t object.
Definition: mode.F90:840
mam_mode::geometric_mean_diameter_of_number_distribution__m
elemental real(kind=musica_dk) function geometric_mean_diameter_of_number_distribution__m(this, mode_state)
Returns the geometric mean diameter [m] of the number distribution of the mode.
Definition: mode.F90:493
mam_mode::wet_volume_to_mass_mixing_ratio__m3_kg
elemental real(kind=musica_dk) function wet_volume_to_mass_mixing_ratio__m3_kg(this, mode_state)
Returns the wet volume to mass mixing ratio for the mode [m3 kg_dryair-1].
Definition: mode.F90:522
mam_mode::longwave_optics_array
subroutine longwave_optics_array(this, environmental_state, aerosol_state, optics)
Returns longwave optical properties.
Definition: mode.F90:246
mam_optics_lookup
The optics_lookup_t type and related functions.
Definition: optics_lookup.F90:8
mam_optics_util
Constants used to calculate MAM optical properties.
Definition: optics_util.F90:8
mam_optics_util::create_optics
class(optics_t) function, pointer, public create_optics(property, native_shortwave_grid, native_longwave_grid, output_grid, interpolation_strategy)
Returns an optics_t object for a given property.
Definition: optics_util.F90:35
mam_optics_util::add_shortwave_property
subroutine, public add_shortwave_property(extinction_optical_depth, single_scatter_albedo, asymmetry_factor, optics)
Adds to shortwave property values.
Definition: optics_util.F90:165
mam_optics_util::add_longwave_property
subroutine, public add_longwave_property(extinction_optical_depth, optics)
Adds to longwave property values.
Definition: optics_util.F90:213
mam_species
The species_t type and related functions.
Definition: species.F90:8
asymmetry_factor
real(kind=musica_dk), dimension(:,:), intent(out), optional asymmetry_factor
Definition: optics_lookup.F90:173
absorption
real(kind=musica_dk), dimension(:,:), intent(out), optional absorption
Definition: optics_lookup.F90:171
i_band
integer i_band
Definition: optics_lookup.F90:176
extinction
real(kind=musica_dk), dimension(:,:), intent(out), optional extinction
Definition: optics_lookup.F90:172
mam_mode::mode_state_t
Aerosol mode state.
Definition: mode.F90:66
mam_mode::mode_t
An aerosol mode.
Definition: mode.F90:22
mam_optics_lookup::optics_lookup_t
Definition: optics_lookup.F90:47
mam_species::species_t
Definition: species.F90:20
Generated on Mon Aug 15 2022 15:01:05 for mam by doxygen 1.9.4