Line data Source code
1 : ! cam_thermo module provides interfaces to compute thermodynamic quantities
2 : module cam_thermo
3 :
4 : use shr_kind_mod, only: r8 => shr_kind_r8
5 : use cam_abortutils, only: endrun
6 : use air_composition, only: thermodynamic_active_species_num
7 : use air_composition, only: thermodynamic_active_species_idx
8 : use air_composition, only: thermodynamic_active_species_idx_dycore
9 : use air_composition, only: thermodynamic_active_species_cp
10 : use air_composition, only: thermodynamic_active_species_R
11 : use air_composition, only: thermodynamic_active_species_mwi
12 : use air_composition, only: thermodynamic_active_species_kv
13 : use air_composition, only: thermodynamic_active_species_kc
14 : use air_composition, only: thermodynamic_active_species_liq_num
15 : use air_composition, only: thermodynamic_active_species_ice_num
16 : use air_composition, only: thermodynamic_active_species_liq_idx
17 : use air_composition, only: thermodynamic_active_species_liq_idx_dycore
18 : use air_composition, only: thermodynamic_active_species_ice_idx
19 : use air_composition, only: thermodynamic_active_species_ice_idx_dycore
20 : use air_composition, only: dry_air_species_num
21 : use air_composition, only: enthalpy_reference_state
22 : use air_composition, only: mmro2, mmrn2, o2_mwi, n2_mwi, mbar
23 :
24 : implicit none
25 : private
26 : save
27 :
28 : ! subroutines to compute thermodynamic quantities
29 : !
30 : ! See Lauritzen et al. (2018) for formulae
31 : ! DOI: 10.1029/2017MS001257
32 : ! https://opensky.ucar.edu/islandora/object/articles:21929
33 :
34 : ! cam_thermo_init: Initialize constituent dependent properties
35 : public :: cam_thermo_init
36 : ! cam_thermo_dry_air_update: Update dry air composition dependent properties
37 : public :: cam_thermo_dry_air_update
38 : ! cam_thermo_water_update: Update water dependent properties
39 : public :: cam_thermo_water_update
40 : ! get_enthalpy: enthalpy quantity = dp*cp*T
41 : public :: get_enthalpy
42 : ! get_virtual_temp: virtual temperature
43 : public :: get_virtual_temp
44 : ! get_sum_species: sum of thermodynamically active species:
45 : ! Note: dp = dp_dry * sum_species
46 : public :: get_sum_species
47 : ! get_virtual_theta: virtual potential temperature
48 : public :: get_virtual_theta
49 : ! cam_thermo_calc_kappav: update species dependent kappa for FV dycore
50 : public :: cam_thermo_calc_kappav
51 : ! get_dp: pressure level thickness from dry dp and dry mixing ratios
52 : public :: get_dp
53 : ! get_pmid_from_dp: full level pressure from dp (approximation depends on dycore)
54 : public :: get_pmid_from_dp
55 : ! get_ps: surface pressure
56 : public :: get_ps
57 : ! get_gz: geopotential
58 : public :: get_gz
59 : ! get_Richardson_number: Richardson number at layer interfaces
60 : public :: get_Richardson_number
61 : ! get_kappa_dry: (generalized) dry kappa = R_dry/cp_dry
62 : public :: get_kappa_dry
63 : ! get_dp_ref: reference pressure layer thickness (include topography)
64 : public :: get_dp_ref
65 : ! get_molecular_diff_coef: molecular diffusion and thermal conductivity
66 : public :: get_molecular_diff_coef
67 : ! get_molecular_diff_coef_reference: reference vertical profile of density,
68 : ! molecular diffusion and thermal conductivity
69 : public :: get_molecular_diff_coef_reference
70 : ! get_rho_dry: dry density from temperature (temp) and
71 : ! pressure (dp_dry and tracer)
72 : public :: get_rho_dry
73 : ! get_exner: Exner pressure
74 : public :: get_exner
75 : ! get_hydrostatic_energy: Vertically integrated total energy
76 : public :: get_hydrostatic_energy
77 :
78 : ! Public variables
79 : ! mixing_ratio options
80 : integer, public, parameter :: DRY_MIXING_RATIO = 1
81 : integer, public, parameter :: MASS_MIXING_RATIO = 2
82 : !--------------- Variables below here are for WACCM-X ---------------------
83 : ! kmvis: molecular viscosity kg/m/s
84 : real(r8), public, protected, allocatable :: kmvis(:,:,:)
85 : ! kmcnd: molecular conductivity J/m/s/K
86 : real(r8), public, protected, allocatable :: kmcnd(:,:,:)
87 :
88 : !------------- Variables for consistent themodynamics --------------------
89 : !
90 :
91 : !
92 : ! Interfaces for public routines
93 : interface get_gz
94 : ! get_gz_geopotential (with dp_dry, ptop, temp, and phis as input)
95 : module procedure get_gz_from_dp_dry_ptop_temp_1hd
96 : ! get_gz_given_dp_Tv_Rdry: geopotential (with dp,dry R and Tv as input)
97 : module procedure get_gz_given_dp_Tv_Rdry_1hd
98 : module procedure get_gz_given_dp_Tv_Rdry_2hd
99 : end interface get_gz
100 :
101 : interface get_enthalpy
102 : module procedure get_enthalpy_1hd
103 : module procedure get_enthalpy_2hd
104 : end interface get_enthalpy
105 :
106 : interface get_virtual_temp
107 : module procedure get_virtual_temp_1hd
108 : module procedure get_virtual_temp_2hd
109 : end interface get_virtual_temp
110 :
111 : interface get_sum_species
112 : module procedure get_sum_species_1hd
113 : module procedure get_sum_species_2hd
114 : end interface get_sum_species
115 :
116 : interface get_dp
117 : module procedure get_dp_1hd
118 : module procedure get_dp_2hd
119 : end interface get_dp
120 :
121 : interface get_pmid_from_dp
122 : module procedure get_pmid_from_dpdry_1hd
123 : module procedure get_pmid_from_dp_1hd
124 : end interface get_pmid_from_dp
125 :
126 : interface get_exner
127 : module procedure get_exner_1hd
128 : end interface get_exner
129 :
130 : interface get_virtual_theta
131 : module procedure get_virtual_theta_1hd
132 : end interface get_virtual_theta
133 :
134 : interface get_Richardson_number
135 : module procedure get_Richardson_number_1hd
136 : end interface get_Richardson_number
137 :
138 : interface get_ps
139 : module procedure get_ps_1hd
140 : module procedure get_ps_2hd
141 : end interface get_ps
142 :
143 : interface get_kappa_dry
144 : module procedure get_kappa_dry_1hd
145 : module procedure get_kappa_dry_2hd
146 : end interface get_kappa_dry
147 :
148 : interface get_dp_ref
149 : module procedure get_dp_ref_1hd
150 : module procedure get_dp_ref_2hd
151 : end interface get_dp_ref
152 :
153 : interface get_rho_dry
154 : module procedure get_rho_dry_1hd
155 : module procedure get_rho_dry_2hd
156 : end interface get_rho_dry
157 :
158 : interface get_molecular_diff_coef
159 : module procedure get_molecular_diff_coef_1hd
160 : module procedure get_molecular_diff_coef_2hd
161 : end interface get_molecular_diff_coef
162 :
163 : interface cam_thermo_calc_kappav
164 : ! Since this routine is currently only used by the FV dycore,
165 : ! a 1-d interface is not needed (but can easily be added)
166 : module procedure cam_thermo_calc_kappav_2hd
167 : end interface cam_thermo_calc_kappav
168 :
169 : interface get_hydrostatic_energy
170 : module procedure get_hydrostatic_energy_1hd
171 : ! This routine is currently only called from the physics so a
172 : ! 2-d interface is not needed (but can easily be added)
173 : end interface get_hydrostatic_energy
174 :
175 : integer, public, parameter :: thermo_budget_num_vars = 10
176 : integer, public, parameter :: wvidx = 1
177 : integer, public, parameter :: wlidx = 2
178 : integer, public, parameter :: wiidx = 3
179 : integer, public, parameter :: seidx = 4 ! enthalpy or internal energy (W/m2) index
180 : integer, public, parameter :: poidx = 5 ! surface potential or potential energy index
181 : integer, public, parameter :: keidx = 6 ! kinetic energy index
182 : integer, public, parameter :: mridx = 7
183 : integer, public, parameter :: moidx = 8
184 : integer, public, parameter :: ttidx = 9
185 : integer, public, parameter :: teidx = 10
186 : character (len = 2) ,public, dimension(thermo_budget_num_vars) :: thermo_budget_vars = &
187 : (/"WV" ,"WL" ,"WI" ,"SE" ,"PO" ,"KE" ,"MR" ,"MO" ,"TT" ,"TE" /)
188 : character (len = 46) ,public, dimension(thermo_budget_num_vars) :: thermo_budget_vars_descriptor = (/&
189 : "Total column water vapor ",&
190 : "Total column liquid water ",&
191 : "Total column frozen water ",&
192 : "Total column enthalpy or internal energy ",&
193 : "Total column srf potential or potential energy",&
194 : "Total column kinetic energy ",&
195 : "Total column wind axial angular momentum ",&
196 : "Total column mass axial angular momentum ",&
197 : "Total column test_tracer ",&
198 : "Total column energy (ke + se + po) "/)
199 :
200 : character (len = 14), public, dimension(thermo_budget_num_vars) :: &
201 : thermo_budget_vars_unit = (/&
202 : "kg/m2 ","kg/m2 ","kg/m2 ","J/m2 ",&
203 : "J/m2 ","J/m2 ","kg*m2/s*rad2 ","kg*m2/s*rad2 ",&
204 : "kg/m2 ","J/m2 "/)
205 : logical ,public, dimension(thermo_budget_num_vars) :: thermo_budget_vars_massv = (/&
206 : .true.,.true.,.true.,.false.,.false.,.false.,.false.,.false.,.true.,.false./)
207 : CONTAINS
208 :
209 : !===========================================================================
210 :
211 1536 : subroutine cam_thermo_init()
212 : use shr_infnan_mod, only: assignment(=), shr_infnan_qnan
213 : use ppgrid, only: pcols, pver, pverp, begchunk, endchunk
214 :
215 : integer :: ierr
216 : character(len=*), parameter :: subname = "cam_thermo_init"
217 : character(len=*), parameter :: errstr = subname//": failed to allocate "
218 :
219 : !------------------------------------------------------------------------
220 : ! Allocate constituent dependent properties
221 : !------------------------------------------------------------------------
222 4608 : allocate(kmvis(pcols,pverp,begchunk:endchunk), stat=ierr)
223 1536 : if (ierr /= 0) then
224 0 : call endrun(errstr//"kmvis")
225 : end if
226 4608 : allocate(kmcnd(pcols,pverp,begchunk:endchunk), stat=ierr)
227 1536 : if (ierr /= 0) then
228 0 : call endrun(errstr//"kmcnd")
229 : end if
230 :
231 : !------------------------------------------------------------------------
232 : ! Initialize constituent dependent properties
233 : !------------------------------------------------------------------------
234 1536 : kmvis(:pcols, :pver, begchunk:endchunk) = shr_infnan_qnan
235 1536 : kmcnd(:pcols, :pver, begchunk:endchunk) = shr_infnan_qnan
236 :
237 1536 : end subroutine cam_thermo_init
238 : !
239 : !***************************************************************************
240 : !
241 : ! cam_thermo_dry_air_update: update dry air species dependent constants for physics
242 : !
243 : !***************************************************************************
244 : !
245 0 : subroutine cam_thermo_dry_air_update(mmr, T, lchnk, ncol, to_dry_factor)
246 : use air_composition, only: dry_air_composition_update
247 : use string_utils, only: int2str
248 : !------------------------------Arguments----------------------------------
249 : !(mmr = dry mixing ratio, if not use to_dry_factor to convert)
250 : real(r8), intent(in) :: mmr(:,:,:) ! constituents array
251 : real(r8), intent(in) :: T(:,:) ! temperature
252 : integer, intent(in) :: lchnk ! Chunk number
253 : integer, intent(in) :: ncol ! number of columns
254 : real(r8), optional, intent(in) :: to_dry_factor(:,:)!if mmr moist convert
255 : !
256 : !---------------------------Local storage-------------------------------
257 0 : real(r8):: sponge_factor(SIZE(mmr, 2))
258 : character(len=*), parameter :: subname = 'cam_thermo_update: '
259 :
260 0 : if (present(to_dry_factor)) then
261 0 : if (SIZE(to_dry_factor, 1) /= ncol) then
262 0 : call endrun(subname//'DIM 1 of to_dry_factor is'//int2str(SIZE(to_dry_factor,1))//'but should be'//int2str(ncol))
263 : end if
264 : end if
265 :
266 0 : sponge_factor = 1.0_r8
267 0 : call dry_air_composition_update(mmr, lchnk, ncol, to_dry_factor=to_dry_factor)
268 0 : call get_molecular_diff_coef(T(:ncol,:), .true., sponge_factor, kmvis(:ncol,:,lchnk), &
269 0 : kmcnd(:ncol,:,lchnk), tracer=mmr(:ncol,:,:), fact=to_dry_factor, &
270 0 : active_species_idx_dycore=thermodynamic_active_species_idx)
271 0 : end subroutine cam_thermo_dry_air_update
272 : !
273 : !***************************************************************************
274 : !
275 : ! cam_thermo_water+update: update water species dependent constants for physics
276 : !
277 : !***************************************************************************
278 : !
279 2984544 : subroutine cam_thermo_water_update(mmr, lchnk, ncol, vcoord, to_dry_factor)
280 : use air_composition, only: water_composition_update
281 : !-----------------------------------------------------------------------
282 : ! Update the physics "constants" that vary
283 : !-------------------------------------------------------------------------
284 :
285 : !------------------------------Arguments----------------------------------
286 :
287 : real(r8), intent(in) :: mmr(:,:,:) ! constituents array
288 : integer, intent(in) :: lchnk ! Chunk number
289 : integer, intent(in) :: ncol ! number of columns
290 : integer, intent(in) :: vcoord
291 : real(r8), optional, intent(in) :: to_dry_factor(:,:)
292 : !
293 : logical :: lcp
294 :
295 4479912 : call water_composition_update(mmr, lchnk, ncol, vcoord, to_dry_factor=to_dry_factor)
296 2984544 : end subroutine cam_thermo_water_update
297 :
298 : !===========================================================================
299 :
300 : !
301 : !***********************************************************************
302 : !
303 : ! Compute enthalpy = cp*T*dp, where dp is pressure level thickness,
304 : ! cp is generalized cp and T temperature
305 : !
306 : ! Note: tracer is in units of m*dp_dry ("mass")
307 : !
308 : !***********************************************************************
309 : !
310 41558400 : subroutine get_enthalpy_1hd(tracer_mass, temp, dp_dry, &
311 41558400 : enthalpy, active_species_idx_dycore)
312 : use air_composition, only: dry_air_species_num, get_cp_dry
313 : ! Dummy arguments
314 : ! tracer_mass: tracer array (mass weighted)
315 : real(r8), intent(in) :: tracer_mass(:,:,:)
316 : ! temp: temperature
317 : real(r8), intent(in) :: temp(:,:)
318 : ! dp_dry: dry presure level thickness
319 : real(r8), intent(in) :: dp_dry(:,:)
320 : ! enthalpy: enthalpy in each column: sum cp*T*dp
321 : real(r8), intent(out) :: enthalpy(:,:)
322 : !
323 : ! active_species_idx_dycore:
324 : ! array of indicies for index of thermodynamic active species in
325 : ! dycore tracer array (if different from physics index)
326 : !
327 : integer, optional, intent(in) :: active_species_idx_dycore(:)
328 :
329 : ! Local vars
330 : integer :: qdx, itrac
331 : character(len=*), parameter :: subname = 'get_enthalpy: '
332 :
333 : !
334 : ! "mass-weighted" cp (dp must be dry)
335 : !
336 41558400 : if (dry_air_species_num == 0) then
337 0 : enthalpy(:,:) = thermodynamic_active_species_cp(0) * &
338 19366214400 : dp_dry(:,:)
339 : else
340 0 : if (present(active_species_idx_dycore)) then
341 : call get_cp_dry(tracer_mass, active_species_idx_dycore, &
342 0 : enthalpy, fact=1.0_r8/dp_dry(:,:))
343 : else
344 : call get_cp_dry(tracer_mass, thermodynamic_active_species_idx, &
345 0 : enthalpy, fact=1.0_r8/dp_dry(:,:))
346 : end if
347 0 : enthalpy(:,:) = enthalpy(:,:) * dp_dry(:,:)
348 : end if
349 : !
350 : ! tracer is in units of m*dp ("mass"), where:
351 : ! m is the dry mixing ratio
352 : ! dp is the dry pressure level thickness
353 : !
354 290908800 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
355 249350400 : if (present(active_species_idx_dycore)) then
356 249350400 : itrac = active_species_idx_dycore(qdx)
357 : else
358 0 : itrac = thermodynamic_active_species_idx(qdx)
359 : end if
360 : enthalpy(:,:) = enthalpy(:,:) + &
361 >11623*10^7 : (thermodynamic_active_species_cp(qdx) * tracer_mass(:,:,itrac))
362 : end do
363 19407772800 : enthalpy(:,:) = enthalpy(:,:) * temp(:,:)
364 :
365 41558400 : end subroutine get_enthalpy_1hd
366 :
367 : !===========================================================================
368 :
369 10389600 : subroutine get_enthalpy_2hd(tracer_mass, temp, dp_dry, &
370 10389600 : enthalpy, active_species_idx_dycore)
371 : ! Dummy arguments
372 : ! tracer_mass: tracer array (mass weighted)
373 : real(r8), intent(in) :: tracer_mass(:,:,:,:)
374 : ! temp: temperature
375 : real(r8), intent(in) :: temp(:,:,:)
376 : ! dp_dry: dry presure level thickness
377 : real(r8), intent(in) :: dp_dry(:,:,:)
378 : ! enthalpy: enthalpy in each column: sum cp*T*dp
379 : real(r8), intent(out) :: enthalpy(:,:,:)
380 : !
381 : ! active_species_idx_dycore:
382 : ! array of indicies for index of thermodynamic active species in
383 : ! dycore tracer array (if different from physics index)
384 : !
385 : integer, optional, intent(in) :: active_species_idx_dycore(:)
386 :
387 : ! Local variables
388 : integer :: jdx
389 : character(len=*), parameter :: subname = 'get_enthalpy_2hd: '
390 :
391 51948000 : do jdx = 1, SIZE(tracer_mass, 2)
392 : call get_enthalpy(tracer_mass(:, jdx, :, :), temp(:, jdx, :), &
393 : dp_dry(:, jdx, :), enthalpy(:, jdx, :), &
394 51948000 : active_species_idx_dycore=active_species_idx_dycore)
395 : end do
396 :
397 10389600 : end subroutine get_enthalpy_2hd
398 :
399 : !===========================================================================
400 :
401 : !**************************************************************************
402 : !
403 : ! get_virtual_temp: Compute virtual temperature T_v
404 : !
405 : ! tracer is in units of dry mixing ratio unless optional argument
406 : ! dp_dry is present in which case tracer is in units of "mass" (=m*dp)
407 : !
408 : ! If temperature is not supplied then just return factor that T
409 : ! needs to be multiplied by to get T_v
410 : !
411 : !**************************************************************************
412 : !
413 623376000 : subroutine get_virtual_temp_1hd(tracer, T_v, temp, dp_dry, sum_q, &
414 311688000 : active_species_idx_dycore)
415 : use cam_abortutils, only: endrun
416 : use string_utils, only: int2str
417 : use air_composition, only: dry_air_species_num, get_R_dry
418 :
419 : ! Dummy Arguments
420 : ! tracer: tracer array
421 : real(r8), intent(in) :: tracer(:, :, :)
422 : ! T_v: virtual temperature
423 : real(r8), intent(out) :: T_v(:, :)
424 : ! temp: temperature
425 : real(r8), optional, intent(in) :: temp(:, :)
426 : ! dp_dry: dry pressure level thickness
427 : real(r8), optional, intent(in) :: dp_dry(:, :)
428 : ! sum_q: sum tracer
429 : real(r8), optional, intent(out) :: sum_q(:, :)
430 : !
431 : ! array of indicies for index of thermodynamic active species in
432 : ! dycore tracer array (if different from physics index)
433 : !
434 : integer, optional, intent(in) :: active_species_idx_dycore(:)
435 :
436 : ! Local Variables
437 : integer :: itrac, qdx
438 623376000 : real(r8) :: sum_species(SIZE(tracer, 1), SIZE(tracer, 2))
439 623376000 : real(r8) :: factor(SIZE(tracer, 1), SIZE(tracer, 2))
440 623376000 : real(r8) :: Rd(SIZE(tracer, 1), SIZE(tracer, 2))
441 623376000 : integer :: idx_local(thermodynamic_active_species_num)
442 : character(len=*), parameter :: subname = 'get_virtual_temp_1hd: '
443 :
444 311688000 : if (present(active_species_idx_dycore)) then
445 311688000 : if (SIZE(active_species_idx_dycore) /= &
446 : thermodynamic_active_species_num) then
447 : call endrun(subname//"SIZE mismatch "// &
448 : int2str(SIZE(active_species_idx_dycore))//' /= '// &
449 0 : int2str(thermodynamic_active_species_num))
450 : end if
451 2493504000 : idx_local = active_species_idx_dycore
452 : else
453 0 : idx_local = thermodynamic_active_species_idx
454 : end if
455 :
456 623376000 : call get_sum_species(tracer, idx_local, sum_species, dp_dry=dp_dry, factor=factor)
457 :
458 311688000 : call get_R_dry(tracer, idx_local, Rd, fact=factor)
459 >14555*10^7 : t_v(:, :) = Rd(:, :)
460 2181816000 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
461 1870128000 : itrac = idx_local(qdx)
462 0 : t_v(:, :) = t_v(:, :) + (thermodynamic_active_species_R(qdx) * &
463 >87179*10^7 : tracer(:, :, itrac) * factor(:, :))
464 : end do
465 311688000 : if (present(temp)) then
466 >14555*10^7 : t_v(:, :) = t_v(:, :) * temp(:, :) / (Rd(:, :) * sum_species)
467 : else
468 0 : t_v(:, :) = t_v(:, :) / (Rd(:, :) * sum_species)
469 : end if
470 311688000 : if (present(sum_q)) then
471 >14555*10^7 : sum_q = sum_species
472 : end if
473 :
474 311688000 : end subroutine get_virtual_temp_1hd
475 :
476 : !===========================================================================
477 :
478 77922000 : subroutine get_virtual_temp_2hd(tracer, T_v, temp, dp_dry, sum_q, &
479 77922000 : active_species_idx_dycore)
480 :
481 : ! Dummy Arguments
482 : ! tracer: tracer array
483 : real(r8), intent(in) :: tracer(:, :, :, :)
484 : ! T_v: virtual temperature
485 : real(r8), intent(out) :: T_v(:, :, :)
486 : ! temp: temperature
487 : real(r8), optional, intent(in) :: temp(:, :, :)
488 : ! dp_dry: dry pressure level thickness
489 : real(r8), optional, intent(in) :: dp_dry(:, :, :)
490 : ! sum_q: sum tracer
491 : real(r8), optional, intent(out) :: sum_q(:, :, :)
492 : !
493 : ! array of indicies for index of thermodynamic active species in
494 : ! dycore tracer array (if different from physics index)
495 : !
496 : integer, optional, intent(in) :: active_species_idx_dycore(:)
497 :
498 : ! Local vars
499 : integer :: jdx
500 : character(len=*), parameter :: subname = 'get_virtual_temp_2hd: '
501 :
502 : ! Rather than do a bunch of copying into temp variables, do the
503 : ! combinatorics
504 389610000 : do jdx = 1, SIZE(tracer, 2)
505 389610000 : if (present(temp) .and. present(dp_dry) .and. present(sum_q)) then
506 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
507 : temp=temp(:, jdx, :), dp_dry=dp_dry(:, jdx, :), &
508 : sum_q=sum_q(:, jdx, :), &
509 0 : active_species_idx_dycore=active_species_idx_dycore)
510 311688000 : else if (present(temp) .and. present(dp_dry)) then
511 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
512 : temp=temp(:, jdx, :), dp_dry=dp_dry(:, jdx, :), &
513 0 : active_species_idx_dycore=active_species_idx_dycore)
514 311688000 : else if (present(temp) .and. present(sum_q)) then
515 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
516 : temp=temp(:, jdx, :), sum_q=sum_q(:, jdx, :), &
517 311688000 : active_species_idx_dycore=active_species_idx_dycore)
518 0 : else if (present(dp_dry) .and. present(sum_q)) then
519 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
520 : dp_dry=dp_dry(:, jdx, :), sum_q=sum_q(:, jdx, :), &
521 0 : active_species_idx_dycore=active_species_idx_dycore)
522 0 : else if (present(temp)) then
523 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
524 : temp=temp(:, jdx, :), &
525 0 : active_species_idx_dycore=active_species_idx_dycore)
526 0 : else if (present(dp_dry)) then
527 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
528 : dp_dry=dp_dry(:, jdx, :), &
529 0 : active_species_idx_dycore=active_species_idx_dycore)
530 0 : else if (present(sum_q)) then
531 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
532 : sum_q=sum_q(:, jdx, :), &
533 0 : active_species_idx_dycore=active_species_idx_dycore)
534 : else
535 : call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :), &
536 0 : active_species_idx_dycore=active_species_idx_dycore)
537 : end if
538 : end do
539 :
540 77922000 : end subroutine get_virtual_temp_2hd
541 :
542 : !===========================================================================
543 :
544 : !
545 : !***************************************************************************
546 : !
547 : ! get_sum_species:
548 : !
549 : ! Compute sum of thermodynamically active species
550 : !
551 : ! tracer is in units of dry mixing ratio unless optional argument
552 : ! dp_dry is present in which case tracer is in units of "mass" (=m*dp)
553 : !
554 : !***************************************************************************
555 : !
556 311688000 : subroutine get_sum_species_1hd(tracer, active_species_idx, &
557 311688000 : sum_species, dp_dry, factor)
558 : use air_composition, only: dry_air_species_num
559 :
560 : ! Dummy arguments
561 : ! tracer: Tracer array
562 : real(r8), intent(in) :: tracer(:, :, :)
563 : ! active_species_idx: Index for thermodynamic active tracers
564 : integer, intent(in) :: active_species_idx(:)
565 : ! dp_dry: Dry pressure level thickness.
566 : ! If present, then tracer is in units of mass
567 : real(r8), optional, intent(in) :: dp_dry(:, :)
568 : ! sum_species: sum species
569 : real(r8), intent(out) :: sum_species(:, :)
570 : ! factor: to moist factor
571 : real(r8), optional, intent(out) :: factor(:, :)
572 : ! Local variables
573 623376000 : real(r8) :: factor_loc(SIZE(tracer, 1), SIZE(tracer, 2))
574 : integer :: qdx, itrac
575 311688000 : if (present(dp_dry)) then
576 0 : factor_loc = 1.0_r8 / dp_dry(:,:)
577 : else
578 >14524*10^7 : factor_loc = 1.0_r8
579 : end if
580 >14524*10^7 : sum_species = 1.0_r8 ! all dry air species sum to 1
581 2181816000 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
582 1870128000 : itrac = active_species_idx(qdx)
583 >87179*10^7 : sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_loc(:,:))
584 : end do
585 311688000 : if (present(factor)) then
586 >14555*10^7 : factor = factor_loc
587 : end if
588 311688000 : end subroutine get_sum_species_1hd
589 :
590 : !===========================================================================
591 :
592 0 : subroutine get_sum_species_2hd(tracer, active_species_idx, &
593 0 : sum_species,dp_dry, factor)
594 :
595 : ! Dummy arguments
596 : ! tracer: Tracer array
597 : real(r8), intent(in) :: tracer(:, :, :, :)
598 : ! active_species_idx: Index for thermodynamic active tracers
599 : integer, intent(in) :: active_species_idx(:)
600 : ! dp_dry: Dry pressure level thickness.
601 : ! If present, then tracer is in units of mass
602 : real(r8), optional, intent(in) :: dp_dry(:, :, :)
603 : ! sum_species: sum species
604 : real(r8), intent(out) :: sum_species(:, :, :)
605 : ! factor: to moist factor
606 : real(r8), optional, intent(out) :: factor(:, :, :)
607 : ! Local variable
608 : integer :: jdx
609 :
610 0 : do jdx = 1, SIZE(tracer, 2)
611 0 : if (present(dp_dry) .and. present(factor)) then
612 : call get_sum_species(tracer(:, jdx, :, :), active_species_idx, &
613 0 : sum_species(:, jdx, :), dp_dry=dp_dry(:, jdx, :), factor=factor(:, jdx, :))
614 0 : else if (present(dp_dry)) then
615 : call get_sum_species(tracer(:, jdx, :, :), active_species_idx, &
616 0 : sum_species(:, jdx, :), dp_dry=dp_dry(:, jdx, :))
617 0 : else if (present(factor)) then
618 : call get_sum_species(tracer(:, jdx, :, :), active_species_idx, &
619 0 : sum_species(:, jdx, :), factor=factor(:, jdx, :))
620 : else
621 : call get_sum_species(tracer(:, jdx, :, :), active_species_idx, &
622 0 : sum_species(:, jdx, :))
623 : end if
624 : end do
625 :
626 0 : end subroutine get_sum_species_2hd
627 :
628 : !===========================================================================
629 :
630 : !***************************************************************************
631 : !
632 : ! get_dp: Compute pressure level thickness from dry pressure and
633 : ! thermodynamic active species mixing ratios
634 : !
635 : ! Tracer can either be in units of dry mixing ratio (mixing_ratio=1) or
636 : ! "mass" (=m*dp_dry) (mixing_ratio=2)
637 : !
638 : !***************************************************************************
639 : !
640 51991200 : subroutine get_dp_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, dp, ps, ptop)
641 : use air_composition, only: dry_air_species_num
642 : use string_utils, only: int2str
643 :
644 : real(r8), intent(in) :: tracer(:, :, :) ! tracers; quantity specified by mixing_ratio arg
645 : integer, intent(in) :: mixing_ratio ! 1 => tracer is dry mixing ratio
646 : ! 2 => tracer is mass (q*dp)
647 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
648 : real(r8), intent(in) :: dp_dry(:, :) ! dry pressure level thickness
649 : real(r8), intent(out) :: dp(:, :) ! pressure level thickness
650 : real(r8), optional,intent(out) :: ps(:) ! surface pressure (if ps present then ptop
651 : ! must be present)
652 : real(r8), optional,intent(in) :: ptop ! pressure at model top
653 :
654 : integer :: idx, kdx, m_cnst, qdx
655 :
656 : character(len=*), parameter :: subname = 'get_dp_1hd: '
657 :
658 24279890400 : dp = dp_dry
659 51991200 : if (mixing_ratio == DRY_MIXING_RATIO) then
660 0 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
661 0 : m_cnst = active_species_idx(qdx)
662 0 : do kdx = 1, SIZE(tracer, 2)
663 0 : do idx = 1, SIZE(tracer, 1)
664 0 : dp(idx, kdx) = dp(idx, kdx) + dp_dry(idx, kdx)*tracer(idx, kdx, m_cnst)
665 : end do
666 : end do
667 : end do
668 51991200 : else if (mixing_ratio == MASS_MIXING_RATIO) then
669 363938400 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
670 311947200 : m_cnst = active_species_idx(qdx)
671 29375028000 : do kdx = 1, SIZE(tracer, 2)
672 >14536*10^7 : do idx = 1, SIZE(tracer, 1)
673 >14505*10^7 : dp(idx, kdx) = dp(idx, kdx) + tracer(idx, kdx, m_cnst)
674 : end do
675 : end do
676 : end do
677 : else
678 0 : call endrun(subname//'unrecognized input ('//int2str(mixing_ratio)//') for mixing_ratio')
679 : end if
680 51991200 : if (present(ps)) then
681 0 : if (present(ptop)) then
682 0 : ps = ptop
683 0 : do kdx = 1, SIZE(tracer, 2)
684 0 : do idx = 1, SIZE(tracer, 1)
685 0 : ps(idx) = ps(idx) + dp(idx, kdx)
686 : end do
687 : end do
688 : else
689 0 : call endrun(subname//'if ps is present ptop must be present')
690 : end if
691 : end if
692 51991200 : end subroutine get_dp_1hd
693 :
694 12997800 : subroutine get_dp_2hd(tracer, mixing_ratio, active_species_idx, dp_dry, dp, ps, ptop)
695 : ! Version of get_dp for arrays that have a second horizontal index
696 : real(r8), intent(in) :: tracer(:,:,:,:) ! tracers; quantity specified by mixing_ratio arg
697 : integer, intent(in) :: mixing_ratio ! 1 => tracer is dry mixing ratio
698 : ! 2 => tracer is mass (q*dp)
699 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
700 : real(r8), intent(in) :: dp_dry(:,:,:) ! dry pressure level thickness
701 : real(r8), intent(out) :: dp(:,:,:) ! pressure level thickness
702 : real(r8), optional,intent(out) :: ps(:,:) ! surface pressure
703 : real(r8), optional,intent(in) :: ptop ! pressure at model top
704 :
705 : integer :: jdx
706 :
707 64989000 : do jdx = 1, SIZE(tracer, 2)
708 64989000 : if (present(ps)) then
709 : call get_dp(tracer(:, jdx, :, :), mixing_ratio, active_species_idx, &
710 0 : dp_dry(:, jdx, :), dp(:, jdx, :), ps=ps(:,jdx), ptop=ptop)
711 : else
712 : call get_dp(tracer(:, jdx, :, :), mixing_ratio, active_species_idx, &
713 51991200 : dp_dry(:, jdx, :), dp(:, jdx, :), ptop=ptop)
714 : end if
715 : end do
716 :
717 12997800 : end subroutine get_dp_2hd
718 : !===========================================================================
719 :
720 : !*************************************************************************************************************************
721 : !
722 : ! compute mid-level (full level) pressure from dry pressure and water tracers
723 : !
724 : !*************************************************************************************************************************
725 : !
726 0 : subroutine get_pmid_from_dpdry_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, pmid, pint, dp)
727 :
728 : real(r8), intent(in) :: tracer(:,:,:) ! tracers; quantity specified by mixing_ratio arg
729 : integer, intent(in) :: mixing_ratio ! 1 => tracer is mixing ratio
730 : ! 2 => tracer is mass (q*dp)
731 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
732 : real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness
733 : real(r8), intent(in) :: ptop ! model top pressure
734 : real(r8), intent(out) :: pmid(:,:) ! mid-level pressure
735 : real(r8), optional, intent(out) :: pint(:,:) ! half-level pressure
736 : real(r8), optional, intent(out) :: dp(:,:) ! presure level thickness
737 :
738 0 : real(r8) :: dp_local(SIZE(tracer, 1), SIZE(tracer, 2)) ! local pressure level thickness
739 0 : real(r8) :: pint_local(SIZE(tracer, 1), SIZE(tracer, 2) + 1)! local interface pressure
740 :
741 0 : call get_dp(tracer, mixing_ratio, active_species_idx, dp_dry, dp_local)
742 :
743 0 : call get_pmid_from_dp(dp_local, ptop, pmid, pint_local)
744 :
745 0 : if (present(pint)) pint=pint_local
746 0 : if (present(dp)) dp=dp_local
747 0 : end subroutine get_pmid_from_dpdry_1hd
748 :
749 : !===========================================================================
750 :
751 : !*************************************************************************************************************************
752 : !
753 : ! compute mid-level (full level) pressure
754 : !
755 : !*************************************************************************************************************************
756 : !
757 311688000 : subroutine get_pmid_from_dp_1hd(dp, ptop, pmid, pint)
758 : use dycore, only: dycore_is
759 : real(r8), intent(in) :: dp(:,:) ! pressure level thickness
760 : real(r8), intent(in) :: ptop ! pressure at model top
761 : real(r8), intent(out) :: pmid(:,:) ! mid (full) level pressure
762 : real(r8), optional, intent(out) :: pint(:,:) ! pressure at interfaces (half levels)
763 :
764 623376000 : real(r8) :: pint_local(SIZE(dp, 1), SIZE(dp,2) + 1)
765 : integer :: kdx
766 :
767 1558440000 : pint_local(:, 1) = ptop
768 29298672000 : do kdx = 2, SIZE(dp, 2) + 1
769 >14524*10^7 : pint_local(:, kdx) = dp(:, kdx - 1) + pint_local(:, kdx - 1)
770 : end do
771 :
772 311688000 : if (dycore_is('LR') .or. dycore_is('FV3')) then
773 0 : do kdx = 1, SIZE(dp, 2)
774 0 : pmid(:, kdx) = dp(:, kdx) / (log(pint_local(:, kdx + 1)) - log(pint_local(:, kdx)))
775 : end do
776 : else
777 29298672000 : do kdx = 1, SIZE(dp, 2)
778 >14524*10^7 : pmid(:, kdx) = 0.5_r8 * (pint_local(:, kdx) + pint_local(:, kdx + 1))
779 : end do
780 : end if
781 >14680*10^7 : if (present(pint)) pint=pint_local
782 311688000 : end subroutine get_pmid_from_dp_1hd
783 :
784 : !===========================================================================
785 :
786 : !****************************************************************************************************************
787 : !
788 : ! Compute Exner pressure
789 : !
790 : !****************************************************************************************************************
791 : !
792 0 : subroutine get_exner_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, inv_exner, exner, poverp0)
793 : use string_utils, only: int2str
794 : real(r8), intent(in) :: tracer(:,:,:) ! tracers; quantity specified by mixing_ratio arg
795 : integer, intent(in) :: mixing_ratio ! 1 => tracer is mixing ratio
796 : ! 2 => tracer is mass (q*dp)
797 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
798 : real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness
799 : real(r8), intent(in) :: ptop ! pressure at model top
800 : real(r8), intent(in) :: p00 ! reference pressure for Exner pressure (usually 1000hPa)
801 : logical , intent(in) :: inv_exner ! logical for outputting inverse Exner or Exner pressure
802 : real(r8), intent(out) :: exner(:,:)
803 : real(r8), optional, intent(out) :: poverp0(:,:) ! for efficiency when a routine needs this variable
804 :
805 0 : real(r8) :: pmid(SIZE(tracer, 1), SIZE(tracer, 2))
806 0 : real(r8) :: kappa_dry(SIZE(tracer, 1), SIZE(tracer, 2))
807 : character(len=*), parameter :: subname = 'get_exner_1hd: '
808 : !
809 : ! compute mid level pressure
810 : !
811 0 : call get_pmid_from_dp(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, pmid)
812 : !
813 : ! compute kappa = Rd / cpd
814 : !
815 0 : if (mixing_ratio == DRY_MIXING_RATIO) then
816 0 : call get_kappa_dry(tracer, active_species_idx, kappa_dry)
817 0 : else if (mixing_ratio == MASS_MIXING_RATIO) then
818 0 : call get_kappa_dry(tracer, active_species_idx, kappa_dry, 1.0_r8 / dp_dry)
819 : else
820 0 : call endrun(subname//'unrecognized input ('//int2str(mixing_ratio)//') for mixing_ratio')
821 : end if
822 0 : if (inv_exner) then
823 0 : exner(:,:) = (p00 / pmid(:,:)) ** kappa_dry(:,:)
824 : else
825 0 : exner(:,:) = (pmid(:,:) / p00) ** kappa_dry(:,:)
826 : end if
827 0 : if (present(poverp0)) poverp0 = pmid(:,:) / p00
828 0 : end subroutine get_exner_1hd
829 :
830 : !===========================================================================
831 :
832 : !****************************************************************************************************************
833 : !
834 : ! Compute virtual potential temperature from dp_dry, m, T and ptop.
835 : !
836 : !****************************************************************************************************************
837 : !
838 0 : subroutine get_virtual_theta_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, temp, theta_v)
839 : real(r8), intent(in) :: tracer(:,:,:) ! tracers; quantity specified by mixing_ratio arg
840 : integer, intent(in) :: mixing_ratio ! 1 => tracer is dry mixing ratio
841 : ! 2 => tracer is mass (q*dp)
842 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
843 : real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness
844 : real(r8), intent(in) :: ptop ! pressure at model top
845 : real(r8), intent(in) :: p00 ! reference pressure for Exner pressure (usually 1000hPa)
846 : real(r8), intent(in) :: temp(:,:) ! temperature
847 : real(r8), intent(out) :: theta_v(:,:) ! virtual potential temperature
848 :
849 0 : real(r8) :: iexner(SIZE(tracer, 1), SIZE(tracer, 2))
850 :
851 0 : call get_exner(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, .true., iexner)
852 :
853 0 : theta_v(:,:) = temp(:,:) * iexner(:,:)
854 :
855 0 : end subroutine get_virtual_theta_1hd
856 :
857 : !===========================================================================
858 :
859 : !****************************************************************************************************************
860 : !
861 : ! Compute geopotential from dry pressure level thichkness, water tracers, model top pressure and temperature
862 : !
863 : !****************************************************************************************************************
864 : !
865 0 : subroutine get_gz_from_dp_dry_ptop_temp_1hd(tracer, mixing_ratio, active_species_idx, &
866 0 : dp_dry, ptop, temp, phis, gz, pmid, dp, T_v)
867 : use air_composition, only: get_R_dry
868 : use string_utils, only: int2str
869 : real(r8), intent(in) :: tracer(:,:,:) ! tracer; quantity specified by mixing_ratio arg
870 : integer, intent(in) :: mixing_ratio ! 1 => tracer is dry mixing ratio
871 : ! 2 => tracer is mass (q*dp)
872 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
873 : real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness
874 : real(r8), intent(in) :: ptop ! pressure at model top
875 : real(r8), intent(in) :: temp(:,:) ! temperature
876 : real(r8), intent(in) :: phis(:) ! surface geopotential
877 : real(r8), intent(out) :: gz(:,:) ! geopotential
878 : real(r8), optional, intent(out) :: pmid(:,:) ! mid-level pressure
879 : real(r8), optional, intent(out) :: dp(:,:) ! pressure level thickness
880 : real(r8), optional, intent(out) :: t_v(:,:) ! virtual temperature
881 :
882 :
883 0 : real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2)) :: pmid_local, t_v_local, dp_local, R_dry
884 0 : real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2) + 1) :: pint
885 : character(len=*), parameter :: subname = 'get_gz_from_dp_dry_ptop_temp_1hd: '
886 :
887 :
888 : call get_pmid_from_dp(tracer, mixing_ratio, active_species_idx, &
889 0 : dp_dry, ptop, pmid_local, pint=pint, dp=dp_local)
890 0 : if (mixing_ratio == DRY_MIXING_RATIO) then
891 0 : call get_virtual_temp(tracer, t_v_local, temp=temp, active_species_idx_dycore=active_species_idx)
892 0 : call get_R_dry(tracer, active_species_idx, R_dry)
893 0 : else if (mixing_ratio == MASS_MIXING_RATIO) then
894 0 : call get_virtual_temp(tracer, t_v_local, temp=temp, dp_dry=dp_dry, active_species_idx_dycore=active_species_idx)
895 0 : call get_R_dry(tracer,active_species_idx, R_dry, fact=1.0_r8 / dp_dry)
896 : else
897 0 : call endrun(subname//'unrecognized input ('//int2str(mixing_ratio)//') for mixing_ratio')
898 : end if
899 0 : call get_gz(dp_local, T_v_local, R_dry, phis, ptop, gz, pmid_local)
900 :
901 0 : if (present(pmid)) pmid=pmid_local
902 0 : if (present(T_v)) T_v=T_v_local
903 0 : if (present(dp)) dp=dp_local
904 0 : end subroutine get_gz_from_dp_dry_ptop_temp_1hd
905 :
906 : !===========================================================================
907 :
908 : !***************************************************************************
909 : !
910 : ! Compute geopotential from pressure level thickness and virtual temperature
911 : !
912 : !***************************************************************************
913 : !
914 311688000 : subroutine get_gz_given_dp_Tv_Rdry_1hd(dp, T_v, R_dry, phis, ptop, gz, pmid)
915 : use dycore, only: dycore_is
916 : real(r8), intent(in) :: dp (:,:) ! pressure level thickness
917 : real(r8), intent(in) :: T_v (:,:) ! virtual temperature
918 : real(r8), intent(in) :: R_dry(:,:) ! R dry
919 : real(r8), intent(in) :: phis (:) ! surface geopotential
920 : real(r8), intent(in) :: ptop ! model top presure
921 : real(r8), intent(out) :: gz(:,:) ! geopotential
922 : real(r8), optional, intent(out) :: pmid(:,:) ! mid-level pressure
923 :
924 :
925 623376000 : real(r8), dimension(SIZE(dp, 1), SIZE(dp, 2)) :: pmid_local
926 623376000 : real(r8), dimension(SIZE(dp, 1), SIZE(dp, 2) + 1) :: pint
927 623376000 : real(r8), dimension(SIZE(dp, 1)) :: gzh, Rdry_tv
928 : integer :: kdx
929 :
930 311688000 : call get_pmid_from_dp(dp, ptop, pmid_local, pint)
931 :
932 : !
933 : ! integrate hydrostatic eqn
934 : !
935 1558440000 : gzh = phis
936 311688000 : if (dycore_is('LR') .or. dycore_is('FV3')) then
937 0 : do kdx = SIZE(dp, 2), 1, -1
938 0 : Rdry_tv(:) = R_dry(:, kdx) * T_v(:, kdx)
939 0 : gz(:, kdx) = gzh(:) + Rdry_tv(:) * (1.0_r8 - pint(:, kdx) / pmid_local(:, kdx))
940 0 : gzh(:) = gzh(:) + Rdry_tv(:) * (log(pint(:, kdx + 1)) - log(pint(:, kdx)))
941 : end do
942 : else
943 29298672000 : do kdx = SIZE(dp,2), 1, -1
944 >14493*10^7 : Rdry_tv(:) = R_dry(:,kdx) * T_v(:, kdx)
945 >14493*10^7 : gz(:,kdx) = gzh(:) + Rdry_tv(:) * 0.5_r8 * dp(:, kdx) / pmid_local(:, kdx)
946 >14524*10^7 : gzh(:) = gzh(:) + Rdry_tv(:) * dp(:, kdx) / pmid_local(:, kdx)
947 : end do
948 : end if
949 >14524*10^7 : if (present(pmid)) pmid=pmid_local
950 311688000 : end subroutine get_gz_given_dp_Tv_Rdry_1hd
951 :
952 77922000 : subroutine get_gz_given_dp_Tv_Rdry_2hd(dp, T_v, R_dry, phis, ptop, gz, pmid)
953 : ! Version of get_gz_given_dp_Tv_Rdry for arrays that have a second horizontal index
954 : real(r8), intent(in) :: dp (:,:,:) ! pressure level thickness
955 : real(r8), intent(in) :: T_v (:,:,:) ! virtual temperature
956 : real(r8), intent(in) :: R_dry(:,:,:) ! R dry
957 : real(r8), intent(in) :: phis (:,:) ! surface geopotential
958 : real(r8), intent(in) :: ptop ! model top presure
959 : real(r8), intent(out) :: gz(:,:,:) ! geopotential
960 : real(r8), optional, intent(out) :: pmid(:,:,:) ! mid-level pressure
961 :
962 : integer :: jdx
963 :
964 389610000 : do jdx = 1, SIZE(dp, 2)
965 389610000 : if (present(pmid)) then
966 : call get_gz(dp(:, jdx, :), T_v(:, jdx, :), R_dry(:, jdx, :), phis(:, jdx), &
967 311688000 : ptop, gz(:, jdx, :), pmid=pmid(:, jdx, :))
968 : else
969 0 : call get_gz(dp(:, jdx, :), T_v(:, jdx, :), R_dry(:, jdx, :), phis(:, jdx), ptop, gz(:, jdx, :))
970 : end if
971 : end do
972 :
973 :
974 77922000 : end subroutine get_gz_given_dp_Tv_Rdry_2hd
975 :
976 : !===========================================================================
977 :
978 : !***************************************************************************
979 : !
980 : ! Compute Richardson number at cell interfaces (half levels)
981 : !
982 : !***************************************************************************
983 : !
984 0 : subroutine get_Richardson_number_1hd(tracer,mixing_ratio, active_species_idx, dp_dry, ptop, &
985 0 : p00, temp, v, Richardson_number, pmid, dp)
986 : real(r8), intent(in) :: tracer(:,:,:) ! tracer; quantity specified by mixing_ratio arg
987 : integer, intent(in) :: mixing_ratio ! 1 => tracer is dry mixing ratio
988 : ! 2 => tracer is mass (q*dp)
989 : integer, intent(in) :: active_species_idx(:) ! index for thermodynamic species in tracer array
990 : real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness
991 : real(r8), intent(in) :: ptop ! pressure at model top
992 : real(r8), intent(in) :: p00 ! reference pressure for Exner pressure (usually 1000hPa)
993 : real(r8), intent(in) :: temp(:,:) ! temperature
994 : real(r8), intent(in) :: v(:,:,:) ! velocity components
995 : real(r8), intent(out) :: Richardson_number(:,:)
996 : real(r8), optional, intent(out) :: pmid(:,:)
997 : real(r8), optional, intent(out) :: dp(:,:)
998 :
999 0 : real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2)) :: gz, theta_v
1000 0 : real(r8), dimension(SIZE(tracer, 1)) :: pt1, pt2, phis
1001 : integer :: kdx, kdxm1
1002 : real(r8), parameter:: ustar2 = 1.E-4_r8
1003 :
1004 0 : phis = 0.0_r8
1005 0 : call get_gz(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, temp, phis, gz, pmid=pmid, dp=dp)
1006 0 : call get_virtual_theta(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, temp, theta_v)
1007 0 : Richardson_number(:, 1) = 0.0_r8
1008 0 : Richardson_number(:, SIZE(tracer, 2) + 1) = 0.0_r8
1009 0 : do kdx = SIZE(tracer, 2), 2, -1
1010 0 : kdxm1 = kdx - 1
1011 0 : pt1(:) = theta_v(:, kdxm1)
1012 0 : pt2(:) = theta_v(:, kdx)
1013 0 : Richardson_number(:, kdx) = (gz(:, kdxm1) - gz(:, kdx)) * (pt1 - pt2) / ( 0.5_r8*(pt1 + pt2) * &
1014 0 : ((v(:, 1, kdxm1) - v(:, 1, kdx)) ** 2 + (v(:, 2, kdxm1) - v(:, 2, kdx)) ** 2 + ustar2) )
1015 : end do
1016 0 : end subroutine get_Richardson_number_1hd
1017 :
1018 : !
1019 : !****************************************************************************************************************
1020 : !
1021 : ! get surface pressure from dry pressure and thermodynamic active species (e.g., forms of water: water vapor, cldliq, etc.)
1022 : !
1023 : !****************************************************************************************************************
1024 : !
1025 0 : subroutine get_ps_1hd(tracer_mass, active_species_idx, dp_dry, ps, ptop)
1026 : use air_composition, only: dry_air_species_num
1027 :
1028 : real(r8), intent(in) :: tracer_mass(:,:,:) ! Tracer array (q*dp)
1029 : real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness
1030 : real(r8), intent(out) :: ps(:) ! surface pressure
1031 : real(r8), intent(in) :: ptop
1032 : integer, intent(in) :: active_species_idx(:)
1033 :
1034 : integer :: idx, kdx, m_cnst, qdx
1035 0 : real(r8) :: dp(SIZE(tracer_mass, 1), SIZE(tracer_mass, 2)) ! dry pressure level thickness
1036 :
1037 0 : dp = dp_dry
1038 0 : do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
1039 0 : m_cnst = active_species_idx(qdx)
1040 0 : do kdx = 1, SIZE(tracer_mass, 2)
1041 0 : do idx = 1, SIZE(tracer_mass, 1)
1042 0 : dp(idx, kdx) = dp(idx, kdx) + tracer_mass(idx, kdx, m_cnst)
1043 : end do
1044 : end do
1045 : end do
1046 0 : ps = ptop
1047 0 : do kdx = 1, SIZE(tracer_mass, 2)
1048 0 : do idx = 1, SIZE(tracer_mass, 1)
1049 0 : ps(idx) = ps(idx) + dp(idx, kdx)
1050 : end do
1051 : end do
1052 0 : end subroutine get_ps_1hd
1053 :
1054 0 : subroutine get_ps_2hd(tracer_mass, active_species_idx, dp_dry, ps, ptop)
1055 : ! Version of get_ps for arrays that have a second horizontal index
1056 : real(r8), intent(in) :: tracer_mass(:,:,:,:) ! Tracer array (q*dp)
1057 : real(r8), intent(in) :: dp_dry(:,:,:) ! dry pressure level thickness
1058 : real(r8), intent(out) :: ps(:,:) ! surface pressure
1059 : real(r8), intent(in) :: ptop
1060 : integer, intent(in) :: active_species_idx(:)
1061 :
1062 : integer :: jdx
1063 :
1064 0 : do jdx = 1, SIZE(tracer_mass, 2)
1065 0 : call get_ps(tracer_mass(:, jdx, :, :), active_species_idx, dp_dry(:, jdx, :), ps(:, jdx), ptop)
1066 : end do
1067 :
1068 0 : end subroutine get_ps_2hd
1069 :
1070 : !===========================================================================
1071 :
1072 : !*************************************************************************************************************************
1073 : !
1074 : ! compute generalized kappa =Rdry/cpdry
1075 : !
1076 : !*************************************************************************************************************************
1077 : !
1078 62337600 : subroutine get_kappa_dry_1hd(tracer, active_species_idx, kappa_dry, fact)
1079 : use air_composition, only: dry_air_species_num, get_R_dry, get_cp_dry
1080 : use physconst, only: rair, cpair
1081 :
1082 : real(r8), intent(in) :: tracer(:,:,:) !tracer array
1083 : integer, intent(in) :: active_species_idx(:) !index of thermodynamic active tracers
1084 : real(r8), intent(out) :: kappa_dry(:,:) !kappa dry
1085 : real(r8), optional, intent(in) :: fact(:,:) !factor for converting tracer to dry mixing ratio
1086 : !
1087 62337600 : real(r8), allocatable, dimension(:,:) :: cp_dry,R_dry
1088 : integer :: ierr
1089 : character(len=*), parameter :: subname = "get_kappa_dry_1hd"
1090 : character(len=*), parameter :: errstr = subname//": failed to allocate "
1091 : !
1092 : ! dry air not species dependent
1093 62337600 : if (dry_air_species_num==0) then
1094 29049321600 : kappa_dry = rair / cpair
1095 : else
1096 0 : allocate(R_dry(SIZE(kappa_dry, 1), SIZE(kappa_dry, 2)), stat=ierr)
1097 0 : if (ierr /= 0) then
1098 0 : call endrun(errstr//"R_dry")
1099 : end if
1100 0 : allocate(cp_dry(SIZE(kappa_dry, 1), SIZE(kappa_dry, 2)), stat=ierr)
1101 0 : if (ierr /= 0) then
1102 0 : call endrun(errstr//"cp_dry")
1103 : end if
1104 0 : call get_cp_dry(tracer, active_species_idx, cp_dry, fact=fact)
1105 0 : call get_R_dry( tracer, active_species_idx, R_dry, fact=fact)
1106 0 : kappa_dry = R_dry / cp_dry
1107 0 : deallocate(R_dry, cp_dry)
1108 : end if
1109 62337600 : end subroutine get_kappa_dry_1hd
1110 :
1111 15584400 : subroutine get_kappa_dry_2hd(tracer, active_species_idx, kappa_dry, fact)
1112 : ! Version of get_kappa_dry for arrays that have a second horizontal index
1113 : real(r8), intent(in) :: tracer(:,:,:,:) !tracer array
1114 : integer, intent(in) :: active_species_idx(:) !index of thermodynamic active tracers
1115 : real(r8), intent(out) :: kappa_dry(:,:,:) !kappa dry
1116 : real(r8), optional, intent(in) :: fact(:,:,:) !factor for converting tracer to dry mixing ratio
1117 :
1118 : integer :: jdx
1119 :
1120 77922000 : do jdx = 1, SIZE(tracer, 2)
1121 77922000 : if (present(fact)) then
1122 0 : call get_kappa_dry(tracer(:, jdx, :, :), active_species_idx, kappa_dry(:, jdx, :), fact=fact(:, jdx, :))
1123 : else
1124 62337600 : call get_kappa_dry(tracer(:, jdx, :, :), active_species_idx, kappa_dry(:, jdx, :))
1125 : end if
1126 : end do
1127 :
1128 15584400 : end subroutine get_kappa_dry_2hd
1129 :
1130 : !===========================================================================
1131 :
1132 : !*************************************************************************************************************************
1133 : !
1134 : ! compute reference pressure levels
1135 : !
1136 : !*************************************************************************************************************************
1137 : !
1138 43200 : subroutine get_dp_ref_1hd(hyai, hybi, ps0, phis, dp_ref, ps_ref)
1139 : use physconst, only: tref, rair
1140 : real(r8), intent(in) :: hyai(:)
1141 : real(r8), intent(in) :: hybi(:)
1142 : real(r8), intent(in) :: ps0
1143 : real(r8), intent(in) :: phis(:)
1144 : real(r8), intent(out) :: dp_ref(:,:)
1145 : real(r8), intent(out) :: ps_ref(:)
1146 : integer :: kdx
1147 : !
1148 : ! use static reference pressure (hydrostatic balance incl. effect of topography)
1149 : !
1150 216000 : ps_ref(:) = ps0 * exp(-phis(:) / (rair * tref))
1151 4060800 : do kdx = 1, SIZE(dp_ref, 2)
1152 20131200 : dp_ref(:,kdx) = ((hyai(kdx + 1) - hyai(kdx)) * ps0 + (hybi(kdx + 1) - hybi(kdx)) * ps_ref(:))
1153 : end do
1154 43200 : end subroutine get_dp_ref_1hd
1155 :
1156 10800 : subroutine get_dp_ref_2hd(hyai, hybi, ps0, phis, dp_ref, ps_ref)
1157 : ! Version of get_dp_ref for arrays that have a second horizontal index
1158 : real(r8), intent(in) :: hyai(:)
1159 : real(r8), intent(in) :: hybi(:)
1160 : real(r8), intent(in) :: ps0
1161 : real(r8), intent(in) :: phis(:,:)
1162 : real(r8), intent(out) :: dp_ref(:,:,:)
1163 : real(r8), intent(out) :: ps_ref(:,:)
1164 : integer :: jdx
1165 :
1166 54000 : do jdx = 1, SIZE(dp_ref, 2)
1167 54000 : call get_dp_ref(hyai, hybi, ps0, phis(:, jdx), dp_ref(:, jdx, :), ps_ref(:, jdx))
1168 : end do
1169 :
1170 10800 : end subroutine get_dp_ref_2hd
1171 :
1172 : !===========================================================================
1173 :
1174 : !*************************************************************************************************************************
1175 : !
1176 : ! compute dry densisty from temperature (temp) and pressure (dp_dry and tracer)
1177 : !
1178 : !*************************************************************************************************************************
1179 : !
1180 0 : subroutine get_rho_dry_1hd(tracer, temp, ptop, dp_dry, tracer_mass, rho_dry, rhoi_dry, &
1181 0 : active_species_idx_dycore)
1182 : use air_composition, only: get_R_dry
1183 : ! args
1184 : real(r8), intent(in) :: tracer(:,:,:) ! Tracer array
1185 : real(r8), intent(in) :: temp(:,:) ! Temperature
1186 : real(r8), intent(in) :: ptop
1187 : real(r8), intent(in) :: dp_dry(:,:)
1188 : logical, intent(in) :: tracer_mass
1189 : real(r8), optional,intent(out) :: rho_dry(:,:)
1190 : real(r8), optional,intent(out) :: rhoi_dry(:,:)
1191 : !
1192 : ! array of indicies for index of thermodynamic active species in dycore tracer array
1193 : ! (if different from physics index)
1194 : !
1195 : integer, optional, intent(in) :: active_species_idx_dycore(:)
1196 :
1197 : ! local vars
1198 : integer :: idx, kdx
1199 0 : real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2)) :: pmid
1200 0 : real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2) + 1) :: pint
1201 0 : real(r8), allocatable :: R_dry(:,:)
1202 0 : integer, dimension(thermodynamic_active_species_num) :: idx_local
1203 : integer :: ierr
1204 : character(len=*), parameter :: subname = "get_rho_dry_1hd"
1205 : character(len=*), parameter :: errstr = subname//": failed to allocate "
1206 :
1207 0 : if (present(active_species_idx_dycore)) then
1208 0 : idx_local = active_species_idx_dycore
1209 : else
1210 0 : idx_local = thermodynamic_active_species_idx
1211 : end if
1212 : !
1213 : ! we assume that air is dry where molecular viscosity may be significant
1214 : !
1215 0 : call get_pmid_from_dp(dp_dry, ptop, pmid, pint=pint)
1216 0 : if (present(rhoi_dry)) then
1217 0 : allocate(R_dry(SIZE(tracer, 1), SIZE(tracer, 2) + 1), stat=ierr)
1218 0 : if (ierr /= 0) then
1219 0 : call endrun(errstr//"R_dry")
1220 : end if
1221 0 : if (tracer_mass) then
1222 0 : call get_R_dry(tracer, idx_local, R_dry, fact=1.0_r8 / dp_dry)
1223 : else
1224 0 : call get_R_dry(tracer, idx_local, R_dry)
1225 : end if
1226 0 : do kdx = 2, SIZE(tracer, 2) + 1
1227 0 : rhoi_dry(:, kdx) = 0.5_r8 * (temp(:, kdx) + temp(:, kdx - 1))!could be more accurate!
1228 0 : rhoi_dry(:, kdx) = pint(:,kdx) / (rhoi_dry(:, kdx) * R_dry(:, kdx)) !ideal gas law for dry air
1229 : end do
1230 : !
1231 : ! extrapolate top level value
1232 : !
1233 0 : kdx=1
1234 0 : rhoi_dry(:, kdx) = 1.5_r8 * (temp(:, kdx) - 0.5_r8 * temp(:, kdx + 1))
1235 0 : rhoi_dry(:, kdx) = pint(:, kdx) / (rhoi_dry(:, kdx) * R_dry(:, kdx)) !ideal gas law for dry air
1236 0 : deallocate(R_dry)
1237 : end if
1238 0 : if (present(rho_dry)) then
1239 0 : allocate(R_dry(SIZE(tracer, 1), size(rho_dry, 2)), stat=ierr)
1240 0 : if (ierr /= 0) then
1241 0 : call endrun(errstr//"R_dry")
1242 : end if
1243 0 : if (tracer_mass) then
1244 0 : call get_R_dry(tracer, idx_local, R_dry, fact=1.0_r8 / dp_dry)
1245 : else
1246 0 : call get_R_dry(tracer, idx_local, R_dry)
1247 : end if
1248 0 : do kdx = 1, SIZE(rho_dry, 2)
1249 0 : do idx = 1, SIZE(rho_dry, 1)
1250 0 : rho_dry(idx, kdx) = pmid(idx, kdx) / (temp(idx, kdx) * R_dry(idx, kdx)) !ideal gas law for dry air
1251 : end do
1252 : end do
1253 0 : deallocate(R_dry)
1254 : end if
1255 0 : end subroutine get_rho_dry_1hd
1256 :
1257 0 : subroutine get_rho_dry_2hd(tracer, temp, ptop, dp_dry, tracer_mass, rho_dry, rhoi_dry, &
1258 0 : active_species_idx_dycore)
1259 : ! Version of get_rho_dry for arrays that have a second horizontal index
1260 : real(r8), intent(in) :: tracer(:,:,:,:) ! Tracer array
1261 : real(r8), intent(in) :: temp(:,:,:) ! Temperature
1262 : real(r8), intent(in) :: ptop
1263 : real(r8), intent(in) :: dp_dry(:,:,:)
1264 : logical, intent(in) :: tracer_mass
1265 : real(r8), optional,intent(out) :: rho_dry(:,:,:)
1266 : real(r8), optional,intent(out) :: rhoi_dry(:,:,:)
1267 : !
1268 : ! array of indicies for index of thermodynamic active species in dycore tracer array
1269 : ! (if different from physics index)
1270 : !
1271 : integer, optional, intent(in) :: active_species_idx_dycore(:)
1272 :
1273 : integer :: jdx
1274 :
1275 0 : do jdx = 1, SIZE(tracer, 2)
1276 0 : if (present(rho_dry) .and. present(rhoi_dry)) then
1277 : call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), &
1278 : tracer_mass, rho_dry=rho_dry(:, jdx, :), rhoi_dry=rhoi_dry(:, jdx, :), &
1279 0 : active_species_idx_dycore=active_species_idx_dycore)
1280 0 : else if (present(rho_dry)) then
1281 : call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), &
1282 0 : tracer_mass, rho_dry=rho_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore)
1283 0 : else if (present(rhoi_dry)) then
1284 : call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), &
1285 0 : tracer_mass, rhoi_dry=rhoi_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore)
1286 : else
1287 : call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), tracer_mass, &
1288 0 : active_species_idx_dycore=active_species_idx_dycore)
1289 : end if
1290 : end do
1291 :
1292 0 : end subroutine get_rho_dry_2hd
1293 : !===========================================================================
1294 :
1295 : !*************************************************************************************************************************
1296 : !
1297 : ! compute 3D molecular diffusion and thermal conductivity
1298 : !
1299 : !*************************************************************************************************************************
1300 : !
1301 0 : subroutine get_molecular_diff_coef_1hd(temp, get_at_interfaces, sponge_factor, kmvis, kmcnd, &
1302 0 : tracer, fact, active_species_idx_dycore, mbarv_in)
1303 : use air_composition, only: dry_air_species_num, get_mbarv
1304 : use air_composition, only: kv1, kc1, kv2, kc2, kv_temp_exp, kc_temp_exp
1305 :
1306 : ! args
1307 : real(r8), intent(in) :: temp(:,:) ! temperature
1308 : logical, intent(in) :: get_at_interfaces ! true: compute kmvis and kmcnd at interfaces
1309 : ! false: compute kmvis and kmcnd at mid-levels
1310 : real(r8), intent(in) :: sponge_factor(:) ! multiply kmvis and kmcnd with sponge_factor
1311 : ! (for sponge layer)
1312 : real(r8), intent(out) :: kmvis(:,:)
1313 : real(r8), intent(out) :: kmcnd(:,:)
1314 : real(r8), intent(in) :: tracer(:,:,:) ! tracer array
1315 : integer, intent(in), optional :: active_species_idx_dycore(:) ! index of active species in tracer
1316 : real(r8), intent(in), optional :: fact(:,:) ! if tracer is in units of mass or moist
1317 : ! fact converts to dry mixing ratio: tracer/fact
1318 : real(r8), intent(in), optional :: mbarv_in(:,:) ! composition dependent atmosphere mean mass
1319 : !
1320 : ! local vars
1321 : !
1322 : integer :: idx, kdx, icnst, ispecies
1323 : real(r8):: mbarvi, mm, residual ! Mean mass at mid level
1324 : real(r8):: cnst_vis, cnst_cnd, temp_local
1325 0 : real(r8), dimension(SIZE(tracer,1), SIZE(sponge_factor, 1)) :: factor, mbarv
1326 0 : integer, dimension(thermodynamic_active_species_num) :: idx_local
1327 : character(len=*), parameter :: subname = 'get_molecular_diff_coef_1hd: '
1328 :
1329 : !--------------------------------------------
1330 : ! Set constants needed for updates
1331 : !--------------------------------------------
1332 :
1333 0 : if (dry_air_species_num==0) then
1334 :
1335 0 : cnst_vis = (kv1 * mmro2 * o2_mwi + kv2 * mmrn2 * n2_mwi) * mbar
1336 0 : cnst_cnd = (kc1 * mmro2 * o2_mwi + kc2 * mmrn2 * n2_mwi) * mbar
1337 0 : if (get_at_interfaces) then
1338 0 : do kdx = 2, SIZE(sponge_factor, 1)
1339 0 : do idx = 1, SIZE(tracer, 1)
1340 0 : temp_local = 0.5_r8 * (temp(idx, kdx) + temp(idx, kdx - 1))
1341 0 : kmvis(idx, kdx) = sponge_factor(kdx) * cnst_vis * temp_local ** kv_temp_exp
1342 0 : kmcnd(idx, kdx) = sponge_factor(kdx) * cnst_cnd * temp_local ** kc_temp_exp
1343 : end do
1344 : end do
1345 : !
1346 : ! extrapolate top level value
1347 : !
1348 0 : kmvis(1:SIZE(tracer, 1), 1) = 1.5_r8 * kmvis(1:SIZE(tracer, 1), 2) - 0.5_r8 * kmvis(1:SIZE(tracer, 1), 3)
1349 0 : kmcnd(1:SIZE(tracer, 1), 1) = 1.5_r8 * kmcnd(1:SIZE(tracer, 1), 2) - 0.5_r8 * kmcnd(1:SIZE(tracer, 1), 3)
1350 : else if (.not. get_at_interfaces) then
1351 0 : do kdx = 1, SIZE(sponge_factor, 1)
1352 0 : do idx = 1, SIZE(tracer, 1)
1353 0 : kmvis(idx, kdx) = sponge_factor(kdx) * cnst_vis * temp(idx, kdx) ** kv_temp_exp
1354 0 : kmcnd(idx, kdx) = sponge_factor(kdx) * cnst_cnd * temp(idx, kdx) ** kc_temp_exp
1355 : end do
1356 : end do
1357 : else
1358 : call endrun(subname//'get_at_interfaces must be .true. or .false.')
1359 : end if
1360 : else
1361 0 : if (present(active_species_idx_dycore)) then
1362 0 : idx_local = active_species_idx_dycore
1363 : else
1364 0 : idx_local = thermodynamic_active_species_idx
1365 : end if
1366 0 : if (present(fact)) then
1367 0 : factor = fact(:,:)
1368 : else
1369 0 : factor = 1.0_r8
1370 : endif
1371 0 : if (present(mbarv_in)) then
1372 0 : mbarv = mbarv_in
1373 : else
1374 0 : call get_mbarv(tracer, idx_local, mbarv, fact=factor)
1375 : end if
1376 : !
1377 : ! major species dependent code
1378 : !
1379 0 : if (get_at_interfaces) then
1380 0 : do kdx = 2, SIZE(sponge_factor, 1)
1381 0 : do idx = 1, SIZE(tracer, 1)
1382 0 : kmvis(idx, kdx) = 0.0_r8
1383 0 : kmcnd(idx, kdx) = 0.0_r8
1384 0 : residual = 1.0_r8
1385 0 : do icnst = 1, dry_air_species_num
1386 0 : ispecies = idx_local(icnst)
1387 0 : mm = 0.5_r8 * (tracer(idx, kdx, ispecies) * factor(idx, kdx) + &
1388 0 : tracer(idx, kdx - 1, ispecies) * factor(idx, kdx-1))
1389 0 : kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
1390 0 : thermodynamic_active_species_mwi(icnst) * mm
1391 0 : kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
1392 0 : thermodynamic_active_species_mwi(icnst) * mm
1393 0 : residual = residual - mm
1394 : end do
1395 0 : icnst = 0 ! N2
1396 0 : kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
1397 0 : thermodynamic_active_species_mwi(icnst) * residual
1398 0 : kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
1399 0 : thermodynamic_active_species_mwi(icnst) * residual
1400 :
1401 0 : temp_local = 0.5_r8 * (temp(idx, kdx - 1) + temp(idx, kdx))
1402 0 : mbarvi = 0.5_r8 * (mbarv(idx, kdx - 1) + mbarv(idx, kdx))
1403 0 : kmvis(idx, kdx) = kmvis(idx, kdx) * mbarvi * temp_local ** kv_temp_exp
1404 0 : kmcnd(idx, kdx) = kmcnd(idx, kdx) * mbarvi * temp_local ** kc_temp_exp
1405 : enddo
1406 : end do
1407 0 : do idx = 1, SIZE(tracer, 1)
1408 0 : kmvis(idx, 1) = 1.5_r8 * kmvis(idx, 2) - .5_r8 * kmvis(idx, 3)
1409 0 : kmcnd(idx, 1) = 1.5_r8 * kmcnd(idx, 2) - .5_r8 * kmcnd(idx, 3)
1410 0 : kmvis(idx, SIZE(sponge_factor, 1) + 1) = kmvis(idx, SIZE(sponge_factor, 1))
1411 0 : kmcnd(idx, SIZE(sponge_factor, 1) + 1) = kmcnd(idx, SIZE(sponge_factor, 1))
1412 : end do
1413 : else if (.not. get_at_interfaces) then
1414 0 : do kdx = 1, SIZE(sponge_factor, 1)
1415 0 : do idx = 1, SIZE(tracer, 1)
1416 0 : kmvis(idx, kdx) = 0.0_r8
1417 0 : kmcnd(idx, kdx) = 0.0_r8
1418 0 : residual = 1.0_r8
1419 0 : do icnst = 1, dry_air_species_num - 1
1420 0 : ispecies = idx_local(icnst)
1421 0 : mm = tracer(idx, kdx, ispecies) * factor(idx, kdx)
1422 0 : kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
1423 0 : thermodynamic_active_species_mwi(icnst) * mm
1424 0 : kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
1425 0 : thermodynamic_active_species_mwi(icnst) * mm
1426 0 : residual = residual - mm
1427 : end do
1428 0 : icnst = dry_air_species_num
1429 0 : kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
1430 0 : thermodynamic_active_species_mwi(icnst) * residual
1431 0 : kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
1432 0 : thermodynamic_active_species_mwi(icnst) * residual
1433 :
1434 0 : kmvis(idx, kdx) = kmvis(idx, kdx) * mbarv(idx, kdx) * temp(idx, kdx) ** kv_temp_exp
1435 0 : kmcnd(idx, kdx) = kmcnd(idx, kdx) * mbarv(idx, kdx) * temp(idx, kdx) ** kc_temp_exp
1436 : end do
1437 : end do
1438 : else
1439 : call endrun(subname//'get_at_interfaces must be .true. or .false.')
1440 : end if
1441 : end if
1442 0 : end subroutine get_molecular_diff_coef_1hd
1443 :
1444 0 : subroutine get_molecular_diff_coef_2hd(temp, get_at_interfaces, sponge_factor, kmvis, kmcnd, &
1445 0 : tracer, fact, active_species_idx_dycore, mbarv_in)
1446 : ! Version of get_molecular_diff_coef for arrays that have a second horizontal index
1447 : real(r8), intent(in) :: temp(:,:,:) ! temperature
1448 : logical, intent(in) :: get_at_interfaces ! true: compute kmvis and kmcnd at interfaces
1449 : ! false: compute kmvis and kmcnd at mid-levels
1450 : real(r8), intent(in) :: sponge_factor(:) ! multiply kmvis and kmcnd with sponge_factor
1451 : ! (for sponge layer)
1452 : real(r8), intent(out) :: kmvis(:,:,:)
1453 : real(r8), intent(out) :: kmcnd(:,:,:)
1454 : real(r8), intent(in) :: tracer(:,:,:,:) ! tracer array
1455 : integer, intent(in), optional :: active_species_idx_dycore(:) ! index of active species in tracer
1456 : real(r8), intent(in), optional :: fact(:,:,:) ! if tracer is in units of mass or moist
1457 : ! fact converts to dry mixing ratio: tracer/fact
1458 : real(r8), intent(in), optional :: mbarv_in(:,:,:) ! composition dependent atmosphere mean mass
1459 : integer :: jdx
1460 :
1461 0 : do jdx = 1, SIZE(tracer, 2)
1462 0 : if (present(fact) .and. present(mbarv_in)) then
1463 : call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
1464 : kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), fact=fact(:, jdx, :), &
1465 0 : active_species_idx_dycore=active_species_idx_dycore, mbarv_in=mbarv_in(:, jdx, :))
1466 0 : else if (present(fact)) then
1467 : call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
1468 : kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), fact=fact(:, jdx, :), &
1469 0 : active_species_idx_dycore=active_species_idx_dycore)
1470 0 : else if (present(mbarv_in)) then
1471 : call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
1472 : kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), &
1473 0 : active_species_idx_dycore=active_species_idx_dycore, mbarv_in=mbarv_in(:, jdx, :))
1474 : else
1475 : call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
1476 : kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), &
1477 0 : active_species_idx_dycore=active_species_idx_dycore)
1478 : end if
1479 : end do
1480 :
1481 0 : end subroutine get_molecular_diff_coef_2hd
1482 : !===========================================================================
1483 :
1484 : !***************************************************************************
1485 : !
1486 : ! compute reference vertical profile of density, molecular diffusion and thermal conductivity
1487 : !
1488 : !***************************************************************************
1489 : !
1490 0 : subroutine get_molecular_diff_coef_reference(tref,press,sponge_factor,kmvis_ref,kmcnd_ref,rho_ref)
1491 : use physconst, only: rair
1492 : use air_composition, only: kv1, kv2, kc1, kc2, kv_temp_exp, kc_temp_exp
1493 : ! args
1494 : real(r8), intent(in) :: tref !reference temperature
1495 : real(r8), intent(in) :: press(:) !pressure
1496 : real(r8), intent(in) :: sponge_factor(:) !multiply kmvis and kmcnd with sponge_factor (for sponge layer)
1497 : real(r8), intent(out) :: kmvis_ref(:) !reference molecular diffusion coefficient
1498 : real(r8), intent(out) :: kmcnd_ref(:) !reference thermal conductivity coefficient
1499 : real(r8), intent(out) :: rho_ref(:) !reference density
1500 :
1501 : ! local vars
1502 : integer :: kdx
1503 :
1504 : !--------------------------------------------
1505 : ! Set constants needed for updates
1506 : !--------------------------------------------
1507 :
1508 0 : do kdx = 1, SIZE(press, 1)
1509 0 : rho_ref(kdx) = press(kdx) / (tref * rair) !ideal gas law for dry air
1510 0 : kmvis_ref(kdx) = sponge_factor(kdx) * &
1511 : (kv1 * mmro2 * o2_mwi + &
1512 : kv2 * mmrn2 * n2_mwi) * mbar * &
1513 0 : tref ** kv_temp_exp
1514 0 : kmcnd_ref(kdx) = sponge_factor(kdx) * &
1515 : (kc1 * mmro2 * o2_mwi + &
1516 : kc2 * mmrn2 * n2_mwi) * mbar * &
1517 0 : tref ** kc_temp_exp
1518 : end do
1519 0 : end subroutine get_molecular_diff_coef_reference
1520 :
1521 : !==========================================================================
1522 :
1523 : !
1524 : !***************************************************************************
1525 : !
1526 : ! cam_thermo_calc_kappav: update species dependent kappa for FV dycore
1527 : !
1528 : !***************************************************************************
1529 : !
1530 0 : subroutine cam_thermo_calc_kappav_2hd(tracer, kappav, cpv)
1531 : use air_composition, only: get_R_dry, get_cp_dry
1532 : ! assumes moist MMRs
1533 :
1534 : ! Dummy arguments
1535 : real(r8), intent(in) :: tracer(:, :, :, :)
1536 : real(r8), intent(out) :: kappav(:, :, :)
1537 : real(r8), optional, intent(out) :: cpv(:, :, :)
1538 :
1539 : ! Local variables
1540 0 : real(r8) :: rgas_var(SIZE(tracer, 1), SIZE(tracer, 2), SIZE(tracer, 3))
1541 0 : real(r8) :: cp_var(SIZE(tracer, 1), SIZE(tracer, 2), SIZE(tracer, 3))
1542 : integer :: ind, jnd, knd
1543 :
1544 : !-----------------------------------------------------------------------
1545 : ! Calculate constituent dependent specific heat, gas constant and cappa
1546 : !-----------------------------------------------------------------------
1547 0 : call get_R_dry(tracer, thermodynamic_active_species_idx, rgas_var)
1548 0 : call get_cp_dry(tracer, thermodynamic_active_species_idx, cp_var)
1549 : !$omp parallel do private(ind,jnd,knd)
1550 0 : do knd = 1, SIZE(tracer, 3)
1551 0 : do jnd = 1, SIZE(tracer, 2)
1552 0 : do ind = 1, SIZE(tracer, 1)
1553 0 : kappav(ind,jnd,knd) = rgas_var(ind,jnd,knd) / cp_var(ind,jnd,knd)
1554 : end do
1555 : end do
1556 : end do
1557 :
1558 0 : if (present(cpv)) then
1559 0 : cpv(:,:,:) = cp_var(:,:,:)
1560 : end if
1561 :
1562 0 : end subroutine cam_thermo_calc_kappav_2hd
1563 :
1564 : !===========================================================================
1565 : !
1566 : !***************************************************************************
1567 : !
1568 : ! compute column integrated total energy consistent with vertical
1569 : ! coordinate as well as vertical integrals of water mass (H2O,wv,liq,ice)
1570 : !
1571 : ! if subroutine is asked to compute "te" then the latent heat terms are
1572 : ! added to the kinetic (ke), internal + geopotential (se) energy terms
1573 : !
1574 : ! subroutine assumes that enthalpy term (rho*cp*T) uses dry air heat capacity
1575 : !
1576 : !***************************************************************************
1577 : !
1578 50669136 : subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, &
1579 50669136 : cp_or_cv, U, V, T, vcoord, ptop, phis, z_mid, dycore_idx, qidx, &
1580 50669136 : te, se, po, ke, wv, H2O, liq, ice)
1581 :
1582 : use cam_logfile, only: iulog
1583 : use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure
1584 : use air_composition, only: wv_idx
1585 : use physconst, only: rga, latvap, latice
1586 :
1587 : ! Dummy arguments
1588 : ! tracer: tracer mixing ratio
1589 : !
1590 : ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry
1591 : real(r8), intent(in) :: tracer(:,:,:)
1592 : logical, intent(in) :: moist_mixing_ratio
1593 : ! pdel: pressure level thickness
1594 : real(r8), intent(in) :: pdel_in(:,:)
1595 : ! cp_or_cv: dry air heat capacity under constant pressure or
1596 : ! constant volume (depends on vcoord)
1597 : real(r8), intent(in) :: cp_or_cv(:,:)
1598 : real(r8), intent(in) :: U(:,:)
1599 : real(r8), intent(in) :: V(:,:)
1600 : real(r8), intent(in) :: T(:,:)
1601 : integer, intent(in) :: vcoord ! vertical coordinate
1602 : real(r8), intent(in), optional :: ptop(:)
1603 : real(r8), intent(in), optional :: phis(:)
1604 : real(r8), intent(in), optional :: z_mid(:,:)
1605 : ! dycore_idx: use dycore index for thermodynamic active species
1606 : logical, intent(in), optional :: dycore_idx
1607 : ! qidx: Index of water vapor
1608 : integer, intent(in), optional :: qidx
1609 : ! H2O: vertically integrated total water
1610 : real(r8), intent(out), optional :: H2O(:)
1611 : ! TE: vertically integrated total energy
1612 : real(r8), intent(out), optional :: te (:)
1613 : ! KE: vertically integrated kinetic energy
1614 : real(r8), intent(out), optional :: ke (:)
1615 : ! SE: vertically integrated enthalpy (pressure coordinate)
1616 : ! or internal energy (z coordinate)
1617 : real(r8), intent(out), optional :: se (:)
1618 : ! PO: vertically integrated PHIS term (pressure coordinate)
1619 : ! or potential energy (z coordinate)
1620 : real(r8), intent(out), optional :: po (:)
1621 : ! WV: vertically integrated water vapor
1622 : real(r8), intent(out), optional :: wv (:)
1623 : ! liq: vertically integrated liquid
1624 : real(r8), intent(out), optional :: liq(:)
1625 : ! ice: vertically integrated ice
1626 : real(r8), intent(out), optional :: ice(:)
1627 :
1628 : ! Local variables
1629 101338272 : real(r8) :: ke_vint(SIZE(tracer, 1)) ! Vertical integral of KE
1630 101338272 : real(r8) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of enthalpy or internal energy
1631 101338272 : real(r8) :: po_vint(SIZE(tracer, 1)) ! Vertical integral of PHIS or potential energy
1632 101338272 : real(r8) :: wv_vint(SIZE(tracer, 1)) ! Vertical integral of wv
1633 101338272 : real(r8) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq
1634 101338272 : real(r8) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice
1635 101338272 : real(r8) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) !moist pressure level thickness
1636 : real(r8) :: latsub ! latent heat of sublimation
1637 :
1638 : integer :: ierr
1639 : integer :: kdx, idx ! coord indices
1640 : integer :: qdx ! tracer index
1641 : integer :: wvidx ! water vapor index
1642 50669136 : integer, allocatable :: species_idx(:)
1643 50669136 : integer, allocatable :: species_liq_idx(:)
1644 50669136 : integer, allocatable :: species_ice_idx(:)
1645 : character(len=*), parameter :: subname = 'get_hydrostatic_energy'
1646 :
1647 152007408 : allocate(species_idx(thermodynamic_active_species_num), stat=ierr)
1648 50669136 : if ( ierr /= 0 ) then
1649 0 : call endrun(subname//': allocation error for species_idx array')
1650 : end if
1651 152007408 : allocate(species_liq_idx(thermodynamic_active_species_liq_num), stat=ierr)
1652 50669136 : if ( ierr /= 0 ) then
1653 0 : call endrun(subname//': allocation error for species_liq_idx array')
1654 : end if
1655 152007408 : allocate(species_ice_idx(thermodynamic_active_species_ice_num), stat=ierr)
1656 50669136 : if ( ierr /= 0 ) then
1657 0 : call endrun(subname//': allocation error for species_ice_idx array')
1658 : end if
1659 :
1660 50669136 : if (present(dycore_idx))then
1661 0 : if (dycore_idx) then
1662 0 : species_idx(:) = thermodynamic_active_species_idx_dycore(:)
1663 0 : species_liq_idx(:) = thermodynamic_active_species_liq_idx_dycore(:)
1664 0 : species_ice_idx(:) = thermodynamic_active_species_ice_idx_dycore(:)
1665 : else
1666 0 : species_idx(:) = thermodynamic_active_species_idx(:)
1667 0 : species_liq_idx(:) = thermodynamic_active_species_liq_idx(:)
1668 0 : species_ice_idx(:) = thermodynamic_active_species_ice_idx(:)
1669 : end if
1670 : else
1671 354683952 : species_idx(:) = thermodynamic_active_species_idx(:)
1672 152007408 : species_liq_idx(:) = thermodynamic_active_species_liq_idx(:)
1673 202676544 : species_ice_idx(:) = thermodynamic_active_species_ice_idx(:)
1674 : end if
1675 :
1676 50669136 : if (present(qidx)) then
1677 0 : wvidx = qidx
1678 : else
1679 50669136 : wvidx = wv_idx
1680 : end if
1681 :
1682 50669136 : if (moist_mixing_ratio) then
1683 78733945584 : pdel = pdel_in
1684 : else
1685 0 : pdel = pdel_in
1686 0 : do qdx = dry_air_species_num+1, thermodynamic_active_species_num
1687 0 : pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx))
1688 : end do
1689 : end if
1690 :
1691 846056736 : ke_vint = 0._r8
1692 846056736 : se_vint = 0._r8
1693 101338272 : select case (vcoord)
1694 : case(vc_moist_pressure, vc_dry_pressure)
1695 50669136 : if (.not. present(ptop).or. (.not. present(phis))) then
1696 0 : write(iulog, *) subname, ' ptop and phis must be present for ', &
1697 0 : 'moist/dry pressure vertical coordinate'
1698 : call endrun(subname//': ptop and phis must be present for '// &
1699 0 : 'moist/dry pressure vertical coordinate')
1700 : end if
1701 846056736 : po_vint = ptop
1702 4762898784 : do kdx = 1, SIZE(tracer, 2)
1703 78733945584 : do idx = 1, SIZE(tracer, 1)
1704 >14794*10^7 : ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * &
1705 >22191*10^7 : 0.5_r8 * (U(idx, kdx)**2 + V(idx, kdx)**2)) * rga
1706 73971046800 : se_vint(idx) = se_vint(idx) + (T(idx, kdx) * &
1707 73971046800 : cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga)
1708 78683276448 : po_vint(idx) = po_vint(idx)+pdel(idx, kdx)
1709 :
1710 : end do
1711 : end do
1712 896725872 : do idx = 1, SIZE(tracer, 1)
1713 846056736 : po_vint(idx) = (phis(idx) * po_vint(idx) * rga)
1714 : end do
1715 : case(vc_height)
1716 0 : if (.not. present(phis)) then
1717 0 : write(iulog, *) subname, ' phis must be present for ', &
1718 0 : 'heigt-based vertical coordinate'
1719 : call endrun(subname//': phis must be present for '// &
1720 0 : 'height-based vertical coordinate')
1721 : end if
1722 0 : po_vint = 0._r8
1723 0 : do kdx = 1, SIZE(tracer, 2)
1724 0 : do idx = 1, SIZE(tracer, 1)
1725 0 : ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * &
1726 0 : 0.5_r8 * (U(idx, kdx)**2 + V(idx, kdx)**2) * rga)
1727 0 : se_vint(idx) = se_vint(idx) + (T(idx, kdx) * &
1728 0 : cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga)
1729 : ! z_mid is height above ground
1730 0 : po_vint(idx) = po_vint(idx) + (z_mid(idx, kdx) + &
1731 0 : phis(idx) * rga) * pdel(idx, kdx)
1732 : end do
1733 : end do
1734 : case default
1735 0 : write(iulog, *) subname, ' vertical coordinate not supported: ', vcoord
1736 50669136 : call endrun(subname//': vertical coordinate not supported')
1737 : end select
1738 50669136 : if (present(te)) then
1739 846056736 : te = se_vint + po_vint+ ke_vint
1740 : end if
1741 50669136 : if (present(se)) then
1742 398059200 : se = se_vint
1743 : end if
1744 50669136 : if (present(po)) then
1745 398059200 : po = po_vint
1746 : end if
1747 50669136 : if (present(ke)) then
1748 398059200 : ke = ke_vint
1749 : end if
1750 : !
1751 : ! vertical integral of total liquid water
1752 : !
1753 50669136 : if (.not.moist_mixing_ratio) then
1754 0 : pdel = pdel_in! set pseudo density to dry
1755 : end if
1756 :
1757 846056736 : wv_vint = 0._r8
1758 4762898784 : do kdx = 1, SIZE(tracer, 2)
1759 78733945584 : do idx = 1, SIZE(tracer, 1)
1760 >22191*10^7 : wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * &
1761 >30059*10^7 : pdel(idx, kdx) * rga)
1762 : end do
1763 : end do
1764 424889136 : if (present(wv)) wv = wv_vint
1765 :
1766 846056736 : liq_vint = 0._r8
1767 152007408 : do qdx = 1, thermodynamic_active_species_liq_num
1768 9576466704 : do kdx = 1, SIZE(tracer, 2)
1769 >15746*10^7 : do idx = 1, SIZE(tracer, 1)
1770 >29588*10^7 : liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * &
1771 >45325*10^7 : tracer(idx, kdx, species_liq_idx(qdx)) * rga)
1772 : end do
1773 : end do
1774 : end do
1775 424889136 : if (present(liq)) liq = liq_vint
1776 :
1777 : !
1778 : ! vertical integral of total frozen (ice) water
1779 : !
1780 846056736 : ice_vint = 0._r8
1781 202676544 : do qdx = 1, thermodynamic_active_species_ice_num
1782 14339365488 : do kdx = 1, SIZE(tracer, 2)
1783 >23620*10^7 : do idx = 1, SIZE(tracer, 1)
1784 >44382*10^7 : ice_vint(idx) = ice_vint(idx) + (pdel(idx, kdx) * &
1785 >67987*10^7 : tracer(idx, kdx, species_ice_idx(qdx)) * rga)
1786 : end do
1787 : end do
1788 : end do
1789 424889136 : if (present(ice)) ice = ice_vint
1790 : ! Compute vertical integrals of total water.
1791 50669136 : if (present(H2O)) then
1792 846056736 : H2O = wv_vint + liq_vint + ice_vint
1793 : end if
1794 : !
1795 : ! latent heat terms depend on enthalpy reference state
1796 : !
1797 50669136 : latsub = latvap + latice
1798 50669136 : if (present(te)) then
1799 101338272 : select case (TRIM(enthalpy_reference_state))
1800 : case('ice')
1801 846056736 : te = te + (latsub * wv_vint) + (latice * liq_vint)
1802 : case('liq')
1803 0 : te = te + (latvap * wv_vint) - (latice * ice_vint)
1804 : case('wv')
1805 0 : te = te - (latvap * liq_vint) - (latsub * ice_vint)
1806 : case default
1807 0 : write(iulog, *) subname, ' enthalpy reference state not ', &
1808 0 : 'supported: ', TRIM(enthalpy_reference_state)
1809 101338272 : call endrun(subname//': enthalpy reference state not supported')
1810 : end select
1811 : end if
1812 50669136 : deallocate(species_idx, species_liq_idx, species_ice_idx)
1813 50669136 : end subroutine get_hydrostatic_energy_1hd
1814 :
1815 : end module cam_thermo
|
