Skip to content

Commit

Permalink
(*)Correct memory declarations in MOM_regridding
Browse files Browse the repository at this point in the history 
  Corrected memory declarations in MOM_regridding.F90 in cases where the
vertical size of the input and output grids are not the same.  Although this is
not know to have caused any particular problems, these inconsistencies could
lead to segmentation faults in cases where the target grid (e.g., diagnostic
output) is larger than the input grid (e.g., the model's native grid).  In some
cases, certain grid generation options were only written to work with the same
size of input and output grids, and error handling has been added to these cases
to gracefully bring down the model if they are used with different grid sizes.
All answers are bitwise identical in the MOM6-examples test suite, but it is
conceivable that this could correct subtle (memory-related) issues in some
configurations.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Feb 16, 2022
1 parent 75bf521 commit fc5253f
Showing 1 changed file with 67 additions and 41 deletions.
108 changes: 67 additions & 41 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ module MOM_regridding

! This file is part of MOM6. See LICENSE.md for the license.

use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_error_handler, only : MOM_error, FATAL, WARNING, assert
use MOM_file_parser, only : param_file_type, get_param, log_param
use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data
use MOM_io, only : verify_variable_units, slasher
Expand Down Expand Up @@ -776,11 +776,11 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_
type(regridding_CS), intent(in) :: CS !< Regridding 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(inout) :: h !< Current 3D grid obtained after
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after
!! the last time step
type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...)
real, dimension(SZI_(G),SZJ_(G), CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate
real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface
real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage
logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment
! Local variables
Expand Down Expand Up @@ -827,7 +827,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_
end select ! type of grid

#ifdef __DO_SAFETY_CHECKS__
call check_remapping_grid(G, GV, h, dzInterface,'in regridding_main')
if (CS%nk == GV%ke) call check_remapping_grid(G, GV, h, dzInterface,'in regridding_main')
#endif

end subroutine regridding_main
Expand Down Expand Up @@ -1086,13 +1086,14 @@ 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 [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
!! [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: frac_shelf_h !< Fractional
!! ice shelf coverage [nondim].
! Local variables
real :: nominalDepth, minThickness, totalThickness, dh ! Depths and thicknesses [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zOld, zNew ! Coordinate interface heights [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
integer :: i, j, k, nz
logical :: ice_shelf

Expand Down Expand Up @@ -1144,17 +1145,23 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)
call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) )

#ifdef __DO_SAFETY_CHECKS__
dh=max(nominalDepth,totalThickness)
if (abs(zNew(1)-zOld(1))>(nz-1)*0.5*epsilon(dh)*dh) then
dh = max(nominalDepth,totalThickness)
if (abs(zNew(1)-zOld(1)) > (nz-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
do k=1,nz+1
write(0,*) 'dzInterface(1) = ', dzInterface(i,j,1), epsilon(dh), nz, CS%nk
do k=1,min(nz,CS%nk)+1
write(0,*) k,zOld(k),zNew(k)
enddo
do k=1,nz
do k=min(nz,CS%nk)+2,CS%nk+1
write(0,*) k,zOld(nz+1),zNew(k)
enddo
do k=1,min(nz,CS%nk)
write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),CS%coordinateResolution(k)
enddo
do k=min(nz,CS%nk)+1,CS%nk
write(0,*) k, 0.0, zNew(k)-zNew(k+1), CS%coordinateResolution(k)
enddo
call MOM_error( FATAL, &
'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
endif
Expand Down Expand Up @@ -1183,14 +1190,15 @@ 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 [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
!! [H ~> m or kg m-2]

! Local variables
integer :: i, j, k
integer :: nz
real :: nominalDepth, totalThickness, dh
real, dimension(SZK_(GV)+1) :: zOld, zNew
real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]

nz = GV%ke

Expand Down Expand Up @@ -1227,12 +1235,18 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface )
write(0,*) 'min_thickness=',CS%min_thickness
write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness
write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz,CS%nk
do k=1,nz+1
do k=1,min(nz,CS%nk)+1
write(0,*) k,zOld(k),zNew(k)
enddo
do k=1,CS%nk
do k=min(nz,CS%nk)+2,CS%nk+1
write(0,*) k,zOld(nz+1),zNew(k)
enddo
do k=1,min(nz,CS%nk)
write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k),CS%coordinateResolution(k)
enddo
do k=min(nz,CS%nk)+1,CS%nk
write(0,*) k,0.0,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
Expand Down Expand Up @@ -1266,22 +1280,23 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel
!------------------------------------------------------------------------------

! Arguments
type(regridding_CS), intent(in) :: CS !< Regridding control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
!! [H ~> m or kg m-2]
type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
type(regridding_CS), intent(in) :: CS !< Regridding control structure
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice
!! shelf coverage [nondim]
! Local variables
integer :: nz
integer :: nz ! The number of layers in the input grid
integer :: i, j, k
real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zOld, zNew ! Old and new interface heights [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
real :: totalThickness ! Total thicknesses [H ~> m or kg m-2]
#ifdef __DO_SAFETY_CHECKS__
Expand Down Expand Up @@ -1355,7 +1370,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel
call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) )

#ifdef __DO_SAFETY_CHECKS__
do k = 2,nz
do k=2,CS%nk
if (zNew(k) > zOld(1)) then
write(0,*) 'zOld=',zOld
write(0,*) 'zNew=',zNew
Expand All @@ -1380,12 +1395,18 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel
write(0,*) 'min_thickness=',CS%min_thickness
write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness
write(0,*) 'zNew(1)-zOld(1) = ',zNew(1)-zOld(1),epsilon(dh),nz
do k=1,nz+1
do k=1,min(nz,CS%nk)+1
write(0,*) k,zOld(k),zNew(k)
enddo
do k=min(nz,CS%nk)+2,CS%nk+1
write(0,*) k,zOld(nz+1),zNew(k)
enddo
do k=1,nz
write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1)
enddo
do k=min(nz,CS%nk)+1,CS%nk
write(0,*) k, 0.0, zNew(k)-zNew(k+1), CS%coordinateResolution(k)
enddo
call MOM_error( FATAL, &
'MOM_regridding, build_rho_grid: top surface has moved!!!' )
endif
Expand Down Expand Up @@ -1416,10 +1437,10 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_she
!! coverage [nondim]

! Local variables
real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(SZK_(GV)) :: p_col ! Layer center pressure in the input column [R L2 T-2 ~> Pa]
real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa]
real, dimension(CS%nk+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
integer :: i, j, k, nki
real :: depth, nominalDepth
real :: h_neglect, h_neglect_edge
Expand Down Expand Up @@ -1496,10 +1517,10 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
!! thermodynamic variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
type(regridding_CS), intent(in) :: CS !< Regridding control structure
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
!! [H ~> m or kg m-2]
type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
type(regridding_CS), intent(in) :: CS !< Regridding control structure

! local variables
integer :: i, j, k, nz ! indices and dimension lengths
Expand All @@ -1512,6 +1533,9 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS)

nz = GV%ke

call assert((GV%ke == CS%nk), "build_grid_adaptive is only written to work "//&
"with the same number of input and target layers.")

! position surface at z = 0.
zInt(:,:,1) = 0.

Expand Down Expand Up @@ -1562,13 +1586,13 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS)
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position
type(regridding_CS), intent(in) :: CS !< Regridding control structure
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position

real, dimension(SZK_(GV)+1) :: z_col ! Interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa]
real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(SZK_(GV)) :: p_col ! Layer center pressure in the input column [R L2 T-2 ~> Pa]
real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(CS%nk+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]

! Local variables
real :: depth ! Depth of the ocean relative to the mean sea surface height in thickness units [H ~> m or kg m-2]
Expand All @@ -1585,8 +1609,10 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS)

nz = GV%ke

if (.not.CS%target_density_set) call MOM_error(FATAL, "build_grid_SLight : "//&
"Target densities must be set before build_grid_SLight is called.")
call assert((GV%ke == CS%nk), "build_grid_SLight is only written to work "//&
"with the same number of input and target layers.")
call assert(CS%target_density_set, "build_grid_SLight : "//&
"Target densities must be set before build_grid_SLight is called.")

! Build grid based on target interface densities
do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
Expand Down Expand Up @@ -1691,13 +1717,13 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS )
!------------------------------------------------------------------------------

! Arguments
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 !< Original layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface
!! depth [H ~> m or kg m-2]
real, intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
type(regridding_CS), intent(in) :: CS !< Regridding 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 !< Original layer thicknesses [H ~> m or kg m-2]
type(regridding_CS), intent(in) :: CS !< Regridding control structure
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface
!! depth [H ~> m or kg m-2]
real, intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]

! Local variables
integer :: i, j, k
Expand Down

0 comments on commit fc5253f

Please sign in to comment.