Skip to content

Compute, rather than read from input, deformation_coef_* fields at model startup #1419

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
mgduda wants to merge 2 commits into
base: develop
Choose a base branch
from
+250 −330
Open

Compute, rather than read from input, deformation_coef_* fields at model startup #1419

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d62480f
Compute, rather than read from input, deformation_coef_* fields at mo…
mgduda Mar 6, 2026
580e004
Remove deformation_coef_* fields and their computation from init_atmo…
mgduda Mar 6, 2026
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: 0 additions & 5 deletions src/core_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -545,11 +545,6 @@
<var name="deriv_two"/>
<var name="defc_a"/>
<var name="defc_b"/>
<var name="deformation_coef_c2"/>
<var name="deformation_coef_s2"/>
<var name="deformation_coef_cs"/>
<var name="deformation_coef_c"/>
<var name="deformation_coef_s"/>
<var name="coeffs_reconstruct"/>
#ifdef MPAS_CAM_DYCORE
<var name="cell_gradient_coef_x"/>
Expand Down
241 changes: 241 additions & 0 deletions src/core_atmosphere/mpas_atm_core.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -419,6 +419,10 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
logical, pointer :: config_do_restart, config_do_DAcycling
logical, pointer :: on_a_sphere
real (kind=RKIND), pointer :: sphere_radius
Copy link
Collaborator

@abishekg7 abishekg7 Mar 19, 2026

Choose a reason for hiding this comment

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

Is the extra newline here required?

Copy link
Contributor Author

@mgduda mgduda Mar 19, 2026

Choose a reason for hiding this comment

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

Using two blank lines after local variable declarations is something that we do in other parts of the code. See, e.g., lines 1658 - 1659 in this same file.

Copy link
Contributor Author

@mgduda mgduda Mar 19, 2026

Choose a reason for hiding this comment

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

For other examples, see mpas_halo.F and the subroutines therein.

call atm_compute_signs(mesh)
call mpas_pool_get_subpool(block % structs, 'diag', diag)
Expand Down Expand Up @@ -470,6 +474,10 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
!!!!! End compute inverses
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call mpas_pool_get_config(mesh, 'on_a_sphere', on_a_sphere)
call mpas_pool_get_config(mesh, 'sphere_radius', sphere_radius)
call atm_initialize_deformation_weights(mesh, nCells, on_a_sphere, sphere_radius)
call atm_adv_coef_compression(mesh)
call atm_couple_coef_3rd_order(mesh, block % configs)
Expand Down Expand Up @@ -1600,5 +1608,238 @@ subroutine mpas_atm_run_compatibility(dminfo, blockList, streamManager, ierr)

end subroutine mpas_atm_run_compatibility


subroutine atm_initialize_deformation_weights(mesh, nCells, on_a_sphere, sphere_radius)

!
! compute the cell coefficients for the deformation calculations
! WCS, 13 July 2010
!

use mpas_vector_operations, only : mpas_fix_periodicity
use mpas_timer, only : mpas_timer_start, mpas_timer_stop
use mpas_geometry_utils, only : mpas_sphere_angle, mpas_plane_angle, mpas_arc_length

implicit none

type (mpas_pool_type), intent(inout) :: mesh
integer, intent(in) :: nCells
logical, intent(in) :: on_a_sphere
real (kind=RKIND), intent(in) :: sphere_radius

! local variables

real (kind=RKIND), dimension(:,:), pointer :: deformation_coef_c2, deformation_coef_s2, deformation_coef_cs
real (kind=RKIND), dimension(:,:), pointer :: deformation_coef_c, deformation_coef_s
integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, cellsOnCell, verticesOnCell
integer, dimension(:), pointer :: nEdgesOnCell
real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
real (kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
real (kind=RKIND), dimension(:), pointer :: xEdge, yEdge, zEdge, angleEdge

real (kind=RKIND), dimension(nCells) :: theta_abs

real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
real (kind=RKIND), dimension(25) :: thetav, thetat, dl_sphere
real (kind=RKIND) :: dl
integer :: i, ip1, ip2, n
integer :: iCell
real (kind=RKIND) :: pii
real (kind=RKIND), dimension(25) :: xp, yp
real (kind=RKIND) :: xe, ye

integer, dimension(25) :: cell_list

integer :: iv, ie
logical :: do_the_cell
real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, dx, dy

real(kind=RKIND), pointer :: x_period, y_period


call mpas_pool_get_config(mesh, 'x_period', x_period)
call mpas_pool_get_config(mesh, 'y_period', y_period)

call mpas_pool_get_array(mesh, 'deformation_coef_c2', deformation_coef_c2)
call mpas_pool_get_array(mesh, 'deformation_coef_s2', deformation_coef_s2)
call mpas_pool_get_array(mesh, 'deformation_coef_cs', deformation_coef_cs)
call mpas_pool_get_array(mesh, 'deformation_coef_c', deformation_coef_c)
call mpas_pool_get_array(mesh, 'deformation_coef_s', deformation_coef_s)
call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(mesh, 'verticesOnCell', verticesOnCell)
call mpas_pool_get_array(mesh, 'xCell', xCell)
call mpas_pool_get_array(mesh, 'yCell', yCell)
call mpas_pool_get_array(mesh, 'zCell', zCell)
call mpas_pool_get_array(mesh, 'xVertex', xVertex)
call mpas_pool_get_array(mesh, 'yVertex', yVertex)
call mpas_pool_get_array(mesh, 'zVertex', zVertex)
call mpas_pool_get_array(mesh, 'xEdge', xEdge)
call mpas_pool_get_array(mesh, 'yEdge', yEdge)
call mpas_pool_get_array(mesh, 'zEdge', zEdge)
call mpas_pool_get_array(mesh, 'angleEdge', angleEdge)

deformation_coef_c2(:,:) = 0.
deformation_coef_s2(:,:) = 0.
deformation_coef_cs(:,:) = 0.
deformation_coef_c(:,:) = 0.
deformation_coef_s(:,:) = 0.

pii = 2.*asin(1.0)

do iCell = 1, nCells

cell_list(1) = iCell
do i=2,nEdgesOnCell(iCell)+1
cell_list(i) = cellsOnCell(i-1,iCell)
end do
n = nEdgesOnCell(iCell) + 1

! check to see if we are reaching outside the halo

do_the_cell = .true.
do i=1,n
if (cell_list(i) > nCells) do_the_cell = .false.
end do


if (.not. do_the_cell) cycle

! compute poynomial fit for this cell if all needed neighbors exist

if (on_a_sphere) then

! xc holds the center point and the vertex points of the cell,
! normalized to a sphere or radius 1.

xc(1) = xCell(iCell)/sphere_radius
yc(1) = yCell(iCell)/sphere_radius
zc(1) = zCell(iCell)/sphere_radius

do i=2,n
iv = verticesOnCell(i-1,iCell)
xc(i) = xVertex(iv)/sphere_radius
yc(i) = yVertex(iv)/sphere_radius
zc(i) = zVertex(iv)/sphere_radius
end do

!
! In case the current cell center lies at exactly z=1.0, the sphere_angle() routine
! may generate an FPE since the triangle it is given will have a zero side length
! adjacent to the vertex whose angle we are trying to find; in this case, simply
! set the value of theta_abs directly
!
if (zc(1) == 1.0) then
theta_abs(iCell) = pii/2.
else
! theta_abs is the angle to the first vertex from the center, normalized so that
! an eastward pointing vector has a angle of 0.
theta_abs(iCell) = pii/2. - mpas_sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
end if

! here we are constructing the tangent-plane cell.
! thetat is the angle in the (x,y) tangent-plane coordinate from
! the cell center to each vertex, normalized so that an
! eastward pointing vector has a angle of 0.

! dl_sphere is the spherical distance from the cell center
! to the sphere vertex points for the cell.

thetat(1) = theta_abs(iCell)
do i=1,n-1

ip2 = i+2
if (ip2 > n) ip2 = 2

thetav(i) = mpas_sphere_angle( xc(1), yc(1), zc(1), &
xc(i+1), yc(i+1), zc(i+1), &
xc(ip2), yc(ip2), zc(ip2) )
dl_sphere(i) = sphere_radius*mpas_arc_length( xc(1), yc(1), zc(1), &
xc(i+1), yc(i+1), zc(i+1) )
if(i.gt.1) thetat(i) = thetat(i-1)+thetav(i-1)
end do

! xp and yp are the tangent-plane vertex points with the cell center at (0,0)

do i=1,n-1
xp(i) = cos(thetat(i)) * dl_sphere(i)
yp(i) = sin(thetat(i)) * dl_sphere(i)
end do

else ! On an x-y plane

do i=1,n-1
iv = verticesOnCell(i,iCell)
xp(i) = mpas_fix_periodicity(xVertex(iv),xCell(iCell),x_period) - xCell(iCell)
yp(i) = mpas_fix_periodicity(yVertex(iv),yCell(iCell),y_period) - yCell(iCell)
end do

do i=1,n-1
ie = edgesOnCell(i,iCell)
xe = mpas_fix_periodicity(xEdge(ie),xCell(iCell),x_period) - xCell(iCell)
ye = mpas_fix_periodicity(yEdge(ie),yCell(iCell),y_period) - yCell(iCell)
thetat(i) = atan2(ye,xe)
end do

theta_abs(iCell) = thetat(1)

end if

! (1) compute cell area on the tangent plane used in the integrals
! (2) compute angle of cell edge normal vector. here we are repurposing thetat
thetat(1) = theta_abs(iCell)

do i=2,n-1
ip1 = i+1
if (ip1 == n) ip1 = 1
thetat(i) = mpas_plane_angle( 0.0_RKIND, 0.0_RKIND, 0.0_RKIND, &
xp(i)-xp(i-1), yp(i)-yp(i-1), 0.0_RKIND, &
xp(ip1)-xp(i), yp(ip1)-yp(i), 0.0_RKIND, &
0.0_RKIND, 0.0_RKIND, 1.0_RKIND)
thetat(i) = thetat(i) + thetat(i-1)
end do

area_cell = 0.
do i=1,n-1
ip1 = i+1
if (ip1 == n) ip1 = 1
dx = xp(ip1)-xp(i)
dy = yp(ip1)-yp(i)
area_cell = area_cell + 0.25*(xp(i)+xp(ip1))*(yp(ip1)-yp(i)) - 0.25*(yp(i)+yp(ip1))*(xp(ip1)-xp(i))
thetat(i) = atan2(dy,dx)-pii/2.
end do

! coefficients - see documentation for the formulas.

do i=1,n-1
ip1 = i+1
if (ip1 == n) ip1 = 1
dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
sint2 = (sin(thetat(i)))**2
cost2 = (cos(thetat(i)))**2
sint_cost = sin(thetat(i))*cos(thetat(i))
deformation_coef_c2(i,iCell) = dl*cost2/area_cell
deformation_coef_s2(i,iCell) = dl*sint2/area_cell
deformation_coef_cs(i,iCell) = dl*sint_cost/area_cell
deformation_coef_c(i,iCell) = dl*cos(thetat(i))/area_cell
deformation_coef_s(i,iCell) = dl*sin(thetat(i))/area_cell
if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
deformation_coef_c2(i,iCell) = - deformation_coef_c2(i,iCell)
deformation_coef_s2(i,iCell) = - deformation_coef_s2(i,iCell)
deformation_coef_cs(i,iCell) = - deformation_coef_cs(i,iCell)
! deformation_coef_c(i,iCell) = - deformation_coef_c(i,iCell)
! deformation_coef_s(i,iCell) = - deformation_coef_s(i,iCell)
end if

end do

end do

end subroutine atm_initialize_deformation_weights

end module atm_core

25 changes: 0 additions & 25 deletions src/core_init_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -478,11 +478,6 @@
<var name="defc_b" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="cell_gradient_coef_x" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="cell_gradient_coef_y" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="deformation_coef_c2" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="deformation_coef_s2" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="deformation_coef_cs" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="deformation_coef_c" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="deformation_coef_s" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="coeffs_reconstruct" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="cf1" packages="met_stage_in"/>
<var name="cf2" packages="met_stage_in"/>
Expand Down Expand Up @@ -588,11 +583,6 @@
<var name="defc_b"/>
<var name="cell_gradient_coef_x"/>
<var name="cell_gradient_coef_y"/>
<var name="deformation_coef_c2"/>
<var name="deformation_coef_s2"/>
<var name="deformation_coef_cs"/>
<var name="deformation_coef_c"/>
<var name="deformation_coef_s"/>
<var name="coeffs_reconstruct"/>
<var name="cf1" packages="vertical_stage_out;met_stage_out"/>
<var name="cf2" packages="vertical_stage_out;met_stage_out"/>
Expand Down Expand Up @@ -1134,21 +1124,6 @@
<var name="cell_gradient_coef_y" type="real" dimensions="maxEdges nCells" units="m^-1"
description="Coefficients for computing the y (meridional) derivative of a cell-centered variable"/>

<var name="deformation_coef_c2" type="real" dimensions="maxEdges nCells" units="unitless"
description="Coefficients for computing the cos-squared terms of the 3d deformation"/>

<var name="deformation_coef_s2" type="real" dimensions="maxEdges nCells" units="unitless"
description="Coefficients for computing the sine-squared terms of the 3d deformation"/>

<var name="deformation_coef_cs" type="real" dimensions="maxEdges nCells" units="unitless"
description="Coefficients for computing the cos-sine terms of the 3d deformation"/>

<var name="deformation_coef_c" type="real" dimensions="maxEdges nCells" units="unitless"
description="Coefficients for computing the cos-only terms of the 3d deformation"/>

<var name="deformation_coef_s" type="real" dimensions="maxEdges nCells" units="unitless"
description="Coefficients for computing the sine-only terms of the 3d deformation"/>

<!-- arrays required for reconstruction of velocity field -->
<var name="coeffs_reconstruct" type="real" dimensions="R3 maxEdges nCells" units="unitless"
description="Coefficients to reconstruct velocity vectors at cell centers"/>
Expand Down
Loading