Skip to content

Use ESMF conservative regridding #1477

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

Open
fvitt wants to merge 1 commit into
base: cam_development
Choose a base branch
from
Open

Use ESMF conservative regridding #1477

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
85 changes: 65 additions & 20 deletions src/utils/ctem_diags_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,6 +243,12 @@ subroutine ctem_diags_calc(phys_state)
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), target :: uv_phys(pver,pcols,begchunk:endchunk)
real(r8), target :: uw_phys(pver,pcols,begchunk:endchunk)
real(r8), target :: vt_phys(pver,pcols,begchunk:endchunk)
real(r8), target :: wt_phys(pver,pcols,begchunk:endchunk)

real(r8) :: ps_phys(pcols,begchunk:endchunk)

real(r8), target :: u_lonlat(beglon:endlon,beglat:endlat,pver)
Expand All @@ -258,16 +264,28 @@ subroutine ctem_diags_calc(phys_state)
real(r8) :: wi_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8) :: ti_lonlat(beglon:endlon,beglat:endlat,pver)

real(r8), target :: uv_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8), target :: uw_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8), target :: vt_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8), target :: wt_lonlat(beglon:endlon,beglat:endlat,pver)

real(r8) :: uvi_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8) :: uwi_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8) :: vti_lonlat(beglon:endlon,beglat:endlat,pver)
real(r8) :: wti_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) :: uv_zm(beglat:endlat,pver)
real(r8) :: uw_zm(beglat:endlat,pver)
real(r8) :: vt_zm(beglat:endlat,pver)
real(r8) :: wt_zm(beglat:endlat,pver)

real(r8) :: ps_zm(beglat:endlat)

real(r8) :: vtp(beglon:endlon,beglat:endlat,pver)
real(r8) :: wtp(beglon:endlon,beglat:endlat,pver)
Expand Down Expand Up @@ -316,8 +334,14 @@ subroutine ctem_diags_calc(phys_state)
! mid point press
p_phys(:,i,lchnk) = phys_state(lchnk)%pmid(i,:)

! surface pressure
ps_phys(i,lchnk) = phys_state(lchnk)%ps(i)

uv_phys(:,i,lchnk) = u_phys(:,i,lchnk)*v_phys(:,i,lchnk)
uw_phys(:,i,lchnk) = u_phys(:,i,lchnk)*w_phys(:,i,lchnk)
vt_phys(:,i,lchnk) = v_phys(:,i,lchnk)*t_phys(:,i,lchnk)
wt_phys(:,i,lchnk) = w_phys(:,i,lchnk)*t_phys(:,i,lchnk)

end do
end do

Expand All @@ -332,12 +356,20 @@ subroutine ctem_diags_calc(phys_state)
physflds(3)%fld => w_phys
physflds(4)%fld => t_phys
physflds(5)%fld => p_phys
physflds(6)%fld => uv_phys
physflds(7)%fld => uw_phys
physflds(8)%fld => vt_phys
physflds(9)%fld => wt_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
lonlatflds(6)%fld => uv_lonlat
lonlatflds(7)%fld => uw_lonlat
lonlatflds(8)%fld => vt_lonlat
lonlatflds(9)%fld => wt_lonlat

call esmf_phys2lonlat_regrid(physflds, lonlatflds)

Expand Down Expand Up @@ -378,6 +410,18 @@ subroutine ctem_diags_calc(phys_state)
call lininterp( t_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
ti_lonlat(i,j,:), pref_mid(:), pver )

call lininterp( uv_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
uvi_lonlat(i,j,:), pref_mid(:), pver )

call lininterp( uw_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
uwi_lonlat(i,j,:), pref_mid(:), pver )

call lininterp( vt_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
vti_lonlat(i,j,:), pref_mid(:), pver )

call lininterp( wt_lonlat(i,j,:), p_lonlat(i,j,:), pver, &
wti_lonlat(i,j,:), pref_mid(:), pver )

end do
end do

Expand All @@ -394,32 +438,33 @@ subroutine ctem_diags_calc(phys_state)
call esmf_zonal_mean_masked(wi_lonlat, wght, wsums, w_zm)
call esmf_zonal_mean_masked(ti_lonlat, wght, wsums, t_zm)

call esmf_zonal_mean_masked(uvi_lonlat, wght, wsums, uv_zm)
call esmf_zonal_mean_masked(uwi_lonlat, wght, wsums, uw_zm)
call esmf_zonal_mean_masked(vti_lonlat, wght, wsums, vt_zm)
call esmf_zonal_mean_masked(wti_lonlat, wght, wsums, wt_zm)

call t_stopf('ctem_diags_calc-zonal_mean-uvwt')

call t_startf('ctem_diags_calc-calc_dev_flx')

! Calculate zonal deviations and fluxes
! Calculate 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)

! u'v' = (uv)zm - uzm*vzm

uvp(i,j,k) = uv_zm(j,k) - (u_zm(j,k)*v_zm(j,k))
uwp(i,j,k) = uw_zm(j,k) - (u_zm(j,k)*w_zm(j,k))
vtp(i,j,k) = vt_zm(j,k) - (v_zm(j,k)*t_zm(j,k))
wtp(i,j,k) = wt_zm(j,k) - (w_zm(j,k)*t_zm(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
uvp(i,j,k) = fillvalue
uwp(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
Expand Down
35 changes: 34 additions & 1 deletion src/utils/esmf_lonlat_grid_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ module esmf_lonlat_grid_mod

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, only: ESMF_KIND_R8, ESMF_INDEX_GLOBAL, ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER
use esmf_check_error_mod, only: check_esmf_error

implicit none
Expand Down Expand Up @@ -287,6 +287,7 @@ subroutine esmf_lonlat_grid_init(nlats_in)

! Set coordinates:

! cell centers
call ESMF_GridAddCoord(lonlat_grid, staggerloc=ESMF_STAGGERLOC_CENTER, rc=ierr)
call check_esmf_error(ierr, subname//'ESMF_GridAddCoord ERROR')

Expand Down Expand Up @@ -314,6 +315,38 @@ subroutine esmf_lonlat_grid_init(nlats_in)
end do
end if

! cell corners
call ESMF_GridAddCoord(lonlat_grid, staggerloc=ESMF_STAGGERLOC_CORNER, rc=ierr)
call check_esmf_error(ierr, subname//'ESMF_GridAddCoord CORNER ERROR')

if (mytid<npes) then
call ESMF_GridGetCoord(lonlat_grid, coordDim=1, &
computationalLBound=lbnd, computationalUBound=ubnd, &
farrayPtr=coordX, staggerloc=ESMF_STAGGERLOC_CORNER, rc=ierr)
call check_esmf_error(ierr, subname//'ESMF_GridGetCoord for longitude corner coords ERROR')

lbnd_lon = lbnd(1)
ubnd_lon = ubnd(1)
do i = lbnd_lon, ubnd_lon
coordX(i) = glons(i) - (0.5_r8 * delx)
end do

call ESMF_GridGetCoord(lonlat_grid, coordDim=2, &
computationalLBound=lbnd, computationalUBound=ubnd, &
farrayPtr=coordY, staggerloc=ESMF_STAGGERLOC_CORNER, rc=ierr)
call check_esmf_error(ierr, subname//'ESMF_GridGetCoord for latitude corner coords ERROR')

lbnd_lat = lbnd(1)
ubnd_lat = min( ubnd(1), nlat)
do i = lbnd_lat, ubnd_lat
coordY(i) = glats(i) - (0.5_r8 * dely)
end do
if (ubnd(1) == nlat+1) then
coordY(nlat+1) = 90._r8
end if

end if

deallocate(nlons_task)
deallocate(nlats_task)

Expand Down
26 changes: 12 additions & 14 deletions src/utils/esmf_phys2lonlat_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ module esmf_phys2lonlat_mod
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_REGRIDMETHOD_CONSERVE
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
Expand Down Expand Up @@ -46,7 +46,7 @@ module esmf_phys2lonlat_mod
real(r8), pointer :: fld(:,:,:) => null()
end type fields_bundle_t

integer, parameter :: nflds = 5
integer, parameter :: nflds = 9

contains

Expand Down Expand Up @@ -101,21 +101,19 @@ subroutine esmf_phys2lonlat_init()
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)
regridMethod=ESMF_REGRIDMETHOD_CONSERVE, &
routeHandle=rh_phys2lonlat_3d, &
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)
regridMethod=ESMF_REGRIDMETHOD_CONSERVE, &
routeHandle=rh_phys2lonlat_2d, &
srcTermProcessing=smm_srctermproc, &
pipelineDepth=smm_pipelinedep, &
rc=rc)
call check_esmf_error(rc, subname//'ESMF_FieldRegridStore 3D routehandle ERROR')

end subroutine esmf_phys2lonlat_init
Expand Down