LCOV - code coverage report
Current view: | top level - physics/cam - check_energy.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 190 | 348 | 54.6 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 11 | 12 | 91.7 %

              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           2 :          read(unitn, check_energy_nl, iostat=ierr)
         110           2 :          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  2365672176 :        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  2323627992 :     state%temp_ini(:ncol,:) = state%T(:ncol,:)
         270  2323627992 :     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  2323627992 :         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    25328376 :   subroutine check_energy_chng(state, tend, name, nstep, ztodt,        &
         326    25328376 :        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    50656752 :     real(r8) :: te_xpd(state%ncol)                 ! expected value (f0 + dt*boundary_flux)
         355    50656752 :     real(r8) :: te_dif(state%ncol)                 ! energy of input state - original energy
         356    50656752 :     real(r8) :: te_tnd(state%ncol)                 ! tendency from last process
         357    50656752 :     real(r8) :: te_rer(state%ncol)                 ! relative error in energy column
         358             :
         359    50656752 :     real(r8) :: tw_xpd(state%ncol)                 ! expected value (w0 + dt*boundary_flux)
         360    50656752 :     real(r8) :: tw_dif(state%ncol)                 ! tw_inp - original water
         361    50656752 :     real(r8) :: tw_tnd(state%ncol)                 ! tendency from last process
         362    50656752 :     real(r8) :: tw_rer(state%ncol)                 ! relative error in water column
         363             :
         364    50656752 :     real(r8) :: te(state%ncol)                     ! vertical integral of total energy
         365    50656752 :     real(r8) :: tw(state%ncol)                     ! vertical integral of total water
         366    50656752 :     real(r8) :: cp_or_cv(state%psetcols,pver)      ! cp or cv depending on vcoord
         367    50656752 :     real(r8) :: scaling(state%psetcols,pver)       ! scaling for conversion of temperature increment
         368    50656752 :     real(r8) :: temp(state%ncol,pver)              ! temperature
         369             :
         370    50656752 :     real(r8) :: se(state%ncol)                     ! enthalpy or internal energy (J/m2)
         371    50656752 :     real(r8) :: po(state%ncol)                     ! surface potential or potential energy (J/m2)
         372    50656752 :     real(r8) :: ke(state%ncol)                     ! kinetic energy    (J/m2)
         373    50656752 :     real(r8) :: wv(state%ncol)                     ! column integrated vapor       (kg/m2)
         374    50656752 :     real(r8) :: liq(state%ncol)                    ! column integrated liquid      (kg/m2)
         375    50656752 :     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    25328376 :     lchnk = state%lchnk
         383    25328376 :     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    25328376 :     if (state%psetcols == pcols) then
         389 40069490832 :        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    75985128 :     call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true.,               &
         397    75985128 :          state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver),                         &
         398   151970256 :          state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver),     &
         399    50656752 :          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   379925640 :          ke=ke(1:ncol),wv=wv(1:ncol),liq=liq(1:ncol),ice=ice(1:ncol))
         402             :     ! compute expected values and tendencies
         403   422924976 :     do i = 1, ncol
         404             :        ! change in static energy and total water
         405   397596600 :        te_dif(i) = te(i) - state%te_cur(i,phys_te_idx)
         406   397596600 :        tw_dif(i) = tw(i) - state%tw_cur(i)
         407             :
         408             :        ! expected tendencies from boundary fluxes for last process
         409   397596600 :        te_tnd(i) = flx_vap(i)*(latvap+latice) - (flx_cnd(i) - flx_ice(i))*1000._r8*latice + flx_sen(i)
         410   397596600 :        tw_tnd(i) = flx_vap(i) - flx_cnd(i) *1000._r8
         411             :
         412             :        ! cummulative tendencies from boundary fluxes
         413   397596600 :        tend%te_tnd(i) = tend%te_tnd(i) + te_tnd(i)
         414   397596600 :        tend%tw_tnd(i) = tend%tw_tnd(i) + tw_tnd(i)
         415             :
         416             :        ! expected new values from previous state plus boundary fluxes
         417   397596600 :        te_xpd(i) = state%te_cur(i,phys_te_idx) + te_tnd(i)*ztodt
         418   397596600 :        tw_xpd(i) = state%tw_cur(i)             + tw_tnd(i)*ztodt
         419             :
         420             :        ! relative error, expected value - input state / previous state
         421   422924976 :        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   422924976 :     tw_rer = 0._r8
         426   422924976 :     where (state%tw_cur(:ncol) > 0._r8)
         427    25328376 :        tw_rer(:ncol) = (tw_xpd(:ncol) - tw(:ncol)) / state%tw_cur(:ncol)
         428             :     end where
         429             :
         430             :     ! error checking
         431    25328376 :     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   422924976 :     do i = 1, ncol
         458   397596600 :       state%te_cur(i,phys_te_idx) = te(i)
         459   422924976 :       state%tw_cur(i)             = tw(i)
         460             :     end do
         461             :
         462             :     !
         463             :     ! Dynamical core total energy
         464             :     !
         465    25328376 :     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    25328376 :     else if (vc_dycore == vc_dry_pressure) then
         484             :       !
         485             :       ! SE specific hydrostatic energy
         486             :       !
         487    25328376 :       if (state%psetcols == pcols) then
         488 39357351144 :         cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk)
         489 39357351144 :         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 39357351144 :       temp(1:ncol,:)   = state%temp_ini(1:ncol,:)+scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:))
         498    75985128 :       call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true.,               &
         499    50656752 :            state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver),                         &
         500   101313504 :            state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), temp(1:ncol,1:pver),        &
         501    50656752 :            vc_dry_pressure, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol),   &
         502   303940512 :            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    25328376 :   end subroutine check_energy_chng
         507             :
         508             :
         509      741888 :   subroutine check_energy_gmean(state, pbuf2d, dtime, nstep)
         510             :
         511    25328376 :     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    62805456 :   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  2323627992 :     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     5956704 :   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     5956704 :     lchnk = state%lchnk
         743     5956704 :     ncol  = state%ncol
         744     5956704 :     call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
         745     5956704 :     call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
         746     5956704 :     call cnst_get_ind('RAINQM', ixrain,   abort=.false.)
         747     5956704 :     call cnst_get_ind('SNOWQM', ixsnow,   abort=.false.)
         748     5956704 :     call cnst_get_ind('GRAUQM', ixgrau,   abort=.false.)
         749             :
         750     5956704 :     do m = 1,pcnst
         751             :
         752     5956704 :        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     5956704 :     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

Generated by: LCOV version 1.14