Merge pull request mom-ocean#730 from adcroft/fix-z-init-many-layers

Fix initialization from data with more levels than the model

src/ALE/coord_hycom.F90

@@ -13,7 +13,7 @@ module coord_hycom
 type, public :: hycom_CS
   private
 
-  !> Number of layers/levels
+  !> Number of layers/levels in generated grid
   integer :: nk
 
   !> Nominal near-surface resolution
@@ -40,8 +40,8 @@ module coord_hycom
 subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS)
   type(hycom_CS),       pointer    :: CS !< Unassociated pointer to hold the control structure
   integer,              intent(in) :: nk !< Number of layers in generated grid
-  real, dimension(:),   intent(in) :: coordinateResolution !< Z-space thicknesses (m)
-  real, dimension(:),   intent(in) :: target_density !< Interface target densities (kg/m3)
+  real, dimension(nk),  intent(in) :: coordinateResolution !< Z-space thicknesses (m)
+  real, dimension(nk+1),intent(in) :: target_density !< Interface target densities (kg/m3)
   type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation
 
   if (associated(CS)) call MOM_error(FATAL, "init_coord_hycom: CS already associated!")
@@ -99,11 +99,12 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
   type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   integer,               intent(in)    :: nz !< Number of levels
   real,                  intent(in)    :: depth !< Depth of ocean bottom (positive in H)
-  real, dimension(nz),   intent(in)    :: T, S !< T and S for column
+  real, dimension(nz),   intent(in)    :: T !< Temperature of column (degC)
+  real, dimension(nz),   intent(in)    :: S !< Salinity of column (psu)
   real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses, (in m or H)
   real, dimension(nz),   intent(in)    :: p_col !< Layer pressure in Pa
-  real, dimension(nz+1), intent(in)    :: z_col ! Interface positions relative to the surface in H units (m or kg m-2)
-  real, dimension(nz+1), intent(inout) :: z_col_new !< Absolute positions of interfaces
+  real, dimension(nz+1), intent(in)    :: z_col !< Interface positions relative to the surface in H units (m or kg m-2)
+  real, dimension(CS%nk+1), intent(inout) :: z_col_new !< Absolute positions of interfaces
   real, optional,        intent(in)    :: zScale !< Scaling factor from the input thicknesses in m
                                                  !! to desired units for zInterface, perhaps m_to_H.
   real,        optional, intent(in)    :: h_neglect !< A negligibly small width for the
@@ -115,7 +116,8 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
 
   ! Local variables
   integer   :: k
-  real, dimension(nz) :: rho_col, h_col_new ! Layer quantities
+  real, dimension(nz) :: rho_col ! Layer quantities
+  real, dimension(CS%nk) :: h_col_new ! New layer thicknesses
   real :: z_scale
   real :: stretching ! z* stretching, converts z* to z.
   real :: nominal_z ! Nominal depth of interface is using z* (m or Pa)
@@ -139,26 +141,26 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
   ! Interpolates for the target interface position with the rho_col profile
   ! Based on global density profile, interpolate to generate a new grid
   call build_and_interpolate_grid(CS%interp_CS, rho_col, nz, h(:), z_col, &
-           CS%target_density, nz, h_col_new, z_col_new, h_neglect, h_neglect_edge)
+           CS%target_density, CS%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge)
 
   ! Sweep down the interfaces and make sure that the interface is at least
   ! as deep as a nominal target z* grid
   nominal_z = 0.
   stretching = z_col(nz+1) / depth ! Stretches z* to z
-  do k = 2, nz+1
+  do k = 2, CS%nk+1
     nominal_z = nominal_z + (z_scale * CS%coordinateResolution(k-1)) * stretching
     z_col_new(k) = max( z_col_new(k), nominal_z )
     z_col_new(k) = min( z_col_new(k), z_col(nz+1) )
   enddo
 
-  if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz
+  if (maximum_depths_set .and. maximum_h_set) then ; do k=2,CS%nk
     ! The loop bounds are 2 & nz so the top and bottom interfaces do not move.
     ! Recall that z_col_new is positive downward.
     z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K), &
                        z_col_new(K-1) + CS%max_layer_thickness(k-1))
-  enddo ; elseif (maximum_depths_set) then ; do K=2,nz
+  enddo ; elseif (maximum_depths_set) then ; do K=2,CS%nk
     z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K))
-  enddo ; elseif (maximum_h_set) then ; do k=2,nz
+  enddo ; elseif (maximum_h_set) then ; do k=2,CS%nk
     z_col_new(K) = min(z_col_new(K), z_col_new(K-1) + CS%max_layer_thickness(k-1))
   enddo ; endif
 end subroutine build_hycom1_column

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

src/ALE/coord_zlike.F90

@@ -10,7 +10,7 @@ module coord_zlike
 !> Control structure containing required parameters for a z-like coordinate
 type, public :: zlike_CS ; private
 
-  !> Number of levels
+  !> Number of levels to be generated
   integer :: nk
 
   !> Minimum thickness allowed for layers, in the same thickness units that will
@@ -39,36 +39,37 @@ subroutine init_coord_zlike(CS, nk, coordinateResolution)
   CS%coordinateResolution = coordinateResolution
 end subroutine init_coord_zlike
 
+!> Deallocates the zlike control structure
 subroutine end_coord_zlike(CS)
-  type(zlike_CS), pointer :: CS
+  type(zlike_CS), pointer :: CS !< Coordinate control structure
 
-  ! nothing to do
+  ! Nothing to do
   if (.not. associated(CS)) return
   deallocate(CS%coordinateResolution)
   deallocate(CS)
 end subroutine end_coord_zlike
 
+!> Set parameters in the zlike structure
 subroutine set_zlike_params(CS, min_thickness)
-  type(zlike_CS), pointer    :: CS
-  real, optional, intent(in) :: min_thickness
+  type(zlike_CS), pointer    :: CS !< Coordinate control structure
+  real, optional, intent(in) :: min_thickness !< Minimum allowed thickness
 
   if (.not. associated(CS)) call MOM_error(FATAL, "set_zlike_params: CS not associated")
 
   if (present(min_thickness)) CS%min_thickness = min_thickness
 end subroutine set_zlike_params
 
 !> Builds a z* coordinate with a minimum thickness
-subroutine build_zstar_column(CS, nz, depth, total_thickness, zInterface, &
+subroutine build_zstar_column(CS, depth, total_thickness, zInterface, &
                               z_rigid_top, eta_orig, zScale)
-  type(zlike_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 or H)
-  real,                  intent(in)    :: total_thickness !< Column thickness (positive in the same units as depth)
-  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 the same units as depth)
-  real, optional,        intent(in)    :: eta_orig !< The actual original height of the top in the same units as depth
-  real, optional,        intent(in)    :: zScale !< Scaling factor from the target coordinate resolution
-                                                 !! in m to desired units for zInterface, perhaps m_to_H
+  type(zlike_CS),           intent(in)    :: CS !< Coordinate control structure
+  real,                     intent(in)    :: depth !< Depth of ocean bottom (positive in m or H)
+  real,                     intent(in)    :: total_thickness !< Column thickness (positive in the same units as depth)
+  real, dimension(CS%nk+1), intent(inout) :: zInterface !< Absolute positions of interfaces
+  real, optional,           intent(in)    :: z_rigid_top !< The height of a rigid top (negative in the same units as depth)
+  real, optional,           intent(in)    :: eta_orig !< The actual original height of the top in the same units as depth
+  real, optional,           intent(in)    :: zScale !< Scaling factor from the target coordinate resolution
+                                                    !! in m to desired units for zInterface, perhaps m_to_H
   ! Local variables
   real :: eta, stretching, dh, min_thickness, z0_top, z_star, z_scale
   integer :: k
@@ -77,7 +78,7 @@ subroutine build_zstar_column(CS, nz, depth, total_thickness, zInterface, &
   z_scale = 1.0 ; if (present(zScale)) z_scale = zScale
 
   new_zstar_def = .false.
-  min_thickness = min( CS%min_thickness, total_thickness/real(nz) )
+  min_thickness = min( CS%min_thickness, total_thickness/real(CS%nk) )
   z0_top = 0.
   if (present(z_rigid_top)) then
     z0_top = z_rigid_top
@@ -101,31 +102,31 @@ subroutine build_zstar_column(CS, nz, depth, total_thickness, zInterface, &
     ! 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
+    do k = 2,CS%nk
       z_star = z_star - CS%coordinateResolution(k-1)*z_scale
       ! 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 )
+      zInterface(k) = max( zInterface(k), -depth + real(CS%nk+1-k) * min_thickness )
     enddo
-    zInterface(nz+1) = -depth
+    zInterface(CS%nk+1) = -depth
 
   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
+    do k = 1,CS%nk
       dh = stretching * CS%coordinateResolution(k)*z_scale ! 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
+    zInterface(CS%nk+1) = -depth
+    do k = CS%nk,1,-1
       if ( zInterface(k) < (zInterface(k+1) + min_thickness) ) then
         zInterface(k) = zInterface(k+1) + min_thickness
       endif

src/framework/MOM_diag_remap.F90

@@ -264,11 +264,11 @@ subroutine diag_remap_update(remap_cs, G, GV, h, T, S, eqn_of_state)
     endif
 
     if (remap_cs%vertical_coord == coordinateMode('ZSTAR')) then
-      call build_zstar_column(get_zlike_CS(remap_cs%regrid_cs), nz, &
+      call build_zstar_column(get_zlike_CS(remap_cs%regrid_cs), &
                               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, &

src/framework/MOM_horizontal_regridding.F90

@@ -132,14 +132,14 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug
   real, parameter :: relc_default = 0.25, crit_default = 1.e-3
 
   integer :: npass
-  integer :: is, ie, js, je, nz
+  integer :: is, ie, js, je
   real    :: relax_coeff, acrit, ares
   logical :: debug_it
 
   debug_it=.false.
   if (PRESENT(debug)) debug_it=debug
 
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
   npass = num_pass_default
   if (PRESENT(num_pass)) npass = num_pass
@@ -298,7 +298,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
   integer :: isc,iec,jsc,jec    ! global compute domain indices
   integer :: isg, ieg, jsg, jeg ! global extent
   integer :: isd, ied, jsd, jed ! data domain indices
-  integer :: ni, nj, nz         ! global grid size
   integer :: id_clock_read
   character(len=12)  :: dim_name(4)
   logical :: debug=.false.
@@ -309,7 +308,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
   real, dimension(SZI_(G),SZJ_(G))  :: good2,fill2
   real, dimension(SZI_(G),SZJ_(G))  :: nlevs
 
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
   isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg
 
@@ -616,7 +615,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
   integer :: isc,iec,jsc,jec    ! global compute domain indices
   integer :: isg, ieg, jsg, jeg ! global extent
   integer :: isd, ied, jsd, jed ! data domain indices
-  integer :: ni, nj, nz         ! global grid size
   integer :: id_clock_read
   integer, dimension(4) :: fld_sz
   character(len=12)  :: dim_name(4)
@@ -628,7 +626,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
   real, dimension(SZI_(G),SZJ_(G))  :: good2,fill2
   real, dimension(SZI_(G),SZJ_(G))  :: nlevs
 
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
   isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg