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!> Tests for the mam_mode module
6
7!> Test program for the mode_t type and related functions
8program test_mode
9
10 implicit none
11
12 call test_mode_t( )
13
14contains
15
16!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
17
18 subroutine test_mode_t( )
19
20 use ai_aerosol_state, only : aerosol_state_t
21 use ai_environmental_state, only : environmental_state_t
22 use ai_optics, only : optics_t, optics_ptr
23 use ai_wavelength_grid, only : wavelength_grid_t
24 use mam_constants, only : kaccellerationbygravity
25 use mam_mode, only : mode_t, mode_state_t
26 use mam_optics_lookup, only : optics_lookup_t
27 use musica_assert, only : assert, die, almost_equal
28 use musica_config, only : config_t
29 use musica_constants, only : dk => musica_dk
30 use musica_property, only : property_t
31
32 character(len=*), parameter :: my_name = "mode_t tests"
33 type(mode_t) :: mode
34 class(aerosol_state_t), pointer :: state_a, state_b
35 type(config_t) :: mode_config, optics_config
36 class(optics_t), pointer :: optics
37 type(optics_ptr), allocatable :: optics_set(:)
38 type(environmental_state_t) :: env_state
39 type(wavelength_grid_t) :: grid
40 type(property_t) :: prop
41 type(optics_lookup_t) :: mock_optics_lookup
42 real(kind=dk), allocatable :: raw_state_a(:), raw_state_b(:)
43 real(kind=dk), allocatable :: optics_values(:)
44 real(kind=dk) :: temp_value
45 complex(kind=dk), allocatable :: complex_array_a(:), complex_array_b(:)
46 real(kind=dk), allocatable :: absorption_short(:), absorption_long(:)
47 real(kind=dk), allocatable :: extinction_short(:), extinction_long(:)
48 real(kind=dk), allocatable :: asymmetry_short(:), asymmetry_long(:)
49 real(kind=dk), allocatable :: ssa_short(:), aod_short(:), aod_long(:)
50 real(kind=dk), allocatable :: compare_values(:)
51 integer :: i_elem, i_band
52 logical :: flag
53
54 call mode_config%from_file( "mode_config.json" )
55 mode = mode_t( mode_config )
56
57 state_a => mode%new_state( )
58 state_b => mode%new_state( )
59
60 ! test mode_state_t functions
61 call state_a%randomize( )
62 call state_b%randomize( )
63 call assert( 119422573, state_a%raw_size( ) .eq. state_b%raw_size( ) )
64 allocate( raw_state_a( state_a%raw_size( ) ) )
65 allocate( raw_state_b( state_b%raw_size( ) ) )
66 call state_a%dump_state( raw_state_a )
67 call state_b%dump_state( raw_state_b )
68 flag = .true.
69 do i_elem = 1, size( raw_state_a )
70 flag = flag .and. raw_state_a( i_elem ) .eq. raw_state_b( i_elem )
71 end do
72 call assert( 904197749, .not. flag )
73 call state_b%load_state( raw_state_a )
74 call state_b%dump_state( raw_state_b )
75 flag = .true.
76 do i_elem = 1, size( raw_state_a )
77 flag = flag .and. raw_state_a( i_elem ) .eq. raw_state_b( i_elem )
78 end do
79 call assert( 167637108, flag )
80
81 deallocate( raw_state_a )
82 deallocate( raw_state_b )
83 allocate( raw_state_a( state_a%raw_size( ) + 15 ) )
84 allocate( raw_state_b( state_b%raw_size( ) + 25 ) )
85 raw_state_a(:) = 12.5_dk
86 raw_state_b(:) = 43.1_dk
87 call state_a%randomize( )
88 call state_b%randomize( )
89 i_elem = 11
90 call state_a%dump_state( raw_state_a, i_elem )
91 call assert( 800815204, i_elem .eq. 11 + state_a%raw_size( ) )
92 i_elem = 16
93 call state_b%dump_state( raw_state_b, i_elem )
94 call assert( 850030683, i_elem .eq. 16 + state_b%raw_size( ) )
95 flag = .true.
96 do i_elem = 1, state_a%raw_size( )
97 flag = flag .and. raw_state_a( 10 + i_elem ) .eq. raw_state_b( 15 + i_elem )
98 end do
99 call assert( 173214135, .not. flag )
100 i_elem = 11
101 call state_b%load_state( raw_state_a, i_elem )
102 call assert( 161232282, i_elem .eq. 11 + state_b%raw_size( ) )
103 i_elem = 16
104 call state_a%dump_state( raw_state_b, i_elem )
105 call assert( 157873509, i_elem .eq. 16 + state_b%raw_size( ) )
106 flag = .true.
107 do i_elem = 1, state_a%raw_size( )
108 flag = flag .and. raw_state_a( 10 + i_elem ) .eq. raw_state_b( 15 + i_elem )
109 end do
110 call assert( 272097388, flag )
111 do i_elem = 1, 10
112 call assert( 661626675, raw_state_a( i_elem ) .eq. 12.5_dk )
113 end do
114 do i_elem = 11 + state_a%raw_size( ), 15 + state_a%raw_size( )
115 call assert( 642927276, raw_state_a( i_elem ) .eq. 12.5_dk )
116 end do
117 do i_elem = 1, 15
118 call assert( 472770372, raw_state_b( i_elem ) .eq. 43.1_dk )
119 end do
120 do i_elem = 16 + state_a%raw_size( ), 25 + state_a%raw_size( )
121 call assert( 867563966, raw_state_b( i_elem ) .eq. 43.1_dk )
122 end do
123
124 deallocate( raw_state_a )
125 deallocate( raw_state_b )
126 deallocate( state_a )
127 deallocate( state_b )
128
129 ! grid accessors
130 grid = wavelength_grid_t( (/ 12.3_dk, 100.0_dk /), &
131 (/ 92.3_dk, 1145.0_dk /) )
132 call assert( 588709033, grid .eq. mode%shortwave_grid( ) )
133 grid = wavelength_grid_t( (/ 12.3_dk, 200.0_dk /), &
134 (/ 92.3_dk, 2290.0_dk /) )
135 call assert( 179027108, grid .eq. mode%longwave_grid( ) )
136
137 ! simple accessors/calculators
138 state_a => mode%new_state( )
139 call state_a%randomize( )
140 allocate( raw_state_a( state_a%raw_size( ) ) )
141 call state_a%dump_state( raw_state_a )
142 ! IMPORTANT: this is just for testing - never rely on the raw state data
143 ! outside of tests. It will change often
144 associate( wet_diameter => raw_state_a( 1 ), &
145 number_mr => raw_state_a( 2 ), &
146 foo_mass => raw_state_a( 3 ), &
147 bar_mass => raw_state_a( 4 ) )
148 select type( state_a )
149 type is( mode_state_t )
150 call assert( 680874740, &
151 mode%geometric_mean_diameter_of_number_distribution__m( state_a ) &
152 .eq. 42.53_dk )
153 call assert( 961896750, &
154 mode%geometric_standard_deviation_of_number_distribution( state_a ) &
155 .eq. 2.4_dk )
156 temp_value = foo_mass * 6.0_dk + bar_mass * 3.0_dk
157 call assert( 267834042, &
158 mode%wet_volume_to_mass_mixing_ratio__m3_kg( state_a ) &
159 .eq. temp_value )
160 call assert( 315596282, mode%wet_number_mode_radius__m( state_a ) &
161 .eq. wet_diameter * 0.5_dk )
162 call assert( 371981602, mode%wet_number_mode_diameter__m( state_a ) &
163 .eq. wet_diameter )
164 temp_value = wet_diameter * exp( 2.0_dk * ( log( 2.4_dk ) )**2 )
165 call assert( 102514747, mode%wet_surface_mode_diameter__m( state_a ) &
166 .eq. temp_value )
167 temp_value = temp_value * 0.5_dk
168 call assert( 543165047, mode%wet_surface_mode_radius__m( state_a ) &
169 .eq. temp_value )
170 call assert( 769707271, mode%number_mixing_ratio__num_mol( state_a ) &
171 .eq. number_mr )
172
173 ! net refractive indices
174 allocate( complex_array_a( 2 ) )
175 allocate( complex_array_b( 2 ) )
176 complex_array_a = mode%net_shortwave_refractive_index( state_a, 2 )
177 complex_array_b = net_refractive_index( state_a, 2 )
178 call assert( 847572900, complex_array_a(1) .eq. complex_array_b(1) )
179 call assert( 934499897, complex_array_a(2) .eq. complex_array_b(2) )
180 complex_array_a = mode%net_longwave_refractive_index( state_a, 2 )
181 complex_array_b = net_refractive_index( state_a, 1 )
182 call assert( 764342993, complex_array_a(1) .eq. complex_array_b(1) )
183 call assert( 594186089, complex_array_a(2) .eq. complex_array_b(2) )
184 deallocate( complex_array_a )
185 deallocate( complex_array_b )
186
187 ! optics functions
188 call env_state%randomize( )
189 allocate( absorption_long( 2 ) )
190 allocate( absorption_short( 2 ) )
191 allocate( extinction_long( 2 ) )
192 allocate( extinction_short( 2 ) )
193 allocate( asymmetry_long( 2 ) )
194 allocate( asymmetry_short( 2 ) )
195 allocate( ssa_short( 2 ) )
196 allocate( aod_short( 2 ) )
197 allocate( aod_long( 2 ) )
198 allocate( compare_values( 2 ) )
199 absorption_long = specific_absorption( mode, state_a, 1 )
200 absorption_short = specific_absorption( mode, state_a, 2 )
201 extinction_long = specific_extinction( mode, state_a, 1 )
202 extinction_short = specific_extinction( mode, state_a, 2 )
203 asymmetry_long = asymmetry_factor( mode, state_a, 1 )
204 asymmetry_short = asymmetry_factor( mode, state_a, 2 )
205
206 absorption_short(:) = min( absorption_short(:), extinction_short(:) )
207 ssa_short(:) = 1.0_dk - absorption_short(:) / &
208 ( max( extinction_short(:), 1.0e-40_dk ) )
209 aod_short(:) = extinction_short(:) * env_state%layer_thickness__Pa( ) / &
210 kaccellerationbygravity
211 aod_long(:) = absorption_long(:) * env_state%layer_thickness__Pa( ) / &
212 kaccellerationbygravity
213
214 ! get optics shortwave
215 call optics_config%empty( )
216 call optics_config%add( "type", 1, my_name )
217 mock_optics_lookup = optics_lookup_t( optics_config )
218 grid = mode%shortwave_grid( )
219 allocate( optics_values( grid%number_of_sections( ) ) )
220
221 ! scalar get optics
222 optics => mode%new_optics( property_t( my_name, "asymmetry factor" ), &
223 grid )
224 call assert( 771601198, optics%native_grid( ) .eq. &
225 mock_optics_lookup%grid( ) )
226 call assert( 560399600, optics%output_grid( ) .eq. grid )
227 call mode%shortwave_optics( env_state, state_a, optics )
228 call optics%get_values( optics_values )
229 compare_values(:) = aod_short(:) * ssa_short(:) * asymmetry_short(:)
230 call assert( 187282691, size( optics_values ) .eq. size( compare_values ) )
231 do i_band = 1, size( optics_values )
232 call assert( 803354202, optics_values( i_band ) .gt. 0.0_dk )
233 call assert( 175300838, optics_values( i_band ) .lt. 1.0e200_dk )
234 call assert( 229780624, almost_equal( optics_values( i_band ), &
235 compare_values( i_band ) ) )
236 end do
237 deallocate( optics )
238
239 ! array get optics
240 allocate( optics_set( 3 ) )
241 prop = property_t( my_name, "layer extinction optical depth" )
242 optics_set( 1 )%ptr_ => mode%new_optics( prop, grid )
243 prop = property_t( my_name, "layer single-scatter albedo" )
244 optics_set( 2 )%ptr_ => mode%new_optics( prop, grid )
245 prop = property_t( my_name, "forward scattered fraction" )
246 optics_set( 3 )%ptr_ => mode%new_optics( prop, grid )
247 call mode%shortwave_optics( env_state, state_a, optics_set )
248 call optics_set( 1 )%ptr_%get_values( optics_values )
249 compare_values(:) = aod_short(:)
250 do i_band = 1, size( optics_values )
251 call assert( 986456716, optics_values( i_band ) .gt. 0.0_dk )
252 call assert( 198775062, optics_values( i_band ) .lt. 1.0e200_dk )
253 call assert( 646142908, almost_equal( optics_values( i_band ), &
254 compare_values( i_band ) ) )
255 end do
256 call optics_set( 2 )%ptr_%get_values( optics_values )
257 compare_values(:) = aod_short(:) * ssa_short(:)
258 do i_band = 1, size( optics_values )
259 call assert( 312998941, optics_values( i_band ) .gt. 0.0_dk )
260 call assert( 360308886, optics_values( i_band ) .lt. 1.0e200_dk )
261 call assert( 255160382, almost_equal( optics_values( i_band ), &
262 compare_values( i_band ) ) )
263 end do
264 call optics_set( 3 )%ptr_%get_values( optics_values )
265 compare_values(:) = aod_short(:) * ssa_short(:) * asymmetry_short(:) &
266 * asymmetry_short(:)
267 do i_band = 1, size( optics_values )
268 call assert( 414788672, optics_values( i_band ) .gt. 0.0_dk )
269 call assert( 244631768, optics_values( i_band ) .lt. 1.0e200_dk )
270 call assert( 139483264, almost_equal( optics_values( i_band ), &
271 compare_values( i_band ) ) )
272 end do
273 deallocate( optics_set )
274 deallocate( optics_values )
275
276 ! get optics longwave
277 call optics_config%empty( )
278 call optics_config%add( "type", 2, my_name )
279 mock_optics_lookup = optics_lookup_t( optics_config )
280 grid = mode%longwave_grid( )
281 allocate( optics_values( grid%number_of_sections( ) ) )
282
283 ! scalar get optics
284 prop = property_t( my_name, "layer absorption optical depth" )
285 optics => mode%new_optics( prop, grid )
286 call assert( 691311390, optics%native_grid( ) .eq. &
287 mock_optics_lookup%grid( ) )
288 call assert( 121096585, optics%output_grid( ) .eq. grid )
289 call mode%longwave_optics( env_state, state_a, optics )
290 call optics%get_values( optics_values )
291 compare_values(:) = aod_long(:)
292 call assert( 517795713, size( optics_values ) .eq. size( compare_values ) )
293 do i_band = 1, size( optics_values )
294 call assert( 295064557, optics_values( i_band ) .gt. 0.0_dk )
295 call assert( 742432403, optics_values( i_band ) .lt. 1.0e200_dk )
296 call assert( 572275499, almost_equal( optics_values( i_band ), &
297 compare_values( i_band ) ) )
298 end do
299
300 ! array get optics
301 allocate( optics_set( 1 ) )
302 optics_set( 1 )%ptr_ => optics
303 call mode%longwave_optics( env_state, state_a, optics_set )
304 call optics_set( 1 )%ptr_%get_values( optics_values )
305 compare_values(:) = aod_long(:)
306 call assert( 836051349, size( optics_values ) .eq. size( compare_values ) )
307 do i_band = 1, size( optics_values )
308 call assert( 383419196, optics_values( i_band ) .gt. 0.0_dk )
309 call assert( 213262292, optics_values( i_band ) .lt. 1.0e200_dk )
310 call assert( 108113788, almost_equal( optics_values( i_band ), &
311 compare_values( i_band ) ) )
312 end do
313
314 deallocate( optics )
315 deallocate( optics_set )
316
317
318 class default
319 call die( 889482580 )
320 end select
321 end associate
322 deallocate( state_a )
323
324 end subroutine test_mode_t
325
326!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
327
328 function net_refractive_index( mode_state, long_or_short )
329
330 use mam_mode, only : mode_state_t
331 use musica_assert, only : die
332 use musica_constants, only : dk => musica_dk
333
334 complex(kind=dk) :: net_refractive_index(2)
335 class(mode_state_t), intent(in) :: mode_state
336 integer, intent(in) :: long_or_short ! 1=long; 2=short
337
338 real(kind=dk), allocatable :: raw_state(:)
339 real(kind=dk) :: real_comp, imag_comp
340
341 ! select coefficients for long/short wave refractive index in mock_species
342 select case( long_or_short )
343 case( 1 )
344 real_comp = 15.6_dk
345 imag_comp = 1.32_dk
346 case( 2 )
347 real_comp = 12.5_dk
348 imag_comp = 2.53_dk
349 case default
350 call die( 530770503 )
351 end select
352
353 allocate( raw_state( mode_state%raw_size( ) ) )
354 call mode_state%dump_state( raw_state )
355 ! IMPORTANT: this is just for testing - never rely on the raw state data
356 ! outside of tests. It will change often
357 associate( foo_mmr => raw_state(3), bar_mmr => raw_state(4) )
358 net_refractive_index(1) = &
359 cmplx( real_comp, 6 + imag_comp, kind=dk ) * ( foo_mmr * 6 ) + &
360 cmplx( real_comp, 3 + imag_comp, kind=dk ) * ( bar_mmr * 3 )
361 net_refractive_index(2) = &
362 cmplx( 2 * real_comp, 6 + imag_comp, kind=dk ) * ( foo_mmr * 6 ) + &
363 cmplx( 2 * real_comp, 3 + imag_comp, kind=dk ) * ( bar_mmr * 3 )
364 net_refractive_index(:) = net_refractive_index(:) &
365 / max( foo_mmr * 6 + bar_mmr * 3, 1.0e-60_dk )
366 end associate
367
368 end function net_refractive_index
369
370!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
371
372 function specific_absorption( mode, mode_state, long_or_short )
373
374 use mam_constants, only : kwaterdensitystp
375 use mam_mode, only : mode_t, mode_state_t
376 use mam_optics_lookup, only : optics_lookup_t
377 use musica_assert, only : assert, die, almost_equal
378 use musica_config, only : config_t
379 use musica_constants, only : dk => musica_dk
380 use musica_math, only : chebyshev_function, &
381 weighted_chebyshev
382
383 real(kind=dk) :: specific_absorption(2)
384 class(mode_t), intent(in) :: mode
385 class(mode_state_t), intent(in) :: mode_state
386 integer, intent(in) :: long_or_short ! 1=long, 2=short
387
388 character(len=*), parameter :: my_name = "mode absorption test"
389 type(config_t) :: lookup_config
390 type(optics_lookup_t) :: mock_optics_lookup
391 real(kind=dk) :: cheby_coeff(3,2)
392 real(kind=dk) :: cheby_func(3)
393 real(kind=dk) :: normal_radius, wet_surface_radius
394 real(kind=dk) :: compare_values(2)
395 real(kind=dk), allocatable :: raw_state(:)
396
397 call lookup_config%empty( )
398 select case( long_or_short )
399 case( 1 )
400 call lookup_config%add( "type", 2, my_name )
401 case( 2 )
402 call lookup_config%add( "type", 1, my_name )
403 case default
404 call die( 373295283 )
405 end select
406 mock_optics_lookup = optics_lookup_t( lookup_config )
407 call mock_optics_lookup%get_optics( (/ ( 0.0_dk, 0.0_dk ), &
408 ( 0.0_dk, 0.0_dk ) /), &
409 absorption = cheby_coeff )
410 allocate( raw_state( mode_state%raw_size( ) ) )
411 call mode_state%dump_state( raw_state )
412 ! IMPORTANT: this is just for testing - never rely on the raw state data
413 ! outside of tests. It will change often
414 associate( wet_diameter => raw_state(1), &
415 foo_mass => raw_state(3), &
416 bar_mass => raw_state(4) )
417 wet_surface_radius = ( 0.5_dk * wet_diameter ) &
418 * exp( 2.0_dk * ( log( 2.4_dk ) )**2 )
419 normal_radius = mock_optics_lookup%normalize_radius( wet_surface_radius )
420 call chebyshev_function( normal_radius, cheby_func )
421 specific_absorption(1) = &
422 weighted_chebyshev( 3, cheby_coeff(:,1), cheby_func ) &
423 ! wet volume to mass mixing ratio
424 * ( foo_mass * 6 + bar_mass * 3 ) * kwaterdensitystp
425 specific_absorption(1) = max( specific_absorption(1), 0.0_dk )
426 specific_absorption(2) = &
427 weighted_chebyshev( 3, cheby_coeff(:,2), cheby_func ) &
428 ! wet volume to mass mixing ratio
429 * ( foo_mass * 6 + bar_mass * 3 ) * kwaterdensitystp
430 specific_absorption(2) = max( specific_absorption(2), 0.0_dk )
431 end associate
432
433 compare_values = mode%specific_absorption__m2_kg( mode_state, 2, 3, &
434 cheby_coeff, cheby_func )
435 call assert( 715576209, compare_values(1) .gt. 0.0_dk )
436 call assert( 770055995, compare_values(1) .lt. 1e200_dk )
437 call assert( 446787188, &
438 almost_equal( compare_values(1), specific_absorption(1) ) )
439 call assert( 654378877, compare_values(2) .gt. 0.0_dk )
440 call assert( 201746724, compare_values(2) .lt. 1e200_dk )
441 call assert( 501266974, &
442 almost_equal( compare_values(2), specific_absorption(2) ) )
443 compare_values = &
444 mode%specific_absorption__m2_kg( mode_state, 2, 3,cheby_coeff, &
445 cheby_func, &
446 max_absorption = (/ 0.0_dk, 0.0_dk /) )
447 call assert( 656083252, compare_values(1) .eq. 0.0_dk )
448 call assert( 485926348, compare_values(2) .eq. 0.0_dk )
449
450 end function specific_absorption
451
452!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
453
454 function specific_extinction( mode, mode_state, long_or_short )
455
456 use mam_constants, only : kwaterdensitystp
457 use mam_mode, only : mode_t, mode_state_t
458 use mam_optics_lookup, only : optics_lookup_t
459 use musica_assert, only : assert, die, almost_equal
460 use musica_config, only : config_t
461 use musica_constants, only : dk => musica_dk
462 use musica_math, only : chebyshev_function, &
463 weighted_chebyshev
464
465 real(kind=dk) :: specific_extinction(2)
466 class(mode_t), intent(in) :: mode
467 class(mode_state_t), intent(in) :: mode_state
468 integer, intent(in) :: long_or_short ! 1=long, 2=short
469
470 character(len=*), parameter :: my_name = "mode extinction test"
471 type(config_t) :: lookup_config
472 type(optics_lookup_t) :: mock_optics_lookup
473 real(kind=dk) :: cheby_coeff(3,2)
474 real(kind=dk) :: cheby_func(3)
475 real(kind=dk) :: normal_radius, wet_surface_radius
476 real(kind=dk) :: compare_values(2)
477 real(kind=dk), allocatable :: raw_state(:)
478
479 call lookup_config%empty( )
480 select case( long_or_short )
481 case( 1 )
482 call lookup_config%add( "type", 2, my_name )
483 case( 2 )
484 call lookup_config%add( "type", 1, my_name )
485 case default
486 call die( 995833204 )
487 end select
488 mock_optics_lookup = optics_lookup_t( lookup_config )
489 call mock_optics_lookup%get_optics( (/ ( 0.0_dk, 0.0_dk ), &
490 ( 0.0_dk, 0.0_dk ) /), &
491 extinction = cheby_coeff )
492 allocate( raw_state( mode_state%raw_size( ) ) )
493 call mode_state%dump_state( raw_state )
494 ! IMPORTANT: this is just for testing - never rely on the raw state data
495 ! outside of tests. It will change often
496 associate( wet_diameter => raw_state(1), &
497 foo_mass => raw_state(3), &
498 bar_mass => raw_state(4) )
499 wet_surface_radius = ( 0.5_dk * wet_diameter ) &
500 * exp( 2.0_dk * ( log( 2.4_dk ) )**2 )
501 normal_radius = mock_optics_lookup%normalize_radius( wet_surface_radius )
502 call chebyshev_function( normal_radius, cheby_func )
503 if( wet_surface_radius .le. mock_optics_lookup%maximum_radius__m( ) ) then
504 specific_extinction(1) = &
505 exp( weighted_chebyshev( 3, cheby_coeff(:,1), cheby_func ) )
506 specific_extinction(2) = &
507 exp( weighted_chebyshev( 3, cheby_coeff(:,2), cheby_func ) )
508 else
509 specific_extinction(:) = &
510 1.5_dk / ( wet_surface_radius * kwaterdensitystp )
511 end if
512 specific_extinction(:) = specific_extinction(:) &
513 ! wet volume to mass mixing ratio
514 * ( foo_mass * 6 + bar_mass * 3 ) * kwaterdensitystp
515 end associate
516
517 compare_values = mode%specific_extinction__m2_kg( mode_state, 2, 3, &
518 cheby_coeff, cheby_func,&
519 mock_optics_lookup )
520 call assert( 933495353, compare_values(1) .gt. 0.0_dk )
521 call assert( 480863200, compare_values(1) .lt. 1e200_dk )
522 call assert( 108102181, &
523 almost_equal( compare_values(1), specific_extinction(1) ) )
524 call assert( 310706296, compare_values(2) .gt. 0.0_dk )
525 call assert( 140549392, compare_values(2) .lt. 1e200_dk )
526 call assert( 837945276, &
527 almost_equal( compare_values(2), specific_extinction(2) ) )
528
529 end function specific_extinction
530
531!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
532
533 function asymmetry_factor( mode, mode_state, long_or_short )
534
535 use mam_constants, only : kwaterdensitystp
536 use mam_mode, only : mode_t, mode_state_t
537 use mam_optics_lookup, only : optics_lookup_t
538 use musica_assert, only : assert, die, almost_equal
539 use musica_config, only : config_t
540 use musica_constants, only : dk => musica_dk
541 use musica_math, only : chebyshev_function, &
542 weighted_chebyshev
543
544 real(kind=dk) :: asymmetry_factor(2)
545 class(mode_t), intent(in) :: mode
546 class(mode_state_t), intent(in) :: mode_state
547 integer, intent(in) :: long_or_short ! 1=long, 2=short
548
549 character(len=*), parameter :: my_name = "mode asymmetry factor test"
550 type(config_t) :: lookup_config
551 type(optics_lookup_t) :: mock_optics_lookup
552 real(kind=dk) :: cheby_coeff(3,2)
553 real(kind=dk) :: cheby_func(3)
554 real(kind=dk) :: normal_radius, wet_surface_radius
555 real(kind=dk) :: compare_values(2)
556 real(kind=dk), allocatable :: raw_state(:)
557
558 call lookup_config%empty( )
559 select case( long_or_short )
560 case( 1 )
561 call lookup_config%add( "type", 2, my_name )
562 case( 2 )
563 call lookup_config%add( "type", 1, my_name )
564 case default
565 call die( 567537861 )
566 end select
567 mock_optics_lookup = optics_lookup_t( lookup_config )
568 call mock_optics_lookup%get_optics( (/ ( 0.0_dk, 0.0_dk ), &
569 ( 0.0_dk, 0.0_dk ) /), &
570 asymmetry_factor = cheby_coeff )
571 allocate( raw_state( mode_state%raw_size( ) ) )
572 call mode_state%dump_state( raw_state )
573 ! IMPORTANT: this is just for testing - never rely on the raw state data
574 ! outside of tests. It will change often
575 associate( wet_diameter => raw_state(1) )
576 wet_surface_radius = ( 0.5_dk * wet_diameter ) &
577 * exp( 2.0_dk * ( log( 2.4_dk ) )**2 )
578 normal_radius = mock_optics_lookup%normalize_radius( wet_surface_radius )
579 call chebyshev_function( normal_radius, cheby_func )
580 asymmetry_factor(1) = weighted_chebyshev( 3, cheby_coeff(:,1), cheby_func )
581 asymmetry_factor(2) = weighted_chebyshev( 3, cheby_coeff(:,2), cheby_func )
582 end associate
583
584 compare_values = mode%asymmetry_factor( mode_state, 2, 3, cheby_coeff, &
585 cheby_func )
586 call assert( 174649801, compare_values(1) .gt. 0.0_dk )
587 call assert( 904492896, compare_values(1) .lt. 1e200_dk )
588 call assert( 734335992, &
589 almost_equal( compare_values(1), asymmetry_factor(1) ) )
590 call assert( 281703839, compare_values(2) .gt. 0.0_dk )
591 call assert( 394022184, compare_values(2) .lt. 1e200_dk )
592 call assert( 223865280, &
593 almost_equal( compare_values(2), asymmetry_factor(2) ) )
594
595 end function asymmetry_factor
596
597!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
598
599end program test_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_optics_lookup
The optics_lookup_t type and related functions.
Definition: optics_lookup.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
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
specific_absorption
real(kind=dk) function, dimension(2) specific_absorption(mode, mode_state, long_or_short)
Definition: mode.F90:373
test_mode_t
subroutine test_mode_t()
Definition: mode.F90:19
net_refractive_index
complex(kind=dk) function, dimension(2) net_refractive_index(mode_state, long_or_short)
Definition: mode.F90:329
specific_extinction
real(kind=dk) function, dimension(2) specific_extinction(mode, mode_state, long_or_short)
Definition: mode.F90:455
test_mode
program test_mode
Test program for the mode_t type and related functions.
Definition: mode.F90:8
Generated on Mon Aug 15 2022 15:01:05 for mam by doxygen 1.9.4