calc_thk_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  c a l c _ t h k _ m
!
!> @file
!!
!! Computation of the ice thickness.
!!
!! @section Copyright
!!
!! Copyright 2009-2017 Ralf Greve, Reinhard Calov, Tatsuru Sato
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Computation of the ice thickness.
!<------------------------------------------------------------------------------
module calc_thk_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m

  implicit none

  real(dp), dimension(0:JMAX,0:IMAX), save :: h, h_neu
  real(dp), dimension(0:JMAX,0:IMAX), save :: mb_source, mb_adjust
  logical,                            save :: flag_solver_explicit

  private
  public :: calc_thk_init
  public :: calc_thk_sia_expl, calc_thk_sia_impl
  public :: calc_thk_expl, calc_thk_impl
  public :: calc_thk_mask_update

contains

!-------------------------------------------------------------------------------
!> Initialisations for the ice thickness computation.
!<------------------------------------------------------------------------------
subroutine calc_thk_init()

implicit none

#if (defined(CALCZS))
  stop ' >>> calc_thk_init: Replace CALCZS by CALCTHK in header file!'
#elif (!defined(CALCTHK))
  stop ' >>> calc_thk_init: Define CALCTHK in header file!'
#endif

!-------- Computation/initialisation of the ice base topography
!                                       and its time derivative --------

#if (MARGIN==1 || MARGIN==2)   /* only grounded ice */

zb_neu   = zl_neu
dzb_dtau = dzl_dtau

#elif (MARGIN==3)   /* grounded and floating ice */

where (maske <= 1_i2b)   ! grounded ice or ice-free land
   zb_neu   = zl_neu
   dzb_dtau = dzl_dtau
elsewhere   ! (maske >= 2_i2b; ocean or floating ice)
   zb_neu   = zb       ! initialisation,
   dzb_dtau = 0.0_dp   ! will be overwritten later
end where

#else
stop ' >>> calc_thk_init: MARGIN must be either 1, 2 or 3!'
#endif

!-------- Initialisation of the ice thickness
!                               and surface topography --------

h = h_c + h_t

zs_neu = zs   ! initialisation,
h_neu  = h    ! will be overwritten later

!-------- Solver type --------

#if (CALCTHK==1 || CALCTHK==4)   /* explicit solver */

  flag_solver_explicit = .true.

#elif (CALCTHK==2 || CALCTHK==3 \
       || calcthk==5 || calcthk==6)   /* implicit solver */

  flag_solver_explicit = .false.

#else
stop ' >>> calc_thk_init: CALCTHK must be between 1 and 6!'
#endif

!-------- Source term for the ice thickness equation --------

#if (MB_ACCOUNT==0)

  mb_source = as_perp - q_b_tot - calv_grounded

#elif (MB_ACCOUNT==1)

  if (flag_solver_explicit) then
     mb_source = 0.0_dp   ! source term will be applied separately
  else
     mb_source = as_perp - q_b_tot - calv_grounded
  end if

#else
  stop ' >>> calc_thk_init: MB_ACCOUNT must be either 0 or 1!'
#endif

end subroutine calc_thk_init

!-------------------------------------------------------------------------------
!> Explicit solver for the diffusive SIA ice surface equation.
!<------------------------------------------------------------------------------
subroutine calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)

implicit none

real(dp), intent(in) :: time, dtime, dxi, deta
real(dp), intent(in) :: z_mar

integer(i4b)                       :: i, j
real(dp)                           :: azs2, azs3
real(dp), dimension(0:JMAX,0:IMAX) :: czs2, czs3

!-------- Abbreviations --------

azs2 = dtime/(dxi*dxi)
azs3 = dtime/(deta*deta)

czs2 = 0.0_dp
czs3 = 0.0_dp

do i=0, imax-1
do j=0, jmax
   czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
               *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
end do
end do

do i=0, imax
do j=0, jmax-1
   czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
               *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
end do
end do

!-------- Solution of the explicit scheme --------

do i=0, imax
do j=0, jmax

   if ( flag_inner_point(j,i) &   ! inner point
#if (MARGIN==1 && MB_ACCOUNT==0)
        .and.(maske(j,i) <= 1_i2b) &   ! grounded ice or ice-free land
#endif
      ) then

      zs_neu(j,i) = zs(j,i) &
                    + dtime*(mb_source(j,i)+dzb_dtau(j,i)) &
                    + ( czs2(j,i)  *(zs(j,i+1)-zs(j,i)  ) &
                       -czs2(j,i-1)*(zs(j,i)  -zs(j,i-1)) &
                       +czs3(j,i)  *(zs(j+1,i)-zs(j,i)  ) &
                       -czs3(j-1,i)*(zs(j,i)  -zs(j-1,i)) ) &
                      *insq_g11_g(j,i)*insq_g22_g(j,i)

   else
      zs_neu(j,i) = zb_neu(j,i)   ! zero-thickness boundary condition
   end if

end do
end do

!-------- Ice thickness --------

h_neu = zs_neu - zb_neu

!-------- Applying the so-far neglected source term --------

#if (MB_ACCOUNT==1)
call apply_mb_source(dtime, z_mar)
#endif

!-------- Adjusting the ice thickness, if needed --------

call thk_adjust(time, dtime)

end subroutine calc_thk_sia_expl

!-------------------------------------------------------------------------------
!> Over-implicit solver for the diffusive SIA ice surface equation.
!<------------------------------------------------------------------------------
subroutine calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, &
                             ii, jj, nn)

use sico_maths_m, only : sor_sprs

implicit none

integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
integer(i4b), intent(in) :: nn(0:jmax,0:imax)
real(dp),     intent(in) :: time, dtime, dxi, deta
real(dp),     intent(in) :: z_mar
real(dp),     intent(in) :: mean_accum

integer(i4b)                       :: i, j
integer(i4b)                       :: k, nnz
real(dp)                           :: azs2, azs3
real(dp), dimension(0:JMAX,0:IMAX) :: czs2, czs3

#if (CALCTHK==2)

integer(i4b)                            :: ierr
integer(i4b)                            :: iter
integer(i4b)                            :: nc, nr
integer(i4b), parameter                 :: nmax   =    (imax+1)*(jmax+1)
integer(i4b), parameter                 :: n_sprs = 10*(imax+1)*(jmax+1)
integer(i4b), allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
integer(i4b), allocatable, dimension(:) :: lgs_a_diag_index
integer(i4b), allocatable, dimension(:) :: lgs_a_index_trim
real(dp),     allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
real(dp),     allocatable, dimension(:) :: lgs_a_value_trim
real(dp)                                :: eps_sor

#elif (CALCTHK==3)

lis_integer                            :: ierr
lis_integer                            :: iter
lis_integer                            :: nc, nr
lis_integer, parameter                 :: nmax   =    (imax+1)*(jmax+1)
lis_integer, parameter                 :: n_sprs = 10*(imax+1)*(jmax+1)
lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
lis_integer, allocatable, dimension(:) :: lgs_a_diag_index
lis_matrix                             :: lgs_a
lis_vector                             :: lgs_b, lgs_x
lis_scalar,  allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
lis_solver                             :: solver
character(len=256)                     :: ch_solver_set_option

#endif

!-------- Abbreviations --------

azs2 = dtime/(dxi*dxi)
azs3 = dtime/(deta*deta)

czs2 = 0.0_dp
czs3 = 0.0_dp

do i=0, imax-1
do j=0, jmax
   czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
               *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
end do
end do

do i=0, imax
do j=0, jmax-1
   czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
               *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
end do
end do

!-------- Assembly of the system of linear equations
!                     (matrix storage: compressed sparse row CSR) --------

#if (CALCTHK==2 || CALCTHK==3)

allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))

lgs_a_value = 0.0_dp
lgs_a_index = 0
lgs_a_ptr   = 0

lgs_b_value = 0.0_dp
lgs_x_value = 0.0_dp

lgs_a_ptr(1) = 1

k = 0

do nr=1, nmax   ! loop over rows

   i = ii(nr)
   j = jj(nr)

   if ( flag_inner_point(j,i) &   ! inner point
#if (MARGIN==1 && MB_ACCOUNT==0)
        .and.(maske(j,i) <= 1_i2b) &   ! grounded ice or ice-free land
#endif
      ) then

      nc = nn(j-1,i)   ! smallest nc (column counter), for zs(j-1,i)
      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
      lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
                        *insq_g11_g(j,i)*insq_g22_g(j,i)
      lgs_a_index(k) = nc

      nc = nn(j,i-1)   ! next nc (column counter), for zs(j,i-1)
      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
      lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
                        *insq_g11_g(j,i)*insq_g22_g(j,i)
      lgs_a_index(k) = nc

      nc = nn(j,i)     ! next nc (column counter), for zs(j,i)
      ! if (nc /= nr) &                     (diagonal element)
      !    stop ' >>> calc_thk_sia_impl: Check for diagonal element failed!'
      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
      lgs_a_value(k) = 1.0_dp &
                       + (czs2(j,i)+czs2(j,i-1)+czs3(j,i)+czs3(j-1,i)) &
                         *ovi_weight &
                         *insq_g11_g(j,i)*insq_g22_g(j,i)
      lgs_a_diag_index(nr) = k
      lgs_a_index(k) = nc

      nc = nn(j,i+1)   ! next nc (column counter), for zs(j,i+1)
      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
      lgs_a_value(k) = -czs2(j,i)*ovi_weight &
                        *insq_g11_g(j,i)*insq_g22_g(j,i)
      lgs_a_index(k) = nc

      nc = nn(j+1,i)   ! largest nc (column counter), for zs(j+1,i)
      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
      lgs_a_value(k) = -czs3(j,i)*ovi_weight &
                        *insq_g11_g(j,i)*insq_g22_g(j,i)
      lgs_a_index(k) = nc

      lgs_b_value(nr) = zs(j,i) &
                          + dtime*(mb_source(j,i)+dzb_dtau(j,i)) &
                          + ( czs2(j,i)*(zs(j,i+1)-zs(j,i)) &
                             -czs2(j,i-1)*(zs(j,i)-zs(j,i-1)) &
                             +czs3(j,i)*(zs(j+1,i)-zs(j,i)) &
                             -czs3(j-1,i)*(zs(j,i)-zs(j-1,i)) ) &
                            *(1.0_dp-ovi_weight) &
                            *insq_g11_g(j,i)*insq_g22_g(j,i)
                                                          ! right-hand side

   else   !  zero-thickness boundary condition

      k = k+1
      ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
      lgs_a_value(k)       = 1.0_dp   ! diagonal element only
      lgs_a_diag_index(nr) = k
      lgs_a_index(k)       = nr
      lgs_b_value(nr)      = zb_neu(j,i)

   end if

   lgs_x_value(nr) = zs(j,i)   ! old surface topography,
                               ! initial guess for solution vector

   lgs_a_ptr(nr+1) = k+1   ! row is completed, store index to next row

end do

nnz = k   ! number of non-zero elements of the matrix

#else
stop ' >>> calc_thk_sia_impl: CALCTHK must be either 2 or 3!'
#endif

!-------- Solution of the system of linear equations --------

#if (CALCTHK==2)

!  ------ Solution with the built-in SOR solver

allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))

do k=1, nnz   ! relocate matrix to trimmed arrays
   lgs_a_value_trim(k) = lgs_a_value(k)
   lgs_a_index_trim(k) = lgs_a_index(k)
end do

deallocate(lgs_a_value, lgs_a_index)

eps_sor = 1.0e-05_dp*mean_accum*dtime   ! convergence parameter

call sor_sprs(lgs_a_value_trim, &
              lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
              lgs_b_value, &
              nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)

do nr=1, nmax
   i = ii(nr)
   j = jj(nr)
   zs_neu(j,i) = lgs_x_value(nr)
end do

deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)

#elif (CALCTHK==3)

!  ------ Settings for Lis

call lis_matrix_create(lis_comm_world, lgs_a, ierr)
call lis_vector_create(lis_comm_world, lgs_b, ierr)
call lis_vector_create(lis_comm_world, lgs_x, ierr)

call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
call lis_vector_set_size(lgs_b, 0, nmax, ierr)
call lis_vector_set_size(lgs_x, 0, nmax, ierr)

do nr=1, nmax

   do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
      call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
                                               lgs_a_value(nc), lgs_a, ierr)
   end do

   call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
   call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)

end do

call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
call lis_matrix_assemble(lgs_a, ierr)

!  ------ Solution with Lis

call lis_solver_create(solver, ierr)

ch_solver_set_option = '-i bicg -p ilu '// &
                       '-maxiter 1000 -tol 1.0e-04 -initx_zeros false'

call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
call chkerr(ierr)

call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
call chkerr(ierr)

call lis_solver_get_iter(solver, iter, ierr)
write(6,'(10x,a,i0)') 'calc_thk_sia_impl: iter = ', iter

lgs_x_value = 0.0_dp
call lis_vector_gather(lgs_x, lgs_x_value, ierr)
call lis_matrix_destroy(lgs_a, ierr)
call lis_vector_destroy(lgs_b, ierr)
call lis_vector_destroy(lgs_x, ierr)
call lis_solver_destroy(solver, ierr)

do nr=1, nmax
   i = ii(nr)
   j = jj(nr)
   zs_neu(j,i) = lgs_x_value(nr)
end do

deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)

#endif

!-------- Ice thickness --------

h_neu = zs_neu - zb_neu

!-------- Applying the so-far neglected source term --------

#if (MB_ACCOUNT==1)
call apply_mb_source(dtime, z_mar)
#endif

!-------- Adjusting the ice thickness, if needed --------

call thk_adjust(time, dtime)

end subroutine calc_thk_sia_impl

!-------------------------------------------------------------------------------
!> Explicit solver for the general ice thickness equation.
!<------------------------------------------------------------------------------
subroutine calc_thk_expl(time, dtime, dxi, deta, z_mar)

implicit none

real(dp), intent(in) :: time, dtime, dxi, deta
real(dp), intent(in) :: z_mar

integer(i4b)                       :: i, j
real(dp), dimension(0:JMAX,0:IMAX) :: dt_darea
real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
real(dp), dimension(0:JMAX,0:IMAX) :: upH_x_1, upH_x_2, upH_y_1, upH_y_2
real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_x_1, sq_g22_x_2
real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_y_1, sq_g11_y_2

!-------- Abbreviations --------

dt_darea = dtime/area

do i=0, imax
do j=0, jmax

   if (flag_inner_point(j,i)) then

      vx_m_1(j,i) = vx_m(j,i-1)
      vx_m_2(j,i) = vx_m(j,i)
      vy_m_1(j,i) = vy_m(j-1,i)
      vy_m_2(j,i) = vy_m(j,i)

      if (vx_m_1(j,i) >= 0.0_dp) then
         uph_x_1(j,i) = h(j,i-1)
      else
         uph_x_1(j,i) = h(j,i)
      end if

      if (vx_m_2(j,i) >= 0.0_dp) then
         uph_x_2(j,i) = h(j,i)
      else
         uph_x_2(j,i) = h(j,i+1)
      end if

      if (vy_m_1(j,i) >= 0.0_dp) then
         uph_y_1(j,i) = h(j-1,i)
      else
         uph_y_1(j,i) = h(j,i)
      end if

      if (vy_m_2(j,i) >= 0.0_dp) then
         uph_y_2(j,i) = h(j,i)
      else
         uph_y_2(j,i) = h(j+1,i)
      end if

      sq_g22_x_1(j,i) = sq_g22_sgx(j,i-1)
      sq_g22_x_2(j,i) = sq_g22_sgx(j,i)
      sq_g11_y_1(j,i) = sq_g11_sgy(j-1,i)
      sq_g11_y_2(j,i) = sq_g11_sgy(j,i)

   else   ! .not.(flag_inner_point(j,i))

      vx_m_1(j,i) = 0.0_dp
      vx_m_2(j,i) = 0.0_dp
      vy_m_1(j,i) = 0.0_dp
      vy_m_2(j,i) = 0.0_dp

      uph_x_1(j,i) = 0.0_dp
      uph_x_2(j,i) = 0.0_dp
      uph_y_1(j,i) = 0.0_dp
      uph_y_2(j,i) = 0.0_dp

      sq_g22_x_1(j,i) = 0.0_dp
      sq_g22_x_2(j,i) = 0.0_dp
      sq_g11_y_1(j,i) = 0.0_dp
      sq_g11_y_2(j,i) = 0.0_dp

   end if

end do
end do

!-------- Solution of the explicit scheme --------

where ( flag_inner_point &   ! inner point
#if (MARGIN==1 && MB_ACCOUNT==0)
        .and.(maske <= 1_i2b) &   ! grounded ice or ice-free land
#endif
      )

   h_neu = h + dtime*mb_source &
             - dt_darea &
               * (  ( vx_m_2*uph_x_2*sq_g22_x_2*deta   &
                     -vx_m_1*uph_x_1*sq_g22_x_1*deta ) &
                  + ( vy_m_2*uph_y_2*sq_g11_y_2*dxi    &
                     -vy_m_1*uph_y_1*sq_g11_y_1*dxi  ) )

elsewhere
   h_neu = 0.0_dp   ! zero-thickness boundary condition
end where

!-------- Applying the so-far neglected source term --------

#if (MB_ACCOUNT==1)
call apply_mb_source(dtime, z_mar)
#endif

!-------- Adjusting the ice thickness, if needed --------

call thk_adjust(time, dtime)

end subroutine calc_thk_expl

!-------------------------------------------------------------------------------
!> Over-implicit solver for the general ice thickness equation.
!<------------------------------------------------------------------------------
subroutine calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, &
                         ii, jj, nn)

use sico_maths_m, only : sor_sprs

implicit none

integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
integer(i4b), intent(in) :: nn(0:jmax,0:imax)
real(dp),     intent(in) :: time, dtime, dxi, deta
real(dp),     intent(in) :: z_mar
real(dp),     intent(in) :: mean_accum

integer(i4b)                       :: i, j
integer(i4b)                       :: k, nnz
real(dp), dimension(0:JMAX,0:IMAX) :: dt_darea
real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
real(dp), dimension(0:JMAX,0:IMAX) :: upH_x_1, upH_x_2, upH_y_1, upH_y_2
real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_x_1, sq_g22_x_2
real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_y_1, sq_g11_y_2

#if (CALCTHK==5)

integer(i4b)                            :: ierr
integer(i4b)                            :: iter
integer(i4b)                            :: nc, nr
integer(i4b), parameter                 :: nmax   =    (imax+1)*(jmax+1)
integer(i4b), parameter                 :: n_sprs = 10*(imax+1)*(jmax+1)
integer(i4b), allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
integer(i4b), allocatable, dimension(:) :: lgs_a_diag_index
integer(i4b), allocatable, dimension(:) :: lgs_a_index_trim
real(dp),     allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
real(dp),     allocatable, dimension(:) :: lgs_a_value_trim
real(dp)                                :: eps_sor

#elif (CALCTHK==6)

lis_integer                            :: ierr
lis_integer                            :: iter
lis_integer                            :: nc, nr
lis_integer, parameter                 :: nmax   =    (imax+1)*(jmax+1)
lis_integer, parameter                 :: n_sprs = 10*(imax+1)*(jmax+1)
lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
lis_integer, allocatable, dimension(:) :: lgs_a_diag_index
lis_matrix                             :: lgs_a
lis_vector                             :: lgs_b, lgs_x
lis_scalar,  allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
lis_solver                             :: solver
character(len=256)                     :: ch_solver_set_option

#endif

!-------- Abbreviations --------

dt_darea = dtime/area

do i=0, imax
do j=0, jmax

   if (flag_inner_point(j,i)) then

      vx_m_1(j,i) = vx_m(j,i-1)
      vx_m_2(j,i) = vx_m(j,i)
      vy_m_1(j,i) = vy_m(j-1,i)
      vy_m_2(j,i) = vy_m(j,i)

      if (vx_m_1(j,i) >= 0.0_dp) then
         uph_x_1(j,i) = h(j,i-1)
      else
         uph_x_1(j,i) = h(j,i)
      end if

      if (vx_m_2(j,i) >= 0.0_dp) then
         uph_x_2(j,i) = h(j,i)
      else
         uph_x_2(j,i) = h(j,i+1)
      end if

      if (vy_m_1(j,i) >= 0.0_dp) then
         uph_y_1(j,i) = h(j-1,i)
      else
         uph_y_1(j,i) = h(j,i)
      end if

      if (vy_m_2(j,i) >= 0.0_dp) then
         uph_y_2(j,i) = h(j,i)
      else
         uph_y_2(j,i) = h(j+1,i)
      end if

      sq_g22_x_1(j,i) = sq_g22_sgx(j,i-1)
      sq_g22_x_2(j,i) = sq_g22_sgx(j,i)
      sq_g11_y_1(j,i) = sq_g11_sgy(j-1,i)
      sq_g11_y_2(j,i) = sq_g11_sgy(j,i)

   else   ! .not.(flag_inner_point(j,i))

      vx_m_1(j,i) = 0.0_dp
      vx_m_2(j,i) = 0.0_dp
      vy_m_1(j,i) = 0.0_dp
      vy_m_2(j,i) = 0.0_dp

      uph_x_1(j,i) = 0.0_dp
      uph_x_2(j,i) = 0.0_dp
      uph_y_1(j,i) = 0.0_dp
      uph_y_2(j,i) = 0.0_dp

      sq_g22_x_1(j,i) = 0.0_dp
      sq_g22_x_2(j,i) = 0.0_dp
      sq_g11_y_1(j,i) = 0.0_dp
      sq_g11_y_2(j,i) = 0.0_dp

   end if

end do
end do

!-------- Assembly of the system of linear equations
!                     (matrix storage: compressed sparse row CSR) --------

#if (CALCTHK==5 || CALCTHK==6)

allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))

lgs_a_value = 0.0_dp
lgs_a_index = 0
lgs_a_ptr   = 0

lgs_b_value = 0.0_dp
lgs_x_value = 0.0_dp

lgs_a_ptr(1) = 1

k = 0

do nr=1, nmax   ! loop over rows

   i = ii(nr)
   j = jj(nr)

   if ( flag_inner_point(j,i) &   ! inner point
#if (MARGIN==1 && MB_ACCOUNT==0)
        .and.(maske(j,i) <= 1_i2b) &   ! grounded ice or ice-free land
#endif
      ) then

      k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc   ! for H(j-1,i)
      if (vy_m_1(j,i) > 0.0_dp) &
         lgs_a_value(k) = -dt_darea(j,i)*vy_m_1(j,i) &
                                        *sq_g11_y_1(j,i)*dxi*ovi_weight

      k=k+1 ; nc=nn(j,i-1) ; lgs_a_index(k)=nc   ! for H(j,i-1)
      if (vx_m_1(j,i) > 0.0_dp) &
         lgs_a_value(k) = -dt_darea(j,i)*vx_m_1(j,i) &
                                        *sq_g22_x_1(j,i)*deta*ovi_weight

      k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k  ! for H(j,i)
      lgs_a_value(k) = 1.0_dp                             ! (diagonal element)
      if (vy_m_1(j,i) < 0.0_dp) &
         lgs_a_value(k) = lgs_a_value(k) &
                          - dt_darea(j,i)*vy_m_1(j,i) &
                                         *sq_g11_y_1(j,i)*dxi*ovi_weight
      if (vx_m_1(j,i) < 0.0_dp) &
         lgs_a_value(k) = lgs_a_value(k) &
                          - dt_darea(j,i)*vx_m_1(j,i) &
                                         *sq_g22_x_1(j,i)*deta*ovi_weight
      if (vx_m_2(j,i) > 0.0_dp) &
         lgs_a_value(k) = lgs_a_value(k) &
                          + dt_darea(j,i)*vx_m_2(j,i) &
                                         *sq_g22_x_2(j,i)*deta*ovi_weight
      if (vy_m_2(j,i) > 0.0_dp) &
         lgs_a_value(k) = lgs_a_value(k) &
                          + dt_darea(j,i)*vy_m_2(j,i) &
                                         *sq_g11_y_2(j,i)*dxi*ovi_weight

      k=k+1 ; nc=nn(j,i+1) ; lgs_a_index(k)=nc   ! for H(j,i+1)
      if (vx_m_2(j,i) < 0.0_dp) &
         lgs_a_value(k) = dt_darea(j,i)*vx_m_2(j,i) &
                                       *sq_g22_x_2(j,i)*deta*ovi_weight

      k=k+1 ; nc=nn(j+1,i) ; lgs_a_index(k)=nc   ! for H(j+1,i)
      if (vy_m_2(j,i) < 0.0_dp) &
         lgs_a_value(k) = dt_darea(j,i)*vy_m_2(j,i) &
                                       *sq_g11_y_2(j,i)*dxi*ovi_weight

      lgs_b_value(nr) = h(j,i) &
                        +dtime*mb_source(j,i) &
                        -(1.0_dp-ovi_weight) &
                           * dt_darea(j,i) &
                             * (  ( vx_m_2(j,i)*uph_x_2(j,i) &
                                               *sq_g22_x_2(j,i)*deta   &
                                   -vx_m_1(j,i)*uph_x_1(j,i) &
                                               *sq_g22_x_1(j,i)*deta ) &
                                + ( vy_m_2(j,i)*uph_y_2(j,i) &
                                               *sq_g11_y_2(j,i)*dxi    &
                                   -vy_m_1(j,i)*uph_y_1(j,i) &
                                               *sq_g11_y_1(j,i)*dxi  ) )
                                                          ! right-hand side

   else   ! zero-thickness boundary condition

      k = k+1
      lgs_a_value(k)       = 1.0_dp   ! diagonal element only
      lgs_a_diag_index(nr) = k
      lgs_a_index(k)       = nr
      lgs_b_value(nr)      = 0.0_dp

   end if

   lgs_x_value(nr) = h(j,i)   ! old ice thickness,
                              ! initial guess for solution vector

   lgs_a_ptr(nr+1) = k+1   ! row is completed, store index to next row

end do

nnz = k   ! number of non-zero elements of the matrix

#else
stop ' >>> calc_thk_impl: CALCTHK must be either 5 or 6!'
#endif

!-------- Solution of the system of linear equations --------

#if (CALCTHK==5)

!  ------ Solution with the built-in SOR solver

allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))

do k=1, nnz   ! relocate matrix to trimmed arrays
   lgs_a_value_trim(k) = lgs_a_value(k)
   lgs_a_index_trim(k) = lgs_a_index(k)
end do

deallocate(lgs_a_value, lgs_a_index)

eps_sor = 1.0e-05_dp*mean_accum*dtime   ! convergence parameter

call sor_sprs(lgs_a_value_trim, &
              lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
              lgs_b_value, &
              nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)

do nr=1, nmax
   i = ii(nr)
   j = jj(nr)
   h_neu(j,i) = lgs_x_value(nr)
end do

deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)

#elif (CALCTHK==6)

!  ------ Settings for Lis

call lis_matrix_create(lis_comm_world, lgs_a, ierr)
call lis_vector_create(lis_comm_world, lgs_b, ierr)
call lis_vector_create(lis_comm_world, lgs_x, ierr)

call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
call lis_vector_set_size(lgs_b, 0, nmax, ierr)
call lis_vector_set_size(lgs_x, 0, nmax, ierr)

do nr=1, nmax

   do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
      call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
                                               lgs_a_value(nc), lgs_a, ierr)
   end do

   call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
   call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)

end do

call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
call lis_matrix_assemble(lgs_a, ierr)

!  ------ Solution with Lis

call lis_solver_create(solver, ierr)

ch_solver_set_option = '-i bicg -p ilu '// &
                       '-maxiter 1000 -tol 1.0e-04 -initx_zeros false'

call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
call chkerr(ierr)

call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
call chkerr(ierr)

call lis_solver_get_iter(solver, iter, ierr)
write(6,'(10x,a,i0)') 'calc_thk_impl: iter = ', iter

lgs_x_value = 0.0_dp
call lis_vector_gather(lgs_x, lgs_x_value, ierr)
call lis_matrix_destroy(lgs_a, ierr)
call lis_vector_destroy(lgs_b, ierr)
call lis_vector_destroy(lgs_x, ierr)
call lis_solver_destroy(solver, ierr)

do nr=1, nmax
   i = ii(nr)
   j = jj(nr)
   h_neu(j,i) = lgs_x_value(nr)
end do

deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)

#endif

!-------- Applying the so-far neglected source term --------

#if (MB_ACCOUNT==1)
call apply_mb_source(dtime, z_mar)
#endif

!-------- Adjusting the ice thickness, if needed --------

call thk_adjust(time, dtime)

end subroutine calc_thk_impl

!-------------------------------------------------------------------------------
!> Ice thickness evolution due to the source term (surface mass balance,
!! basal mass balance and calving).
!! Determination of the several components of the mass balance.
!<------------------------------------------------------------------------------
subroutine apply_mb_source(dtime, z_mar)

  implicit none

  real(dp), intent(in) :: dtime
  real(dp), intent(in) :: z_mar

  real(dp), dimension(0:JMAX,0:IMAX) :: H_smb, H_adv
  logical,  dimension(0:JMAX,0:IMAX) :: flag_ocean_point

  if (flag_solver_explicit) then
     mb_source = as_perp - q_b_tot - calv_grounded
                 ! define source term for the ice thickness equation
  else
     where (flag_inner_point) h_neu = h_neu - dtime*mb_source
                 ! subtract source term from the ice thickness solution
  end if

! Store thickness as local variable for updating
  h_smb = h_neu
  
! ### At the end, forbid appearance of ice past two outer grid points
! ### (this relates to the one outer points in calc_top and, therefore, applies
! ### for non-cyclic RB for domain boundary of topography only)
  h_smb(0,:)      = 0.0_dp;   h_smb(1,:)    = 0.0_dp
  h_smb(jmax-1,:) = 0.0_dp;   h_smb(jmax,:) = 0.0_dp
  h_smb(:,0)      = 0.0_dp;   h_smb(:,1)    = 0.0_dp
  h_smb(:,imax-1) = 0.0_dp;   h_smb(:,imax) = 0.0_dp

! Now, for inner points *two* steps from the boundaries
! ### Apply mass balance source term ###
  h_smb(2:jmax-2,2:imax-2) &
     = max( h_neu(2:jmax-2,2:imax-2) + dtime*mb_source(2:jmax-2,2:imax-2), &
            0.0_dp)

! Note: Still the ice mass which went past the continental shelf and to the second 
! outer points will by counted later in this scheme by using H_adv and H_smb. This 
! is why the next loop runs from 1 to IMAX and from 1 to JMAX.

#if (MARGIN==1)
  ! Ocean kill: if ocean point and ice exists, remove the ice 
  where (maske == 2_i2b .and. h_smb > eps_dp) h_smb = 0.0_dp
#else   /* margin = 2 or 3  */
  ! If outside continental shelf, remove the ice 
  where (zl0 < z_mar) h_smb = 0.0_dp   !!! Is this OK? ... !!!
#endif
 
  ! Forbid negative thickness
  where (h_smb < 0.0_dp) h_smb = 0.0_dp
  
  ! Set mass balance and melt quatities zero
  mb_source_apl = 0.0_dp
  as_perp_apl   = 0.0_dp
  runoff_top    = 0.0_dp
  disc_top      = 0.0_dp
  runoff_bot    = 0.0_dp
  mask_run      = 0_i2b

  ! For inner points *one* step from the boundaries. All mass which went here by
  ! ice flow will be counted.

  h_adv = h_neu

  where (flag_inner_point)

    where (maske == 2_i2b)
       flag_ocean_point = .true.
    elsewhere
       flag_ocean_point = .false.
    end where
  
    where (h_smb >= eps_dp)   ! ## CASE 1: glaciated point ##
        
      ! Actual mass balance, based on temporal change in ice sheet
      mb_source_apl = (h_smb-h_adv)/dtime
        
      ! Store melt quantities here for later accounting
      runoff_top  = runoff
      disc_top    = calv_grounded
      runoff_bot  = q_b_tot

      ! Actual SMB
      as_perp_apl = mb_source_apl + (disc_top + runoff_bot)

      ! Store melting mask value. This became an ice point. Regarded
      ! as land always and the melt will be attributed to land.
      mask_run = 1_i2b ! grounded ice   

    elsewhere (h_smb < eps_dp .and. h_adv >= eps_dp)
                              ! ## CASE 2: ice disappeared ##
        
      ! Actual mass balance, based on temporal change in ice sheet
      mb_source_apl = (h_smb-h_adv)/dtime
        
      ! Assume all hidden melt comes from the surface
      ! (since no basal calculations were made at these points)
      where (runoff+calv_grounded > 0.0_dp)
      ! share by the respective melting continents by their fractions
        runoff_top = -mb_source_apl*runoff/(runoff+calv_grounded)
        disc_top   = -mb_source_apl*calv_grounded/(runoff+calv_grounded)
      elsewhere (runoff > 0.0_dp .and. calv_grounded <= 0.0_dp)
      ! only ordinary surface melt applied 
        runoff_top = -mb_source_apl
      elsewhere (calv_grounded > 0.0_dp .and. runoff <= 0.0_dp)
      ! only surface discharge applied
        disc_top   = -mb_source_apl
      elsewhere
      ! safety measure (hidden would be small here anyhow)
        runoff_top = -0.5_dp*mb_source_apl
        disc_top   = -0.5_dp*mb_source_apl
      end where

      ! Actual SMB
      as_perp_apl = mb_source_apl + (disc_top + runoff_bot)

      ! Store melting mask value. This grid point is ice free now, therefore it
      ! can be either ocean or land (using actual maske value). 
      where (flag_ocean_point)   ! Ocean point (hidden melt)
        mask_run = -2_i2b
      elsewhere                  ! Land (hidden melt)
        mask_run = -1_i2b
      end where

    end where

  end where

! Update the ice thickness
  h_neu = h_smb

end subroutine apply_mb_source

!-------------------------------------------------------------------------------
!> Adjustment of the newly computed ice thickness distribution.
!<------------------------------------------------------------------------------
  subroutine thk_adjust(time, dtime)

  implicit none

  real(dp), intent(in) :: time, dtime

  real(dp), dimension(0:JMAX,0:IMAX) :: H_neu_tmp
  real(dp)                           :: dtime_inv

!-------- Saving computed H_neu before any adjustments --------

  h_neu_tmp = h_neu

!-------- Correct negative thickness values --------

  where (h_neu < 0.0_dp) h_neu = 0.0_dp

!-------- Further adjustments --------

#if (THK_EVOL==0)

!  ------ No evolution of the ice thickness

  h_neu  = h   ! newly computed ice thickness is discarded

#elif (THK_EVOL==1)

!  ------ No adjustment, ice thickness evolves freely;
!         thus nothing to be done

#elif (THK_EVOL==2)

!  ------ Adjustment due to prescribed target topography

  if (time >= time_target_topo_final) then

     h_neu = h_target

  else if (time >= time_target_topo_init) then

     target_topo_tau = (time_target_topo_final-time)/3.0_dp

     h_neu =   ( target_topo_tau*h_neu + dtime*h_target ) &
             / ( target_topo_tau       + dtime )

  end if

#elif (THK_EVOL==3)

!  ------ Adjustment due to prescribed target topography with relaxation time

  ! relaxation time target_topo_tau defined in routine sico_init

  h_neu =   ( target_topo_tau*h_neu + dtime*h_target ) &
          / ( target_topo_tau       + dtime )

#elif (THK_EVOL==4)

!  ------ Maximum ice extent constrained by prescribed mask

  where (maske_maxextent == 0_i2b) &   ! not allowed to glaciate
     h_neu = 0.0_dp

#else
  stop ' >>> thk_adjust: THK_EVOL must be between 0 and 4!'
#endif

!-------- Computation of the mass balance adjustment --------

  dtime_inv = 1.0_dp/dtime

  mb_adjust = (h_neu-h_neu_tmp)*dtime_inv

end subroutine thk_adjust

!-------------------------------------------------------------------------------
!> Update of the ice-land-ocean mask etc.
!<------------------------------------------------------------------------------
subroutine calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, &
                                n_calc_thk_mask_update_aux)

use topograd_m

implicit none

real(dp),     intent(in) :: time
real(dp),     intent(in) :: dtime, dxi, deta
real(dp),     intent(in) :: z_sl, z_mar
integer(i2b), intent(in) :: n_calc_thk_mask_update_aux

integer(i4b)                       :: i, j
real(dp)                           :: year_sec_inv
real(dp), dimension(0:JMAX,0:IMAX) :: H_neu_tmp
real(dp)                           :: dtime_inv

!-------- Term abbreviations --------

year_sec_inv = 1.0_dp/year_sec

!-------- Saving computed H_neu before any modifications --------

h_neu_tmp = h_neu

!-------- Update of the mask --------

dtime_inv = 1.0_dp/dtime

if (n_calc_thk_mask_update_aux == 1_i2b) then
   call calc_thk_mask_update_aux1(time, dtime, dtime_inv, z_sl, z_mar)
else if (n_calc_thk_mask_update_aux == 2_i2b) then
   call calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)
else if (n_calc_thk_mask_update_aux == 3_i2b) then
   call calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)
else
   write(6, fmt='(a)') ' >>> calc_thk_mask_update:'
   write(6, fmt='(a)') '        n_calc_thk_mask_update_aux has no valid value!'
   stop
end if

!-------- Time derivatives --------

dzs_dtau  = (zs_neu-zs)*dtime_inv
dzb_dtau  = (zb_neu-zb)*dtime_inv
dzm_dtau  = dh_t_dtau+dzb_dtau
dh_c_dtau = dzs_dtau-dzm_dtau

#if (THK_EVOL==2)
if ( abs((time-time_target_topo_final)*year_sec_inv) < eps ) then
   dzs_dtau  = 0.0_dp   ! Introduced
   dzb_dtau  = 0.0_dp   ! by
   dzm_dtau  = 0.0_dp   ! Tatsuru Sato
   dh_c_dtau = 0.0_dp   ! for
   dh_t_dtau = 0.0_dp   ! stability reasons
end if
#endif

!-------- New gradients --------

#if (TOPOGRAD==0)
call topograd_1(dxi, deta, 2)
#elif (TOPOGRAD==1)
call topograd_2(dxi, deta, 2)
#endif

!-------- Compute the volume flux across the CTS, am_perp --------

#if (CALCMOD==1)

do i=1, imax-1
do j=1, jmax-1

   if ( ((maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b)) &
        .and.(n_cts(j,i) == 1_i2b)) then

      am_perp_st(j,i) = &
                0.5_dp*(vx_c(0,j,i)+vx_c(0,j,i-1))*dzm_dxi_g(j,i) &
              + 0.5_dp*(vy_c(0,j,i)+vy_c(0,j-1,i))*dzm_deta_g(j,i) &
              - 0.5_dp*(vz_c(0,j,i)+vz_t(ktmax-1,j,i))
      am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)

   else
      am_perp_st(j,i) = 0.0_dp
      am_perp(j,i)    = 0.0_dp
   end if

end do
end do

#endif

!-------- Compute the applied surface mass balance as_perp_apl --------

#if (MB_ACCOUNT==0)

where (     (maske==0_i2b).or.(maske_old==0_i2b) &
        .or.(maske==3_i2b).or.(maske_old==3_i2b) )
                  ! grounded or floating ice before or after the time step

   as_perp_apl = as_perp + mb_adjust

elsewhere   ! maske==1_i2b or 2_i2b, ice-free land or sea

   as_perp_apl = 0.0_dp

end where

#endif

end subroutine calc_thk_mask_update

!-------------------------------------------------------------------------------
!> Update of the ice-land-ocean mask for SIA-only dynamics of ice sheets
!! without ice shelves.
!<------------------------------------------------------------------------------
subroutine calc_thk_mask_update_aux1(time, dtime, dtime_inv, z_sl, z_mar)

implicit none

real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar

integer(i4b) :: i, j, kc, kt
real(dp)     :: rhosw_rho_ratio, rho_rhosw_ratio
logical      :: flag_calving_event

real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance

real(dp), parameter :: H_isol_max = 1.0e+03_dp
                       ! Thickness limit of isolated glaciated points

!-------- Term abbreviations --------

rhosw_rho_ratio = rho_sw/rho
rho_rhosw_ratio = rho/rho_sw

!-------- Update of the mask --------

zs_neu    = zb_neu + h_neu   ! ice surface topography
h_sea_neu = z_sl - zl_neu    ! sea depth
h_balance = h_sea_neu*rhosw_rho_ratio

#if (THK_EVOL>=1)

do i=0, imax
do j=0, jmax

   if (maske(j,i) <= 1_i2b) then   ! grounded ice or ice-free land

      if (h_neu(j,i) >= eps_dp) then
         maske(j,i) = 0_i2b   ! grounded ice
      else
         maske(j,i) = 1_i2b   ! ice-free land
      end if

#if (MARGIN==2)

   else   ! (maske(j,i) == 2_i2b); sea

      if (h_neu(j,i) >= eps_dp) then

         if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then

#if (MARINE_ICE_FORMATION==1)
            maske(j,i) = 2_i2b   ! floating ice cut off -> sea
#elif (MARINE_ICE_FORMATION==2)
            maske(j,i) = 0_i2b   ! "underwater ice"
#endif
         else
            maske(j,i) = 0_i2b   ! grounded ice
         end if

#if (MARINE_ICE_CALVING==2 \
      || marine_ice_calving==4 \
      || marine_ice_calving==6)
         if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b   ! sea
#elif (MARINE_ICE_CALVING==3 \
        || marine_ice_calving==5 \
        || marine_ice_calving==7)
         if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b   ! sea
#endif

      else
         maske(j,i) = 2_i2b   ! sea
      end if

#endif

   end if

end do
end do

#endif

!  ------ Adjustment due to prescribed target topography

#if (THK_EVOL==2)
if (time >= time_target_topo_final) maske = maske_target
#endif

!-------- Detection of the grounding line and the calving front
!                          (not relevant for SIA-only dynamics) --------

flag_grounding_line_1 = .false.
flag_grounding_line_2 = .false.

flag_calving_front_1  = .false.
flag_calving_front_2  = .false.

!-------- Correction of zs_neu and zb_neu for ice-free land and sea --------

do i=0, imax
do j=0, jmax

   if (maske(j,i) == 1_i2b) then   ! ice-free land

      zs_neu(j,i) = zb_neu(j,i)   ! this prevents zs_neu(j,i)
      h_neu(j,i)  = 0.0_dp        ! from being below zb_neu(j,i)

   else if (maske(j,i) == 2_i2b) then   ! sea

      zs_neu(j,i) = zb_neu(j,i)   ! this prevents zs_neu(j,i)
      h_neu(j,i)  = 0.0_dp        ! from being below zb_neu(j,i)

   else if (maske(j,i) == 3_i2b) then   ! floating ice

      write(6, fmt='(a)') ' >>> calc_thk_mask_update_aux1:'
      write(6, fmt='(a)') '          maske(j,i)==3 not allowed for'
      write(6, fmt='(a)') '          SIA-only dynamics!'
      stop

   end if

end do
end do

!-------- Computation of further quantities --------

do i=0, imax
do j=0, jmax

   if (maske(j,i) == 0_i2b) then   ! grounded ice

!  ------ Limit the thickness of isolated glaciated points to H_isol_max

      if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
           .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
           .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
           .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
         ) then   ! all nearest neighbours ice-free
         h_neu(j,i)   = min(h_neu(j,i), h_isol_max)
         zs_neu(j,i)  = zb_neu(j,i)+h_neu(j,i)
      end if

!  ------ Further stuff

      if (n_cts(j,i) == 1) then
         h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
         zm_neu(j,i)  = zb_neu(j,i)+h_t_neu(j,i)
         h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
      else
         h_t_neu(j,i) = 0.0_dp
         zm_neu(j,i)  = zb_neu(j,i)
         h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
      end if

   else   ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea

      zs_neu(j,i) = zb_neu(j,i)
      zm_neu(j,i) = zb_neu(j,i)
      h_c_neu(j,i) = 0.0_dp
      h_t_neu(j,i) = 0.0_dp
      h_neu(j,i)   = 0.0_dp

   end if

end do
end do

end subroutine calc_thk_mask_update_aux1

!-------------------------------------------------------------------------------
!> Update of the ice-land-ocean mask for hybrid SIA/SStA dynamics of ice sheets
!! without ice shelves.
!<------------------------------------------------------------------------------
subroutine calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)

implicit none

real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar

integer(i4b) :: i, j, kc, kt
real(dp)     :: rhosw_rho_ratio, rho_rhosw_ratio
logical      :: flag_calving_event

real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance

real(dp), parameter :: H_isol_max = 1.0e+03_dp
                       ! Thickness limit of isolated glaciated points

!-------- Term abbreviations --------

rhosw_rho_ratio = rho_sw/rho
rho_rhosw_ratio = rho/rho_sw

!-------- Update of the mask --------

zs_neu    = zb_neu + h_neu   ! ice surface topography
h_sea_neu = z_sl - zl_neu    ! sea depth
h_balance = h_sea_neu*rhosw_rho_ratio

#if (THK_EVOL>=1)

do i=0, imax
do j=0, jmax

   if (maske(j,i) <= 1_i2b) then   ! grounded ice or ice-free land

      if (h_neu(j,i) >= eps_dp) then
         maske(j,i) = 0_i2b   ! grounded ice
      else
         maske(j,i) = 1_i2b   ! ice-free land
      end if

#if (MARGIN==2)

   else   ! (maske(j,i) == 2_i2b); sea

      if (h_neu(j,i) >= eps_dp) then

         if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then

#if (MARINE_ICE_FORMATION==1)
            maske(j,i) = 2_i2b   ! floating ice cut off -> sea
#elif (MARINE_ICE_FORMATION==2)
            maske(j,i) = 0_i2b   ! "underwater ice"
#endif
         else
            maske(j,i) = 0_i2b   ! grounded ice
         end if

#if (MARINE_ICE_CALVING==2 \
      || marine_ice_calving==4 \
      || marine_ice_calving==6)
         if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b   ! sea
#elif (MARINE_ICE_CALVING==3 \
        || marine_ice_calving==5 \
        || marine_ice_calving==7)
         if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b   ! sea
#endif

      else
         maske(j,i) = 2_i2b   ! sea
      end if

#endif

   end if

end do
end do

#endif

!  ------ Adjustment due to prescribed target topography

#if (THK_EVOL==2)
if (time >= time_target_topo_final) maske = maske_target
#endif

!-------- Detection of the grounding line and the calving front
!                   (not relevant for hybrid SIA/SStA dynamics
!                                 of ice sheets without ice shelves) --------

flag_grounding_line_1 = .false.
flag_grounding_line_2 = .false.

flag_calving_front_1  = .false.
flag_calving_front_2  = .false.

!-------- Correction of zs_neu and zb_neu for ice-free land and sea --------

do i=0, imax
do j=0, jmax

   if (maske(j,i) == 1_i2b) then   ! ice-free land

      zs_neu(j,i) = zb_neu(j,i)   ! this prevents zs_neu(j,i)
      h_neu(j,i)  = 0.0_dp        ! from being below zb_neu(j,i)

   else if (maske(j,i) == 2_i2b) then   ! sea

      zs_neu(j,i) = zb_neu(j,i)   ! this prevents zs_neu(j,i)
      h_neu(j,i)  = 0.0_dp        ! from being below zb_neu(j,i)

   else if (maske(j,i) == 3_i2b) then   ! floating ice

      write(6, fmt='(a)') ' >>> calc_thk_mask_update_aux2:'
      write(6, fmt='(a)') '          maske(j,i)==3 not allowed for'
      write(6, fmt='(a)') '          hybrid SIA/SStA dynamics of ice sheets'
      write(6, fmt='(a)') '          without ice shelves!'
      stop

   end if

end do
end do

!-------- Computation of further quantities --------

do i=0, imax
do j=0, jmax

   if (maske(j,i) == 0_i2b) then   ! grounded ice

!  ------ Limit the thickness of isolated glaciated points to H_isol_max

      if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
           .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
           .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
           .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
         ) then   ! all nearest neighbours ice-free
         h_neu(j,i)   = min(h_neu(j,i), h_isol_max)
         zs_neu(j,i)  = zb_neu(j,i)+h_neu(j,i)
      end if

!  ------ Further stuff

      if (n_cts(j,i) == 1) then
         h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
         zm_neu(j,i)  = zb_neu(j,i)+h_t_neu(j,i)
         h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
      else
         h_t_neu(j,i) = 0.0_dp
         zm_neu(j,i)  = zb_neu(j,i)
         h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
      end if

   else   ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea

      zs_neu(j,i) = zb_neu(j,i)
      zm_neu(j,i) = zb_neu(j,i)
      h_c_neu(j,i) = 0.0_dp
      h_t_neu(j,i) = 0.0_dp
      h_neu(j,i)   = 0.0_dp

   end if

end do
end do

end subroutine calc_thk_mask_update_aux2

!-------------------------------------------------------------------------------
!> Update of the ice-land-ocean mask for coupled SIA/SSA or
!! SIA/SStA/SSA dynamics of ice sheets with ice shelves.
!<------------------------------------------------------------------------------
subroutine calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)

implicit none

real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar

integer(i4b) :: i, j, kc, kt
real(dp)     :: rhosw_rho_ratio, rho_rhosw_ratio
logical      :: flag_calving_event

real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance

real(dp), parameter :: H_isol_max = 1.0e+03_dp
                       ! Thickness limit of isolated glaciated points

!-------- Term abbreviations --------

rhosw_rho_ratio = rho_sw/rho
rho_rhosw_ratio = rho/rho_sw

!-------- Update of the mask --------

zs_neu    = zb_neu + h_neu   ! ice surface topography
h_sea_neu = z_sl - zl_neu    ! sea depth
h_balance = h_sea_neu*rhosw_rho_ratio

#if (THK_EVOL>=1)

do i=1, imax-1
do j=1, jmax-1

   ! grounding_line migration check

   if ( ( maske(j,i) <= 1_i2b &
         .and.    (maske(j,i+1)>1_i2b.or.maske(j,i-1)>1_i2b &
               .or.maske(j+1,i)>1_i2b.or.maske(j-1,i)>1_i2b) ) &
        .or. &
        ( maske(j,i)>=2_i2b &
         .and.    (maske(j,i+1)==0_i2b.or.maske(j,i-1)==0_i2b  &
               .or.maske(j+1,i)==0_i2b.or.maske(j-1,i)==0_i2b) ) ) then

      if (h_neu(j,i) >= eps_dp) then

         if (h_neu(j,i)<h_balance(j,i).and.zl_neu(j,i)<z_sl) then
            maske(j,i)    = 3_i2b
            zb_neu(j,i)   = z_sl-rho_rhosw_ratio*h_neu(j,i)
            dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
            zs_neu(j,i)   = zb_neu(j,i)+h_neu(j,i)
         else
            maske(j,i)    = 0_i2b
            zb_neu(j,i)   = zl_neu(j,i)
            dzb_dtau(j,i) = dzl_dtau(j,i)
            zs_neu(j,i)   = zb_neu(j,i)+h_neu(j,i)
         end if

      else   ! if (H_neu(j,i) <= eps_dp) then

         if (zl_neu(j,i)>z_sl) then
            maske(j,i)    = 1_i2b
            zb_neu(j,i)   = zl_neu(j,i)
            dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
            zs_neu(j,i)   = zb_neu(j,i)
         else
            maske(j,i)    = 2_i2b
            zb_neu(j,i)   = z_sl
            dzb_dtau(j,i) = 0.0_dp
            zs_neu(j,i)   = z_sl
         end if

      end if

   else if (maske(j,i) <= 1_i2b) then   ! grounded ice or ice-free land

      if (h_neu(j,i) >= eps_dp) then   ! can change maske 0 or 1
         maske(j,i)    = 0_i2b
         zb_neu(j,i)   = zl_neu(j,i)
         dzb_dtau(j,i) = dzl_dtau(j,i)
         zs_neu(j,i)   = zb_neu(j,i)+h_neu(j,i)
      else
         maske(j,i)    = 1_i2b
         zb_neu(j,i)   = zl_neu(j,i)
         dzb_dtau(j,i) = dzl_dtau(j,i)
         zs_neu(j,i)   = zl_neu(j,i)
      end if

   else   ! if (maske(j,i)==2.or.maske(j,i)==3) then

      if (h_neu(j,i) > eps_dp) then

         if (h_neu(j,i)<h_balance(j,i).and.zl_neu(j,i)<z_sl) then
            maske(j,i)    = 3_i2b
            zb_neu(j,i)   = z_sl-rho_rhosw_ratio*h_neu(j,i)
            dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
            zs_neu(j,i)   = zb_neu(j,i)+h_neu(j,i)
         else
            maske(j,i)    = 0_i2b
            zb_neu(j,i)   = zl_neu(j,i)
            dzb_dtau(j,i) = dzl_dtau(j,i)
            zs_neu(j,i)   = zb_neu(j,i)+h_neu(j,i)
         end if

      else   ! if (H_neu(j,i) <= eps_dp) then

         if (zl_neu(j,i)>z_sl) then
            maske(j,i)    = 1_i2b
            zb_neu(j,i)   = zl_neu(j,i)
            dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
            zs_neu(j,i)   = zb_neu(j,i)
         else
            maske(j,i)    = 2_i2b
            zb_neu(j,i)   = z_sl
            dzb_dtau(j,i) = 0.0_dp
            zs_neu(j,i)   = z_sl
         end if

      end if

   end if

end do
end do

#if (ICE_SHELF_CALVING==2)

do

   flag_calving_front_1 = .false.
   flag_calving_event   = .false.

   do i=1, imax-1
   do j=1, jmax-1

      if ( (maske(j,i)==3_i2b) &   ! floating ice
           .and. &
             (    (maske(j,i+1)==2_i2b)   &   ! with
              .or.(maske(j,i-1)==2_i2b)   &   ! one
              .or.(maske(j+1,i)==2_i2b)   &   ! neighbouring
              .or.(maske(j-1,i)==2_i2b) ) &   ! sea point
         ) &
         flag_calving_front_1(j,i) = .true.   ! preliminary detection
                                              ! of the calving front

      if ( flag_calving_front_1(j,i).and.(h_neu(j,i) < h_calv) ) then
         flag_calving_event = .true.  ! calving event,
         maske(j,i)         = 2_i2b   ! floating ice point changes to sea point
      end if

   end do
   end do

   if (.not.flag_calving_event) exit

end do

#elif (ICE_SHELF_CALVING==9)

#if (defined(MISMIPP))

do i=0, imax
do j=0, jmax

   if ((maske(j,i)==3_i2b).and.(xi(i) > lx)) then
      maske(j,i) = 2_i2b   ! floating ice point changes to sea point
   end if

end do
end do

#else
write(6, fmt='(a)') ' >>> calc_thk_mask_update_aux3:'
write(6, fmt='(a)') '          Option ICE_SHELF_CALVING==9'
write(6, fmt='(a)') '          only defined for MISMIP+!'
stop
#endif

#endif

#endif

!  ------ Adjustment due to prescribed target topography

#if (THK_EVOL==2)
if (time >= time_target_topo_final) maske = maske_target
#endif

!-------- Detection of the grounding line and the calving front --------

flag_grounding_line_1 = .false.
flag_grounding_line_2 = .false.

flag_calving_front_1  = .false.
flag_calving_front_2  = .false.

do i=1, imax-1
do j=1, jmax-1

   if ( (maske(j,i)==0_i2b) &   ! grounded ice
        .and. &
          (    (maske(j,i+1)==3_i2b)   &   ! with
           .or.(maske(j,i-1)==3_i2b)   &   ! one
           .or.(maske(j+1,i)==3_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==3_i2b) ) &   ! floating ice point
      ) &
   flag_grounding_line_1(j,i) = .true.

   if ( (maske(j,i)==3_i2b) &   ! floating ice
        .and. &
          (    (maske(j,i+1)==0_i2b)   &   ! with
           .or.(maske(j,i-1)==0_i2b)   &   ! one
           .or.(maske(j+1,i)==0_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==0_i2b) ) &   ! grounded ice point
      ) &
   flag_grounding_line_2(j,i) = .true.

   if ( (maske(j,i)==3_i2b) &   ! floating ice
        .and. &
          (    (maske(j,i+1)==2_i2b)   &   ! with
           .or.(maske(j,i-1)==2_i2b)   &   ! one
           .or.(maske(j+1,i)==2_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==2_i2b) ) &   ! sea point
      ) &
   flag_calving_front_1(j,i) = .true.

   if ( (maske(j,i)==2_i2b) &   ! sea
        .and. &
          (    (maske(j,i+1)==3_i2b)   &   ! with
           .or.(maske(j,i-1)==3_i2b)   &   ! one
           .or.(maske(j+1,i)==3_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==3_i2b) ) &   ! floating ice point
      ) &
   flag_calving_front_2(j,i) = .true.

end do
end do

do i=1, imax-1

   j=0
   if ( (maske(j,i)==2_i2b) &   ! sea
        .and. (maske(j+1,i)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

   j=jmax
   if ( (maske(j,i)==2_i2b) &   ! sea
        .and. (maske(j-1,i)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

end do

do j=1, jmax-1

   i=0
   if ( (maske(j,i)==2_i2b) &   ! sea
        .and. (maske(j,i+1)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

   i=imax
   if ( (maske(j,i)==2_i2b) &   ! sea
        .and. (maske(j,i-1)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

end do

!-------- Correction of zs_neu and zb_neu
!         for ice-free land, sea and floating ice --------

do i=0, imax
do j=0, jmax

   if (maske(j,i) == 1_i2b) then   ! ice-free land

      zs_neu(j,i) = zb_neu(j,i)   ! this prevents zs_neu(j,i)
      h_neu(j,i)  = 0.0_dp        ! from being below zb_neu(j,i)

   else if (maske(j,i) == 2_i2b) then   ! sea

      zs_neu(j,i) = z_sl          ! both zs_neu(j,i) and zb_neu(j,i)
      zb_neu(j,i) = z_sl          ! set to the current sea level stand
      h_neu(j,i)  = 0.0_dp

   else if (maske(j,i) == 3_i2b) then   ! floating ice

      h_neu(j,i) = zs_neu(j,i)-zb_neu(j,i)   ! ice thickness

      zb_neu(j,i) = z_sl - rho_rhosw_ratio*h_neu(j,i)   ! floating condition
      zs_neu(j,i) = zb_neu(j,i) + h_neu(j,i)

   end if

end do
end do

!-------- Computation of further quantities --------

do i=0, imax
do j=0, jmax

   if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
                                        ! grounded or floating ice

!  ------ Limit the thickness of isolated glaciated points to H_isol_max

      if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
           .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
           .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
           .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
         ) then   ! all nearest neighbours ice-free
         h_neu(j,i)   = min(h_neu(j,i), h_isol_max)
         zs_neu(j,i)  = zb_neu(j,i)+h_neu(j,i)
      end if

!  ------ Further stuff

      if (n_cts(j,i) == 1) then
         h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
         zm_neu(j,i)  = zb_neu(j,i)+h_t_neu(j,i)
         h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
      else
         h_t_neu(j,i) = 0.0_dp
         zm_neu(j,i)  = zb_neu(j,i)
         h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
      end if

   else   ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea

      zs_neu(j,i) = zb_neu(j,i)
      zm_neu(j,i) = zb_neu(j,i)
      h_c_neu(j,i) = 0.0_dp
      h_t_neu(j,i) = 0.0_dp
      h_neu(j,i)   = 0.0_dp

   end if

end do
end do

end subroutine calc_thk_mask_update_aux3

!-------------------------------------------------------------------------------

end module calc_thk_m
!