Line data Source code
1 :
2 : module check_energy
3 :
4 : !---------------------------------------------------------------------------------
5 : ! Purpose:
6 : !
7 : ! Module to check
8 : ! 1. vertically integrated total energy and water conservation for each
9 : ! column within the physical parameterizations
10 : !
11 : ! 2. global mean total energy conservation between the physics output state
12 : ! and the input state on the next time step.
13 : !
14 : ! 3. add a globally uniform heating term to account for any change of total energy in 2.
15 : !
16 : ! Author: Byron Boville Oct 31, 2002
17 : !
18 : ! Modifications:
19 : ! 03.03.29 Boville Add global energy check and fixer.
20 : !
21 : !---------------------------------------------------------------------------------
22 :
23 : use shr_kind_mod, only: r8 => shr_kind_r8
24 : use ppgrid, only: pcols, pver, begchunk, endchunk
25 : use spmd_utils, only: masterproc
26 :
27 : use gmean_mod, only: gmean
28 : use physconst, only: gravit, rga, latvap, latice, cpair, rair
29 : use air_composition, only: cpairv, rairv, cp_or_cv_dycore
30 : use physics_types, only: physics_state, physics_tend, physics_ptend, physics_ptend_init
31 : use constituents, only: cnst_get_ind, pcnst, cnst_name, cnst_get_type_byind
32 : use time_manager, only: is_first_step
33 : use cam_logfile, only: iulog
34 : use scamMod, only: single_column, use_camiop, heat_glob_scm
35 : use cam_history, only: outfld, write_camiop
36 :
37 : implicit none
38 : private
39 :
40 : ! Public types:
41 : public check_tracers_data
42 :
43 : ! Public methods
44 : public :: check_energy_readnl ! read namelist values
45 : public :: check_energy_register ! register fields in physics buffer
46 : public :: check_energy_get_integrals ! get energy integrals computed in check_energy_gmean
47 : public :: check_energy_init ! initialization of module
48 : public :: check_energy_timestep_init ! timestep initialization of energy integrals and cumulative boundary fluxes
49 : public :: check_energy_chng ! check changes in integrals against cumulative boundary fluxes
50 : public :: check_energy_gmean ! global means of physics input and output total energy
51 : public :: check_energy_fix ! add global mean energy difference as a heating
52 : public :: check_tracers_init ! initialize tracer integrals and cumulative boundary fluxes
53 : public :: check_tracers_chng ! check changes in integrals against cumulative boundary fluxes
54 :
55 : public :: tot_energy_phys ! calculate and output total energy and axial angular momentum diagnostics
56 :
57 : ! Private module data
58 :
59 : logical :: print_energy_errors = .false.
60 :
61 : real(r8) :: teout_glob ! global mean energy of output state
62 : real(r8) :: teinp_glob ! global mean energy of input state
63 : real(r8) :: tedif_glob ! global mean energy difference
64 : real(r8) :: psurf_glob ! global mean surface pressure
65 : real(r8) :: ptopb_glob ! global mean top boundary pressure
66 : real(r8) :: heat_glob ! global mean heating rate
67 :
68 : ! Physics buffer indices
69 :
70 : integer :: teout_idx = 0 ! teout index in physics buffer
71 : integer :: dtcore_idx = 0 ! dtcore index in physics buffer
72 : integer :: dqcore_idx = 0 ! dqcore index in physics buffer
73 : integer :: ducore_idx = 0 ! ducore index in physics buffer
74 : integer :: dvcore_idx = 0 ! dvcore index in physics buffer
75 :
76 : type check_tracers_data
77 : real(r8) :: tracer(pcols,pcnst) ! initial vertically integrated total (kinetic + static) energy
78 : real(r8) :: tracer_tnd(pcols,pcnst) ! cumulative boundary flux of total energy
79 : integer :: count(pcnst) ! count of values with significant imbalances
80 : end type check_tracers_data
81 :
82 :
83 : !===============================================================================
84 : contains
85 : !===============================================================================
86 :
87 1536 : subroutine check_energy_readnl(nlfile)
88 :
89 : use namelist_utils, only: find_group_name
90 : use units, only: getunit, freeunit
91 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_logical
92 : use cam_abortutils, only: endrun
93 :
94 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
95 :
96 : ! Local variables
97 : integer :: unitn, ierr
98 : character(len=*), parameter :: sub = 'check_energy_readnl'
99 :
100 : namelist /check_energy_nl/ print_energy_errors
101 : !-----------------------------------------------------------------------------
102 :
103 : ! Read namelist
104 1536 : if (masterproc) then
105 2 : unitn = getunit()
106 2 : open( unitn, file=trim(nlfile), status='old' )
107 2 : call find_group_name(unitn, 'check_energy_nl', status=ierr)
108 2 : if (ierr == 0) then
109 0 : read(unitn, check_energy_nl, iostat=ierr)
110 0 : if (ierr /= 0) then
111 0 : call endrun(sub//': FATAL: reading namelist')
112 : end if
113 : end if
114 2 : close(unitn)
115 2 : call freeunit(unitn)
116 : end if
117 :
118 1536 : call mpi_bcast(print_energy_errors, 1, mpi_logical, mstrid, mpicom, ierr)
119 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: print_energy_errors")
120 :
121 1536 : if (masterproc) then
122 2 : write(iulog,*) 'check_energy options:'
123 2 : write(iulog,*) ' print_energy_errors =', print_energy_errors
124 : end if
125 :
126 1536 : end subroutine check_energy_readnl
127 :
128 : !===============================================================================
129 :
130 1536 : subroutine check_energy_register()
131 : !
132 : ! Register fields in the physics buffer.
133 : !
134 : !-----------------------------------------------------------------------
135 :
136 : use physics_buffer, only : pbuf_add_field, dtype_r8, dyn_time_lvls
137 : use physics_buffer, only : pbuf_register_subcol
138 : use subcol_utils, only : is_subcol_on
139 :
140 : !-----------------------------------------------------------------------
141 :
142 : ! Request physics buffer space for fields that persist across timesteps.
143 :
144 4608 : call pbuf_add_field('TEOUT', 'global',dtype_r8 , (/pcols,dyn_time_lvls/), teout_idx)
145 6144 : call pbuf_add_field('DTCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),dtcore_idx)
146 : ! DQCORE refers to dycore tendency of water vapor
147 6144 : call pbuf_add_field('DQCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),dqcore_idx)
148 6144 : call pbuf_add_field('DUCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),ducore_idx)
149 6144 : call pbuf_add_field('DVCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),dvcore_idx)
150 1536 : if(is_subcol_on()) then
151 0 : call pbuf_register_subcol('TEOUT', 'phys_register', teout_idx)
152 0 : call pbuf_register_subcol('DTCORE', 'phys_register', dtcore_idx)
153 0 : call pbuf_register_subcol('DQCORE', 'phys_register', dqcore_idx)
154 0 : call pbuf_register_subcol('DUCORE', 'phys_register', ducore_idx)
155 0 : call pbuf_register_subcol('DVCORE', 'phys_register', dvcore_idx)
156 : end if
157 :
158 1536 : end subroutine check_energy_register
159 :
160 : !===============================================================================
161 :
162 1489176 : subroutine check_energy_get_integrals( tedif_glob_out, heat_glob_out )
163 :
164 : !-----------------------------------------------------------------------
165 : ! Purpose: Return energy integrals
166 : !-----------------------------------------------------------------------
167 :
168 : real(r8), intent(out), optional :: tedif_glob_out
169 : real(r8), intent(out), optional :: heat_glob_out
170 :
171 : !-----------------------------------------------------------------------
172 :
173 1489176 : if ( present(tedif_glob_out) ) then
174 0 : tedif_glob_out = tedif_glob
175 : endif
176 1489176 : if ( present(heat_glob_out) ) then
177 1489176 : heat_glob_out = heat_glob
178 : endif
179 :
180 1536 : end subroutine check_energy_get_integrals
181 :
182 : !================================================================================================
183 :
184 1536 : subroutine check_energy_init()
185 : !
186 : ! Initialize the energy conservation module
187 : !
188 : !-----------------------------------------------------------------------
189 : use cam_history, only: addfld, add_default, horiz_only
190 : use phys_control, only: phys_getopts
191 :
192 : implicit none
193 :
194 : logical :: history_budget, history_waccm
195 : integer :: history_budget_histfile_num ! output history file number for budget fields
196 :
197 : !-----------------------------------------------------------------------
198 :
199 : call phys_getopts( history_budget_out = history_budget, &
200 : history_budget_histfile_num_out = history_budget_histfile_num, &
201 1536 : history_waccm_out = history_waccm )
202 :
203 : ! register history variables
204 1536 : call addfld('TEINP', horiz_only, 'A', 'J/m2', 'Total energy of physics input')
205 1536 : call addfld('TEOUT', horiz_only, 'A', 'J/m2', 'Total energy of physics output')
206 1536 : call addfld('TEFIX', horiz_only, 'A', 'J/m2', 'Total energy after fixer')
207 1536 : call addfld('EFIX', horiz_only, 'A', 'W/m2', 'Effective sensible heat flux due to energy fixer')
208 3072 : call addfld('DTCORE', (/ 'lev' /), 'A', 'K/s' , 'T tendency due to dynamical core')
209 3072 : call addfld('DQCORE', (/ 'lev' /), 'A', 'kg/kg/s' , 'Water vapor tendency due to dynamical core')
210 :
211 1536 : if ( history_budget ) then
212 0 : call add_default ('DTCORE', history_budget_histfile_num, ' ')
213 : end if
214 1536 : if ( history_waccm ) then
215 0 : call add_default ('DTCORE', 1, ' ')
216 : end if
217 :
218 1536 : end subroutine check_energy_init
219 :
220 : !===============================================================================
221 :
222 1495368 : subroutine check_energy_timestep_init(state, tend, pbuf, col_type)
223 1536 : use cam_thermo, only: get_hydrostatic_energy
224 : use physics_buffer, only: physics_buffer_desc, pbuf_set_field
225 : use cam_abortutils, only: endrun
226 : use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure
227 : use physics_types, only: phys_te_idx, dyn_te_idx
228 : !-----------------------------------------------------------------------
229 : ! Compute initial values of energy and water integrals,
230 : ! zero cumulative tendencies
231 : !-----------------------------------------------------------------------
232 : !------------------------------Arguments--------------------------------
233 :
234 : type(physics_state), intent(inout) :: state
235 : type(physics_tend ), intent(inout) :: tend
236 : type(physics_buffer_desc), pointer :: pbuf(:)
237 : integer, optional :: col_type ! Flag inidicating whether using grid or subcolumns
238 : !---------------------------Local storage-------------------------------
239 2990736 : real(r8) :: cp_or_cv(state%psetcols,pver)
240 : integer lchnk ! chunk identifier
241 : integer ncol ! number of atmospheric columns
242 : !-----------------------------------------------------------------------
243 :
244 1495368 : lchnk = state%lchnk
245 1495368 : ncol = state%ncol
246 :
247 : ! cp_or_cv needs to be allocated to a size which matches state and ptend
248 : ! If psetcols == pcols, cpairv is the correct size and just copy into cp_or_cv
249 : ! If psetcols > pcols and all cpairv match cpair, then assign the constant cpair
250 :
251 1495368 : if (state%psetcols == pcols) then
252 662448024 : cp_or_cv(:,:) = cpairv(:,:,lchnk)
253 0 : else if (state%psetcols > pcols .and. all(cpairv(:,:,lchnk) == cpair)) then
254 0 : cp_or_cv(1:ncol,:) = cpair
255 : else
256 0 : call endrun('check_energy_timestep_init: cpairv is not allowed to vary when subcolumns are turned on')
257 : end if
258 : !
259 : ! CAM physics total energy
260 : !
261 4486104 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
262 4486104 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
263 8972208 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), &
264 2990736 : vc_physics, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol),&
265 22430520 : te = state%te_ini(1:ncol,phys_te_idx), H2O = state%tw_ini(1:ncol))
266 : !
267 : ! Dynamical core total energy
268 : !
269 650693736 : state%temp_ini(:ncol,:) = state%T(:ncol,:)
270 650693736 : state%z_ini(:ncol,:) = state%zm(:ncol,:)
271 1495368 : if (vc_dycore == vc_height) then
272 : !
273 : ! MPAS specific hydrostatic energy computation (internal energy)
274 : !
275 0 : if (state%psetcols == pcols) then
276 0 : cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk)
277 : else
278 0 : cp_or_cv(:ncol,:) = cpair-rair
279 : endif
280 :
281 0 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
282 0 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
283 0 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), &
284 0 : vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), &
285 0 : z_mid = state%z_ini(1:ncol,:), &
286 0 : te = state%te_ini(1:ncol,dyn_te_idx), H2O = state%tw_ini(1:ncol))
287 1495368 : else if (vc_dycore == vc_dry_pressure) then
288 : !
289 : ! SE specific hydrostatic energy (enthalpy)
290 : !
291 1495368 : if (state%psetcols == pcols) then
292 650693736 : cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk)
293 : else
294 0 : cp_or_cv(:ncol,:) = cpair
295 : endif
296 4486104 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
297 2990736 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
298 8972208 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), &
299 2990736 : vc_dry_pressure, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), &
300 20935152 : te = state%te_ini(1:ncol,dyn_te_idx), H2O = state%tw_ini(1:ncol))
301 : else
302 : !
303 : ! dycore energy is the same as physics
304 : !
305 0 : state%te_ini(1:ncol,dyn_te_idx) = state%te_ini(1:ncol,phys_te_idx)
306 : end if
307 51433704 : state%te_cur(:ncol,:) = state%te_ini(:ncol,:)
308 24969168 : state%tw_cur(:ncol) = state%tw_ini(:ncol)
309 :
310 : ! zero cummulative boundary fluxes
311 26464536 : tend%te_tnd(:ncol) = 0._r8
312 26464536 : tend%tw_tnd(:ncol) = 0._r8
313 :
314 1495368 : state%count = 0
315 :
316 : ! initialize physics buffer
317 1495368 : if (is_first_step()) then
318 3096 : call pbuf_set_field(pbuf, teout_idx, state%te_ini(:,dyn_te_idx), col_type=col_type)
319 : end if
320 :
321 1495368 : end subroutine check_energy_timestep_init
322 :
323 : !===============================================================================
324 :
325 16411896 : subroutine check_energy_chng(state, tend, name, nstep, ztodt, &
326 16411896 : flx_vap, flx_cnd, flx_ice, flx_sen)
327 1495368 : use cam_thermo, only: get_hydrostatic_energy
328 : use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure
329 : use cam_abortutils, only: endrun
330 : use physics_types, only: phys_te_idx, dyn_te_idx
331 : !-----------------------------------------------------------------------
332 : ! Check that the energy and water change matches the boundary fluxes
333 : !-----------------------------------------------------------------------
334 : !------------------------------Arguments--------------------------------
335 :
336 : type(physics_state) , intent(inout) :: state
337 : type(physics_tend ) , intent(inout) :: tend
338 : character*(*),intent(in) :: name ! parameterization name for fluxes
339 : integer , intent(in ) :: nstep ! current timestep number
340 : real(r8), intent(in ) :: ztodt ! 2 delta t (model time increment)
341 : real(r8), intent(in ) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s)
342 : real(r8), intent(in ) :: flx_cnd(:) ! (pcols) -boundary flux of liquid+ice (m/s) (precip?)
343 : real(r8), intent(in ) :: flx_ice(:) ! (pcols) -boundary flux of ice (m/s) (snow?)
344 : real(r8), intent(in ) :: flx_sen(:) ! (pcols) -boundary flux of sensible heat (w/m2)
345 :
346 : !******************** BAB ******************************************************
347 : !******* Note that the precip and ice fluxes are in precip units (m/s). ********
348 : !******* I would prefer to have kg/m2/s. ********
349 : !******* I would also prefer liquid (not total) and ice fluxes ********
350 : !*******************************************************************************
351 :
352 : !---------------------------Local storage-------------------------------
353 :
354 32823792 : real(r8) :: te_xpd(state%ncol) ! expected value (f0 + dt*boundary_flux)
355 32823792 : real(r8) :: te_dif(state%ncol) ! energy of input state - original energy
356 32823792 : real(r8) :: te_tnd(state%ncol) ! tendency from last process
357 32823792 : real(r8) :: te_rer(state%ncol) ! relative error in energy column
358 :
359 32823792 : real(r8) :: tw_xpd(state%ncol) ! expected value (w0 + dt*boundary_flux)
360 32823792 : real(r8) :: tw_dif(state%ncol) ! tw_inp - original water
361 32823792 : real(r8) :: tw_tnd(state%ncol) ! tendency from last process
362 32823792 : real(r8) :: tw_rer(state%ncol) ! relative error in water column
363 :
364 32823792 : real(r8) :: te(state%ncol) ! vertical integral of total energy
365 32823792 : real(r8) :: tw(state%ncol) ! vertical integral of total water
366 32823792 : real(r8) :: cp_or_cv(state%psetcols,pver) ! cp or cv depending on vcoord
367 32823792 : real(r8) :: scaling(state%psetcols,pver) ! scaling for conversion of temperature increment
368 32823792 : real(r8) :: temp(state%ncol,pver) ! temperature
369 :
370 32823792 : real(r8) :: se(state%ncol) ! enthalpy or internal energy (J/m2)
371 32823792 : real(r8) :: po(state%ncol) ! surface potential or potential energy (J/m2)
372 32823792 : real(r8) :: ke(state%ncol) ! kinetic energy (J/m2)
373 32823792 : real(r8) :: wv(state%ncol) ! column integrated vapor (kg/m2)
374 32823792 : real(r8) :: liq(state%ncol) ! column integrated liquid (kg/m2)
375 32823792 : real(r8) :: ice(state%ncol) ! column integrated ice (kg/m2)
376 :
377 : integer lchnk ! chunk identifier
378 : integer ncol ! number of atmospheric columns
379 : integer i ! column index
380 : !-----------------------------------------------------------------------
381 :
382 16411896 : lchnk = state%lchnk
383 16411896 : ncol = state%ncol
384 :
385 : ! If psetcols == pcols, cpairv is the correct size and just copy into cp_or_cv
386 : ! If psetcols > pcols and all cpairv match cpair, then assign the constant cpair
387 :
388 16411896 : if (state%psetcols == pcols) then
389 7270469928 : cp_or_cv(:,:) = cpairv(:,:,lchnk)
390 0 : else if (state%psetcols > pcols .and. all(cpairv(:,:,:) == cpair)) then
391 0 : cp_or_cv(:,:) = cpair
392 : else
393 0 : call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on')
394 : end if
395 :
396 49235688 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
397 49235688 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
398 98471376 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), &
399 32823792 : vc_physics, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), &
400 : te = te(1:ncol), H2O = tw(1:ncol), se=se(1:ncol),po=po(1:ncol), &
401 246178440 : ke=ke(1:ncol),wv=wv(1:ncol),liq=liq(1:ncol),ice=ice(1:ncol))
402 : ! compute expected values and tendencies
403 274040496 : do i = 1, ncol
404 : ! change in static energy and total water
405 257628600 : te_dif(i) = te(i) - state%te_cur(i,phys_te_idx)
406 257628600 : tw_dif(i) = tw(i) - state%tw_cur(i)
407 :
408 : ! expected tendencies from boundary fluxes for last process
409 257628600 : te_tnd(i) = flx_vap(i)*(latvap+latice) - (flx_cnd(i) - flx_ice(i))*1000._r8*latice + flx_sen(i)
410 257628600 : tw_tnd(i) = flx_vap(i) - flx_cnd(i) *1000._r8
411 :
412 : ! cummulative tendencies from boundary fluxes
413 257628600 : tend%te_tnd(i) = tend%te_tnd(i) + te_tnd(i)
414 257628600 : tend%tw_tnd(i) = tend%tw_tnd(i) + tw_tnd(i)
415 :
416 : ! expected new values from previous state plus boundary fluxes
417 257628600 : te_xpd(i) = state%te_cur(i,phys_te_idx) + te_tnd(i)*ztodt
418 257628600 : tw_xpd(i) = state%tw_cur(i) + tw_tnd(i)*ztodt
419 :
420 : ! relative error, expected value - input state / previous state
421 274040496 : te_rer(i) = (te_xpd(i) - te(i)) / state%te_cur(i,phys_te_idx)
422 : end do
423 :
424 : ! relative error for total water (allow for dry atmosphere)
425 274040496 : tw_rer = 0._r8
426 274040496 : where (state%tw_cur(:ncol) > 0._r8)
427 16411896 : tw_rer(:ncol) = (tw_xpd(:ncol) - tw(:ncol)) / state%tw_cur(:ncol)
428 : end where
429 :
430 : ! error checking
431 16411896 : if (print_energy_errors) then
432 0 : if (any(abs(te_rer(1:ncol)) > 1.E-14_r8 .or. abs(tw_rer(1:ncol)) > 1.E-10_r8)) then
433 0 : do i = 1, ncol
434 : ! the relative error threshold for the water budget has been reduced to 1.e-10
435 : ! to avoid messages generated by QNEG3 calls
436 : ! PJR- change to identify if error in energy or water
437 0 : if (abs(te_rer(i)) > 1.E-14_r8 ) then
438 0 : state%count = state%count + 1
439 0 : write(iulog,*) "significant energy conservation error after ", name, &
440 0 : " count", state%count, " nstep", nstep, "chunk", lchnk, "col", i
441 0 : write(iulog,*) te(i),te_xpd(i),te_dif(i),tend%te_tnd(i)*ztodt, &
442 0 : te_tnd(i)*ztodt,te_rer(i)
443 : endif
444 0 : if ( abs(tw_rer(i)) > 1.E-10_r8) then
445 0 : state%count = state%count + 1
446 0 : write(iulog,*) "significant water conservation error after ", name, &
447 0 : " count", state%count, " nstep", nstep, "chunk", lchnk, "col", i
448 0 : write(iulog,*) tw(i),tw_xpd(i),tw_dif(i),tend%tw_tnd(i)*ztodt, &
449 0 : tw_tnd(i)*ztodt,tw_rer(i)
450 : end if
451 : end do
452 : end if
453 : end if
454 :
455 : ! copy new value to state
456 :
457 274040496 : do i = 1, ncol
458 257628600 : state%te_cur(i,phys_te_idx) = te(i)
459 274040496 : state%tw_cur(i) = tw(i)
460 : end do
461 :
462 : !
463 : ! Dynamical core total energy
464 : !
465 16411896 : if (vc_dycore == vc_height) then
466 : !
467 : ! compute cv if vertical coordinate is height: cv = cp - R
468 : !
469 : ! Note: cp_or_cv set above for pressure coordinate
470 0 : if (state%psetcols == pcols) then
471 0 : cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk)
472 : else
473 0 : cp_or_cv(:ncol,:) = cpair-rair
474 : endif
475 0 : scaling(:,:) = cpairv(:,:,lchnk)/cp_or_cv(:,:) !cp/cv scaling
476 0 : temp(1:ncol,:) = state%temp_ini(1:ncol,:)+scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:))
477 0 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
478 0 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
479 0 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), temp(1:ncol,1:pver), &
480 0 : vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), &
481 0 : z_mid = state%z_ini(1:ncol,:), &
482 0 : te = state%te_cur(1:ncol,dyn_te_idx), H2O = state%tw_cur(1:ncol))
483 16411896 : else if (vc_dycore == vc_dry_pressure) then
484 : !
485 : ! SE specific hydrostatic energy
486 : !
487 16411896 : if (state%psetcols == pcols) then
488 7141464792 : cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk)
489 7141464792 : scaling(:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk)
490 : else
491 0 : cp_or_cv(:ncol,:) = cpair
492 0 : scaling(:ncol,:) = 1.0_r8
493 : endif
494 : !
495 : ! enthalpy scaling for energy consistency
496 : !
497 7141464792 : temp(1:ncol,:) = state%temp_ini(1:ncol,:)+scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:))
498 49235688 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
499 32823792 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
500 65647584 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), temp(1:ncol,1:pver), &
501 32823792 : vc_dry_pressure, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), &
502 196942752 : te = state%te_cur(1:ncol,dyn_te_idx), H2O = state%tw_cur(1:ncol))
503 : else
504 0 : state%te_cur(1:ncol,dyn_te_idx) = te(1:ncol)
505 : end if
506 16411896 : end subroutine check_energy_chng
507 :
508 :
509 741888 : subroutine check_energy_gmean(state, pbuf2d, dtime, nstep)
510 :
511 16411896 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk
512 : use physics_types, only: dyn_te_idx
513 : use cam_history, only: write_camiop
514 : !-----------------------------------------------------------------------
515 : ! Compute global mean total energy of physics input and output states
516 : ! computed consistently with dynamical core vertical coordinate
517 : ! (under hydrostatic assumption)
518 : !-----------------------------------------------------------------------
519 : !------------------------------Arguments--------------------------------
520 :
521 : type(physics_state), intent(in ), dimension(begchunk:endchunk) :: state
522 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
523 :
524 : real(r8), intent(in) :: dtime ! physics time step
525 : integer , intent(in) :: nstep ! current timestep number
526 :
527 : !---------------------------Local storage-------------------------------
528 : integer :: ncol ! number of active columns
529 : integer :: lchnk ! chunk index
530 :
531 741888 : real(r8) :: te(pcols,begchunk:endchunk,4)
532 : ! total energy of input/output states (copy)
533 : real(r8) :: te_glob(4) ! global means of total energy
534 370944 : real(r8), pointer :: teout(:)
535 : !-----------------------------------------------------------------------
536 :
537 : ! Copy total energy out of input and output states
538 1866312 : do lchnk = begchunk, endchunk
539 1495368 : ncol = state(lchnk)%ncol
540 : ! input energy using dynamical core energy formula
541 24969168 : te(:ncol,lchnk,1) = state(lchnk)%te_ini(:ncol,dyn_te_idx)
542 : ! output energy
543 1495368 : call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),teout_idx, teout)
544 :
545 24969168 : te(:ncol,lchnk,2) = teout(1:ncol)
546 : ! surface pressure for heating rate
547 24969168 : te(:ncol,lchnk,3) = state(lchnk)%pint(:ncol,pver+1)
548 : ! model top pressure for heating rate (not constant for z-based vertical coordinate!)
549 25340112 : te(:ncol,lchnk,4) = state(lchnk)%pint(:ncol,1)
550 : end do
551 :
552 : ! Compute global means of input and output energies and of
553 : ! surface pressure for heating rate (assume uniform ptop)
554 370944 : call gmean(te, te_glob, 4)
555 :
556 370944 : if (begchunk .le. endchunk) then
557 370944 : teinp_glob = te_glob(1)
558 370944 : teout_glob = te_glob(2)
559 370944 : psurf_glob = te_glob(3)
560 370944 : ptopb_glob = te_glob(4)
561 :
562 : ! Global mean total energy difference
563 370944 : tedif_glob = teinp_glob - teout_glob
564 370944 : heat_glob = -tedif_glob/dtime * gravit / (psurf_glob - ptopb_glob)
565 370944 : if (masterproc) then
566 483 : write(iulog,'(1x,a9,1x,i8,5(1x,e25.17))') "nstep, te", nstep, teinp_glob, teout_glob, &
567 966 : heat_glob, psurf_glob, ptopb_glob
568 : end if
569 : else
570 0 : heat_glob = 0._r8
571 : end if ! (begchunk .le. endchunk)
572 :
573 741888 : end subroutine check_energy_gmean
574 :
575 : !===============================================================================
576 5981472 : subroutine check_energy_fix(state, ptend, nstep, eshflx)
577 :
578 : !-----------------------------------------------------------------------
579 : ! Add heating rate required for global mean total energy conservation
580 : !-----------------------------------------------------------------------
581 : !------------------------------Arguments--------------------------------
582 :
583 : type(physics_state), intent(in ) :: state
584 : type(physics_ptend), intent(out) :: ptend
585 :
586 : integer , intent(in ) :: nstep ! time step number
587 : real(r8), intent(out ) :: eshflx(pcols) ! effective sensible heat flux
588 :
589 : !---------------------------Local storage-------------------------------
590 : integer :: i ! column
591 : integer :: ncol ! number of atmospheric columns in chunk
592 : integer :: lchnk ! chunk number
593 : real(r8) :: heat_out(pcols)
594 : !-----------------------------------------------------------------------
595 1495368 : lchnk = state%lchnk
596 1495368 : ncol = state%ncol
597 :
598 1495368 : call physics_ptend_init(ptend, state%psetcols, 'chkenergyfix', ls=.true.)
599 :
600 : #if ( defined OFFLINE_DYN )
601 : ! disable the energy fix for offline driver
602 : heat_glob = 0._r8
603 : #endif
604 :
605 : ! Special handling of energy fix for SCAM - supplied via CAMIOP - zero's for normal IOPs
606 1495368 : if (single_column) then
607 0 : if ( use_camiop) then
608 0 : heat_glob = heat_glob_scm(1)
609 : else
610 0 : heat_glob = 0._r8
611 : endif
612 : endif
613 650693736 : ptend%s(:ncol,:pver) = heat_glob
614 :
615 1495368 : if (nstep > 0 .and. write_camiop) then
616 0 : heat_out(:ncol) = heat_glob
617 0 : call outfld('heat_glob', heat_out(:ncol), pcols, lchnk)
618 : endif
619 :
620 : ! compute effective sensible heat flux
621 24969168 : do i = 1, ncol
622 24969168 : eshflx(i) = heat_glob * (state%pint(i,pver+1) - state%pint(i,1)) * rga
623 : end do
624 :
625 1495368 : return
626 370944 : end subroutine check_energy_fix
627 :
628 :
629 : !===============================================================================
630 2984544 : subroutine check_tracers_init(state, tracerint)
631 :
632 : !-----------------------------------------------------------------------
633 : ! Compute initial values of tracers integrals,
634 : ! zero cumulative tendencies
635 : !-----------------------------------------------------------------------
636 :
637 : !------------------------------Arguments--------------------------------
638 :
639 : type(physics_state), intent(in) :: state
640 : type(check_tracers_data), intent(out) :: tracerint
641 :
642 : !---------------------------Local storage-------------------------------
643 :
644 : real(r8) :: tr(pcols) ! vertical integral of tracer
645 : real(r8) :: trpdel(pcols, pver) ! pdel for tracer
646 :
647 : integer ncol ! number of atmospheric columns
648 : integer i,k,m ! column, level,constituent indices
649 : integer :: ixcldice, ixcldliq ! CLDICE and CLDLIQ indices
650 : integer :: ixrain, ixsnow ! RAINQM and SNOWQM indices
651 : integer :: ixgrau ! GRAUQM index
652 : !-----------------------------------------------------------------------
653 :
654 2984544 : ncol = state%ncol
655 2984544 : call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
656 2984544 : call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
657 2984544 : call cnst_get_ind('RAINQM', ixrain, abort=.false.)
658 2984544 : call cnst_get_ind('SNOWQM', ixsnow, abort=.false.)
659 2984544 : call cnst_get_ind('GRAUQM', ixgrau, abort=.false.)
660 :
661 :
662 2984544 : do m = 1,pcnst
663 :
664 2984544 : if ( any(m == (/ 1, ixcldliq, ixcldice, &
665 : ixrain, ixsnow, ixgrau /)) ) exit ! dont process water substances
666 : ! they are checked in check_energy
667 :
668 0 : if (cnst_get_type_byind(m).eq.'dry') then
669 0 : trpdel(:ncol,:) = state%pdeldry(:ncol,:)
670 : else
671 0 : trpdel(:ncol,:) = state%pdel(:ncol,:)
672 : endif
673 :
674 : ! Compute vertical integrals of tracer
675 0 : tr = 0._r8
676 0 : do k = 1, pver
677 0 : do i = 1, ncol
678 0 : tr(i) = tr(i) + state%q(i,k,m)*trpdel(i,k)*rga
679 : end do
680 : end do
681 :
682 : ! Compute vertical integrals of frozen static tracers and total water.
683 0 : do i = 1, ncol
684 0 : tracerint%tracer(i,m) = tr(i)
685 : end do
686 :
687 : ! zero cummulative boundary fluxes
688 0 : tracerint%tracer_tnd(:ncol,m) = 0._r8
689 :
690 0 : tracerint%count(m) = 0
691 :
692 : end do
693 :
694 2984544 : return
695 : end subroutine check_tracers_init
696 :
697 : !===============================================================================
698 5969088 : subroutine check_tracers_chng(state, tracerint, name, nstep, ztodt, cflx)
699 :
700 : !-----------------------------------------------------------------------
701 : ! Check that the tracers and water change matches the boundary fluxes
702 : ! these checks are not save when there are tracers transformations, as
703 : ! they only check to see whether a mass change in the column is
704 : ! associated with a flux
705 : !-----------------------------------------------------------------------
706 :
707 : use cam_abortutils, only: endrun
708 :
709 :
710 : implicit none
711 :
712 : !------------------------------Arguments--------------------------------
713 :
714 : type(physics_state) , intent(in ) :: state
715 : type(check_tracers_data), intent(inout) :: tracerint! tracers integrals and boundary fluxes
716 : character*(*),intent(in) :: name ! parameterization name for fluxes
717 : integer , intent(in ) :: nstep ! current timestep number
718 : real(r8), intent(in ) :: ztodt ! 2 delta t (model time increment)
719 : real(r8), intent(in ) :: cflx(pcols,pcnst) ! boundary flux of tracers (kg/m2/s)
720 :
721 : !---------------------------Local storage-------------------------------
722 :
723 : real(r8) :: tracer_inp(pcols,pcnst) ! total tracer of new (input) state
724 : real(r8) :: tracer_xpd(pcols,pcnst) ! expected value (w0 + dt*boundary_flux)
725 : real(r8) :: tracer_dif(pcols,pcnst) ! tracer_inp - original tracer
726 : real(r8) :: tracer_tnd(pcols,pcnst) ! tendency from last process
727 : real(r8) :: tracer_rer(pcols,pcnst) ! relative error in tracer column
728 :
729 : real(r8) :: tr(pcols) ! vertical integral of tracer
730 : real(r8) :: trpdel(pcols, pver) ! pdel for tracer
731 :
732 : integer lchnk ! chunk identifier
733 : integer ncol ! number of atmospheric columns
734 : integer i,k ! column, level indices
735 : integer :: ixcldice, ixcldliq ! CLDICE and CLDLIQ indices
736 : integer :: ixrain, ixsnow ! RAINQM and SNOWQM indices
737 : integer :: ixgrau ! GRAUQM index
738 : integer :: m ! tracer index
739 : character(len=8) :: tracname ! tracername
740 : !-----------------------------------------------------------------------
741 :
742 5969088 : lchnk = state%lchnk
743 5969088 : ncol = state%ncol
744 5969088 : call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
745 5969088 : call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
746 5969088 : call cnst_get_ind('RAINQM', ixrain, abort=.false.)
747 5969088 : call cnst_get_ind('SNOWQM', ixsnow, abort=.false.)
748 5969088 : call cnst_get_ind('GRAUQM', ixgrau, abort=.false.)
749 :
750 5969088 : do m = 1,pcnst
751 :
752 5969088 : if ( any(m == (/ 1, ixcldliq, ixcldice, &
753 : ixrain, ixsnow, ixgrau /)) ) exit ! dont process water substances
754 : ! they are checked in check_energy
755 0 : tracname = cnst_name(m)
756 0 : if (cnst_get_type_byind(m).eq.'dry') then
757 0 : trpdel(:ncol,:) = state%pdeldry(:ncol,:)
758 : else
759 0 : trpdel(:ncol,:) = state%pdel(:ncol,:)
760 : endif
761 :
762 : ! Compute vertical integrals tracers
763 0 : tr = 0._r8
764 0 : do k = 1, pver
765 0 : do i = 1, ncol
766 0 : tr(i) = tr(i) + state%q(i,k,m)*trpdel(i,k)*rga
767 : end do
768 : end do
769 :
770 : ! Compute vertical integrals of tracer
771 0 : do i = 1, ncol
772 0 : tracer_inp(i,m) = tr(i)
773 : end do
774 :
775 : ! compute expected values and tendencies
776 0 : do i = 1, ncol
777 : ! change in tracers
778 0 : tracer_dif(i,m) = tracer_inp(i,m) - tracerint%tracer(i,m)
779 :
780 : ! expected tendencies from boundary fluxes for last process
781 0 : tracer_tnd(i,m) = cflx(i,m)
782 :
783 : ! cummulative tendencies from boundary fluxes
784 0 : tracerint%tracer_tnd(i,m) = tracerint%tracer_tnd(i,m) + tracer_tnd(i,m)
785 :
786 : ! expected new values from original values plus boundary fluxes
787 0 : tracer_xpd(i,m) = tracerint%tracer(i,m) + tracerint%tracer_tnd(i,m)*ztodt
788 :
789 : ! relative error, expected value - input value / original
790 0 : tracer_rer(i,m) = (tracer_xpd(i,m) - tracer_inp(i,m)) / tracerint%tracer(i,m)
791 : end do
792 :
793 : !! final loop for error checking
794 : ! do i = 1, ncol
795 :
796 : !! error messages
797 : ! if (abs(enrgy_rer(i)) > 1.E-14 .or. abs(water_rer(i)) > 1.E-14) then
798 : ! tracerint%count = tracerint%count + 1
799 : ! write(iulog,*) "significant conservations error after ", name, &
800 : ! " count", tracerint%count, " nstep", nstep, "chunk", lchnk, "col", i
801 : ! write(iulog,*) enrgy_inp(i),enrgy_xpd(i),enrgy_dif(i),tracerint%enrgy_tnd(i)*ztodt, &
802 : ! enrgy_tnd(i)*ztodt,enrgy_rer(i)
803 : ! write(iulog,*) water_inp(i),water_xpd(i),water_dif(i),tracerint%water_tnd(i)*ztodt, &
804 : ! water_tnd(i)*ztodt,water_rer(i)
805 : ! end if
806 : ! end do
807 :
808 :
809 : ! final loop for error checking
810 0 : if ( maxval(tracer_rer) > 1.E-14_r8 ) then
811 0 : write(iulog,*) "CHECK_TRACERS TRACER large rel error"
812 0 : write(iulog,*) tracer_rer
813 : endif
814 :
815 0 : do i = 1, ncol
816 : ! error messages
817 0 : if (abs(tracer_rer(i,m)) > 1.E-14_r8 ) then
818 0 : tracerint%count = tracerint%count + 1
819 0 : write(iulog,*) "CHECK_TRACERS TRACER significant conservation error after ", name, &
820 0 : " count", tracerint%count, " nstep", nstep, "chunk", lchnk, "col",i
821 0 : write(iulog,*)' process name, tracname, index ', name, tracname, m
822 0 : write(iulog,*)" input integral ",tracer_inp(i,m)
823 0 : write(iulog,*)" expected integral ", tracer_xpd(i,m)
824 0 : write(iulog,*)" input - inital integral ",tracer_dif(i,m)
825 0 : write(iulog,*)" cumulative tend ",tracerint%tracer_tnd(i,m)*ztodt
826 0 : write(iulog,*)" process tend ",tracer_tnd(i,m)*ztodt
827 0 : write(iulog,*)" relative error ",tracer_rer(i,m)
828 0 : call endrun()
829 : end if
830 : end do
831 : end do
832 :
833 5969088 : return
834 : end subroutine check_tracers_chng
835 :
836 : !#######################################################################
837 :
838 8959824 : subroutine tot_energy_phys(state, outfld_name_suffix,vc)
839 : use physconst, only: rga,rearth,omega
840 : use cam_thermo, only: get_hydrostatic_energy,thermo_budget_num_vars,thermo_budget_vars, &
841 : wvidx,wlidx,wiidx,seidx,poidx,keidx,moidx,mridx,ttidx,teidx
842 : use cam_history, only: outfld
843 : use dyn_tests_utils, only: vc_physics, vc_height, vc_dry_pressure
844 :
845 : use cam_abortutils, only: endrun
846 : use cam_history_support, only: max_fieldname_len
847 : use cam_budget, only: thermo_budget_history
848 : !------------------------------Arguments--------------------------------
849 :
850 : type(physics_state), intent(inout) :: state
851 : character(len=*), intent(in) :: outfld_name_suffix ! suffix for "outfld"
852 : integer, optional, intent(in) :: vc ! vertical coordinate
853 :
854 : !---------------------------Local storage-------------------------------
855 : real(r8) :: se(pcols) ! Dry Static energy (J/m2)
856 : real(r8) :: po(pcols) ! surface potential or potential energy (J/m2)
857 : real(r8) :: ke(pcols) ! kinetic energy (J/m2)
858 : real(r8) :: wv(pcols) ! column integrated vapor (kg/m2)
859 : real(r8) :: liq(pcols) ! column integrated liquid (kg/m2)
860 : real(r8) :: ice(pcols) ! column integrated ice (kg/m2)
861 : real(r8) :: tt(pcols) ! column integrated test tracer (kg/m2)
862 : real(r8) :: mr(pcols) ! column integrated wind axial angular momentum (kg*m2/s)
863 : real(r8) :: mo(pcols) ! column integrated mass axial angular momentum (kg*m2/s)
864 : real(r8) :: tt_tmp,mr_tmp,mo_tmp,cos_lat
865 : real(r8) :: mr_cnst, mo_cnst
866 : real(r8) :: cp_or_cv(pcols,pver) ! cp for pressure-based vcoord and cv for height vcoord
867 : real(r8) :: temp(pcols,pver) ! temperature
868 : real(r8) :: scaling(pcols,pver) ! scaling for conversion of temperature increment
869 :
870 : integer :: lchnk ! chunk identifier
871 : integer :: ncol ! number of atmospheric columns
872 : integer :: i,k ! column, level indices
873 : integer :: vc_loc ! local vertical coordinate variable
874 : integer :: ixtt ! test tracer index
875 : character(len=max_fieldname_len) :: name_out(thermo_budget_num_vars)
876 :
877 : !-----------------------------------------------------------------------
878 :
879 8959824 : if (.not.thermo_budget_history) return
880 :
881 0 : do i=1,thermo_budget_num_vars
882 0 : name_out(i)=trim(thermo_budget_vars(i))//'_'//trim(outfld_name_suffix)
883 : end do
884 :
885 0 : lchnk = state%lchnk
886 0 : ncol = state%ncol
887 :
888 0 : if (present(vc)) then
889 0 : vc_loc = vc
890 : else
891 0 : vc_loc = vc_physics
892 : end if
893 :
894 0 : if (state%psetcols == pcols) then
895 0 : if (vc_loc == vc_height .or. vc_loc == vc_dry_pressure) then
896 0 : cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk)
897 : else
898 0 : cp_or_cv(:ncol,:) = cpairv(:ncol,:,lchnk)
899 : end if
900 : else
901 0 : call endrun('tot_energy_phys: energy diagnostics not implemented/tested for subcolumns')
902 : end if
903 :
904 0 : if (vc_loc == vc_height .or. vc_loc == vc_dry_pressure) then
905 0 : scaling(:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv(:ncol,:)!scaling for energy consistency
906 : else
907 0 : scaling(:ncol,:) = 1.0_r8 !internal energy / enthalpy same as CAM physics
908 : end if
909 : ! scale accumulated temperature increment for internal energy / enthalpy consistency
910 0 : temp(1:ncol,:) = state%temp_ini(1:ncol,:)+scaling(1:ncol,:)*(state%T(1:ncol,:)- state%temp_ini(1:ncol,:))
911 0 : call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., &
912 0 : state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), &
913 0 : state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), temp(1:ncol,1:pver), &
914 0 : vc_loc, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), &
915 0 : z_mid = state%z_ini(1:ncol,:), se = se(1:ncol), &
916 : po = po(1:ncol), ke = ke(1:ncol), wv = wv(1:ncol), liq = liq(1:ncol), &
917 0 : ice = ice(1:ncol))
918 :
919 0 : call cnst_get_ind('TT_LW' , ixtt , abort=.false.)
920 0 : tt = 0._r8
921 0 : if (ixtt > 1) then
922 0 : if (name_out(ttidx) == 'TT_pAM'.or.name_out(ttidx) == 'TT_zAM') then
923 : !
924 : ! after dme_adjust mixing ratios are all wet
925 : !
926 0 : do k = 1, pver
927 0 : do i = 1, ncol
928 0 : tt_tmp = state%q(i,k,ixtt)*state%pdel(i,k)*rga
929 0 : tt (i) = tt(i) + tt_tmp
930 : end do
931 : end do
932 : else
933 0 : do k = 1, pver
934 0 : do i = 1, ncol
935 0 : tt_tmp = state%q(i,k,ixtt)*state%pdeldry(i,k)*rga
936 0 : tt (i) = tt(i) + tt_tmp
937 : end do
938 : end do
939 : end if
940 : end if
941 :
942 0 : call outfld(name_out(seidx) ,se , pcols ,lchnk )
943 0 : call outfld(name_out(poidx) ,po , pcols ,lchnk )
944 0 : call outfld(name_out(keidx) ,ke , pcols ,lchnk )
945 0 : call outfld(name_out(wvidx) ,wv , pcols ,lchnk )
946 0 : call outfld(name_out(wlidx) ,liq , pcols ,lchnk )
947 0 : call outfld(name_out(wiidx) ,ice , pcols ,lchnk )
948 0 : call outfld(name_out(ttidx) ,tt , pcols ,lchnk )
949 0 : call outfld(name_out(teidx) ,se+ke+po, pcols ,lchnk )
950 : !
951 : ! Axial angular momentum diagnostics
952 : !
953 : ! Code follows
954 : !
955 : ! Lauritzen et al., (2014): Held-Suarez simulations with the Community Atmosphere Model
956 : ! Spectral Element (CAM-SE) dynamical core: A global axial angularmomentum analysis using Eulerian
957 : ! and floating Lagrangian vertical coordinates. J. Adv. Model. Earth Syst. 6,129-140,
958 : ! doi:10.1002/2013MS000268
959 : !
960 : ! MR is equation (6) without \Delta A and sum over areas (areas are in units of radians**2)
961 : ! MO is equation (7) without \Delta A and sum over areas (areas are in units of radians**2)
962 : !
963 :
964 0 : mr_cnst = rga*rearth**3
965 0 : mo_cnst = rga*omega*rearth**4
966 :
967 0 : mr = 0.0_r8
968 0 : mo = 0.0_r8
969 0 : do k = 1, pver
970 0 : do i = 1, ncol
971 0 : cos_lat = cos(state%lat(i))
972 0 : mr_tmp = mr_cnst*state%u(i,k)*state%pdel(i,k)*cos_lat
973 0 : mo_tmp = mo_cnst*state%pdel(i,k)*cos_lat**2
974 :
975 0 : mr(i) = mr(i) + mr_tmp
976 0 : mo(i) = mo(i) + mo_tmp
977 : end do
978 : end do
979 :
980 0 : call outfld(name_out(mridx) ,mr, pcols,lchnk )
981 0 : call outfld(name_out(moidx) ,mo, pcols,lchnk )
982 :
983 8959824 : end subroutine tot_energy_phys
984 :
985 :
986 8959824 : end module check_energy
|
