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.
nikizadehgfdl · Mar 5, 2018 · dae13f5 · dae13f5
1 parent a46f701
commit dae13f5
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 16 deletions.
diff --git 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
@@ -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
@@ -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, &