Merge branch 'dev/master' of github.com:NOAA-GFDL/MOM6 into dev/master

src/ALE/MOM_ALE.F90

@@ -705,7 +705,7 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
   if (present(old_remap)) use_remapping_core_w = old_remap
   n_points = nk_src
 
-!$OMP parallel default(none) shared(CS,G,h_src,s_src,h_dst,s_dst &
+!$OMP parallel default(none) shared(CS,G,GV,h_src,s_src,h_dst,s_dst &
 !$OMP                               ,ignore_vanished_layers, use_remapping_core_w, nk_src,dx ) &
 !$OMP                        firstprivate(n_points)
 !$OMP do
@@ -917,6 +917,7 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz )
   real :: filt_len, strat_tol, index_scale
   real :: tmpReal, compress_fraction
   real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
+  real :: height_of_rigid_surface
   integer :: nz_fixed_sfc
   real :: rho_target(GV%ke+1) ! Target density used in HYBRID mode
   real, dimension(size(dz))   :: h_max  ! Maximum layer thicknesses, in m.
@@ -951,6 +952,15 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz )
   call initialize_regridding( GV%ke, coordMode, interpScheme, regridCS, &
                               compressibility_fraction=compress_fraction )
 
+  if (coordMode(1:2) == 'Z*') then
+    call get_param(param_file, mod, "ZSTAR_RIGID_SURFACE_THRESHOLD", height_of_rigid_surface, &
+                 "A threshold height used to detect the presence of a rigid-surface\n"//&
+                 "depressing the upper-surface of the model, such as an ice-shelf.\n"//&
+                 "This is a temporary work around for initialization under an ice-shelf.", &
+                 units='m', default=-1.E30)
+    call set_regrid_params( regridCS, height_of_rigid_surface=height_of_rigid_surface*GV%m_to_H)
+  endif
+
   call get_param(param_file, mod, "ALE_COORDINATE_CONFIG", string, &
                  "Determines how to specify the coordinate\n"//&
                  "resolution. Valid options are:\n"//&

src/ALE/MOM_regridding.F90

@@ -133,6 +133,13 @@ module MOM_regridding
   !! If false, integrate from the bottom upward, as does the rest of the model.
   logical :: integrate_downward_for_e = .true.
 
+
+  !> The height to use as a threshold for detection of a "rigid lid" such as an ice shelf
+  !! base. If the free-surface is depressed below this height the z* coordinate generator
+  !! assumes the depression is due to a rigid surface. This is a kludge until we provide
+  !! the position of such a surface to the ALE code.
+  real :: height_of_rigid_surface = -1.E30
+
 end type
 
 ! The following routines are visible to the outside world
@@ -337,7 +344,7 @@ subroutine calc_h_new_by_dz(G, GV, h, dzInterface, h_new)
   ! Local variables
   integer :: i, j, k
 
-!$OMP parallel do default(none) shared(G,h,dzInterface,h_new)
+!$OMP parallel do default(none) shared(G,GV,h,dzInterface,h_new)
   do j = G%jsc-1,G%jec+1
     do i = G%isc-1,G%iec+1
       if (G%mask2dT(i,j)>0.) then
@@ -362,7 +369,7 @@ subroutine check_remapping_grid( G, GV, h, dzInterface, msg )
   ! Local variables
   integer :: i, j
 
-!$OMP parallel do default(none) shared(G,h,dzInterface,msg)
+!$OMP parallel do default(none) shared(G,GV,h,dzInterface,msg)
   do j = G%jsc-1,G%jec+1
     do i = G%isc-1,G%iec+1
       if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, G%bathyT(i,j), h(i,j,:), dzInterface(i,j,:), msg )
@@ -602,13 +609,19 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface )
         totalThickness = totalThickness + h(i,j,k)
       end do
 
-      call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew)
-
       zOld(nz+1) = - nominalDepth
       do k = nz,1,-1
         zOld(k) = zOld(k+1) + h(i,j,k)
       enddo
 
+      if (totalThickness-nominalDepth<CS%height_of_rigid_surface) then
+        call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew, &
+                                z_rigid_top = totalThickness-nominalDepth, &
+                                eta_orig = zOld(1))
+      else
+        call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew)
+      endif
+
       ! Calculate the final change in grid position after blending new and old grids
       call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) )
 
@@ -637,40 +650,74 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface )
 end subroutine build_zstar_grid
 
 !> Builds a z* coordinate with a minimum thickness
-subroutine build_zstar_column( CS, nz, depth, total_thickness, zInterface)
+subroutine build_zstar_column(CS, nz, depth, total_thickness, zInterface, z_rigid_top, eta_orig)
   type(regridding_CS),   intent(in)    :: CS !< Regridding control structure
   integer,               intent(in)    :: nz !< Number of levels
   real,                  intent(in)    :: depth !< Depth of ocean bottom (positive in m)
   real,                  intent(in)    :: total_thickness !< Column thickness (positive in m)
   real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces
+  real, optional,        intent(in)    :: z_rigid_top !< The height of a rigid top (negative in m)
+  real, optional,        intent(in)    :: eta_orig !< The actual original height of the top (m)
   ! Local variables
-  real :: eta, stretching, dh, min_thickness
+  real :: eta, stretching, dh, min_thickness, z0_top, z_star
   integer :: k
+  logical :: new_zstar_def
 
+  new_zstar_def = .false.
   min_thickness = min( CS%min_thickness, total_thickness/real(nz) )
+  z0_top = 0.
+  if (present(z_rigid_top)) then
+    z0_top = z_rigid_top
+    new_zstar_def = .true.
+  endif
 
-  ! Position of free-surface
+  ! Position of free-surface (or the rigid top, for which eta ~ z0_top)
   eta = total_thickness - depth
+  if (present(eta_orig)) eta = eta_orig
+
+  ! Conventional z* coordinate:
+  !   z* = (z-eta) / stretching   where stretching = (H+eta)/H
+  !   z = eta + stretching * z*
+  ! The above gives z*(z=eta) = 0, z*(z=-H) = -H.
+  ! With a rigid top boundary at eta = z0_top then
+  !   z* = z0 + (z-eta) / stretching   where stretching = (H+eta)/(H+z0)
+  !   z = eta + stretching * (z*-z0) * stretching
+  stretching = total_thickness / ( depth + z0_top )
+
+  if (new_zstar_def) then
+    ! z_star is the notional z* coordinate in absence of upper/lower topography
+    z_star = 0. ! z*=0 at the free-surface
+    zInterface(1) = eta ! The actual position of the top of the column
+    do k = 2,nz
+      z_star = z_star - CS%coordinateResolution(k-1)
+      ! This ensures that z is below a rigid upper surface (ice shelf bottom)
+      zInterface(k) = min( eta + stretching * ( z_star - z0_top ), z0_top )
+      ! This ensures that the layer in inflated
+      zInterface(k) = min( zInterface(k), zInterface(k-1) - min_thickness )
+      ! This ensures that z is above or at the topography
+      zInterface(k) = max( zInterface(k), -depth + real(nz+1-k) * min_thickness )
+    enddo
+    zInterface(nz+1) = -depth
 
-  ! z* = (z-eta) / stretching   where stretching = (H+eta)/H
-  ! z = eta + stretching * z*
-  stretching = total_thickness / depth
-
-  ! Integrate down from the top for a notional new grid, ignoring topography
-  zInterface(1) = eta
-  do k = 1,nz
-    dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing
-    zInterface(k+1) = zInterface(k) - dh
-  enddo
+  else
+    ! Integrate down from the top for a notional new grid, ignoring topography
+    ! The starting position is offset by z0_top which, if z0_top<0, will place
+    ! interfaces above the rigid boundary.
+    zInterface(1) = eta
+    do k = 1,nz
+      dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing
+      zInterface(k+1) = zInterface(k) - dh
+    enddo
 
-  ! Integrating up from the bottom adjusting interface position to accommodate
-  ! inflating layers without disturbing the interface above
-  zInterface(nz+1) = -depth
-  do k = nz,1,-1
-    if ( zInterface(k) < (zInterface(k+1) + min_thickness) ) then
-      zInterface(k) = zInterface(k+1) + min_thickness
-    endif
-  enddo
+    ! Integrating up from the bottom adjusting interface position to accommodate
+    ! inflating layers without disturbing the interface above
+    zInterface(nz+1) = -depth
+    do k = nz,1,-1
+      if ( zInterface(k) < (zInterface(k+1) + min_thickness) ) then
+        zInterface(k) = zInterface(k+1) + min_thickness
+      endif
+    enddo
+  endif
 
 end subroutine build_zstar_column
 
@@ -1699,7 +1746,7 @@ subroutine adjust_interface_motion( nk, min_thickness, h_old, dz_int )
                      'implied h<0 is larger than roundoff!')
     endif
   enddo
-  do k = nk,1,-1
+  do k = nk,2,-1
     h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
     if (h_new<min_thickness) dz_int(k) = ( dz_int(k+1) - h_old(k) ) + min_thickness ! Implies next h_new = min_thickness
     h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
@@ -1714,7 +1761,7 @@ subroutine adjust_interface_motion( nk, min_thickness, h_old, dz_int )
                      'Repeated adjustment for roundoff h<0 failed!')
     endif
   enddo
-  if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
+ !if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
 
 end subroutine adjust_interface_motion
 
@@ -2597,7 +2644,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
              depth_of_time_filter_shallow, depth_of_time_filter_deep, &
              compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, &
              nlay_ML_to_interior, fix_haloclines, halocline_filt_len, &
-             halocline_strat_tol, integrate_downward_for_e)
+             halocline_strat_tol, integrate_downward_for_e, height_of_rigid_surface)
   type(regridding_CS), intent(inout) :: CS !< Regridding control structure
   logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
   real,    optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the new grid (m)
@@ -2613,6 +2660,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
   real,    optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for spuriously unstable water mass profiles (m)
   real,    optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic halocline region.
   logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward from the top.
+  real,    optional, intent(in) :: height_of_rigid_surface !< Threshold height for detection of a rigid upper surface.
 
   if (present(boundary_extrapolation)) CS%boundary_extrapolation = boundary_extrapolation
   if (present(min_thickness)) CS%min_thickness = min_Thickness
@@ -2640,6 +2688,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
     CS%halocline_strat_tol = halocline_strat_tol
   endif
   if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e
+  if (present(height_of_rigid_surface)) CS%height_of_rigid_surface = height_of_rigid_surface
 
 end subroutine set_regrid_params
 

src/core/MOM_barotropic.F90

@@ -3541,6 +3541,7 @@ subroutine find_face_areas(Datu, Datv, G, GV, CS, MS, eta, halo, add_max)
   else
 !$OMP do
     do j=js-hs,je+hs ; do I=is-1-hs,ie+hs
+      Datu(I, j) = 0.0
       !Would be "if (G%mask2dCu(I,j)>0.) &" is G was valid on BT domain
       if (CS%bathyT(i+1,j)+CS%bathyT(i,j)>0.) &
         Datu(I,j) = 2.0*CS%dy_Cu(I,j) * GV%m_to_H * &
@@ -3549,6 +3550,7 @@ subroutine find_face_areas(Datu, Datv, G, GV, CS, MS, eta, halo, add_max)
     enddo ; enddo
 !$OMP do
     do J=js-1-hs,je+hs ; do i=is-hs,ie+hs
+      Datv(i, J) = 0.0
       !Would be "if (G%mask2dCv(i,J)>0.) &" is G was valid on BT domain
       if (CS%bathyT(i,j+1)+CS%bathyT(i,j)>0.) &
         Datv(i,J) = 2.0*CS%dx_Cv(i,J) * GV%m_to_H * &

src/initialization/MOM_state_initialization.F90

@@ -77,6 +77,7 @@ module MOM_state_initialization
 use MOM_ALE, only : ALE_remap_scalar, ALE_build_grid
 use MOM_regridding, only : regridding_CS, set_regrid_params
 use MOM_remapping, only : remapping_CS, initialize_remapping
+use MOM_remapping, only : remapping_core_h
 use MOM_tracer_initialization_from_Z, only : horiz_interp_and_extrap_tracer
 
 implicit none ; private
@@ -846,7 +847,7 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h)
   type(ocean_grid_type),                    intent(in)    :: G !< Ocean grid structure
   type(verticalGrid_type),                  intent(in)    :: GV !< Vertical grid structure
   type(ALE_CS),                             pointer       :: ALE_CSp !< ALE control structure
-  type(thermo_var_ptrs),                    intent(in)    :: tv !< Thermodynamics structure
+  type(thermo_var_ptrs),                    intent(inout) :: tv !< Thermodynamics structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness (H units, m or Pa)
   ! Local variables
   character(len=200) :: mod = "trim_for_ice"
@@ -856,6 +857,8 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h)
   character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
   real :: scale_factor, min_thickness
   integer :: i, j, k
+  logical :: use_remapping
+  type(remapping_CS), pointer :: remap_CS => NULL()
 
   call get_param(PF, mod, "SURFACE_PRESSURE_FILE", p_surf_file, &
                  "The initial condition file for the surface height.", &
@@ -875,6 +878,13 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h)
   if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:)
   call get_param(PF, mod, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
                  units='m', default=1.e-3)
+  call get_param(PF, mod, "TRIMMING_USES_REMAPPING", use_remapping, &
+                 'When trimming the column, also remap T and S.', &
+                 default=.false.)
+  if (use_remapping) then
+    allocate(remap_CS)
+    call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true.)
+  endif
 
   ! Find edge values of T and S used in reconstructions
   if ( associated(ALE_CSp) ) then ! This should only be associated if we are in ALE mode
@@ -893,37 +903,42 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h)
 
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
     call cut_off_column_top(GV%ke, tv, GV%Rho0, G%G_earth, G%bathyT(i,j), min_thickness, &
-               tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:))
+               tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), &
+               p_surf(i,j), h(i,j,:), remap_CS)
   enddo ; enddo
 
 end subroutine trim_for_ice
 
 !> Adjust the layer thicknesses by cutting away the top at the depth where the hydrostatic
 !! pressure matches p_surf
-subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, T, T_t, T_b, S, S_t, S_b, p_surf, h)
+subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, &
+                              T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS)
   integer,               intent(in)    :: nk !< Number of layers
   type(thermo_var_ptrs), intent(in)    :: tv !< Thermodynamics structure
   real,                  intent(in)    :: Rho0 !< Reference density (kg/m3)
   real,                  intent(in)    :: G_earth !< Gravitational acceleration (m/s2)
   real,                  intent(in)    :: depth !< Depth of ocean column (m)
   real,                  intent(in)    :: min_thickness !< Smallest thickness allowed (m)
-  real, dimension(nk),   intent(in)    :: T !< 
-  real, dimension(nk),   intent(in)    :: T_t !< 
-  real, dimension(nk),   intent(in)    :: T_b !< 
-  real, dimension(nk),   intent(in)    :: S !< 
-  real, dimension(nk),   intent(in)    :: S_t !< 
-  real, dimension(nk),   intent(in)    :: S_b !< 
+  real, dimension(nk),   intent(inout) :: T !< Layer mean temperature
+  real, dimension(nk),   intent(in)    :: T_t !< Temperature at top of layer
+  real, dimension(nk),   intent(in)    :: T_b !< Temperature at bottom of layer
+  real, dimension(nk),   intent(inout) :: S !< Layer mean salinity
+  real, dimension(nk),   intent(in)    :: S_t !< Salinity at top of layer
+  real, dimension(nk),   intent(in)    :: S_b !< Salinity at bottom of layer
   real,                  intent(in)    :: p_surf !< Imposed pressure on ocean at surface (Pa)
   real, dimension(nk),   intent(inout) :: h !< Layer thickness (H units, m or Pa)
+  type(remapping_CS),    pointer       :: remap_CS ! Remapping structure for remapping T and S, if associated
   ! Local variables
   real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions
+  real, dimension(nk) :: h0, S0, T0, h1, S1, T1
   real :: P_t, P_b, z_out, e_top
   integer :: k
 
   ! Calculate original interface positions
   e(nk+1) = -depth
   do k=nk,1,-1
     e(K) = e(K+1) + h(k)
+    h0(k) = h(nk+1-k) ! Keep a copy to use in remapping
   enddo
 
   P_t = 0.
@@ -958,6 +973,22 @@ subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, T, T_
     enddo
   endif
 
+  ! Now we need to remap but remapping assumes the surface is at the
+  ! same place in the two columns so we turn the column upside down.
+  if (associated(remap_CS)) then
+    do k=1,nk
+      S0(k) = S(nk+1-k)
+      T0(k) = T(nk+1-k)
+      h1(k) = h(nk+1-k)
+    enddo
+    call remapping_core_h(nk, h0, T0, nk, h1, T1, remap_CS )
+    call remapping_core_h(nk, h0, S0, nk, h1, S1, remap_CS )
+    do k=1,nk
+      S(k) = S1(nk+1-k)
+      T(k) = T1(nk+1-k)
+    enddo
+  endif
+
 end subroutine cut_off_column_top
 
 ! -----------------------------------------------------------------------------
@@ -2486,6 +2517,7 @@ subroutine MOM_state_init_tests(G, GV, tv)
   real, dimension(nk+1) :: e
   integer :: k
   real :: P_tot, P_t, P_b, z_out
+  type(remapping_CS), pointer :: remap_CS => NULL()
 
   do k = 1, nk
     h(k) = 100.
@@ -2521,7 +2553,7 @@ subroutine MOM_state_init_tests(G, GV, tv)
   write(0,*) ''
   write(0,*) h
   call cut_off_column_top(nk, tv, GV%Rho0, G%G_earth, -e(nk+1), GV%Angstrom, &
-               T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h)
+               T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS)
   write(0,*) h
 
 end subroutine MOM_state_init_tests