From dae13f5e9c1743ec28406749fa08a2d342a63a46 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 5 Mar 2018 14:06:45 -0500 Subject: [PATCH] Fixed build_sigma_column() to generate CS%nk+1 interfaces - As other coordinate generating kernels, build_sigma_column() now generates CS%nk+1 interfaces rather than the (now removed) dummy argument, nk. --- src/ALE/MOM_regridding.F90 | 14 +++++++------- src/ALE/coord_sigma.F90 | 15 +++++++-------- src/framework/MOM_diag_remap.F90 | 2 +- 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 654ac91097..7bc3f72be4 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1112,7 +1112,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) 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 !< Layer thicknesses, in H - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H. + real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth in H. real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage. ! Local variables integer :: i, j, k @@ -1212,7 +1212,7 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) 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 !< Layer thicknesses, in H - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H. + real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth in H. ! Local variables integer :: i, j, k @@ -1239,7 +1239,7 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) totalThickness = totalThickness + h(i,j,k) end do - call build_sigma_column(CS%sigma_CS, nz, nominalDepth, totalThickness, zNew) + call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew) ! Calculate the final change in grid position after blending new and old grids zOld(nz+1) = -nominalDepth @@ -1251,21 +1251,21 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) #ifdef __DO_SAFETY_CHECKS__ dh=max(nominalDepth,totalThickness) - if (abs(zNew(1)-zOld(1))>(nz-1)*0.5*epsilon(dh)*dh) then + if (abs(zNew(1)-zOld(1))>(CS%nk-1)*0.5*epsilon(dh)*dh) then write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness - write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz + write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz,CS%nk do k=1,nz+1 write(0,*) k,zOld(k),zNew(k) enddo - do k=1,nz + do k=1,CS%nk write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k),CS%coordinateResolution(k) enddo call MOM_error( FATAL, & 'MOM_regridding, build_sigma_grid: top surface has moved!!!' ) endif dzInterface(i,j,1) = 0. - dzInterface(i,j,nz+1) = 0. + dzInterface(i,j,CS%nk+1) = 0. #endif end do diff --git a/src/ALE/coord_sigma.F90 b/src/ALE/coord_sigma.F90 index c353461e27..416ab757e2 100644 --- a/src/ALE/coord_sigma.F90 +++ b/src/ALE/coord_sigma.F90 @@ -59,18 +59,17 @@ end subroutine set_sigma_params !> Build a sigma coordinate column -subroutine build_sigma_column(CS, nz, depth, totalThickness, zInterface) - type(sigma_CS), intent(in) :: CS !< Coordinate control structure - integer, intent(in) :: nz !< Number of levels - real, intent(in) :: depth !< Depth of ocean bottom (positive in m) - real, intent(in) :: totalThickness !< Column thickness (positive in m) - real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces +subroutine build_sigma_column(CS, depth, totalThickness, zInterface) + type(sigma_CS), intent(in) :: CS !< Coordinate control structure + real, intent(in) :: depth !< Depth of ocean bottom (positive in m) + real, intent(in) :: totalThickness !< Column thickness (positive in m) + real, dimension(CS%nk+1), intent(inout) :: zInterface !< Absolute positions of interfaces ! Local variables integer :: k - zInterface(nz+1) = -depth - do k = nz,1,-1 + zInterface(CS%nk+1) = -depth + do k = CS%nk,1,-1 zInterface(k) = zInterface(k+1) + (totalThickness * CS%coordinateResolution(k)) ! Adjust interface position to accomodate inflating layers ! without disturbing the interface above diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 31ce2f08e8..03bb65feb9 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -268,7 +268,7 @@ subroutine diag_remap_update(remap_cs, G, GV, h, T, S, eqn_of_state) G%bathyT(i,j)*GV%m_to_H, sum(h(i,j,:)), & zInterfaces, zScale=GV%m_to_H) elseif (remap_cs%vertical_coord == coordinateMode('SIGMA')) then - call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), nz, & + call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), & GV%m_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then call build_rho_column(get_rho_CS(remap_cs%regrid_cs), G%ke, &