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