diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml
index a2946cd243..8e951e8b5b 100644
--- a/bld/namelist_files/namelist_definition.xml
+++ b/bld/namelist_files/namelist_definition.xml
@@ -538,6 +538,13 @@ Default: .false., unless it is overridden (WACCM with interactive chemistry and
configurations do this)
+
+Number of latitude grid points for zonal average Transformed Eulerian Mean (TEM)
+diagnostics history fields. Values greater than zero turn on the TEM diagostics.
+Default: 0
+
+
Number of zonal mean basis functions (number of m=0 spherical harmonics) used in
@@ -3226,7 +3233,7 @@ Default: .false.
- The fraction of the boundary layer (PBL) depth, over which to mix the initial ZM convective parcel properties (fraction).
+ The fraction of the boundary layer (PBL) depth, over which to mix the initial ZM convective parcel properties (fraction).
Default: 0.5
diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml
index e51f8e05ff..24d6bc3ec8 100644
--- a/cime_config/testdefs/testlist_cam.xml
+++ b/cime_config/testdefs/testlist_cam.xml
@@ -1138,7 +1138,7 @@
-
+
diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/shell_commands
new file mode 100644
index 0000000000..eb40ad83e0
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/shell_commands
@@ -0,0 +1,2 @@
+./xmlchange ROF_NCPL=\$ATM_NCPL
+./xmlchange GLC_NCPL=\$ATM_NCPL
diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_cam
new file mode 100644
index 0000000000..ea040c7c83
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_cam
@@ -0,0 +1,14 @@
+mfilt=1,1,1,1,1,1,1,1,1,1
+ndens=1,1,1,1,1,1,1,1,1,1
+nhtfrq=9,9,9,9,9,9,9,9,9,9
+write_nstep0=.true.
+inithist='ENDOFRUN'
+
+cam_physics_mesh = '$ATM_DOMAIN_MESH'
+ctem_diags_numlats = 90
+
+fincl2 = 'THtem','Utem','Vtem','Wtem','VTHtem','WTHtem','UVtem', 'UWtem', 'MSKtem', 'PStem'
+fincl3 = 'Uzm', 'Vzm','Wzm', 'THzm','VTHzm','WTHzm','UVzm', 'UWzm', 'PSzm'
+
+do_circulation_diags = .false.
+fincl7 = ' '
diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_clm
new file mode 100644
index 0000000000..0d83b5367b
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_clm
@@ -0,0 +1,27 @@
+!----------------------------------------------------------------------------------
+! Users should add all user specific namelist changes below in the form of
+! namelist_var = new_namelist_value
+!
+! Include namelist variables for drv_flds_in ONLY if -megan and/or -drydep options
+! are set in the CLM_NAMELIST_OPTS env variable.
+!
+! EXCEPTIONS:
+! Set use_cndv by the compset you use and the CLM_BLDNML_OPTS -dynamic_vegetation setting
+! Set use_vichydro by the compset you use and the CLM_BLDNML_OPTS -vichydro setting
+! Set use_cn by the compset you use and CLM_BLDNML_OPTS -bgc setting
+! Set use_crop by the compset you use and CLM_BLDNML_OPTS -crop setting
+! Set spinup_state by the CLM_BLDNML_OPTS -bgc_spinup setting
+! Set irrigate by the CLM_BLDNML_OPTS -irrig setting
+! Set dtime with L_NCPL option
+! Set fatmlndfrc with LND_DOMAIN_PATH/LND_DOMAIN_FILE options
+! Set finidat with RUN_REFCASE/RUN_REFDATE/RUN_REFTOD options for hybrid or branch cases
+! (includes $inst_string for multi-ensemble cases)
+! Set glc_grid with CISM_GRID option
+! Set glc_smb with GLC_SMB option
+! Set maxpatch_glcmec with GLC_NEC option
+! Set glc_do_dynglacier with GLC_TWO_WAY_COUPLING env variable
+!----------------------------------------------------------------------------------
+hist_nhtfrq = 9
+hist_mfilt = 1
+hist_ndens = 1
+
diff --git a/doc/ChangeLog b/doc/ChangeLog
index 0adcdda10d..f5e1809e66 100644
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,4 +1,89 @@
+===============================================================
+
+Tag name: cam6_4_097
+Originator(s): fvitt
+Date: 13 Jun 2025
+One-line Summary: An alternative method for zonal-mean TEM diagnostics
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1319
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+Provides an method alternative to of computing zonal-mean terms of the
+Transformed Eulerian Mean (TEM) diagnostics independent of dynamical core.
+
+The method here uses regridding capabilities provided in the ESMF library to map the dynamical
+fields needed in TEM diagnostics from the physics grid to a distributed regular longitude -
+latitude grid. The TEM diagnostics are computed on this regular longitude - latitude grid and
+computing zonal means is straightforward.
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist:
+
+New namelist option ctem_diags_numlats :
+ Number of latitude grid points for zonal average Transformed Eulerian Mean (TEM)
+ diagnostics history fields. Values greater than zero turn on the TEM diagostics.
+
+List any changes to the defaults for the boundary datasets: N/A
+Describe any substantial timing or memory changes:
+ Performance impact is minimal
+
+Code reviewed by: nusbaume
+
+List all files eliminated: N/A
+
+List all files added and what they do:
+
+A src/utils/ctem_diags_modc.F90
+A src/utils/esmf_check_error_mod.F90
+A src/utils/esmf_lonlat_grid_mod.F90
+A src/utils/esmf_phys2lonlat_mod.F90
+A src/utils/esmf_phys_mesh_mod.F90
+A src/utils/esmf_zonal_mean_mod.F90
+ - module files for computing TEM diagnostics and zonal means
+
+A cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/shell_commands
+A cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_cam
+A cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_clm
+ - for testing the new TEM diagnostics
+
+List all existing files that have been modified, and describe the changes:
+
+M bld/namelist_files/namelist_definition.xml
+ - new namelist option described above
+
+M src/control/runtime_opts.F90
+ - read namelist options for ctem_diags_mod module
+
+M src/physics/cam/physpkg.F90
+M src/physics/cam7/physpkg.F90
+ - interface with the ctem_diags_mod module
+
+M cime_config/testdefs/testlist_cam.xml
+ - add test for TEM diags to aux_cam test of FHISTC_WXma
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below. All failed tests must be justified.
+
+derecho/intel/aux_cam:
+ DIFF ERS_Ln9.ne30pg3_ne30pg3_mg17.FHISTC_WXma.derecho_intel.cam-outfrq9s_ctem
+ - new TEM test -- no baseline to compare against
+
+ FAIL SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s
+ - pre-existing failure due to build-namelist error requiring CLM/CTSM external update
+
+derecho/nvhpc/aux_cam: All PASS
+
+izumi/nag/aux_cam: All PASS
+
+izumi/gnu/aux_cam: All PASS
+
+Summarize any changes to answers: bit-for-bit unchanged
+
+===============================================================
===============================================================
Tag name: cam6_4_096
diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90
index cb1a95c8e7..cfafafd303 100644
--- a/src/control/runtime_opts.F90
+++ b/src/control/runtime_opts.F90
@@ -105,6 +105,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon)
use surface_emissions_mod, only: surface_emissions_readnl
use elevated_emissions_mod, only: elevated_emissions_readnl
use atm_stream_ndep, only: stream_ndep_readnl
+ use ctem_diags_mod, only: ctem_diags_readnl
!---------------------------Arguments-----------------------------------
@@ -211,6 +212,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon)
call surface_emissions_readnl(nlfilename)
call elevated_emissions_readnl(nlfilename)
call stream_ndep_readnl(nlfilename)
+ call ctem_diags_readnl(nlfilename)
end subroutine read_namelist
diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90
index 78fe050101..faa956cff8 100644
--- a/src/physics/cam/physpkg.F90
+++ b/src/physics/cam/physpkg.F90
@@ -164,6 +164,7 @@ subroutine phys_register
use upper_bc, only: ubc_fixed_conc
use surface_emissions_mod, only: surface_emissions_reg
use elevated_emissions_mod, only: elevated_emissions_reg
+ use ctem_diags_mod, only: ctem_diags_reg
!---------------------------Local variables-----------------------------
!
@@ -355,6 +356,9 @@ subroutine phys_register
call HCOI_Chunk_Init()
endif
+ ! TEM diagnostics
+ call ctem_diags_reg()
+
! This needs to be last as it requires all pbuf fields to be added
if (cam_snapshot_before_num > 0 .or. cam_snapshot_after_num > 0) then
call pbuf_cam_snapshot_register()
@@ -777,6 +781,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
use elevated_emissions_mod, only: elevated_emissions_init
use ccpp_constituent_prop_mod, only: ccpp_const_props_init
+ use ctem_diags_mod, only: ctem_diags_init
! Input/output arguments
type(physics_state), pointer :: phys_state(:)
@@ -1059,6 +1064,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
dtcore_idx = pbuf_get_index('DTCORE')
dqcore_idx = pbuf_get_index('DQCORE')
+ call ctem_diags_init()
+
end subroutine phys_init
!
@@ -1083,6 +1090,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out)
#if ( defined OFFLINE_DYN )
use metdata, only: get_met_srf1
#endif
+ use ctem_diags_mod, only: ctem_diags_calc
!
! Input arguments
!
@@ -1127,6 +1135,9 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out)
call pbuf_allocate(pbuf2d, 'physpkg')
call diag_allocate()
+ ! TEM diagnostics
+ call ctem_diags_calc(phys_state)
+
!-----------------------------------------------------------------------
! Advance time information
!-----------------------------------------------------------------------
@@ -1316,6 +1327,7 @@ subroutine phys_final( phys_state, phys_tend, pbuf2d )
use microp_aero, only : microp_aero_final
use phys_grid_ctem, only : phys_grid_ctem_final
use nudging, only: Nudge_Model, nudging_final
+ use ctem_diags_mod, only: ctem_diags_final
!-----------------------------------------------------------------------
!
@@ -1346,6 +1358,8 @@ subroutine phys_final( phys_state, phys_tend, pbuf2d )
call HCOI_Chunk_Final
endif
+ call ctem_diags_final()
+
end subroutine phys_final
diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90
index 4efd6c1247..b360e0ce6b 100644
--- a/src/physics/cam7/physpkg.F90
+++ b/src/physics/cam7/physpkg.F90
@@ -158,6 +158,7 @@ subroutine phys_register
use hemco_interface, only: HCOI_Chunk_Init
use surface_emissions_mod, only: surface_emissions_reg
use elevated_emissions_mod, only: elevated_emissions_reg
+ use ctem_diags_mod, only: ctem_diags_reg
!---------------------------Local variables-----------------------------
!
@@ -345,6 +346,9 @@ subroutine phys_register
call HCOI_Chunk_Init()
endif
+ ! TEM diagnostics
+ call ctem_diags_reg()
+
! This needs to be last as it requires all pbuf fields to be added
if (cam_snapshot_before_num > 0 .or. cam_snapshot_after_num > 0) then
call pbuf_cam_snapshot_register()
@@ -772,6 +776,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
use elevated_emissions_mod, only: elevated_emissions_init
use ccpp_constituent_prop_mod, only: ccpp_const_props_init
+ use ctem_diags_mod, only: ctem_diags_init
! Input/output arguments
type(physics_state), pointer :: phys_state(:)
@@ -1053,6 +1058,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
psl_idx = pbuf_get_index('PSL')
+ call ctem_diags_init()
+
end subroutine phys_init
!
@@ -1076,6 +1083,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out)
#if ( defined OFFLINE_DYN )
use metdata, only: get_met_srf1
#endif
+ use ctem_diags_mod, only: ctem_diags_calc
!
! Input arguments
!
@@ -1121,6 +1129,9 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out)
call pbuf_allocate(pbuf2d, 'physpkg')
call diag_allocate()
+ ! TEM diagnostics
+ call ctem_diags_calc(phys_state)
+
!-----------------------------------------------------------------------
! Advance time information
!-----------------------------------------------------------------------
@@ -1310,6 +1321,7 @@ subroutine phys_final( phys_state, phys_tend, pbuf2d )
use phys_grid_ctem, only: phys_grid_ctem_final
use nudging, only: Nudge_Model, nudging_final
use hemco_interface, only: HCOI_Chunk_Final
+ use ctem_diags_mod, only: ctem_diags_final
!-----------------------------------------------------------------------
!
@@ -1336,10 +1348,12 @@ subroutine phys_final( phys_state, phys_tend, pbuf2d )
if(Nudge_Model) call nudging_final()
if(use_hemco) then
- ! cleanup hemco
- call HCOI_Chunk_Final
+ ! cleanup hemco
+ call HCOI_Chunk_Final
endif
+ call ctem_diags_final()
+
end subroutine phys_final
diff --git a/src/utils/ctem_diags_mod.F90 b/src/utils/ctem_diags_mod.F90
new file mode 100644
index 0000000000..e354e7fa64
--- /dev/null
+++ b/src/utils/ctem_diags_mod.F90
@@ -0,0 +1,498 @@
+!-----------------------------------------------------------------------------
+! For physics grid (and dycore) independent circulation diagnostics
+! -- terms of the Transformed Eulerian Mean (TEM) equation
+!
+! This uses ESMF utilities to remap dynamical fields (U,V, etc) from the physics
+! grid to a regular longitude / latitude grid where it is convenient to compute
+! the zonal mean terms of the TEM equation
+!-----------------------------------------------------------------------------
+module ctem_diags_mod
+ use shr_kind_mod, only: r8 => shr_kind_r8, cx => SHR_KIND_CX
+ use ppgrid, only: begchunk, endchunk, pcols, pver, pverp
+ use physics_types, only: physics_state
+ use phys_grid, only: get_ncols_p
+ use spmd_utils, only: masterproc
+ use ref_pres, only: pref_mid
+ use esmf_lonlat_grid_mod, only: beglon=>lon_beg, endlon=>lon_end, beglat=>lat_beg, endlat=>lat_end
+ use cam_history, only: addfld, outfld, horiz_only
+ use cam_history_support, only : fillvalue
+ use perf_mod, only: t_startf, t_stopf
+ use cam_logfile, only: iulog
+ use cam_abortutils, only: endrun
+
+ implicit none
+
+ private
+
+ public :: ctem_diags_readnl
+ public :: ctem_diags_reg
+ public :: ctem_diags_init
+ public :: ctem_diags_calc
+ public :: ctem_diags_final
+
+ integer :: ctem_diags_numlats = 0
+ logical :: ctem_diags_active = .false.
+
+contains
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine ctem_diags_readnl(nlfile)
+ use namelist_utils, only : find_group_name
+ use spmd_utils, only : mpicom, masterprocid, mpi_integer, mpi_success
+ use string_utils, only : to_lower
+
+ character(len=*), intent(in) :: nlfile
+ integer :: unitn, ierr
+ character(len=cx) :: iomsg
+
+ character(len=*), parameter :: prefix = 'ctem_diags_readnl: '
+
+ namelist /ctem_diags_nl/ ctem_diags_numlats
+
+ if (masterproc) then
+ ! read namelist
+ open( newunit=unitn, file=trim(nlfile), status='old' )
+ call find_group_name(unitn, 'ctem_diags_nl', status=ierr)
+ if (ierr == 0) then
+ read(unitn, ctem_diags_nl, iostat=ierr, iomsg=iomsg)
+ if (ierr /= 0) then
+ call endrun(prefix//'ctem_diags_nl: ERROR reading namelist: '//trim(iomsg))
+ end if
+ end if
+ close(unitn)
+ end if
+
+ call mpi_bcast(ctem_diags_numlats, 1, mpi_integer, masterprocid, mpicom, ierr)
+ if (ierr /= mpi_success) call endrun(prefix//'mpi_bcast error : ctem_diags_numlats')
+
+ ctem_diags_active = ctem_diags_numlats > 0
+
+ if (masterproc) then
+ write(iulog,*) prefix//'ctem_diags_numlats: ', ctem_diags_numlats
+ write(iulog,*) prefix//'ctem_diags_active : ', ctem_diags_active
+ end if
+
+ end subroutine ctem_diags_readnl
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine ctem_diags_reg()
+ use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap, cam_grid_register
+ use esmf_lonlat_grid_mod, only: glats, nlat, glons, nlon
+ use esmf_lonlat_grid_mod, only: esmf_lonlat_grid_init
+ use esmf_phys_mesh_mod, only: esmf_phys_mesh_init
+ use esmf_phys2lonlat_mod, only: esmf_phys2lonlat_init
+
+ integer, parameter :: zm_decomp = 331 ! Must be unique within CAM
+ integer, parameter :: reg_decomp = 332
+
+ type(horiz_coord_t), pointer :: zmlon_coord
+ type(horiz_coord_t), pointer :: zmlat_coord
+
+ integer(iMap), pointer :: grid_map(:,:)
+ real(r8) :: zmlons(1)
+
+ integer(iMap), pointer :: coord_map(:) => null()
+ type(horiz_coord_t), pointer :: lon_coord
+ type(horiz_coord_t), pointer :: lat_coord
+ integer :: i, j, ind, astat
+
+ character(len=*), parameter :: subname = 'ctem_diags_reg: '
+
+ if (.not.ctem_diags_active) return
+
+ ! initialize grids and mapping
+ call esmf_lonlat_grid_init(ctem_diags_numlats)
+ call esmf_phys_mesh_init()
+ call esmf_phys2lonlat_init()
+
+ ! Zonal mean grid for history fields
+ zmlons = 0._r8
+
+ zmlat_coord => horiz_coord_create('zmlat', '', nlat, 'latitude', 'degrees_north', 1, nlat, glats)
+ zmlon_coord => horiz_coord_create('zmlon', '', 1, 'longitude', 'degrees_east', 1, 1, zmlons)
+
+ ! grid decomposition map
+ allocate(grid_map(4,endlat-beglat+1), stat=astat)
+ if (astat/=0) then
+ call endrun(subname//'not able to allocate grid_map array')
+ end if
+
+ ind = 0
+ do j = beglat, endlat
+ ind = ind + 1
+ grid_map(1,ind) = 1
+ grid_map(2,ind) = j
+ if (beglon==1) then
+ grid_map(3,ind) = 1
+ grid_map(4,ind) = j
+ else
+ grid_map(3,ind) = 0
+ grid_map(4,ind) = 0
+ end if
+ end do
+
+ ! register the zonal average grid
+ call cam_grid_register('ctem_zm', zm_decomp, zmlat_coord, zmlon_coord, grid_map, &
+ unstruct=.false., zonal_grid=.true.)
+
+ nullify(grid_map)
+
+ ! for the lon-lat grid
+ allocate(grid_map(4, ((endlon - beglon + 1) * (endlat - beglat + 1))), stat=astat)
+ if (astat/=0) then
+ call endrun(subname//'not able to allocate grid_map array')
+ end if
+
+ ind = 0
+ do i = beglat, endlat
+ do j = beglon, endlon
+ ind = ind + 1
+ grid_map(1, ind) = j
+ grid_map(2, ind) = i
+ grid_map(3, ind) = j
+ grid_map(4, ind) = i
+ end do
+ end do
+
+ allocate(coord_map(endlat - beglat + 1), stat=astat)
+ if (astat/=0) then
+ call endrun(subname//'not able to allocate coord_map array')
+ end if
+
+ if (beglon==1) then
+ coord_map = (/ (i, i = beglat, endlat) /)
+ else
+ coord_map = 0
+ end if
+ lat_coord => horiz_coord_create('reglat', '', nlat, 'latitude', 'degrees_north', beglat, endlat, &
+ glats(beglat:endlat), map=coord_map)
+
+ nullify(coord_map)
+
+ allocate(coord_map(endlon - beglon + 1), stat=astat)
+ if (astat/=0) then
+ call endrun(subname//'not able to allocate coord_map array')
+ end if
+
+ if (beglat==1) then
+ coord_map = (/ (i, i = beglon, endlon) /)
+ else
+ coord_map = 0
+ end if
+
+ lon_coord => horiz_coord_create('reglon', '', nlon, 'longitude', 'degrees_east', beglon, endlon, &
+ glons(beglon:endlon), map=coord_map)
+
+ nullify(coord_map)
+
+ call cam_grid_register('ctem_lonlat', reg_decomp, lat_coord, lon_coord, grid_map, unstruct=.false.)
+
+ nullify(grid_map)
+
+ end subroutine ctem_diags_reg
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine ctem_diags_init()
+
+ if (.not.ctem_diags_active) return
+
+ ! fields on reg lon lat grid
+ call addfld ('THtem', (/'lev'/), 'A','K', 'Potential temp', gridname='ctem_lonlat' )
+ call addfld ('Utem', (/'lev'/), 'A','m s-1', 'Zonal wind', gridname='ctem_lonlat' )
+ call addfld ('Vtem', (/'lev'/), 'A','m s-1', 'Meridional wind', gridname='ctem_lonlat' )
+ call addfld ('Wtem', (/'lev'/), 'A','m s-1', 'Vertical wind', gridname='ctem_lonlat' )
+ call addfld ('VTHtem',(/'lev'/), 'A','K m s-1','Meridional Heat Flux:', gridname='ctem_lonlat')
+ call addfld ('WTHtem',(/'lev'/), 'A','K m s-1','Vertical Heat Flux:', gridname='ctem_lonlat')
+ call addfld ('UVtem', (/'lev'/), 'A','m2 s-2', 'Meridional Flux of Zonal Momentum', gridname='ctem_lonlat')
+ call addfld ('UWtem', (/'lev'/), 'A','m2 s-2', 'Vertical Flux of Zonal Momentum', gridname='ctem_lonlat')
+
+ ! fields on zonal mean grid
+ call addfld ('Uzm', (/'lev'/), 'A','m s-1', 'Zonal-Mean zonal wind', gridname='ctem_zm' )
+ call addfld ('Vzm', (/'lev'/), 'A','m s-1', 'Zonal-Mean meridional wind', gridname='ctem_zm' )
+ call addfld ('Wzm', (/'lev'/), 'A','m s-1', 'Zonal-Mean vertical wind', gridname='ctem_zm' )
+ call addfld ('THzm', (/'lev'/), 'A','K', 'Zonal-Mean potential temp', gridname='ctem_zm' )
+ call addfld ('VTHzm',(/'lev'/), 'A','K m s-1','Zonal-Mean meridional Heat Flux:', gridname='ctem_zm')
+ call addfld ('WTHzm',(/'lev'/), 'A','K m s-1','Zonal-Mean vertical Heat Flux:', gridname='ctem_zm')
+ call addfld ('UVzm', (/'lev'/), 'A','m2 s-2', 'Zonal-Mean meridional Flux of Zonal Momentum', gridname='ctem_zm')
+ call addfld ('UWzm', (/'lev'/), 'A','m2 s-2', 'Zonal-Mean vertical Flux of Zonal Momentum', gridname='ctem_zm')
+
+ call addfld ('PSzm', horiz_only, 'A', 'Pa', 'Zonal-Mean surface pressure', gridname='ctem_zm' )
+ call addfld ('PStem', horiz_only, 'A', 'Pa', 'Surface Pressure', gridname='ctem_lonlat')
+ call addfld ('MSKtem',horiz_only, 'A', '1', 'TEM mask', gridname='ctem_lonlat' )
+
+ end subroutine ctem_diags_init
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine ctem_diags_calc(phys_state)
+ use air_composition, only: mbarv ! g/mole
+ use shr_const_mod, only: rgas => shr_const_rgas ! J/K/kmole
+ use shr_const_mod, only: grav => shr_const_g ! m/s2
+ use esmf_phys2lonlat_mod, only: esmf_phys2lonlat_regrid
+ use esmf_zonal_mean_mod, only: esmf_zonal_mean_calc, esmf_zonal_mean_wsums, esmf_zonal_mean_masked
+ use interpolate_data, only: lininterp
+ use esmf_phys2lonlat_mod, only: fields_bundle_t, nflds
+
+ type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
+
+ real(r8), target :: u_phys(pver,pcols,begchunk:endchunk)
+ real(r8), target :: v_phys(pver,pcols,begchunk:endchunk)
+ real(r8), target :: w_phys(pver,pcols,begchunk:endchunk)
+ real(r8), target :: t_phys(pver,pcols,begchunk:endchunk)
+ real(r8), target :: p_phys(pver,pcols,begchunk:endchunk)
+ real(r8) :: ps_phys(pcols,begchunk:endchunk)
+
+ real(r8), target :: u_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8), target :: v_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8), target :: w_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8), target :: t_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8), target :: p_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: ps_lonlat(beglon:endlon,beglat:endlat)
+ real(r8) :: mskind1(beglon:endlon,beglat:endlat) ! vertical index where mountain masking begins
+
+ real(r8) :: ui_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: vi_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: wi_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: ti_lonlat(beglon:endlon,beglat:endlat,pver)
+
+ real(r8) :: u_zm(beglat:endlat,pver)
+ real(r8) :: v_zm(beglat:endlat,pver)
+ real(r8) :: w_zm(beglat:endlat,pver)
+ real(r8) :: t_zm(beglat:endlat,pver)
+ real(r8) :: ps_zm(beglat:endlat)
+
+ real(r8) :: ud_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: vd_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: wd_lonlat(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: td_lonlat(beglon:endlon,beglat:endlat,pver)
+
+ real(r8) :: vtp(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: wtp(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: uwp(beglon:endlon,beglat:endlat,pver)
+ real(r8) :: uvp(beglon:endlon,beglat:endlat,pver)
+
+ real(r8) :: vtp_zm(beglat:endlat,pver)
+ real(r8) :: wtp_zm(beglat:endlat,pver)
+ real(r8) :: uwp_zm(beglat:endlat,pver)
+ real(r8) :: uvp_zm(beglat:endlat,pver)
+ real(r8) :: wsums(beglat:endlat,pver)
+
+ integer :: lchnk, ncol, i, j, k
+ real(r8) :: sheight(pver) ! pressure scale height (m)
+
+ real(r8) :: outtmp(beglon:endlon,pver)
+ integer :: outcnt
+
+ real(r8) :: wght(beglon:endlon,beglat:endlat,pver)
+
+ type(fields_bundle_t) :: physflds(nflds)
+ type(fields_bundle_t) :: lonlatflds(nflds)
+
+ if (.not.ctem_diags_active) return
+
+ call t_startf('ctem_diags_calc')
+
+ call t_startf('ctem_diags_calc-setarrs')
+
+ do lchnk = begchunk,endchunk
+ ncol = phys_state(lchnk)%ncol
+ do i = 1,ncol
+ ! wind components
+ u_phys(:,i,lchnk) = phys_state(lchnk)%u(i,:)
+ v_phys(:,i,lchnk) = phys_state(lchnk)%v(i,:)
+
+ ! scale height
+ sheight(:) = phys_state(lchnk)%t(i,:) * rgas / ( mbarv(i,:,lchnk) * grav ) ! meters
+
+ ! vertical velocity
+ w_phys(:,i,lchnk) = -sheight(:) * phys_state(lchnk)%omega(i,:) / phys_state(lchnk)%pmid(i,:)
+
+ ! potential temperature
+ t_phys(:,i,lchnk) = phys_state(lchnk)%t(i,:) * phys_state(lchnk)%exner(i,:)
+
+ ! mid point press
+ p_phys(:,i,lchnk) = phys_state(lchnk)%pmid(i,:)
+
+ ps_phys(i,lchnk) = phys_state(lchnk)%ps(i)
+
+ end do
+ end do
+
+ call t_stopf('ctem_diags_calc-setarrs')
+
+ call t_startf('ctem_diags_calc-regrid')
+
+ ! regrid to lon/lat grid
+
+ physflds(1)%fld => u_phys
+ physflds(2)%fld => v_phys
+ physflds(3)%fld => w_phys
+ physflds(4)%fld => t_phys
+ physflds(5)%fld => p_phys
+
+ lonlatflds(1)%fld => u_lonlat
+ lonlatflds(2)%fld => v_lonlat
+ lonlatflds(3)%fld => w_lonlat
+ lonlatflds(4)%fld => t_lonlat
+ lonlatflds(5)%fld => p_lonlat
+
+ call esmf_phys2lonlat_regrid(physflds, lonlatflds)
+
+ call esmf_phys2lonlat_regrid(ps_phys, ps_lonlat)
+
+ call t_stopf('ctem_diags_calc-regrid')
+
+ call t_startf('ctem_diags_calc-zonal_mean-ps')
+ call esmf_zonal_mean_calc(ps_lonlat, ps_zm)
+ call t_stopf('ctem_diags_calc-zonal_mean-ps')
+
+ call t_startf('ctem_diags_calc-interp')
+
+ ! vertically intepolate to ref press
+ do i = beglon,endlon
+ do j = beglat,endlat
+
+ mskind1(i,j) = pver
+
+ do k = 1,pver
+ if (ps_lonlat(i,j)>=pref_mid(k)) then
+ wght(i,j,k) = 1._r8
+ mskind1(i,j) = k ! vertical index where masking begins
+ else
+ wght(i,j,k) = 0._r8
+ end if
+ end do
+
+ call lininterp( u_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
+ ui_lonlat(i,j,:), pref_mid(:), pver )
+
+ call lininterp( v_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
+ vi_lonlat(i,j,:), pref_mid(:), pver )
+
+ call lininterp( w_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
+ wi_lonlat(i,j,:), pref_mid(:), pver )
+
+ call lininterp( t_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
+ ti_lonlat(i,j,:), pref_mid(:), pver )
+
+ end do
+ end do
+
+ call t_stopf('ctem_diags_calc-interp')
+
+ call t_startf('ctem_diags_calc-zonal_mean-uvwt')
+
+ ! calculate zonal means from interpolated fields
+ ! mask out mountains from the zonal mean calculations
+ wsums = esmf_zonal_mean_wsums(wght)
+
+ call esmf_zonal_mean_masked(ui_lonlat, wght, wsums, u_zm)
+ call esmf_zonal_mean_masked(vi_lonlat, wght, wsums, v_zm)
+ call esmf_zonal_mean_masked(wi_lonlat, wght, wsums, w_zm)
+ call esmf_zonal_mean_masked(ti_lonlat, wght, wsums, t_zm)
+
+ call t_stopf('ctem_diags_calc-zonal_mean-uvwt')
+
+ call t_startf('ctem_diags_calc-calc_dev_flx')
+
+ ! Calculate zonal deviations and fluxes
+ do k = 1,pver
+ do j = beglat,endlat
+ do i = beglon,endlon
+ if (wght(i,j,k)>0._r8) then
+ ud_lonlat(i,j,k) = ui_lonlat(i,j,k) - u_zm(j,k)
+ vd_lonlat(i,j,k) = vi_lonlat(i,j,k) - v_zm(j,k)
+ wd_lonlat(i,j,k) = wi_lonlat(i,j,k) - w_zm(j,k)
+ td_lonlat(i,j,k) = ti_lonlat(i,j,k) - t_zm(j,k)
+ vtp(i,j,k) = vd_lonlat(i,j,k) * td_lonlat(i,j,k)
+ wtp(i,j,k) = wd_lonlat(i,j,k) * td_lonlat(i,j,k)
+ uwp(i,j,k) = ud_lonlat(i,j,k) * wd_lonlat(i,j,k)
+ uvp(i,j,k) = ud_lonlat(i,j,k) * vd_lonlat(i,j,k)
+ else
+ ud_lonlat(i,j,k) = fillvalue
+ vd_lonlat(i,j,k) = fillvalue
+ wd_lonlat(i,j,k) = fillvalue
+ td_lonlat(i,j,k) = fillvalue
+ vtp(i,j,k) = fillvalue
+ wtp(i,j,k) = fillvalue
+ uwp(i,j,k) = fillvalue
+ uvp(i,j,k) = fillvalue
+ end if
+ end do
+ end do
+ end do
+
+ call t_stopf('ctem_diags_calc-calc_dev_flx')
+
+ call t_startf('ctem_diags_calc-zonal_mean-p')
+
+ call esmf_zonal_mean_masked(vtp, wght, wsums, vtp_zm)
+ call esmf_zonal_mean_masked(wtp, wght, wsums, wtp_zm)
+ call esmf_zonal_mean_masked(uwp, wght, wsums, uwp_zm)
+ call esmf_zonal_mean_masked(uvp, wght, wsums, uvp_zm)
+
+ call t_stopf('ctem_diags_calc-zonal_mean-p')
+
+ call t_startf('ctem_diags_calc-output')
+
+ outcnt = endlon-beglon+1
+
+ ! output diagnostics
+ do j = beglat,endlat
+ outtmp(beglon:endlon,1:pver) = ti_lonlat(beglon:endlon,j,1:pver)
+ call outfld('THtem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = ui_lonlat(beglon:endlon,j,1:pver)
+ call outfld('Utem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = vi_lonlat(beglon:endlon,j,1:pver)
+ call outfld('Vtem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = wi_lonlat(beglon:endlon,j,1:pver)
+ call outfld('Wtem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = vtp(beglon:endlon,j,1:pver)
+ call outfld('VTHtem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = wtp(beglon:endlon,j,1:pver)
+ call outfld('WTHtem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = uvp(beglon:endlon,j,1:pver)
+ call outfld('UVtem',outtmp, outcnt, j)
+ outtmp(beglon:endlon,1:pver) = uwp(beglon:endlon,j,1:pver)
+ call outfld('UWtem',outtmp, outcnt, j)
+
+ call outfld('PStem', ps_lonlat(:,j), outcnt, j)
+ call outfld('MSKtem',mskind1(:,j), outcnt, j)
+
+ call outfld('Uzm', u_zm(j,:), 1,j)
+ call outfld('Vzm', v_zm(j,:), 1,j)
+ call outfld('Wzm', w_zm(j,:), 1,j)
+ call outfld('THzm', t_zm(j,:), 1,j)
+ call outfld('PSzm', ps_zm(j), 1,j)
+
+ call outfld('VTHzm',vtp_zm(j,:),1,j)
+ call outfld('WTHzm',wtp_zm(j,:),1,j)
+ call outfld('UVzm', uvp_zm(j,:),1,j)
+ call outfld('UWzm', uwp_zm(j,:),1,j)
+ end do
+
+ call t_stopf('ctem_diags_calc-output')
+
+ call t_stopf('ctem_diags_calc')
+
+ end subroutine ctem_diags_calc
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine ctem_diags_final()
+ use esmf_phys2lonlat_mod, only: esmf_phys2lonlat_destroy
+ use esmf_lonlat_grid_mod, only: esmf_lonlat_grid_destroy
+ use esmf_phys_mesh_mod, only: esmf_phys_mesh_destroy
+
+ if (.not.ctem_diags_active) return
+
+ call esmf_phys2lonlat_destroy()
+ call esmf_lonlat_grid_destroy()
+ call esmf_phys_mesh_destroy()
+
+ end subroutine ctem_diags_final
+
+end module ctem_diags_mod
diff --git a/src/utils/esmf_check_error_mod.F90 b/src/utils/esmf_check_error_mod.F90
new file mode 100644
index 0000000000..eb9db18a90
--- /dev/null
+++ b/src/utils/esmf_check_error_mod.F90
@@ -0,0 +1,33 @@
+!------------------------------------------------------------------------------
+! ESMF error handler
+!------------------------------------------------------------------------------
+module esmf_check_error_mod
+ use shr_kind_mod, only: cl=>SHR_KIND_CL
+ use spmd_utils, only: masterproc
+ use cam_logfile, only: iulog
+ use cam_abortutils, only: endrun
+ use ESMF, only: ESMF_SUCCESS
+
+ implicit none
+
+ private
+ public :: check_esmf_error
+
+contains
+
+ subroutine check_esmf_error( rc, errmsg )
+
+ integer, intent(in) :: rc
+ character(len=*), intent(in) :: errmsg
+
+ character(len=cl) :: errstr
+
+ if (rc /= ESMF_SUCCESS) then
+ write(errstr,'(a,i6)') 'esmf_zonal_mod::'//trim(errmsg)//' -- ESMF ERROR code: ',rc
+ if (masterproc) write(iulog,*) trim(errstr)
+ call endrun(trim(errstr))
+ end if
+
+ end subroutine check_esmf_error
+
+end module esmf_check_error_mod
diff --git a/src/utils/esmf_lonlat_grid_mod.F90 b/src/utils/esmf_lonlat_grid_mod.F90
new file mode 100644
index 0000000000..23b575ac61
--- /dev/null
+++ b/src/utils/esmf_lonlat_grid_mod.F90
@@ -0,0 +1,338 @@
+!-------------------------------------------------------------------------------
+! Encapsulates an ESMF regular longitude / latitude grid
+!-------------------------------------------------------------------------------
+module esmf_lonlat_grid_mod
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use spmd_utils, only: masterproc, mpicom
+ use cam_logfile, only: iulog
+ use cam_abortutils, only: endrun
+
+ use ESMF, only: ESMF_Grid, ESMF_GridCreate1PeriDim, ESMF_GridAddCoord
+ use ESMF, only: ESMF_GridGetCoord, ESMF_GridDestroy
+ use ESMF, only: ESMF_KIND_R8, ESMF_INDEX_GLOBAL, ESMF_STAGGERLOC_CENTER
+ use esmf_check_error_mod, only: check_esmf_error
+
+ implicit none
+
+ public
+
+ type(ESMF_Grid), protected :: lonlat_grid
+
+ integer, protected :: nlon = 0
+ integer, protected :: nlat = 0
+
+ integer, protected :: lon_beg = -1
+ integer, protected :: lon_end = -1
+ integer, protected :: lat_beg = -1
+ integer, protected :: lat_end = -1
+
+ real(r8), allocatable, protected :: glats(:)
+ real(r8), allocatable, protected :: glons(:)
+
+ integer, protected :: zonal_comm ! zonal direction MPI communicator
+
+contains
+
+ subroutine esmf_lonlat_grid_init(nlats_in)
+ use phys_grid, only: get_grid_dims
+ use mpi, only: mpi_comm_size, mpi_comm_rank, MPI_PROC_NULL, MPI_INTEGER
+
+ integer, intent(in) :: nlats_in
+
+ real(r8) :: delx, dely
+
+ integer :: npes, ierr, mytid, irank, mytidi, mytidj
+ integer :: i,j, n
+ integer :: ntasks_lon, ntasks_lat
+ integer :: lons_per_task, lons_overflow, lats_per_task, lats_overflow
+ integer :: task_cnt
+ integer :: mynlats, mynlons
+
+ integer, allocatable :: mytidi_send(:)
+ integer, allocatable :: mytidj_send(:)
+ integer, allocatable :: mytidi_recv(:)
+ integer, allocatable :: mytidj_recv(:)
+
+ integer, allocatable :: nlons_send(:)
+ integer, allocatable :: nlats_send(:)
+ integer, allocatable :: nlons_recv(:)
+ integer, allocatable :: nlats_recv(:)
+
+ integer, allocatable :: nlons_task(:)
+ integer, allocatable :: nlats_task(:)
+
+ integer, parameter :: minlats_per_pe = 2
+ integer, parameter :: minlons_per_pe = 2
+
+ integer, allocatable :: petmap(:,:,:)
+ integer :: petcnt, astat
+
+ integer :: lbnd_lat, ubnd_lat, lbnd_lon, ubnd_lon
+ integer :: lbnd(1), ubnd(1)
+ real(ESMF_KIND_R8), pointer :: coordX(:), coordY(:)
+
+ character(len=*), parameter :: subname = 'esmf_lonlat_grid_init: '
+
+ ! create reg lon lat grid
+
+ nlat = nlats_in
+ dely = 180._r8/nlat
+
+ nlon = 2*nlat
+ delx = 360._r8/nlon
+
+ allocate(glons(nlon), stat=astat)
+ if (astat/=0) then
+ call endrun(subname//'not able to allocate glons array')
+ end if
+ allocate(glats(nlat), stat=astat)
+ if (astat/=0) then
+ call endrun(subname//'not able to allocate glats array')
+ end if
+
+ glons(1) = 0._r8
+ glats(1) = -90._r8 + 0.5_r8 * dely
+
+ do i = 2,nlon
+ glons(i) = glons(i-1) + delx
+ end do
+ do i = 2,nlat
+ glats(i) = glats(i-1) + dely
+ end do
+
+ ! decompose the grid across mpi tasks ...
+
+ call mpi_comm_size(mpicom, npes, ierr)
+ call mpi_comm_rank(mpicom, mytid, ierr)
+
+ decomp_loop: do ntasks_lon = 1,nlon
+ ntasks_lat = npes/ntasks_lon
+ if ( (minlats_per_pe*ntasks_latmytid) exit jloop
+ end do
+ enddo jloop
+ endif
+
+ mynlats = lat_end-lat_beg+1
+ mynlons = lon_end-lon_beg+1
+
+ if (mynlats shr_kind_r8
+ use cam_logfile, only: iulog
+ use cam_abortutils, only: endrun
+ use spmd_utils, only: masterproc
+ use ppgrid, only: pver
+
+ use ESMF, only: ESMF_RouteHandle, ESMF_Field, ESMF_ArraySpec, ESMF_ArraySpecSet
+ use ESMF, only: ESMF_FieldCreate, ESMF_FieldRegridStore
+ use ESMF, only: ESMF_FieldGet, ESMF_FieldRegrid
+ use ESMF, only: ESMF_KIND_I4, ESMF_KIND_R8, ESMF_TYPEKIND_R8
+ use ESMF, only: ESMF_REGRIDMETHOD_BILINEAR, ESMF_POLEMETHOD_ALLAVG, ESMF_EXTRAPMETHOD_NEAREST_IDAVG
+ use ESMF, only: ESMF_TERMORDER_SRCSEQ, ESMF_MESHLOC_ELEMENT, ESMF_STAGGERLOC_CENTER
+ use ESMF, only: ESMF_FieldDestroy, ESMF_RouteHandleDestroy
+ use esmf_check_error_mod, only: check_esmf_error
+
+ implicit none
+
+ private
+
+ public :: esmf_phys2lonlat_init
+ public :: esmf_phys2lonlat_regrid
+ public :: esmf_phys2lonlat_destroy
+ public :: fields_bundle_t
+ public :: nflds
+
+ type(ESMF_RouteHandle) :: rh_phys2lonlat_3d
+ type(ESMF_RouteHandle) :: rh_phys2lonlat_2d
+
+ type(ESMF_Field) :: physfld_3d
+ type(ESMF_Field) :: lonlatfld_3d
+
+ type(ESMF_Field) :: physfld_2d
+ type(ESMF_Field) :: lonlatfld_2d
+
+ interface esmf_phys2lonlat_regrid
+ module procedure esmf_phys2lonlat_regrid_2d
+ module procedure esmf_phys2lonlat_regrid_3d
+ end interface esmf_phys2lonlat_regrid
+
+ type :: fields_bundle_t
+ real(r8), pointer :: fld(:,:,:) => null()
+ end type fields_bundle_t
+
+ integer, parameter :: nflds = 5
+
+contains
+
+ !------------------------------------------------------------------------------
+ !------------------------------------------------------------------------------
+ subroutine esmf_phys2lonlat_init()
+ use esmf_phys_mesh_mod, only: physics_grid_mesh
+ use esmf_lonlat_grid_mod, only: lonlat_grid
+
+ type(ESMF_ArraySpec) :: arrayspec
+ integer(ESMF_KIND_I4), pointer :: factorIndexList(:,:)
+ real(ESMF_KIND_R8), pointer :: factorList(:)
+ integer :: smm_srctermproc, smm_pipelinedep, rc
+
+ character(len=*), parameter :: subname = 'esmf_phys2lonlat_init: '
+
+ smm_srctermproc = 0
+ smm_pipelinedep = 16
+
+ ! create ESMF fields
+
+ ! 3D phys fld
+ call ESMF_ArraySpecSet(arrayspec, 3, ESMF_TYPEKIND_R8, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_ArraySpecSet 3D phys fld ERROR')
+
+ physfld_3d = ESMF_FieldCreate(physics_grid_mesh, arrayspec, &
+ gridToFieldMap=(/3/), meshloc=ESMF_MESHLOC_ELEMENT, &
+ ungriddedLBound=(/1,1/), ungriddedUBound=(/pver,nflds/), rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldCreate 3D phys fld ERROR')
+
+ ! 3D lon lat grid
+ call ESMF_ArraySpecSet(arrayspec, 4, ESMF_TYPEKIND_R8, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_ArraySpecSet 3D lonlat fld ERROR')
+
+ lonlatfld_3d = ESMF_FieldCreate( lonlat_grid, arrayspec, staggerloc=ESMF_STAGGERLOC_CENTER, &
+ ungriddedLBound=(/1,1/), ungriddedUBound=(/pver,nflds/), rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldCreate 3D lonlat fld ERROR')
+
+ ! 2D phys fld
+ call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_ArraySpecSet 2D phys fld ERROR')
+
+ physfld_2d = ESMF_FieldCreate(physics_grid_mesh, arrayspec, &
+ meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldCreate 2D phys fld ERROR')
+
+ ! 2D lon/lat grid
+ call ESMF_ArraySpecSet(arrayspec, 2, ESMF_TYPEKIND_R8, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_ArraySpecSet 2D lonlat fld ERROR')
+
+ lonlatfld_2d = ESMF_FieldCreate( lonlat_grid, arrayspec, staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldCreate 2D lonlat fld ERROR')
+
+ call ESMF_FieldRegridStore(srcField=physfld_3d, dstField=lonlatfld_3d, &
+ regridMethod=ESMF_REGRIDMETHOD_BILINEAR, &
+ polemethod=ESMF_POLEMETHOD_ALLAVG, &
+ extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_IDAVG, &
+ routeHandle=rh_phys2lonlat_3d, factorIndexList=factorIndexList, &
+ factorList=factorList, srcTermProcessing=smm_srctermproc, &
+ pipelineDepth=smm_pipelinedep, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldRegridStore 3D routehandle ERROR')
+
+ call ESMF_FieldRegridStore(srcField=physfld_2d, dstField=lonlatfld_2d, &
+ regridMethod=ESMF_REGRIDMETHOD_BILINEAR, &
+ polemethod=ESMF_POLEMETHOD_ALLAVG, &
+ extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_IDAVG, &
+ routeHandle=rh_phys2lonlat_2d, factorIndexList=factorIndexList, &
+ factorList=factorList, srcTermProcessing=smm_srctermproc, &
+ pipelineDepth=smm_pipelinedep, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldRegridStore 3D routehandle ERROR')
+
+ end subroutine esmf_phys2lonlat_init
+
+ !------------------------------------------------------------------------------
+ !------------------------------------------------------------------------------
+ subroutine esmf_phys2lonlat_regrid_3d(physflds, lonlatflds)
+ use esmf_lonlat_grid_mod, only: lon_beg,lon_end,lat_beg,lat_end
+ use ppgrid, only: pcols, pver, begchunk, endchunk
+ use phys_grid, only: get_ncols_p
+
+ type(fields_bundle_t), intent(in) :: physflds(nflds)
+ type(fields_bundle_t), intent(inout) :: lonlatflds(nflds)
+
+ integer :: i, ichnk, ncol, ifld, ilev, icol, rc
+ real(ESMF_KIND_R8), pointer :: physptr(:,:,:)
+ real(ESMF_KIND_R8), pointer :: lonlatptr(:,:,:,:)
+
+ character(len=*), parameter :: subname = 'esmf_phys2lonlat_regrid_3d: '
+
+ call ESMF_FieldGet(physfld_3d, localDe=0, farrayPtr=physptr, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldGet physptr')
+
+ i = 0
+ do ichnk = begchunk, endchunk
+ ncol = get_ncols_p(ichnk)
+ do icol = 1,ncol
+ i = i+1
+ do ifld = 1,nflds
+ do ilev = 1,pver
+ physptr(ilev,ifld,i) = physflds(ifld)%fld(ilev,icol,ichnk)
+ end do
+ end do
+ end do
+ end do
+
+ call ESMF_FieldRegrid(physfld_3d, lonlatfld_3d, rh_phys2lonlat_3d, &
+ termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldRegrid physfld_3d->lonlatfld_3d')
+
+ call ESMF_FieldGet(lonlatfld_3d, localDe=0, farrayPtr=lonlatptr, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldGet lonlatptr')
+
+ do ifld = 1,nflds
+ lonlatflds(ifld)%fld(lon_beg:lon_end,lat_beg:lat_end,1:pver) = lonlatptr(lon_beg:lon_end,lat_beg:lat_end,1:pver,ifld)
+ end do
+
+ end subroutine esmf_phys2lonlat_regrid_3d
+
+ !------------------------------------------------------------------------------
+ !------------------------------------------------------------------------------
+ subroutine esmf_phys2lonlat_regrid_2d(physarr, lonlatarr)
+ use esmf_lonlat_grid_mod, only: lon_beg,lon_end,lat_beg,lat_end
+ use ppgrid, only: pcols, pver, begchunk, endchunk
+ use phys_grid, only: get_ncols_p
+
+ real(r8),intent(in) :: physarr(pcols,begchunk:endchunk)
+ real(r8),intent(out) :: lonlatarr(lon_beg:lon_end,lat_beg:lat_end)
+
+ integer :: i, ichnk, ncol, icol, rc
+ real(ESMF_KIND_R8), pointer :: physptr(:)
+ real(ESMF_KIND_R8), pointer :: lonlatptr(:,:)
+
+ character(len=*), parameter :: subname = 'esmf_phys2lonlat_regrid_2d: '
+
+ call ESMF_FieldGet(physfld_2d, localDe=0, farrayPtr=physptr, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldGet physptr')
+
+ i = 0
+ do ichnk = begchunk, endchunk
+ ncol = get_ncols_p(ichnk)
+ do icol = 1,ncol
+ i = i+1
+ physptr(i) = physarr(icol,ichnk)
+ end do
+ end do
+
+ call ESMF_FieldRegrid(physfld_2d, lonlatfld_2d, rh_phys2lonlat_2d, &
+ termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldRegrid physfld_3d->lonlatfld_3d')
+
+ call ESMF_FieldGet(lonlatfld_2d, localDe=0, farrayPtr=lonlatptr, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldGet lonlatptr')
+
+ lonlatarr(lon_beg:lon_end,lat_beg:lat_end) = lonlatptr(lon_beg:lon_end,lat_beg:lat_end)
+
+ end subroutine esmf_phys2lonlat_regrid_2d
+
+ !------------------------------------------------------------------------------
+ !------------------------------------------------------------------------------
+ subroutine esmf_phys2lonlat_destroy()
+
+ integer :: rc
+ character(len=*), parameter :: subname = 'esmf_phys2lonlat_destroy: '
+
+ call ESMF_RouteHandleDestroy(rh_phys2lonlat_3d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldDestroy rh_phys2lonlat_3d')
+
+ call ESMF_RouteHandleDestroy(rh_phys2lonlat_2d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldDestroy rh_phys2lonlat_2d')
+
+ call ESMF_FieldDestroy(lonlatfld_3d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldDestroy lonlatfld_3d')
+
+ call ESMF_FieldDestroy(lonlatfld_2d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldDestroy lonlatfld_2d')
+
+ call ESMF_FieldDestroy(physfld_3d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldDestroy physfld_3d')
+
+ call ESMF_FieldDestroy(physfld_2d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_FieldDestroy physfld_2d')
+
+ end subroutine esmf_phys2lonlat_destroy
+
+end module esmf_phys2lonlat_mod
diff --git a/src/utils/esmf_phys_mesh_mod.F90 b/src/utils/esmf_phys_mesh_mod.F90
new file mode 100644
index 0000000000..4235452dc5
--- /dev/null
+++ b/src/utils/esmf_phys_mesh_mod.F90
@@ -0,0 +1,202 @@
+!-------------------------------------------------------------------------------
+! Encapsulates the CAM physics grid mesh
+!-------------------------------------------------------------------------------
+module esmf_phys_mesh_mod
+ use shr_kind_mod, only: r8 => shr_kind_r8, cs=>shr_kind_cs, cl=>shr_kind_cl
+ use spmd_utils, only: masterproc
+ use cam_logfile, only: iulog
+ use cam_abortutils, only: endrun
+ use ESMF, only: ESMF_DistGrid, ESMF_DistGridCreate, ESMF_MeshCreate
+ use ESMF, only: ESMF_FILEFORMAT_ESMFMESH,ESMF_MeshGet,ESMF_Mesh, ESMF_SUCCESS
+ use ESMF, only: ESMF_MeshDestroy, ESMF_DistGridDestroy
+ use esmf_check_error_mod, only: check_esmf_error
+
+ implicit none
+
+ private
+
+ public :: esmf_phys_mesh_init
+ public :: esmf_phys_mesh_destroy
+ public :: physics_grid_mesh
+
+ ! phys_mesh: Local copy of physics grid
+ type(ESMF_Mesh), protected :: physics_grid_mesh
+
+ ! dist_grid_2d: DistGrid for 2D fields
+ type(ESMF_DistGrid) :: dist_grid_2d
+
+contains
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine esmf_phys_mesh_init()
+ use phys_control, only: phys_getopts
+ use phys_grid, only: get_ncols_p, get_gcol_p, get_rlon_all_p, get_rlat_all_p
+ use ppgrid, only: pcols, begchunk, endchunk
+ use shr_const_mod,only: shr_const_pi
+
+ ! Local variables
+ integer :: ncols
+ integer :: chnk, col, dindex
+ integer, allocatable :: decomp(:)
+ character(len=cl) :: grid_file
+ integer :: spatialDim
+ integer :: numOwnedElements
+ real(r8), pointer :: ownedElemCoords(:)
+ real(r8), pointer :: lat(:), latMesh(:)
+ real(r8), pointer :: lon(:), lonMesh(:)
+ real(r8) :: lats(pcols) ! array of chunk latitudes
+ real(r8) :: lons(pcols) ! array of chunk longitude
+ character(len=cs) :: tempc1,tempc2
+ character(len=300) :: errstr
+
+ integer :: i, c, n, total_cols, rc
+
+ real(r8), parameter :: abstol = 1.e-3_r8
+ real(r8), parameter :: radtodeg = 180.0_r8/shr_const_pi
+ character(len=*), parameter :: subname = 'esmf_phys_mesh_init: '
+
+ ! Find the physics grid file
+ call phys_getopts(physics_grid_out=grid_file)
+
+ ! Compute the local decomp
+ total_cols = 0
+ do chnk = begchunk, endchunk
+ total_cols = total_cols + get_ncols_p(chnk)
+ end do
+ allocate(decomp(total_cols), stat=rc)
+ if (rc/=0) then
+ call endrun(subname//'not able to allocate decomp')
+ end if
+
+ dindex = 0
+ do chnk = begchunk, endchunk
+ ncols = get_ncols_p(chnk)
+ do col = 1, ncols
+ dindex = dindex + 1
+ decomp(dindex) = get_gcol_p(chnk, col)
+ end do
+ end do
+
+ ! Create a DistGrid based on the physics decomp
+ dist_grid_2d = ESMF_DistGridCreate(arbSeqIndexList=decomp, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_DistGridCreate')
+
+ ! Create an ESMF_mesh for the physics decomposition
+ physics_grid_mesh = ESMF_MeshCreate(trim(grid_file), ESMF_FILEFORMAT_ESMFMESH, &
+ elementDistgrid=dist_grid_2d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_MeshCreate')
+
+ ! Check that the mesh coordinates are consistent with the model physics column coordinates
+
+ ! obtain mesh lats and lons
+ call ESMF_MeshGet(physics_grid_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_MeshGet')
+
+ if (numOwnedElements /= total_cols) then
+ write(tempc1,'(i10)') numOwnedElements
+ write(tempc2,'(i10)') total_cols
+ call endrun(subname//"ERROR numOwnedElements "// &
+ trim(tempc1) //" not equal to local size "// trim(tempc2))
+ end if
+
+ allocate(ownedElemCoords(spatialDim*numOwnedElements), stat=rc)
+ if (rc/=0) then
+ call endrun(subname//'not able to allocate ownedElemCoords')
+ end if
+
+ allocate(lonMesh(total_cols), stat=rc)
+ if (rc/=0) then
+ call endrun(subname//'not able to allocate lonMesh')
+ end if
+
+ allocate(latMesh(total_cols), stat=rc)
+ if (rc/=0) then
+ call endrun(subname//'not able to allocate latMesh')
+ end if
+
+ call ESMF_MeshGet(physics_grid_mesh, ownedElemCoords=ownedElemCoords)
+ call check_esmf_error(rc, subname//'ESMF_MeshGet')
+
+ do n = 1,total_cols
+ lonMesh(n) = ownedElemCoords(2*n-1)
+ latMesh(n) = ownedElemCoords(2*n)
+ end do
+
+ ! obtain internally generated cam lats and lons
+ allocate(lon(total_cols), stat=rc);
+ if (rc/=0) then
+ call endrun(subname//'not able to allocate lon')
+ end if
+
+ lon(:) = 0._r8
+
+ allocate(lat(total_cols), stat=rc);
+ if (rc/=0) then
+ call endrun(subname//'not able to allocate lat')
+ end if
+
+ lat(:) = 0._r8
+
+ n=0
+ do c = begchunk, endchunk
+ ncols = get_ncols_p(c)
+ ! latitudes and longitudes returned in radians
+ call get_rlat_all_p(c, ncols, lats)
+ call get_rlon_all_p(c, ncols, lons)
+ do i=1,ncols
+ n = n+1
+ lat(n) = lats(i)*radtodeg
+ lon(n) = lons(i)*radtodeg
+ end do
+ end do
+
+ errstr = ''
+ ! error check differences between internally generated lons and those read in
+ do n = 1,total_cols
+ if (abs(lonMesh(n) - lon(n)) > abstol) then
+ if ( (abs(lonMesh(n)-lon(n)) > 360._r8+abstol) .or. (abs(lonMesh(n)-lon(n)) < 360._r8-abstol) ) then
+ write(errstr,100) n,lon(n),lonMesh(n), abs(lonMesh(n)-lon(n))
+ write(iulog,*) trim(errstr)
+ endif
+ end if
+ if (abs(latMesh(n) - lat(n)) > abstol) then
+ ! poles in the 4x5 SCRIP file seem to be off by 1 degree
+ if (.not.( (abs(lat(n))>88.0_r8) .and. (abs(latMesh(n))>88.0_r8) )) then
+ write(errstr,101) n,lat(n),latMesh(n), abs(latMesh(n)-lat(n))
+ write(iulog,*) trim(errstr)
+ endif
+ end if
+ end do
+
+ if ( len_trim(errstr) > 0 ) then
+ call endrun(subname//'physics mesh coords do not match model coords')
+ end if
+
+ ! deallocate memory
+ deallocate(ownedElemCoords)
+ deallocate(lon, lonMesh)
+ deallocate(lat, latMesh)
+ deallocate(decomp)
+
+100 format('esmf_phys_mesh_init: coord mismatch... n, lon(n), lonmesh(n), diff_lon = ',i6,2(f21.13,3x),d21.5)
+101 format('esmf_phys_mesh_init: coord mismatch... n, lat(n), latmesh(n), diff_lat = ',i6,2(f21.13,3x),d21.5)
+
+ end subroutine esmf_phys_mesh_init
+
+ !-----------------------------------------------------------------------------
+ !-----------------------------------------------------------------------------
+ subroutine esmf_phys_mesh_destroy()
+
+ integer :: rc
+ character(len=*), parameter :: subname = 'esmf_phys_mesh_destroy: '
+
+ call ESMF_MeshDestroy(physics_grid_mesh, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_MeshDestroy phys_mesh')
+
+ call ESMF_DistGridDestroy(dist_grid_2d, rc=rc)
+ call check_esmf_error(rc, subname//'ESMF_DistGridDestroy dist_grid_2d')
+
+ end subroutine esmf_phys_mesh_destroy
+
+end module esmf_phys_mesh_mod
diff --git a/src/utils/esmf_zonal_mean_mod.F90 b/src/utils/esmf_zonal_mean_mod.F90
new file mode 100644
index 0000000000..822c7ae639
--- /dev/null
+++ b/src/utils/esmf_zonal_mean_mod.F90
@@ -0,0 +1,155 @@
+!------------------------------------------------------------------------------
+! Provides methods for calculating zonal means on the ESMF regular longitude
+! / latitude grid
+!------------------------------------------------------------------------------
+module esmf_zonal_mean_mod
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use cam_logfile, only: iulog
+ use cam_abortutils, only: endrun
+ use spmd_utils, only: masterproc
+ use cam_history_support, only : fillvalue
+
+ implicit none
+
+ private
+
+ public :: esmf_zonal_mean_calc
+ public :: esmf_zonal_mean_masked
+ public :: esmf_zonal_mean_wsums
+
+ interface esmf_zonal_mean_calc
+ module procedure esmf_zonal_mean_calc_2d
+ module procedure esmf_zonal_mean_calc_3d
+ end interface esmf_zonal_mean_calc
+
+contains
+
+ !------------------------------------------------------------------------------
+ ! Calculates zonal means of 3D fields. The wght option can be used to mask out
+ ! regions such as mountains.
+ !------------------------------------------------------------------------------
+ subroutine esmf_zonal_mean_calc_3d(lonlatarr, zmarr, wght)
+ use ppgrid, only: pver
+ use esmf_lonlat_grid_mod, only: lon_beg,lon_end,lat_beg,lat_end, nlon
+ use esmf_lonlat_grid_mod, only: zonal_comm
+ use shr_reprosum_mod,only: shr_reprosum_calc
+
+ real(r8), intent(in) :: lonlatarr(lon_beg:lon_end,lat_beg:lat_end,pver)
+ real(r8), intent(out) :: zmarr(lat_beg:lat_end,pver)
+
+ real(r8), optional, intent(in) :: wght(lon_beg:lon_end,lat_beg:lat_end,pver)
+
+ real(r8) :: tmparr(lon_beg:lon_end,pver)
+ real(r8) :: gsum(pver)
+
+ real(r8) :: wsums(lat_beg:lat_end,pver)
+
+ integer :: numlons, ilat, ilev
+
+ numlons = lon_end-lon_beg+1
+
+ ! zonal mean
+ if (present(wght)) then
+
+ wsums = esmf_zonal_mean_wsums(wght)
+ call esmf_zonal_mean_masked(lonlatarr, wght, wsums, zmarr)
+
+ else
+
+ do ilat = lat_beg, lat_end
+ call shr_reprosum_calc(lonlatarr(lon_beg:lon_end,ilat,:), gsum, numlons, numlons, pver, gbl_count=nlon, commid=zonal_comm)
+ zmarr(ilat,:) = gsum(:)/nlon
+ end do
+
+ end if
+
+ end subroutine esmf_zonal_mean_calc_3d
+
+ !------------------------------------------------------------------------------
+ ! Computes zonal mean for 2D lon / lat fields.
+ !------------------------------------------------------------------------------
+ subroutine esmf_zonal_mean_calc_2d(lonlatarr, zmarr)
+ use esmf_lonlat_grid_mod, only: lon_beg,lon_end,lat_beg,lat_end, nlon
+ use esmf_lonlat_grid_mod, only: zonal_comm
+ use shr_reprosum_mod,only: shr_reprosum_calc
+
+ real(r8), intent(in) :: lonlatarr(lon_beg:lon_end,lat_beg:lat_end)
+ real(r8), intent(out) :: zmarr(lat_beg:lat_end)
+
+ real(r8) :: gsum(lat_beg:lat_end)
+
+ integer :: numlons, numlats
+
+ numlons = lon_end-lon_beg+1
+ numlats = lat_end-lat_beg+1
+
+ ! zonal mean
+
+ call shr_reprosum_calc(lonlatarr, gsum, numlons, numlons, numlats, gbl_count=nlon, commid=zonal_comm)
+ zmarr(:) = gsum(:)/nlon
+
+ end subroutine esmf_zonal_mean_calc_2d
+
+ !------------------------------------------------------------------------------
+ ! Computes longitude sums of grid cell weights.
+ !------------------------------------------------------------------------------
+ function esmf_zonal_mean_wsums(wght) result(wsums)
+ use esmf_lonlat_grid_mod, only: lon_beg,lon_end,lat_beg,lat_end, nlon
+ use esmf_lonlat_grid_mod, only: zonal_comm
+ use shr_reprosum_mod,only: shr_reprosum_calc
+ use ppgrid, only: pver
+
+ real(r8), intent(in) :: wght(lon_beg:lon_end,lat_beg:lat_end,pver)
+
+ real(r8) :: wsums(lat_beg:lat_end,pver)
+ integer :: numlons, ilat
+
+ numlons = lon_end-lon_beg+1
+
+ do ilat = lat_beg, lat_end
+
+ call shr_reprosum_calc(wght(lon_beg:lon_end,ilat,:), wsums(ilat,1:pver), &
+ numlons, numlons, pver, gbl_count=nlon, commid=zonal_comm)
+
+ end do
+
+ end function esmf_zonal_mean_wsums
+
+ !------------------------------------------------------------------------------
+ ! Masks out regions (e.g. mountains) from zonal mean calculation.
+ !------------------------------------------------------------------------------
+ subroutine esmf_zonal_mean_masked(lonlatarr, wght, wsums, zmarr)
+ use esmf_lonlat_grid_mod, only: lon_beg,lon_end,lat_beg,lat_end, nlon
+ use esmf_lonlat_grid_mod, only: zonal_comm
+ use shr_reprosum_mod,only: shr_reprosum_calc
+ use ppgrid, only: pver
+
+ real(r8), intent(in) :: lonlatarr(lon_beg:lon_end,lat_beg:lat_end,pver)
+ real(r8), intent(in) :: wght(lon_beg:lon_end,lat_beg:lat_end,pver) ! grid cell weights
+ real(r8), intent(in) :: wsums(lat_beg:lat_end,pver) ! pre-computed sums of grid cell weights
+ real(r8), intent(out) :: zmarr(lat_beg:lat_end,pver) ! zonal means
+
+ real(r8) :: tmparr(lon_beg:lon_end,pver)
+ integer :: numlons, ilat, ilev
+ real(r8) :: gsum(pver)
+
+ numlons = lon_end-lon_beg+1
+
+ do ilat = lat_beg, lat_end
+
+ tmparr(lon_beg:lon_end,:) = wght(lon_beg:lon_end,ilat,:)*lonlatarr(lon_beg:lon_end,ilat,:)
+ call shr_reprosum_calc(tmparr, gsum, numlons, numlons, pver, gbl_count=nlon, commid=zonal_comm)
+
+ do ilev = 1,pver
+ if (wsums(ilat,ilev)>0._r8) then
+ zmarr(ilat,ilev) = gsum(ilev)/wsums(ilat,ilev)
+ else
+ zmarr(ilat,ilev) = fillvalue
+ end if
+ end do
+
+ end do
+
+ end subroutine esmf_zonal_mean_masked
+
+end module esmf_zonal_mean_mod