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