From 27c13b9b50c9b7341c630dcc5bfcff4d33606e36 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 26 Oct 2022 15:56:02 -0400 Subject: [PATCH] +Add ALE_remap_interface_vals and ALE_remap_vertex_vals Added ALE_remap_interface_vals and ALE_remap_vertex_vals to handle the interpolation of variables that are at the interfaces atop tracer cells or above the corners of the tracers cells from one grid to another. Because these are not yet used (but have been tested in calls that will be added with the next commit, all answers are bitwise identical, but there are two new publicly visible routines. --- src/ALE/MOM_ALE.F90 | 85 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 9f8b2336f9..47ebf30bfb 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -129,6 +129,8 @@ module MOM_ALE public ALE_remap_scalar public ALE_remap_tracers public ALE_remap_velocities +public ALE_remap_interface_vals +public ALE_remap_vertex_vals public ALE_PLM_edge_values public TS_PLM_edge_values public TS_PPM_edge_values @@ -974,6 +976,89 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old, h_new, u, v, OBC, dzInterface, end subroutine ALE_remap_velocities +!> Interpolate to find an updated array of values at interfaces after remapping. +subroutine ALE_remap_interface_vals(CS, G, GV, h_old, h_new, int_val) + type(ALE_CS), intent(in) :: CS !< ALE control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & + intent(inout) :: int_val !< The interface values to interpolate [A] + + real :: val_src(GV%ke) ! A column of interface values on the source grid [A] + real :: val_tgt(GV%ke) ! A column of interface values on the target grid [A] + real :: h_src(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] + real :: h_tgt(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] + real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] + logical :: show_call_tree + integer :: i, j, k, nz + + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j)>0.) then + do k=1,nz + h_src(k) = h_old(i,j,k) + h_tgt(k) = h_new(i,j,k) + enddo + + do K=1,nz+1 + val_src(K) = int_val(i,j,K) + enddo + + call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.) + + do K=1,nz+1 + int_val(i,j,K) = val_tgt(K) + enddo + endif ; enddo ; enddo + +end subroutine ALE_remap_interface_vals + +!> Interpolate to find an updated array of values at vertices of tracer cells after remapping. +subroutine ALE_remap_vertex_vals(CS, G, GV, h_old, h_new, vert_val) + type(ALE_CS), intent(in) :: CS !< ALE control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid + !! [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), & + intent(inout) :: vert_val !< The interface values to interpolate [A] + + real :: val_src(GV%ke) ! A column of interface values on the source grid [A] + real :: val_tgt(GV%ke) ! A column of interface values on the target grid [A] + real :: h_src(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] + real :: h_tgt(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] + real :: I_mask_sum ! The inverse of the tracer point masks surrounding a corner [nondim] + real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] + logical :: show_call_tree + integer :: i, j, k, nz + + do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB + if ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) > 0.0 ) then + I_mask_sum = 1.0 / ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1))) + + do k=1,nz + h_src(k) = ((G%mask2dT(i,j) * h_old(i,j,k) + G%mask2dT(i+1,j+1) * h_old(i+1,j+1,k)) + & + (G%mask2dT(i+1,j) * h_old(i+1,j,k) + G%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * I_mask_sum + h_tgt(k) = ((G%mask2dT(i,j) * h_new(i,j,k) + G%mask2dT(i+1,j+1) * h_new(i+1,j+1,k)) + & + (G%mask2dT(i+1,j) * h_new(i+1,j,k) + G%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * I_mask_sum + enddo + + do K=1,nz+1 + val_src(K) = vert_val(I,J,K) + enddo + + call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.) + + do K=1,nz+1 + vert_val(I,J,K) = val_tgt(K) + enddo + endif ; enddo ; enddo + +end subroutine ALE_remap_vertex_vals !> Mask out thicknesses to 0 when their running sum exceeds a specified value. subroutine apply_partial_cell_mask(h1, h_mask)