Skip to content

cam6_4_097: An alternative for TEM diagnostics #1319

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
Jun 13, 2025
Merged

cam6_4_097: An alternative for TEM diagnostics #1319

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8056408
An alternative for TEM diagnostics
fvitt Jun 2, 2025
8b39771
Changes suggested by code reviewer
fvitt Jun 10, 2025
fe4be3d
Fix typos
fvitt Jun 10, 2025
cf57d22
Exercise the TEM diags in an aux_cam test
fvitt Jun 13, 2025
92477aa
Merge remote-tracking branch 'ESCOMP/cam_development' into esmf_tem_d…
fvitt Jun 13, 2025
fd7f708
ChangeLog draft
fvitt Jun 13, 2025
6878d14
ChangeLog update
fvitt Jun 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions src/physics/cam/physpkg.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -356,13 +356,14 @@ 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()
end if

call ctem_diags_reg()

end subroutine phys_register


Expand Down
5 changes: 3 additions & 2 deletions src/physics/cam7/physpkg.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,13 +346,14 @@ 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()
end if

call ctem_diags_reg()

end subroutine phys_register


Expand Down
47 changes: 32 additions & 15 deletions src/utils/ctem_diags_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
! the zonal mean terms of the TEM equation
!-----------------------------------------------------------------------------
module ctem_diags_mod
use shr_kind_mod, only: r8 => shr_kind_r8
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
Expand Down Expand Up @@ -44,6 +44,7 @@ subroutine ctem_diags_readnl(nlfile)

character(len=*), intent(in) :: nlfile
integer :: unitn, ierr
character(len=cx) :: iomsg

character(len=*), parameter :: prefix = 'ctem_diags_readnl: '

Expand All @@ -54,9 +55,9 @@ subroutine ctem_diags_readnl(nlfile)
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)
read(unitn, ctem_diags_nl, iostat=ierr, iomsg=iomsg)
if (ierr /= 0) then
call endrun(prefix//'ctem_diags_nl: ERROR reading namelist')
call endrun(prefix//'ctem_diags_nl: ERROR reading namelist: '//trim(iomsg))
end if
end if
close(unitn)
Expand Down Expand Up @@ -97,6 +98,8 @@ subroutine ctem_diags_reg()
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
Expand All @@ -112,6 +115,9 @@ subroutine ctem_diags_reg()

! 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
Expand All @@ -134,7 +140,11 @@ subroutine ctem_diags_reg()
nullify(grid_map)

! for the lon-lat grid
allocate(grid_map(4, ((endlon - beglon + 1) * (endlat - beglat + 1))))
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
Expand All @@ -146,7 +156,11 @@ subroutine ctem_diags_reg()
end do
end do

allocate(coord_map(endlat - beglat + 1))
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
Expand All @@ -157,7 +171,11 @@ subroutine ctem_diags_reg()

nullify(coord_map)

allocate(coord_map(endlon - beglon + 1))
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
Expand All @@ -183,10 +201,9 @@ subroutine ctem_diags_init()

! 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-Mean zonal wind', gridname='ctem_lonlat' )
call addfld ('Vtem', (/'lev'/), 'A','m s-1', 'Zonal-Mean meridional wind', gridname='ctem_lonlat' )
call addfld ('Wtem', (/'lev'/), 'A','m s-1', 'Zonal-Mean vertical wind', 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')
Expand All @@ -197,14 +214,14 @@ subroutine ctem_diags_init()
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','Meridional Heat Flux:', gridname='ctem_zm')
call addfld ('WTHzm',(/'lev'/), 'A','K m s-1','Vertical Heat Flux:', gridname='ctem_zm')
call addfld ('UVzm', (/'lev'/), 'A','m2 s-2', 'Meridional Flux of Zonal Momentum', gridname='ctem_zm')
call addfld ('UWzm', (/'lev'/), 'A','m2 s-2', 'Vertical Flux of Zonal Momentum', 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 ,eridional Flux of Zonal Momentum', gridname='ctem_zm')
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo here:

Suggested change
call addfld ('UVzm', (/'lev'/), 'A','m2 s-2', 'Zonal-Mean ,eridional Flux of Zonal Momentum', gridname='ctem_zm')
call addfld ('UVzm', (/'lev'/), 'A','m2 s-2', 'Zonal-Mean meridional Flux of Zonal Momentum', gridname='ctem_zm')

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

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' )
call addfld ('MSKtem',horiz_only, 'A', '1', 'TEM mask', gridname='ctem_lonlat' )

end subroutine ctem_diags_init

Expand Down
67 changes: 53 additions & 14 deletions src/utils/esmf_lonlat_grid_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ subroutine esmf_lonlat_grid_init(nlats_in)
integer, parameter :: minlons_per_pe = 2

integer, allocatable :: petmap(:,:,:)
integer :: petcnt
integer :: petcnt, astat

integer :: lbnd_lat, ubnd_lat, lbnd_lon, ubnd_lon
integer :: lbnd(1), ubnd(1)
Expand All @@ -81,8 +81,14 @@ subroutine esmf_lonlat_grid_init(nlats_in)
nlon = 2*nlat
delx = 360._r8/nlon

allocate(glons(nlon))
allocate(glats(nlat))
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
Expand Down Expand Up @@ -171,20 +177,44 @@ subroutine esmf_lonlat_grid_init(nlats_in)

call mpi_comm_split(mpicom,mytidj,mytid,zonal_comm,ierr)

allocate(mytidi_send(npes))
allocate(mytidj_send(npes))
allocate(mytidi_recv(npes))
allocate(mytidj_recv(npes))
allocate(mytidi_send(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate mytidi_send( array')
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo here:

Suggested change
call endrun(subname//'not able to allocate mytidi_send( array')
call endrun(subname//'not able to allocate mytidi_send array')

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

end if
allocate(mytidj_send(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate mytidj_send array')
end if
allocate(mytidi_recv(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate mytidi_recv array')
end if
allocate(mytidj_recv(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate mytidj_recv array')
end if

mytidi_send = mytidi
mytidj_send = mytidj
call mpi_alltoall(mytidi_send, 1, MPI_INTEGER, mytidi_recv, 1, MPI_INTEGER, mpicom, ierr)
call mpi_alltoall(mytidj_send, 1, MPI_INTEGER, mytidj_recv, 1, MPI_INTEGER, mpicom, ierr)

allocate(nlons_send(npes))
allocate(nlats_send(npes))
allocate(nlons_recv(npes))
allocate(nlats_recv(npes))
allocate(nlons_send(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate nlons_send array')
end if
allocate(nlats_send(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate nlats_send array')
end if
allocate(nlons_recv(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate nlons_recv array')
end if
allocate(nlats_recv(npes), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate nlats_recv array')
end if

nlons_send(:) = mynlons
nlats_send(:) = mynlats
Expand All @@ -195,8 +225,14 @@ subroutine esmf_lonlat_grid_init(nlats_in)
deallocate(nlons_send)
deallocate(nlats_send)

allocate(nlons_task(ntasks_lon))
allocate(nlats_task(ntasks_lat))
allocate(nlons_task(ntasks_lon), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate nlons_task array')
end if
allocate(nlats_task(ntasks_lat), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate nlats_task array')
end if

do i = 1, ntasks_lon
loop1: do n = 1, npes
Expand Down Expand Up @@ -227,7 +263,10 @@ subroutine esmf_lonlat_grid_init(nlats_in)

! set up 2D ESMF lon lat grid

allocate(petmap(ntasks_lon,ntasks_lat,1))
allocate(petmap(ntasks_lon,ntasks_lat,1), stat=astat)
if (astat/=0) then
call endrun(subname//'not able to allocate petmap array')
end if

petcnt = 0
do j = 1,ntasks_lat
Expand Down
6 changes: 3 additions & 3 deletions src/utils/esmf_phys2lonlat_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,14 +127,14 @@ subroutine esmf_phys2lonlat_regrid_3d(physflds, lonlatflds)
use ppgrid, only: pcols, pver, begchunk, endchunk
use phys_grid, only: get_ncols_p

type(fields_bundle_t) :: physflds(nflds)
type(fields_bundle_t) :: lonlatflds(nflds)
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: '
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')
Expand Down
12 changes: 6 additions & 6 deletions src/utils/esmf_phys_mesh_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ subroutine esmf_phys_mesh_init()
end do
allocate(decomp(total_cols), stat=rc)
if (rc/=0) then
call endrun(subname//'not able to allacate decomp')
call endrun(subname//'not able to allocate decomp')
end if

dindex = 0
Expand Down Expand Up @@ -102,17 +102,17 @@ subroutine esmf_phys_mesh_init()

allocate(ownedElemCoords(spatialDim*numOwnedElements), stat=rc)
if (rc/=0) then
call endrun(subname//'not able to allacate ownedElemCoords')
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 allacate lonMesh')
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 allacate latMesh')
call endrun(subname//'not able to allocate latMesh')
end if

call ESMF_MeshGet(physics_grid_mesh, ownedElemCoords=ownedElemCoords)
Expand All @@ -126,14 +126,14 @@ subroutine esmf_phys_mesh_init()
! obtain internally generated cam lats and lons
allocate(lon(total_cols), stat=rc);
if (rc/=0) then
call endrun(subname//'not able to allacate lon')
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 allacate lat')
call endrun(subname//'not able to allocate lat')
end if

lat(:) = 0._r8
Expand Down
29 changes: 6 additions & 23 deletions src/utils/esmf_zonal_mean_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ subroutine esmf_zonal_mean_calc_3d(lonlatarr, zmarr, wght)
real(r8) :: tmparr(lon_beg:lon_end,pver)
real(r8) :: gsum(pver)

real(r8) :: wsum(pver)
real(r8) :: wsums(lat_beg:lat_end,pver)

integer :: numlons, ilat, ilev

Expand All @@ -51,29 +51,13 @@ subroutine esmf_zonal_mean_calc_3d(lonlatarr, zmarr, wght)
! zonal mean
if (present(wght)) then

do ilat = lat_beg, lat_end

tmparr(lon_beg:lon_end,:) = wght(lon_beg:lon_end,ilat,:)
call shr_reprosum_calc(tmparr, wsum, numlons, numlons, pver, gbl_count=nlon, commid=zonal_comm)

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 (wsum(ilev)>0._r8) then
zmarr(ilat,ilev) = gsum(ilev)/wsum(ilev)
else
zmarr(ilat,ilev) = fillvalue
end if
end do

end do
wsums = esmf_zonal_mean_wsums(wght)
call esmf_zonal_mean_masked(lonlatarr, wght, wsums, zmarr)

else

do ilat = lat_beg, lat_end
tmparr(lon_beg:lon_end,:) = lonlatarr(lon_beg:lon_end,ilat,:)
call shr_reprosum_calc(tmparr, gsum, numlons, numlons, pver, gbl_count=nlon, commid=zonal_comm)
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

Expand Down Expand Up @@ -118,15 +102,14 @@ function esmf_zonal_mean_wsums(wght) result(wsums)
real(r8), intent(in) :: wght(lon_beg:lon_end,lat_beg:lat_end,pver)

real(r8) :: wsums(lat_beg:lat_end,pver)
real(r8) :: tmparr(lon_beg:lon_end,pver)
integer :: numlons, ilat

numlons = lon_end-lon_beg+1

do ilat = lat_beg, lat_end

tmparr(lon_beg:lon_end,:) = wght(lon_beg:lon_end,ilat,:)
call shr_reprosum_calc(tmparr, wsums(ilat,1:pver), numlons, numlons, pver, gbl_count=nlon, commid=zonal_comm)
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

Expand Down