Line data Source code
1 : ! air_composition module defines major species of the atmosphere and manages
2 : ! the physical properties that are dependent on the composition of air
3 : module air_composition
4 :
5 : use shr_kind_mod, only: r8 => shr_kind_r8
6 : use cam_abortutils, only: endrun
7 :
8 : implicit none
9 : private
10 : save
11 :
12 : public :: air_composition_readnl
13 : public :: air_composition_init
14 : public :: dry_air_composition_update
15 : public :: water_composition_update
16 :
17 : ! get_cp_dry: (generalized) heat capacity for dry air
18 : public :: get_cp_dry
19 : ! get_cp: (generalized) heat capacity
20 : public :: get_cp
21 : ! get_R_dry: (generalized) dry air gas constant
22 : public :: get_R_dry
23 : ! get_R: Compute generalized R
24 : public :: get_R
25 : ! get_mbarv: molecular weight of dry air
26 : public :: get_mbarv
27 :
28 : private :: air_species_info
29 :
30 : integer, parameter :: unseti = -HUGE(1)
31 : real(r8), parameter :: unsetr = HUGE(1.0_r8)
32 :
33 : ! composition of air
34 : !
35 : integer, parameter :: num_names_max = 20 ! Should match namelist definition
36 : character(len=6) :: dry_air_species(num_names_max)
37 : character(len=6) :: water_species_in_air(num_names_max)
38 :
39 : integer, protected, public :: dry_air_species_num
40 : integer, protected, public :: water_species_in_air_num
41 :
42 : ! Thermodynamic variables
43 : integer, protected, public :: thermodynamic_active_species_num = unseti
44 : integer, allocatable, protected, public :: thermodynamic_active_species_idx(:)
45 : integer, allocatable, public :: thermodynamic_active_species_idx_dycore(:)
46 : real(r8), allocatable, protected, public :: thermodynamic_active_species_cp(:)
47 : real(r8), allocatable, protected, public :: thermodynamic_active_species_cv(:)
48 : real(r8), allocatable, protected, public :: thermodynamic_active_species_R(:)
49 : ! thermodynamic_active_species_mwi: inverse molecular weights dry air
50 : real(r8), allocatable, protected, public :: thermodynamic_active_species_mwi(:)
51 : ! thermodynamic_active_species_kv: molecular diffusion
52 : real(r8), allocatable, protected, public :: thermodynamic_active_species_kv(:)
53 : ! thermodynamic_active_species_kc: thermal conductivity
54 : real(r8), allocatable, protected, public :: thermodynamic_active_species_kc(:)
55 : !
56 : ! for energy computations liquid and ice species need to be identified
57 : !
58 : ! thermodynamic_active_species_liq_num: number of liquid water species
59 : integer, protected, public :: thermodynamic_active_species_liq_num = unseti
60 : ! thermodynamic_active_species_ice_num: number of frozen water species
61 : integer, protected, public :: thermodynamic_active_species_ice_num = unseti
62 : ! thermodynamic_active_species_liq_idx: index of liquid water species
63 : integer, allocatable, protected, public :: thermodynamic_active_species_liq_idx(:)
64 : ! thermodynamic_active_species_liq_idx_dycore: index of liquid water species
65 : integer, allocatable, public :: thermodynamic_active_species_liq_idx_dycore(:)
66 : ! thermodynamic_active_species_ice_idx: index of ice water species
67 : integer, allocatable, protected, public :: thermodynamic_active_species_ice_idx(:)
68 : ! thermodynamic_active_species_ice_idx_dycore: index of ice water species
69 : integer, allocatable, public :: thermodynamic_active_species_ice_idx_dycore(:)
70 : ! enthalpy_reference_state: choices: 'ice', 'liq', 'wv'
71 : character(len=3), public, protected :: enthalpy_reference_state = 'xxx'
72 :
73 : integer, protected, public :: wv_idx = -1 ! Water vapor index
74 :
75 : !------------- Variables for consistent themodynamics --------------------
76 : !
77 :
78 : ! standard dry air (constant composition)
79 : real(r8), public, protected :: mmro2 = unsetr ! Mass mixing ratio of O2
80 : real(r8), public, protected :: mmrn2 = unsetr ! Mass mixing ratio of N2
81 : real(r8), public, protected :: o2_mwi = unsetr ! Inverse mol. weight of O2
82 : real(r8), public, protected :: n2_mwi = unsetr ! Inverse mol. weight of N2
83 : real(r8), public, protected :: mbar = unsetr ! Mean mass at mid level
84 :
85 : ! coefficients in expressions for molecular diffusion coefficients
86 : ! kv1,..,kv3 are coefficients for kmvis calculation
87 : ! kc1,..,kc3 are coefficients for kmcnd calculation
88 : ! Liu, H.-L., et al. (2010), Thermosphere extension of the Whole Atmosphere Community Climate Model,
89 : ! J. Geophys. Res., 115, A12302, doi:10.1029/2010JA015586.
90 : real(r8), public, parameter :: kv1 = 4.03_r8 * 1.e-7_r8
91 : real(r8), public, parameter :: kv2 = 3.42_r8 * 1.e-7_r8
92 : real(r8), public, parameter :: kv3 = 3.9_r8 * 1.e-7_r8
93 : real(r8), public, parameter :: kc1 = 56._r8 * 1.e-5_r8
94 : real(r8), public, parameter :: kc2 = 56._r8 * 1.e-5_r8
95 : real(r8), public, parameter :: kc3 = 75.9_r8 * 1.e-5_r8
96 :
97 : real(r8), public, parameter :: kv_temp_exp = 0.69_r8
98 : real(r8), public, parameter :: kc_temp_exp = 0.69_r8
99 :
100 : ! cpairv: composition dependent specific heat at constant pressure
101 : real(r8), public, protected, allocatable :: cpairv(:,:,:)
102 : ! rairv: composition dependent gas "constant"
103 : real(r8), public, protected, allocatable :: rairv(:,:,:)
104 : ! cappav: rairv / cpairv
105 : real(r8), public, protected, allocatable :: cappav(:,:,:)
106 : ! mbarv: composition dependent atmosphere mean mass
107 : real(r8), public, protected, allocatable :: mbarv(:,:,:)
108 : ! cp_or_cv_dycore: enthalpy or internal energy scaling factor for
109 : ! energy consistency
110 : real(r8), public, protected, allocatable :: cp_or_cv_dycore(:,:,:)
111 : !
112 : ! Interfaces for public routines
113 : interface get_cp_dry
114 : module procedure get_cp_dry_1hd
115 : module procedure get_cp_dry_2hd
116 : end interface get_cp_dry
117 :
118 : interface get_cp
119 : module procedure get_cp_1hd
120 : module procedure get_cp_2hd
121 : end interface get_cp
122 :
123 : interface get_R_dry
124 : module procedure get_R_dry_1hd
125 : module procedure get_R_dry_2hd
126 : end interface get_R_dry
127 :
128 : interface get_R
129 : module procedure get_R_1hd
130 : module procedure get_R_2hd
131 : end interface get_R
132 :
133 : interface get_mbarv
134 : module procedure get_mbarv_1hd
135 : end interface get_mbarv
136 :
137 : CONTAINS
138 :
139 : ! Read namelist variables.
140 768 : subroutine air_composition_readnl(nlfile)
141 : use namelist_utils, only: find_group_name
142 : use spmd_utils, only: masterproc, mpicom, masterprocid
143 : use spmd_utils, only: mpi_character
144 : use cam_logfile, only: iulog
145 :
146 : ! Dummy argument: filepath for file containing namelist input
147 : character(len=*), intent(in) :: nlfile
148 :
149 : ! Local variables
150 : integer :: unitn, ierr, indx
151 : integer, parameter :: lsize = 76
152 : character(len=*), parameter :: subname = 'air_composition_readnl :: '
153 : character(len=lsize) :: banner
154 : character(len=lsize) :: bline
155 :
156 : ! Variable components of dry air and water species in air
157 : namelist /air_composition_nl/ dry_air_species, water_species_in_air
158 : !-----------------------------------------------------------------------
159 :
160 768 : banner = repeat('*', lsize)
161 768 : bline = "***"//repeat(' ', lsize - 6)//"***"
162 :
163 : ! Read variable components of dry air and water species in air
164 16128 : dry_air_species = (/ (' ', indx = 1, num_names_max) /)
165 16128 : water_species_in_air = (/ (' ', indx = 1, num_names_max) /)
166 :
167 768 : if (masterproc) then
168 2 : open(newunit=unitn, file=trim(nlfile), status='old')
169 2 : call find_group_name(unitn, 'air_composition_nl', status=ierr)
170 2 : if (ierr == 0) then
171 2 : read(unitn, air_composition_nl, iostat=ierr)
172 2 : if (ierr /= 0) then
173 0 : call endrun(subname//'ERROR reading namelist, air_composition_nl')
174 : end if
175 : end if
176 2 : close(unitn)
177 : end if
178 :
179 : call mpi_bcast(dry_air_species, len(dry_air_species)*num_names_max, &
180 768 : mpi_character, masterprocid, mpicom, ierr)
181 768 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: dry_air_species")
182 : call mpi_bcast(water_species_in_air, &
183 : len(water_species_in_air)*num_names_max, mpi_character, &
184 768 : masterprocid, mpicom, ierr)
185 768 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: water_species_in_air")
186 :
187 768 : dry_air_species_num = 0
188 768 : water_species_in_air_num = 0
189 16128 : do indx = 1, num_names_max
190 15360 : if ( (LEN_TRIM(dry_air_species(indx)) > 0) .and. &
191 : (TRIM(dry_air_species(indx)) /= 'N2')) then
192 2304 : dry_air_species_num = dry_air_species_num + 1
193 : end if
194 16128 : if (LEN_TRIM(water_species_in_air(indx)) > 0) then
195 3840 : water_species_in_air_num = water_species_in_air_num + 1
196 : end if
197 : end do
198 :
199 : ! Initialize number of thermodynamically active species
200 : thermodynamic_active_species_num = &
201 768 : dry_air_species_num + water_species_in_air_num
202 :
203 768 : if (masterproc) then
204 2 : write(iulog, *) banner
205 2 : write(iulog, *) bline
206 :
207 2 : if (dry_air_species_num == 0) then
208 0 : write(iulog, *) " Thermodynamic properties of dry air are ", &
209 0 : "fixed at troposphere values"
210 : else
211 2 : write(iulog, *) " Thermodynamic properties of dry air are ", &
212 4 : "based on variable composition of the following species:"
213 8 : do indx = 1, dry_air_species_num
214 8 : write(iulog, *) ' ', trim(dry_air_species(indx))
215 : end do
216 2 : write(iulog,*) ' '
217 : end if
218 2 : write(iulog,*) " Thermodynamic properties of moist air are ", &
219 4 : "based on variable composition of the following water species:"
220 12 : do indx = 1, water_species_in_air_num
221 12 : write(iulog, *) ' ', trim(water_species_in_air(indx))
222 : end do
223 2 : write(iulog, *) bline
224 2 : write(iulog, *) banner
225 : end if
226 :
227 768 : end subroutine air_composition_readnl
228 :
229 : !===========================================================================
230 :
231 768 : subroutine air_composition_init()
232 : use string_utils, only: int2str
233 : use spmd_utils, only: masterproc
234 : use cam_logfile, only: iulog
235 : use physconst, only: r_universal, cpair, rair, cpwv, rh2o, cpliq, cpice, mwdry
236 : use constituents, only: cnst_get_ind, cnst_mw
237 : use ppgrid, only: pcols, pver, begchunk, endchunk
238 : integer :: icnst, ix, isize, ierr, idx
239 : integer :: liq_num, ice_num
240 1536 : integer :: liq_idx(water_species_in_air_num)
241 1536 : integer :: ice_idx(water_species_in_air_num)
242 : logical :: has_liq, has_ice
243 : real(r8) :: mw
244 :
245 : character(len=*), parameter :: subname = 'composition_init'
246 : character(len=*), parameter :: errstr = subname//": failed to allocate "
247 :
248 : !
249 : ! define cp and R for species in species_name
250 : !
251 : ! Last major species in namelist dry_air_species is derived from the
252 : ! other major species (since the sum of dry mixing ratios for
253 : ! major species of dry air add must add to one)
254 : !
255 : ! cv = R * dofx / 2; cp = R * (1 + (dofx / 2))
256 : ! DOF == Degrees of Freedom
257 : ! dof1 = monatomic ideal gas, 3 translational DOF
258 : real(r8), parameter :: dof1 = 3._r8
259 : real(r8), parameter :: cv1 = 0.5_r8 * r_universal * dof1
260 : real(r8), parameter :: cp1 = 0.5_r8 * r_universal * (2._r8 + dof1)
261 : ! dof2 = diatomic ideal gas, 3 translational + 2 rotational = 5 DOF
262 : real(r8), parameter :: dof2 = 5._r8
263 : real(r8), parameter :: cv2 = 0.5_r8 * r_universal * dof2
264 : real(r8), parameter :: cp2 = 0.5_r8 * r_universal * (2._r8 + dof2)
265 : ! dof3 = polyatomic ideal gas, 3 translational + 3 rotational = 6 DOF
266 : real(r8), parameter :: dof3 = 6._r8
267 : real(r8), parameter :: cv3 = 0.5_r8 * r_universal * dof3
268 : real(r8), parameter :: cp3 = 0.5_r8 * r_universal * (2._r8 + dof3)
269 :
270 768 : liq_num = 0
271 768 : ice_num = 0
272 768 : has_liq = .false.
273 768 : has_ice = .false.
274 : ! standard dry air (constant composition)
275 768 : o2_mwi = 1._r8 / 32._r8
276 768 : n2_mwi = 1._r8 / 28._r8
277 768 : mmro2 = 0.235_r8
278 768 : mmrn2 = 0.765_r8
279 768 : mbar = 1._r8 / ((mmro2 * o2_mwi) + (mmrn2 * n2_mwi))
280 :
281 : ! init for variable composition dry air
282 :
283 768 : isize = dry_air_species_num + water_species_in_air_num
284 2304 : allocate(thermodynamic_active_species_idx(isize), stat=ierr)
285 768 : if (ierr /= 0) then
286 0 : call endrun(errstr//"thermodynamic_active_species_idx")
287 : end if
288 1536 : allocate(thermodynamic_active_species_idx_dycore(isize), stat=ierr)
289 768 : if (ierr /= 0) then
290 0 : call endrun(errstr//"thermodynamic_active_species_idx_dycore")
291 : end if
292 2304 : allocate(thermodynamic_active_species_cp(0:isize), stat=ierr)
293 768 : if (ierr /= 0) then
294 0 : call endrun(errstr//"thermodynamic_active_species_cp")
295 : end if
296 1536 : allocate(thermodynamic_active_species_cv(0:isize), stat=ierr)
297 768 : if (ierr /= 0) then
298 0 : call endrun(errstr//"thermodynamic_active_species_cv")
299 : end if
300 1536 : allocate(thermodynamic_active_species_R(0:isize), stat=ierr)
301 768 : if (ierr /= 0) then
302 0 : call endrun(errstr//"thermodynamic_active_species_R")
303 : end if
304 :
305 768 : isize = dry_air_species_num
306 2304 : allocate(thermodynamic_active_species_mwi(0:isize), stat=ierr)
307 768 : if (ierr /= 0) then
308 0 : call endrun(errstr//"thermodynamic_active_species_mwi")
309 : end if
310 1536 : allocate(thermodynamic_active_species_kv(0:isize), stat=ierr)
311 768 : if (ierr /= 0) then
312 0 : call endrun(errstr//"thermodynamic_active_species_kv")
313 : end if
314 1536 : allocate(thermodynamic_active_species_kc(0:isize), stat=ierr)
315 768 : if (ierr /= 0) then
316 0 : call endrun(errstr//"thermodynamic_active_species_kc")
317 : end if
318 : !------------------------------------------------------------------------
319 : ! Allocate constituent dependent properties
320 : !------------------------------------------------------------------------
321 2304 : allocate(cpairv(pcols,pver,begchunk:endchunk), stat=ierr)
322 768 : if (ierr /= 0) then
323 0 : call endrun(errstr//"cpairv")
324 : end if
325 2304 : allocate(rairv(pcols,pver,begchunk:endchunk), stat=ierr)
326 768 : if (ierr /= 0) then
327 0 : call endrun(errstr//"rairv")
328 : end if
329 2304 : allocate(cappav(pcols,pver,begchunk:endchunk), stat=ierr)
330 768 : if (ierr /= 0) then
331 0 : call endrun(errstr//"cappav")
332 : end if
333 2304 : allocate(mbarv(pcols,pver,begchunk:endchunk), stat=ierr)
334 768 : if (ierr /= 0) then
335 0 : call endrun(errstr//"mbarv")
336 : end if
337 2304 : allocate(cp_or_cv_dycore(pcols,pver,begchunk:endchunk), stat=ierr)
338 768 : if (ierr /= 0) then
339 0 : call endrun(errstr//"cp_or_cv_dycore")
340 : end if
341 :
342 6912 : thermodynamic_active_species_idx = -HUGE(1)
343 6912 : thermodynamic_active_species_idx_dycore = -HUGE(1)
344 7680 : thermodynamic_active_species_cp = 0.0_r8
345 7680 : thermodynamic_active_species_cv = 0.0_r8
346 7680 : thermodynamic_active_species_R = 0.0_r8
347 3840 : thermodynamic_active_species_mwi = 0.0_r8
348 3840 : thermodynamic_active_species_kv = 0.0_r8
349 3840 : thermodynamic_active_species_kc = 0.0_r8
350 : !------------------------------------------------------------------------
351 : ! Initialize constituent dependent properties
352 : !------------------------------------------------------------------------
353 5094912 : cpairv(:pcols, :pver, begchunk:endchunk) = cpair
354 5094912 : rairv(:pcols, :pver, begchunk:endchunk) = rair
355 5094912 : cappav(:pcols, :pver, begchunk:endchunk) = rair / cpair
356 5094912 : mbarv(:pcols, :pver, begchunk:endchunk) = mwdry
357 : !
358 768 : if (dry_air_species_num > 0) then
359 : !
360 : ! The last major species in dry_air_species is derived from the
361 : ! others and constants associated with it are initialized here
362 : !
363 768 : if (TRIM(dry_air_species(dry_air_species_num + 1)) == 'N2') then
364 768 : call air_species_info('N', ix, mw)
365 768 : mw = 2.0_r8 * mw
366 768 : icnst = 0 ! index for the derived tracer N2
367 768 : thermodynamic_active_species_cp(icnst) = cp2 / mw
368 768 : thermodynamic_active_species_cv(icnst) = cv2 / mw !N2
369 768 : thermodynamic_active_species_R (icnst) = r_universal / mw
370 768 : thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
371 768 : thermodynamic_active_species_kv(icnst) = kv2
372 768 : thermodynamic_active_species_kc(icnst) = kc2
373 : !
374 : ! if last major species is not N2 then add code here
375 : !
376 : else
377 0 : write(iulog, *) subname, ' derived major species not found: ', &
378 0 : dry_air_species(dry_air_species_num)
379 0 : call endrun(subname//': derived major species not found')
380 : end if
381 : else
382 : !
383 : ! dry air is not species dependent
384 : !
385 0 : icnst = 0
386 0 : thermodynamic_active_species_cp (icnst) = cpair
387 0 : thermodynamic_active_species_cv (icnst) = cpair - rair
388 0 : thermodynamic_active_species_R (icnst) = rair
389 : end if
390 : !
391 : !************************************************************************
392 : !
393 : ! add prognostic components of dry air
394 : !
395 : !************************************************************************
396 : !
397 768 : icnst = 1
398 3072 : do idx = 1, dry_air_species_num
399 4608 : select case (TRIM(dry_air_species(idx)))
400 : !
401 : ! O
402 : !
403 : case('O')
404 768 : call air_species_info('O', ix, mw)
405 768 : thermodynamic_active_species_idx(icnst) = ix
406 768 : thermodynamic_active_species_cp (icnst) = cp1 / mw
407 768 : thermodynamic_active_species_cv (icnst) = cv1 / mw
408 768 : thermodynamic_active_species_R (icnst) = r_universal / mw
409 768 : thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
410 768 : thermodynamic_active_species_kv(icnst) = kv3
411 768 : thermodynamic_active_species_kc(icnst) = kc3
412 768 : icnst = icnst + 1
413 : !
414 : ! O2
415 : !
416 : case('O2')
417 768 : call air_species_info('O2', ix, mw)
418 768 : thermodynamic_active_species_idx(icnst) = ix
419 768 : thermodynamic_active_species_cp (icnst) = cp2 / mw
420 768 : thermodynamic_active_species_cv (icnst) = cv2 / mw
421 768 : thermodynamic_active_species_R (icnst) = r_universal / mw
422 768 : thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
423 768 : thermodynamic_active_species_kv(icnst) = kv1
424 768 : thermodynamic_active_species_kc(icnst) = kc1
425 768 : icnst = icnst + 1
426 : !
427 : ! H
428 : !
429 : case('H')
430 768 : call air_species_info('H', ix, mw)
431 768 : thermodynamic_active_species_idx(icnst) = ix
432 768 : thermodynamic_active_species_cp (icnst) = cp1 / mw
433 768 : thermodynamic_active_species_cv (icnst) = cv1 / mw
434 768 : thermodynamic_active_species_R (icnst) = r_universal / mw
435 768 : thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
436 : ! Hydrogen not included in calculation of diffusivity and conductivity
437 768 : thermodynamic_active_species_kv(icnst) = 0.0_r8
438 768 : thermodynamic_active_species_kc(icnst) = 0.0_r8
439 768 : icnst = icnst + 1
440 : !
441 : ! If support for more major species is to be included add code here
442 : !
443 : case default
444 0 : write(iulog, *) subname, ' dry air component not found: ', &
445 0 : dry_air_species(idx)
446 4608 : call endrun(subname//': dry air component not found')
447 : end select
448 :
449 3072 : if (masterproc) then
450 6 : write(iulog, *) "Dry air composition ", &
451 6 : TRIM(dry_air_species(idx)), &
452 6 : icnst-1,thermodynamic_active_species_idx(icnst-1), &
453 6 : thermodynamic_active_species_mwi(icnst-1), &
454 6 : thermodynamic_active_species_cp(icnst-1), &
455 12 : thermodynamic_active_species_cv(icnst-1)
456 : end if
457 : end do
458 768 : isize = dry_air_species_num+1
459 768 : icnst = 0 ! N2
460 768 : if(isize > 0) then
461 768 : if(masterproc) then
462 2 : write(iulog, *) "Dry air composition ", &
463 2 : TRIM(dry_air_species(idx)), &
464 2 : icnst, -1, thermodynamic_active_species_mwi(icnst), &
465 2 : thermodynamic_active_species_cp(icnst), &
466 4 : thermodynamic_active_species_cv(icnst)
467 : end if
468 : end if
469 : !
470 : !************************************************************************
471 : !
472 : ! Add non-dry components of moist air (water vapor and condensates)
473 : !
474 : !************************************************************************
475 : !
476 768 : icnst = dry_air_species_num + 1
477 4608 : do idx = 1, water_species_in_air_num
478 7680 : select case (TRIM(water_species_in_air(idx)))
479 : !
480 : ! Q
481 : !
482 : case('Q')
483 768 : call air_species_info('Q', ix, mw)
484 768 : wv_idx = ix
485 768 : thermodynamic_active_species_idx(icnst) = ix
486 768 : thermodynamic_active_species_cp (icnst) = cpwv
487 768 : thermodynamic_active_species_cv (icnst) = cv3 / mw
488 768 : thermodynamic_active_species_R (icnst) = rh2o
489 768 : icnst = icnst + 1
490 : !
491 : ! CLDLIQ
492 : !
493 : case('CLDLIQ')
494 768 : call air_species_info('CLDLIQ', ix, mw)
495 768 : thermodynamic_active_species_idx(icnst) = ix
496 768 : thermodynamic_active_species_cp (icnst) = cpliq
497 768 : thermodynamic_active_species_cv (icnst) = cpliq
498 768 : liq_num = liq_num+1
499 768 : liq_idx (liq_num) = ix
500 768 : icnst = icnst + 1
501 768 : has_liq = .true.
502 : !
503 : ! CLDICE
504 : !
505 : case('CLDICE')
506 768 : call air_species_info('CLDICE', ix, mw)
507 768 : thermodynamic_active_species_idx(icnst) = ix
508 768 : thermodynamic_active_species_cp (icnst) = cpice
509 768 : thermodynamic_active_species_cv (icnst) = cpice
510 768 : ice_num = ice_num+1
511 768 : ice_idx(ice_num) = ix
512 768 : icnst = icnst + 1
513 768 : has_ice = .true.
514 : !
515 : ! RAINQM
516 : !
517 : case('RAINQM')
518 768 : call air_species_info('RAINQM', ix, mw)
519 768 : thermodynamic_active_species_idx(icnst) = ix
520 768 : thermodynamic_active_species_cp (icnst) = cpliq
521 768 : thermodynamic_active_species_cv (icnst) = cpliq
522 768 : liq_num = liq_num+1
523 768 : liq_idx(liq_num) = ix
524 768 : icnst = icnst + 1
525 768 : has_liq = .true.
526 : !
527 : ! SNOWQM
528 : !
529 : case('SNOWQM')
530 768 : call air_species_info('SNOWQM', ix, mw)
531 768 : thermodynamic_active_species_idx(icnst) = ix
532 768 : thermodynamic_active_species_cp (icnst) = cpice
533 768 : thermodynamic_active_species_cv (icnst) = cpice
534 768 : ice_num = ice_num+1
535 768 : ice_idx(ice_num) = ix
536 768 : icnst = icnst + 1
537 768 : has_ice = .true.
538 : !
539 : ! GRAUQM
540 : !
541 : case('GRAUQM')
542 0 : call air_species_info('GRAUQM', ix, mw)
543 0 : thermodynamic_active_species_idx(icnst) = ix
544 0 : thermodynamic_active_species_cp (icnst) = cpice
545 0 : thermodynamic_active_species_cv (icnst) = cpice
546 0 : ice_num = ice_num+1
547 0 : ice_idx(ice_num) = ix
548 0 : icnst = icnst + 1
549 0 : has_ice = .true.
550 : !
551 : ! If support for more major species is to be included add code here
552 : !
553 : case default
554 0 : write(iulog, *) subname, ' moist air component not found: ', &
555 0 : water_species_in_air(idx)
556 7680 : call endrun(subname//': moist air component not found')
557 : end select
558 : !
559 : !
560 : !
561 3840 : if (masterproc) then
562 10 : write(iulog, *) "Thermodynamic active species ", &
563 20 : TRIM(water_species_in_air(idx))
564 10 : write(iulog, *) " global index : ", &
565 20 : icnst-1
566 10 : write(iulog, *) " thermodynamic_active_species_idx : ", &
567 20 : thermodynamic_active_species_idx(icnst-1)
568 10 : write(iulog, *) " cp : ", &
569 20 : thermodynamic_active_species_cp(icnst-1)
570 10 : write(iulog, *) " cv : ", &
571 20 : thermodynamic_active_species_cv(icnst-1)
572 10 : if (has_liq) then
573 4 : write(iulog, *) " register phase (liquid or ice) :", &
574 8 : " liquid"
575 : end if
576 10 : if (has_ice) then
577 4 : write(iulog, *) " register phase (liquid or ice) :", &
578 8 : " ice"
579 : end if
580 10 : write(iulog, *) " "
581 : end if
582 3840 : has_liq = .false.
583 4608 : has_ice = .false.
584 : end do
585 :
586 2304 : allocate(thermodynamic_active_species_liq_idx(liq_num), stat=ierr)
587 768 : if (ierr /= 0) then
588 0 : call endrun(errstr//"thermodynamic_active_species_liq_idx")
589 : end if
590 1536 : allocate(thermodynamic_active_species_liq_idx_dycore(liq_num), stat=ierr)
591 768 : if (ierr /= 0) then
592 0 : call endrun(errstr//"thermodynamic_active_species_liq_idx_dycore")
593 : end if
594 2304 : allocate(thermodynamic_active_species_ice_idx(ice_num), stat=ierr)
595 768 : if (ierr /= 0) then
596 0 : call endrun(errstr//"thermodynamic_active_species_ice_idx")
597 : end if
598 1536 : allocate(thermodynamic_active_species_ice_idx_dycore(ice_num), stat=ierr)
599 768 : if (ierr /= 0) then
600 0 : call endrun(errstr//"thermodynamic_active_species_ice_idx_dycore")
601 : end if
602 :
603 3072 : thermodynamic_active_species_liq_idx = liq_idx(1:liq_num)
604 768 : thermodynamic_active_species_liq_num = liq_num
605 :
606 : ! array initialized by the dycore
607 2304 : thermodynamic_active_species_liq_idx_dycore = -99
608 :
609 3072 : thermodynamic_active_species_ice_idx = ice_idx(1:ice_num)
610 768 : thermodynamic_active_species_ice_num = ice_num
611 :
612 : ! array initialized by the dycore
613 2304 : thermodynamic_active_species_ice_idx_dycore = -99
614 :
615 768 : if (water_species_in_air_num /= 1 + liq_num+ice_num) then
616 0 : write(iulog, '(2a,2(i0,a))') subname, &
617 0 : " water_species_in_air_num = ", &
618 0 : water_species_in_air_num, ", should be ", &
619 0 : (1 + liq_num + ice_num), " (1 + liq_num + ice_num)"
620 0 : call endrun(subname//': water_species_in_air_num /= 1+liq_num+ice_num')
621 : end if
622 768 : enthalpy_reference_state = 'ice'
623 768 : if (masterproc) then
624 2 : write(iulog, *) 'Enthalpy reference state : ', &
625 4 : TRIM(enthalpy_reference_state)
626 : end if
627 768 : end subroutine air_composition_init
628 :
629 : !===========================================================================
630 : !-----------------------------------------------------------------------
631 : ! dry_air_composition_update: Update the physics "constants" that vary
632 : !-------------------------------------------------------------------------
633 : !===========================================================================
634 :
635 683136 : subroutine dry_air_composition_update(mmr, lchnk, ncol, to_dry_factor)
636 : use cam_abortutils, only: endrun
637 : !(mmr = dry mixing ratio, if not, use to_dry_factor to convert!)
638 : real(r8), intent(in) :: mmr(:,:,:) ! mixing ratios for species dependent dry air
639 : integer, intent(in) :: lchnk ! Chunk number
640 : integer, intent(in) :: ncol ! number of columns
641 : real(r8), optional, intent(in) :: to_dry_factor(:,:)
642 :
643 683136 : call get_R_dry(mmr(:ncol, :, :), thermodynamic_active_species_idx, &
644 2049408 : rairv(:ncol, :, lchnk), fact=to_dry_factor)
645 : call get_cp_dry(mmr(:ncol,:,:), thermodynamic_active_species_idx, &
646 1366272 : cpairv(:ncol,:,lchnk), fact=to_dry_factor)
647 : call get_mbarv(mmr(:ncol,:,:), thermodynamic_active_species_idx, &
648 1366272 : mbarv(:ncol,:,lchnk), fact=to_dry_factor)
649 1155182976 : cappav(:ncol,:,lchnk) = rairv(:ncol,:,lchnk) / cpairv(:ncol,:,lchnk)
650 683136 : end subroutine dry_air_composition_update
651 :
652 : !===========================================================================
653 : !---------------------------------------------------------------------------
654 : ! water_composition_update: Update generalized cp or cv depending on dycore
655 : !---------------------------------------------------------------------------
656 : !===========================================================================
657 :
658 21888 : subroutine water_composition_update(mmr, lchnk, ncol, vcoord, to_dry_factor)
659 : use cam_abortutils, only: endrun
660 : use string_utils, only: int2str
661 : use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure
662 : real(r8), intent(in) :: mmr(:,:,:) ! constituents array
663 : integer, intent(in) :: lchnk ! Chunk number
664 : integer, intent(in) :: ncol ! number of columns
665 : integer, intent(in) :: vcoord
666 : real(r8), optional, intent(in) :: to_dry_factor(:,:)
667 :
668 : character(len=*), parameter :: subname = 'water_composition_update'
669 :
670 21888 : if (vcoord==vc_dry_pressure) then
671 0 : call get_cp(mmr(:ncol,:,:),.false.,cp_or_cv_dycore(:ncol,:,lchnk), factor=to_dry_factor, &
672 0 : active_species_idx_dycore=thermodynamic_active_species_idx,cpdry=cpairv(:ncol,:,lchnk))
673 21888 : else if (vcoord==vc_height) then
674 0 : call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, &
675 0 : cp_or_cv_dycore(:ncol,:,lchnk), fact=to_dry_factor, Rdry=rairv(:ncol,:,lchnk))
676 : !
677 : ! internal energy coefficient for MPAS
678 : ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353)
679 : !
680 0 : cp_or_cv_dycore(:ncol,:,lchnk)=cp_or_cv_dycore(:ncol,:,lchnk)*&
681 0 : (cpairv(:ncol,:,lchnk)-rairv(:ncol,:,lchnk)) /rairv(:ncol,:,lchnk)
682 21888 : else if (vcoord==vc_moist_pressure) then
683 : ! no update needed for moist pressure vcoord
684 : else
685 0 : call endrun(subname//" vertical coordinate not supported; vcoord="// int2str(vcoord))
686 : end if
687 21888 : end subroutine water_composition_update
688 :
689 : !===========================================================================
690 : !***************************************************************************
691 : !
692 : ! get_cp_dry: Compute dry air heat capacity under constant pressure
693 : !
694 : !***************************************************************************
695 : !
696 876672 : subroutine get_cp_dry_1hd(tracer, active_species_idx, cp_dry, fact)
697 : use cam_abortutils, only: endrun
698 : use string_utils, only: int2str
699 : use physconst, only: cpair
700 :
701 : ! Dummy arguments
702 : ! tracer: tracer array
703 : real(r8), intent(in) :: tracer(:,:,:)
704 : integer, intent(in) :: active_species_idx(:)
705 : ! fact: optional dry pressure level thickness
706 : real(r8), optional, intent(in) :: fact(:,:)
707 : ! cp_dry: dry air heat capacity under constant pressure
708 : real(r8), intent(out) :: cp_dry(:,:)
709 :
710 : ! Local variables
711 : integer :: idx, kdx , m_cnst, qdx
712 : ! factor: dry pressure level thickness
713 1753344 : real(r8) :: factor(SIZE(cp_dry, 1), SIZE(cp_dry, 2))
714 1753344 : real(r8) :: residual(SIZE(cp_dry, 1), SIZE(cp_dry, 2))
715 : real(r8) :: mmr
716 : character(len=*), parameter :: subname = 'get_cp_dry_1hd: '
717 :
718 876672 : if (dry_air_species_num == 0) then
719 : ! dry air heat capacity not species dependent
720 0 : cp_dry = cpair
721 : else
722 : ! dry air heat capacity is species dependent
723 876672 : if (present(fact)) then
724 0 : if (SIZE(fact, 1) /= SIZE(factor, 1)) then
725 : call endrun(subname//"SIZE mismatch in dimension 1 "// &
726 0 : int2str(SIZE(fact, 1))//' /= '//int2str(SIZE(factor, 1)))
727 : end if
728 0 : if (SIZE(fact, 2) /= SIZE(factor, 2)) then
729 : call endrun(subname//"SIZE mismatch in dimension 2 "// &
730 0 : int2str(SIZE(fact, 2))//' /= '//int2str(SIZE(factor, 2)))
731 : end if
732 0 : factor = fact(:,:)
733 : else
734 1476686592 : factor = 1.0_r8
735 : end if
736 :
737 1476686592 : cp_dry = 0.0_r8
738 1476686592 : residual = 1.0_r8
739 3506688 : do qdx = 1, dry_air_species_num
740 2630016 : m_cnst = active_species_idx(qdx)
741 328111488 : do kdx = 1, SIZE(cp_dry, 2)
742 4430059776 : do idx = 1, SIZE(cp_dry, 1)
743 4102824960 : mmr = tracer(idx, kdx, m_cnst) * factor(idx, kdx)
744 4102824960 : cp_dry(idx, kdx) = cp_dry(idx, kdx) + &
745 4102824960 : (thermodynamic_active_species_cp(qdx) * mmr)
746 4427429760 : residual(idx, kdx) = residual(idx, kdx) - mmr
747 : end do
748 : end do
749 : end do
750 876672 : qdx = 0 ! N2
751 109078272 : do kdx = 1, SIZE(cp_dry, 2)
752 1476686592 : do idx = 1, SIZE(cp_dry, 1)
753 2735216640 : cp_dry(idx, kdx) = cp_dry(idx, kdx) + &
754 4211026560 : (thermodynamic_active_species_cp(qdx) * residual(idx, kdx))
755 : end do
756 : end do
757 : end if
758 876672 : end subroutine get_cp_dry_1hd
759 :
760 : !===========================================================================
761 :
762 64512 : subroutine get_cp_dry_2hd(tracer, active_species_idx, cp_dry, fact)
763 : ! Version of get_cp_dry for arrays that have a second horizontal index
764 :
765 : ! Dummy arguments
766 : ! tracer: tracer array
767 : real(r8), intent(in) :: tracer(:,:,:,:)
768 : integer, intent(in) :: active_species_idx(:)
769 : ! fact: optional dry pressure level thickness
770 : real(r8), optional, intent(in) :: fact(:,:,:)
771 : ! cp_dry: dry air heat capacity under constant pressure
772 : real(r8), intent(out) :: cp_dry(:,:,:)
773 :
774 : ! Local variable
775 : integer :: jdx
776 :
777 258048 : do jdx = 1, SIZE(cp_dry, 2)
778 258048 : if (present(fact)) then
779 : call get_cp_dry(tracer(:,jdx,:,:), active_species_idx, &
780 0 : cp_dry(:,jdx,:), fact=fact(:,jdx,:))
781 : else
782 : call get_cp_dry(tracer(:,jdx,:,:), active_species_idx, &
783 193536 : cp_dry(:,jdx,:))
784 : end if
785 : end do
786 :
787 64512 : end subroutine get_cp_dry_2hd
788 :
789 : !===========================================================================
790 : !
791 : !***************************************************************************
792 : !
793 : ! get_cp: Compute generalized heat capacity at constant pressure
794 : !
795 : !***************************************************************************
796 : !
797 0 : subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry)
798 : use cam_abortutils, only: endrun
799 : use string_utils, only: int2str
800 :
801 : ! Dummy arguments
802 : ! tracer: Tracer array
803 : !
804 : ! factor not present then tracer must be dry mixing ratio
805 : ! if factor present tracer*factor must be dry mixing ratio
806 : !
807 : real(r8), intent(in) :: tracer(:,:,:)
808 : ! inv_cp: output inverse cp instead of cp
809 : logical, intent(in) :: inv_cp
810 : real(r8), intent(out) :: cp(:,:)
811 : ! dp: if provided then tracer is mass not mixing ratio
812 : real(r8), optional, intent(in) :: factor(:,:)
813 : ! active_species_idx_dycore: array of indices for index of
814 : ! thermodynamic active species in dycore tracer array
815 : ! (if different from physics index)
816 : integer, optional, intent(in) :: active_species_idx_dycore(:)
817 : real(r8),optional, intent(in) :: cpdry(:,:)
818 :
819 : ! LOCAL VARIABLES
820 : integer :: qdx, itrac
821 0 : real(r8) :: sum_species(SIZE(cp, 1), SIZE(cp, 2))
822 0 : real(r8) :: sum_cp(SIZE(cp, 1), SIZE(cp, 2))
823 0 : real(r8) :: factor_local(SIZE(cp, 1), SIZE(cp, 2))
824 0 : integer :: idx_local(thermodynamic_active_species_num)
825 : character(LEN=*), parameter :: subname = 'get_cp_1hd: '
826 :
827 0 : if (present(active_species_idx_dycore)) then
828 0 : if (SIZE(active_species_idx_dycore) /= &
829 : thermodynamic_active_species_num) then
830 : call endrun(subname//"SIZE mismatch "// &
831 : int2str(SIZE(active_species_idx_dycore))//' /= '// &
832 0 : int2str(thermodynamic_active_species_num))
833 : end if
834 0 : idx_local = active_species_idx_dycore
835 : else
836 0 : idx_local = thermodynamic_active_species_idx
837 : end if
838 :
839 0 : if (present(factor)) then
840 0 : factor_local = factor
841 : else
842 0 : factor_local = 1.0_r8
843 : end if
844 :
845 0 : sum_species = 1.0_r8 ! all dry air species sum to 1
846 0 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
847 0 : itrac = idx_local(qdx)
848 0 : sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_local(:,:))
849 : end do
850 :
851 0 : if (dry_air_species_num == 0) then
852 0 : sum_cp = thermodynamic_active_species_cp(0)
853 0 : else if (present(cpdry)) then
854 : !
855 : ! if cpdry is known don't recompute
856 : !
857 0 : sum_cp = cpdry
858 : else
859 0 : call get_cp_dry(tracer, idx_local, sum_cp, fact=factor_local)
860 : end if
861 0 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
862 0 : itrac = idx_local(qdx)
863 : sum_cp(:,:) = sum_cp(:,:)+ &
864 0 : thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac)* factor_local(:,:)
865 : end do
866 0 : if (inv_cp) then
867 0 : cp = sum_species / sum_cp
868 : else
869 0 : cp = sum_cp / sum_species
870 : end if
871 0 : end subroutine get_cp_1hd
872 :
873 : !===========================================================================
874 :
875 0 : subroutine get_cp_2hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry)
876 : ! Version of get_cp for arrays that have a second horizontal index
877 : use cam_abortutils, only: endrun
878 : use string_utils, only: int2str
879 :
880 : ! Dummy arguments
881 : ! tracer: Tracer array
882 : !
883 : real(r8), intent(in) :: tracer(:,:,:,:)
884 : ! inv_cp: output inverse cp instead of cp
885 : logical, intent(in) :: inv_cp
886 : real(r8), intent(out) :: cp(:,:,:)
887 : real(r8), optional, intent(in) :: factor(:,:,:)
888 : real(r8), optional, intent(in) :: cpdry(:,:,:)
889 :
890 : ! active_species_idx_dycore: array of indicies for index of
891 : ! thermodynamic active species in dycore tracer array
892 : ! (if different from physics index)
893 : integer, optional, intent(in) :: active_species_idx_dycore(:)
894 :
895 : ! Local variables
896 : integer :: jdx
897 : integer :: idx_local(thermodynamic_active_species_num)
898 : character(len=*), parameter :: subname = 'get_cp_2hd: '
899 :
900 0 : do jdx = 1, SIZE(cp, 2)
901 0 : if (present(factor).and.present(cpdry)) then
902 : call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
903 0 : factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:))
904 0 : else if (present(factor)) then
905 : call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
906 0 : factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore)
907 0 : else if (present(cpdry)) then
908 : call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
909 0 : active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:))
910 : else
911 : call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
912 0 : active_species_idx_dycore=active_species_idx_dycore)
913 : end if
914 : end do
915 :
916 0 : end subroutine get_cp_2hd
917 :
918 : !===========================================================================
919 :
920 : !***************************************************************************
921 : !
922 : ! get_R_dry: Compute generalized dry air gas constant R
923 : !
924 : !***************************************************************************
925 : !
926 876672 : subroutine get_R_dry_1hd(tracer, active_species_idx_dycore, R_dry, fact)
927 : use physconst, only: rair
928 :
929 : ! tracer: tracer array
930 : real(r8), intent(in) :: tracer(:, :, :)
931 : ! active_species_idx_dycore: index of active species in tracer
932 : integer, intent(in) :: active_species_idx_dycore(:)
933 : ! R_dry: dry air R
934 : real(r8), intent(out) :: R_dry(:, :)
935 : ! fact: optional factor for converting tracer to dry mixing ratio
936 : real(r8), optional, intent(in) :: fact(:, :)
937 :
938 : ! Local variables
939 : integer :: idx, kdx, m_cnst, qdx
940 1753344 : real(r8) :: factor(SIZE(tracer, 1), SIZE(tracer, 2))
941 1753344 : real(r8) :: residual(SIZE(R_dry, 1), SIZE(R_dry, 2))
942 : real(r8) :: mmr
943 :
944 876672 : if (dry_air_species_num == 0) then
945 : !
946 : ! dry air not species dependent
947 : !
948 0 : R_dry = rair
949 : else
950 876672 : if (present(fact)) then
951 0 : factor = fact(:,:)
952 : else
953 1476686592 : factor = 1.0_r8
954 : end if
955 :
956 1476686592 : R_dry = 0.0_r8
957 1476686592 : residual = 1.0_r8
958 3506688 : do qdx = 1, dry_air_species_num
959 2630016 : m_cnst = active_species_idx_dycore(qdx)
960 328111488 : do kdx = 1, SIZE(R_dry, 2)
961 4430059776 : do idx = 1, SIZE(R_dry, 1)
962 4102824960 : mmr = tracer(idx, kdx, m_cnst) * factor(idx, kdx)
963 4102824960 : R_dry(idx, kdx) = R_dry(idx, kdx) + &
964 4102824960 : (thermodynamic_active_species_R(qdx) * mmr)
965 4427429760 : residual(idx, kdx) = residual(idx, kdx) - mmr
966 : end do
967 : end do
968 : end do
969 : !
970 : ! N2 derived from the others
971 : !
972 876672 : qdx = 0
973 109078272 : do kdx = 1, SIZE(R_dry, 2)
974 1476686592 : do idx = 1, SIZE(R_dry, 1)
975 2735216640 : R_dry(idx, kdx) = R_dry(idx, kdx) + &
976 4211026560 : (thermodynamic_active_species_R(qdx) * residual(idx, kdx))
977 : end do
978 : end do
979 : end if
980 876672 : end subroutine get_R_dry_1hd
981 :
982 : !===========================================================================
983 :
984 64512 : subroutine get_R_dry_2hd(tracer, active_species_idx_dycore, R_dry, fact)
985 : ! Version of get_R_dry for arrays that have a second horizontal index
986 :
987 : ! tracer: tracer array
988 : real(r8), intent(in) :: tracer(:, :, :, :)
989 : ! active_species_idx_dycore: index of active species in tracer
990 : integer, intent(in) :: active_species_idx_dycore(:)
991 : ! R_dry: dry air R
992 : real(r8), intent(out) :: R_dry(:, :, :)
993 : ! fact: optional factor for converting tracer to dry mixing ratio
994 : real(r8), optional, intent(in) :: fact(:, :, :)
995 :
996 : ! Local variable
997 : integer :: jdx
998 :
999 258048 : do jdx = 1, SIZE(tracer, 2)
1000 258048 : if (present(fact)) then
1001 : call get_R_dry(tracer(:, jdx, :, :), active_species_idx_dycore, &
1002 0 : R_dry(:, jdx, :), fact=fact(:, jdx, :))
1003 : else
1004 : call get_R_dry(tracer(:, jdx, :, :), active_species_idx_dycore, &
1005 193536 : R_dry(:, jdx, :))
1006 : end if
1007 : end do
1008 :
1009 64512 : end subroutine get_R_dry_2hd
1010 :
1011 : !===========================================================================
1012 : !
1013 : !***************************************************************************
1014 : !
1015 : ! get_R: Compute generalized R
1016 : ! This code (both 1hd and 2hd) is currently unused and untested
1017 : !
1018 : !***************************************************************************
1019 : !
1020 0 : subroutine get_R_1hd(tracer, active_species_idx, R, fact, Rdry)
1021 : use cam_abortutils, only: endrun
1022 : use string_utils, only: int2str
1023 : use physconst, only: rair
1024 :
1025 : ! Dummy arguments
1026 : ! tracer: !tracer array
1027 : real(r8), intent(in) :: tracer(:, :, :)
1028 : ! active_species_idx: index of active species in tracer
1029 : integer, intent(in) :: active_species_idx(:)
1030 : ! R: generalized gas constant
1031 : real(r8), intent(out) :: R(:, :)
1032 : ! fact: optional factor for converting tracer to dry mixing ratio
1033 : real(r8), optional, intent(in) :: fact(:, :)
1034 : real(r8), optional, intent(in) :: Rdry(:, :)
1035 :
1036 : ! Local variables
1037 : integer :: qdx, itrac
1038 0 : real(r8) :: factor(SIZE(tracer, 1), SIZE(tracer, 2))
1039 0 : real(r8) :: sum_species(SIZE(R, 1), SIZE(R, 2))
1040 0 : integer :: idx_local(thermodynamic_active_species_num)
1041 :
1042 : character(len=*), parameter :: subname = 'get_R_1hd: '
1043 :
1044 0 : if (present(fact)) then
1045 0 : if (SIZE(fact, 1) /= SIZE(factor, 1)) then
1046 : call endrun(subname//"SIZE mismatch in dimension 1 "// &
1047 0 : int2str(SIZE(fact, 1))//' /= '//int2str(SIZE(factor, 1)))
1048 : end if
1049 0 : if (SIZE(fact, 2) /= SIZE(factor, 2)) then
1050 : call endrun(subname//"SIZE mismatch in dimension 2 "// &
1051 0 : int2str(SIZE(fact, 2))//' /= '//int2str(SIZE(factor, 2)))
1052 : end if
1053 0 : factor = fact(:,:)
1054 : else
1055 0 : factor = 1.0_r8
1056 : end if
1057 :
1058 0 : if (dry_air_species_num == 0) then
1059 0 : R = rair
1060 0 : else if (present(Rdry)) then
1061 0 : R = Rdry
1062 : else
1063 0 : call get_R_dry(tracer, active_species_idx, R, fact=factor)
1064 : end if
1065 :
1066 0 : idx_local = active_species_idx
1067 0 : sum_species = 1.0_r8 ! all dry air species sum to 1
1068 0 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
1069 0 : itrac = idx_local(qdx)
1070 : sum_species(:,:) = sum_species(:,:) + &
1071 0 : (tracer(:,:,itrac) * factor(:,:))
1072 : end do
1073 0 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
1074 0 : itrac = idx_local(qdx)
1075 : R(:,:) = R(:,:) + &
1076 0 : (thermodynamic_active_species_R(qdx) * tracer(:,:,itrac) * &
1077 0 : factor(:,:))
1078 : end do
1079 0 : R = R / sum_species
1080 0 : end subroutine get_R_1hd
1081 :
1082 : !===========================================================================
1083 :
1084 0 : subroutine get_R_2hd(tracer, active_species_idx, R, fact)
1085 :
1086 : ! Dummy arguments
1087 : ! tracer: !tracer array
1088 : real(r8), intent(in) :: tracer(:, :, :, :)
1089 : ! active_species_idx: index of active species in tracer
1090 : integer, intent(in) :: active_species_idx(:)
1091 : ! R: generalized gas constant
1092 : real(r8), intent(out) :: R(:, :, :)
1093 : ! fact: optional factor for converting tracer to dry mixing ratio
1094 : real(r8), optional, intent(in) :: fact(:, :, :)
1095 :
1096 : ! Local variable
1097 : integer :: jdx
1098 :
1099 0 : do jdx = 1, SIZE(tracer, 2)
1100 0 : if (present(fact)) then
1101 : call get_R(tracer(:, jdx, :, :), active_species_idx, &
1102 0 : R(:, jdx, :), fact=fact(:, jdx, :))
1103 : else
1104 : call get_R(tracer(:, jdx, :, :), active_species_idx, &
1105 0 : R(:, jdx, :))
1106 : end if
1107 : end do
1108 :
1109 0 : end subroutine get_R_2hd
1110 :
1111 : !===========================================================================
1112 :
1113 : !*************************************************************************************************************************
1114 : !
1115 : ! compute molecular weight dry air
1116 : !
1117 : !*************************************************************************************************************************
1118 : !
1119 1366272 : subroutine get_mbarv_1hd(tracer, active_species_idx, mbarv_in, fact)
1120 : use physconst, only: mwdry
1121 : real(r8), intent(in) :: tracer(:,:,:) !tracer array
1122 : integer, intent(in) :: active_species_idx(:) !index of active species in tracer
1123 : real(r8), intent(out) :: mbarv_in(:,:) !molecular weight of dry air
1124 : real(r8), optional, intent(in) :: fact(:,:) !factor for converting tracer to dry mixing ratio
1125 :
1126 : integer :: idx, kdx, m_cnst, qdx
1127 2732544 : real(r8):: factor(SIZE(mbarv_in, 1), SIZE(mbarv_in, 2))
1128 2732544 : real(r8):: residual(SIZE(tracer, 1), SIZE(mbarv_in, 2))
1129 : real(r8):: mm
1130 : !
1131 : ! dry air not species dependent
1132 : !
1133 1366272 : if (dry_air_species_num==0) then
1134 0 : mbarv_in = mwdry
1135 : else
1136 1366272 : if (present(fact)) then
1137 1155866112 : factor(:,:) = fact(:,:)
1138 : else
1139 1155182976 : factor(:,:) = 1.0_r8
1140 : endif
1141 :
1142 2310365952 : mbarv_in = 0.0_r8
1143 2310365952 : residual = 1.0_r8
1144 5465088 : do qdx = 1, dry_air_species_num
1145 4098816 : m_cnst = active_species_idx(qdx)
1146 538311168 : do kdx = 1, SIZE(mbarv_in, 2)
1147 6931097856 : do idx = 1, SIZE(mbarv_in, 1)
1148 6394152960 : mm = tracer(idx, kdx, m_cnst) * factor(idx, kdx)
1149 6394152960 : mbarv_in(idx, kdx) = mbarv_in(idx, kdx) + thermodynamic_active_species_mwi(qdx) * mm
1150 6926999040 : residual(idx, kdx) = residual(idx, kdx) - mm
1151 : end do
1152 : end do
1153 : end do
1154 1366272 : qdx = 0 ! N2
1155 178981632 : do kdx = 1, SIZE(mbarv_in, 2)
1156 2310365952 : do idx = 1, SIZE(mbarv_in, 1)
1157 2308999680 : mbarv_in(idx, kdx) = mbarv_in(idx, kdx) + thermodynamic_active_species_mwi(qdx) * residual(idx, kdx)
1158 : end do
1159 : end do
1160 2310365952 : mbarv_in(:,:) = 1.0_r8 / mbarv_in(:,:)
1161 : end if
1162 1366272 : end subroutine get_mbarv_1hd
1163 :
1164 : !===========================================================================
1165 :
1166 6912 : subroutine air_species_info(name, index, molec_weight, caller)
1167 : use cam_abortutils, only: endrun
1168 : use cam_logfile, only: iulog
1169 : use constituents, only: cnst_get_ind, cnst_mw
1170 : ! Find the constituent index of <name> and return it in
1171 : ! <index>. Return the constituent molecular weight in
1172 : ! <molec_weight>
1173 :
1174 : ! Dummy arguments
1175 : character(len=*), intent(in) :: name
1176 : integer, intent(out) :: index
1177 : real(r8), intent(out) :: molec_weight
1178 : character(len=*), optional, intent(in) :: caller
1179 : ! Local parameter
1180 : character(len=*), parameter :: subname = 'air_species_info: '
1181 :
1182 6912 : call cnst_get_ind(trim(name), index, abort=.false.)
1183 6912 : if (index < 1) then
1184 0 : if (present(caller)) then
1185 0 : write(iulog, *) trim(caller), ": air component not found, '", &
1186 0 : trim(name), "'"
1187 : call endrun(trim(caller)//": air component not found, '"// &
1188 0 : trim(name)//"'")
1189 : else
1190 0 : write(iulog, *) subname, "air component not found, '", &
1191 0 : trim(name), "'"
1192 : call endrun(subname//"air component not found, '"// &
1193 0 : trim(name)//"'")
1194 : end if
1195 : else
1196 6912 : molec_weight = cnst_mw(index)
1197 : end if
1198 :
1199 6912 : end subroutine air_species_info
1200 :
1201 :
1202 : end module air_composition
|
