ExtrudedHybridLevels Subroutine

  subroutine ExtrudedHybridLevels(h,depths,sigma,zLevels)
  !! Fills the level heights for ExtrudedMesh with nzFlat = size(depths) flat
  !! interfaces near the surface over nzS = size(sigma)-1 terrain-following layers.
  !! With d_0 = 0 and d_F = depths(nzFlat) (d_F = 0 when nzFlat = 0), the nz = nzS+nzFlat
  !! levels at every geometry node are
  !!
  !!   z_k        = -( (1+sigma_k) d_F - sigma_k h ) ,  k = 0..nzS
  !!   z_{nzS+m}  = -d_{nzFlat-m} ,                    m = 1..nzFlat
  !!
  !! so z_0 = -h, z_nzS = -d_F, z_nz = 0, and the layers between z = -d_F and the
  !! bottom follow the terrain with the stretching sigma. The form used for z_k makes
  !! z_0 = -h and z_nzS = -d_F exact in floating point. With nzFlat = 0 this is the
  !! sigma coordinate z_k = sigma_k h (see ExtrudedSigmaLevels); with a constant h it
  !! gives flat z-levels throughout.
  !!
  !! The levels are evaluated pointwise from h, so they are continuous across lateral
  !! faces exactly when h is: two elements that share a geometry node must be given the
  !! same h there, as they are when h is sampled from one bathymetry function at the
  !! mesh's node coordinates.
  !!
  !!  Input
  !!    - h(1:nGeo+1,1:nGeo+1,1:nElem2D) : Water depth at the surface geometry nodes, in
  !!      the length unit of the mesh, with h > d_F (> 0) everywhere
  !!    - depths(1:nzFlat) : Depths of the flat interfaces below the surface, same unit,
  !!      0 < d_1 < ... < d_nzFlat < min(h). May be empty.
  !!    - sigma(0:nzS) : Terrain-following stretching, nzS >= 1, sigma(0) = -1,
  !!      sigma(nzS) = 0, strictly increasing
  !!
  !!  Output
  !!    - zLevels(1:nGeo+1,1:nGeo+1,0:nz,1:nElem2D) : Level heights, nz = nzS+nzFlat
    implicit none
    real(prec),intent(in) :: h(1:,1:,1:)
    real(prec),intent(in) :: depths(1:)
    real(prec),intent(in) :: sigma(0:)
    real(prec),intent(out) :: zLevels(1:,1:,0:,1:)
    ! Local
    integer :: nzFlat,nzS,nz,i,j,k,m,e
    real(prec) :: dF

    nzFlat = size(depths)
    nzS = size(sigma)-1
    nz = nzS+nzFlat

    if(nzS < 1) then
      error stop "ExtrudedHybridLevels : sigma must carry at least two values (one terrain-following layer)"
    endif
    if(sigma(0) /= -1.0_prec .or. sigma(nzS) /= 0.0_prec) then
      error stop "ExtrudedHybridLevels : sigma must run from sigma(0) = -1 to sigma(nzS) = 0"
    endif
    do k = 1,nzS
      if(.not.(sigma(k) > sigma(k-1))) then
        error stop "ExtrudedHybridLevels : sigma must be strictly increasing"
      endif
    enddo
    if(nzFlat > 0) then
      if(.not.(depths(1) > 0.0_prec)) then
        error stop "ExtrudedHybridLevels : the flat interface depths must be positive"
      endif
      do m = 2,nzFlat
        if(.not.(depths(m) > depths(m-1))) then
          error stop "ExtrudedHybridLevels : the flat interface depths must be strictly increasing"
        endif
      enddo
      dF = depths(nzFlat)
    else
      dF = 0.0_prec
    endif
    if(size(zLevels,1) /= size(h,1) .or. size(zLevels,2) /= size(h,2) .or. &
       size(zLevels,3) /= nz+1 .or. size(zLevels,4) /= size(h,3)) then
      error stop "ExtrudedHybridLevels : zLevels must be dimensioned (1:nGeo+1,1:nGeo+1,0:nz,1:nElem2D)"
    endif
    if(size(h) > 0) then
      if(.not.(minval(h) > dF)) then
        error stop "ExtrudedHybridLevels : the depth h must exceed the deepest flat interface (and 0)"
      endif
    endif

    do e = 1,size(h,3)
      do j = 1,size(h,2)
        do i = 1,size(h,1)
          do k = 0,nzS
            zLevels(i,j,k,e) = -((1.0_prec+sigma(k))*dF-sigma(k)*h(i,j,e))
          enddo
          do m = 1,nzFlat-1
            zLevels(i,j,nzS+m,e) = -depths(nzFlat-m)
          enddo
          zLevels(i,j,nz,e) = 0.0_prec
        enddo
      enddo
    enddo

  endsubroutine ExtrudedHybridLevels