calc_vxy_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  c a l c _ v x y _ m
!
!> @file
!!
!! Computation of the horizontal velocity vx, vy.
!!
!! @section Copyright
!!
!! Copyright 2009-2017 Ralf Greve, Tatsuru Sato, Thomas Goelles, Jorge Bernales
!!
!! @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 horizontal velocity vx, vy.
!<------------------------------------------------------------------------------
module calc_vxy_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m

  implicit none

  private
  public :: calc_vxy_b_sia, calc_vxy_sia, calc_vxy_static, calc_vxy_ssa

contains

!-------------------------------------------------------------------------------
!> Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice
!! approximation.
!<------------------------------------------------------------------------------
subroutine calc_vxy_b_sia(time, z_sl)

implicit none

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

integer(i4b) :: i, j
integer(i4b) :: r_smw_1, s_smw_1, r_smw_2, s_smw_2
real(dp), dimension(0:JMAX,0:IMAX) :: gamma_slide_inv
real(dp) :: gamma_slide_inv_1, gamma_slide_inv_2
real(dp) :: smw_coeff_1, smw_coeff_2
real(dp), dimension(0:JMAX,0:IMAX) :: p_b, p_b_red, tau_b
real(dp) :: cvxy1, cvxy1a, cvxy1b, ctau1, ctau1a, ctau1b
real(dp) :: temp_diff
real(dp) :: c_Hw_slide, Hw0_slide, Hw0_slide_inv
real(dp) :: vh_max
real(dp) :: year_sec_inv, time_in_years
real(dp) :: ramp_up_factor
logical, dimension(0:JMAX,0:IMAX) :: sub_melt_flag

year_sec_inv  = 1.0_dp/year_sec
time_in_years = time*year_sec_inv

!-------- Sliding-law coefficients --------

#if (defined(ANT) && (!defined(SEDI_SLIDE) || SEDI_SLIDE==1))
      ! Antarctica without soft sediment,
      ! one sliding law for the entire modelling domain

gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)

do i=0, imax
do j=0, jmax
   c_slide(j,i)         = c_slide*year_sec_inv
   p_weert(j,i)         = p_weert
   q_weert(j,i)         = q_weert
   gamma_slide_inv(j,i) = gamma_slide_inv_1
   sub_melt_flag(j,i)   = (gamma_slide >= eps)
end do
end do

#elif (defined(ANT) && SEDI_SLIDE==2)
        ! Antarctica with soft sediment,
        ! different sliding laws for hard rock and soft sediment

gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)

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

   if (maske_sedi(j,i) /= 7) then   ! no soft sediment

#if (defined(TRANS_LENGTH_SL))
      c_slide(j,i)         = ( (1.0_dp-r_mask_sedi(j,i))*c_slide &
                                      +r_mask_sedi(j,i) *c_slide_sedi ) &
                             *year_sec_inv
#else
      c_slide(j,i)         = c_slide*year_sec_inv
#endif

      p_weert(j,i)         = p_weert
      q_weert(j,i)         = q_weert
      gamma_slide_inv(j,i) = gamma_slide_inv_1
      sub_melt_flag(j,i)   = (gamma_slide >= eps)

   else   ! maske_sedi(j,i) == 7, soft sediment

#if (defined(TRANS_LENGTH_SL))
      c_slide(j,i)         = ( (1.0_dp-r_mask_sedi(j,i))*c_slide &
                                      +r_mask_sedi(j,i) *c_slide_sedi ) &
                             *year_sec_inv
#else
      c_slide(j,i)         = c_slide_sedi*year_sec_inv
#endif

      p_weert(j,i)         = p_weert_sedi
      q_weert(j,i)         = q_weert_sedi
      gamma_slide_inv(j,i) = gamma_slide_inv_2
      sub_melt_flag(j,i)   = (gamma_slide_sedi >= eps)

   end if

end do
end do

#elif (defined(GRL) && ICE_STREAM==1 \
        || defined(asf) && ice_stream==1)
        ! Greenland or Austfonna without ice streams,
        ! one sliding law for the entire modelling domain,
        ! with surface meltwater effect

gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)

r_smw_1     = r_smw
s_smw_1     = s_smw
smw_coeff_1 = smw_coeff*year_sec**s_smw_1

do i=0, imax
do j=0, jmax
   c_slide(j,i)         = c_slide*year_sec_inv &
                                 *(1.0_dp+smw_coeff_1*runoff(j,i)**s_smw_1 &
                                   /max((h_c(j,i)+h_t(j,i))**r_smw_1, eps))
   p_weert(j,i)         = p_weert
   q_weert(j,i)         = q_weert
   gamma_slide_inv(j,i) = gamma_slide_inv_1
   sub_melt_flag(j,i)   = (gamma_slide >= eps)
end do
end do

#elif (defined(GRL) && ICE_STREAM==2 \
        || defined(asf) && ice_stream==2)
        ! Greenland or Austfonna with ice streams,
        ! different sliding laws for hard rock and ice streams

gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)

r_smw_1     = r_smw
s_smw_1     = s_smw
smw_coeff_1 = smw_coeff     *year_sec**s_smw_1
r_smw_2     = r_smw_sedi
s_smw_2     = s_smw_sedi
smw_coeff_2 = smw_coeff_sedi*year_sec**s_smw_2

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

   if (maske_sedi(j,i) /= 7) then   ! no ice stream
      c_slide(j,i)         = c_slide*year_sec_inv &
                                    *(1.0_dp+smw_coeff_1*runoff(j,i)**s_smw_1 &
                                      /max((h_c(j,i)+h_t(j,i))**r_smw_1, eps))
      p_weert(j,i)         = p_weert
      q_weert(j,i)         = q_weert
      gamma_slide_inv(j,i) = gamma_slide_inv_1
      sub_melt_flag(j,i)   = (gamma_slide >= eps)
   else   ! maske_sedi(j,i) == 7, ice stream
      c_slide(j,i)         = c_slide_sedi*year_sec_inv &
                                         *(1.0_dp+smw_coeff_2*runoff(j,i)**s_smw_2 &
                                           /max((h_c(j,i)+h_t(j,i))**r_smw_2, eps))
      p_weert(j,i)         = p_weert_sedi
      q_weert(j,i)         = q_weert_sedi
      gamma_slide_inv(j,i) = gamma_slide_inv_2
      sub_melt_flag(j,i)   = (gamma_slide_sedi >= eps)
   end if

end do
end do

#elif (defined(XYZ) && defined(HEINO))
        ! ISMIP HEINO,
        ! different sliding laws for hard rock and soft sediment

gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)

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

   if (maske_sedi(j,i) /= 7) then   ! no soft sediment
      c_slide(j,i)         = c_slide*year_sec_inv
      p_weert(j,i)         = p_weert
      q_weert(j,i)         = q_weert
      gamma_slide_inv(j,i) = gamma_slide_inv_1
      sub_melt_flag(j,i)   = (gamma_slide >= eps)
   else   ! maske_sedi(j,i) == 7, soft sediment
      c_slide(j,i)         = c_slide_sedi*year_sec_inv
      p_weert(j,i)         = p_weert_sedi
      q_weert(j,i)         = q_weert_sedi
      gamma_slide_inv(j,i) = gamma_slide_inv_2
      sub_melt_flag(j,i)   = (gamma_slide_sedi >= eps)
   end if

end do
end do

#elif (defined(SCAND) \
        || defined(nhem) \
        || defined(tibet) \
        || defined(nmars) \
        || defined(smars) \
        || defined(emtp2sge) \
        || defined(xyz))
           ! one sliding law for the entire modelling domain

gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)

do i=0, imax
do j=0, jmax
   c_slide(j,i)         = c_slide*year_sec_inv
   p_weert(j,i)         = p_weert
   q_weert(j,i)         = q_weert
   gamma_slide_inv(j,i) = gamma_slide_inv_1
   sub_melt_flag(j,i)   = (gamma_slide >= eps)
end do
end do

#else

stop ' >>> calc_vxy_b_sia: No valid domain specified!'

#endif

do i=0, imax
do j=0, jmax
   p_weert_inv(j,i) = 1.0_dp/max(real(p_weert(j,i),dp), eps)
end do
end do

!  ------ Ramping up basal sliding

ramp_up_factor = 1.0_dp

#if (defined(TIME_RAMP_UP_SLIDE))

if ( (time_ramp_up_slide > eps_dp) &
     .and. &
     ((time_in_years-(time_init0)) < (time_ramp_up_slide)) ) then

   ramp_up_factor = (time_in_years-(time_init0))/(time_ramp_up_slide)

   ramp_up_factor = max(min(ramp_up_factor, 1.0_dp), 0.0_dp)
                                  ! constrain to interval [0,1]

   ramp_up_factor = ramp_up_factor*ramp_up_factor*ramp_up_factor &
                                  *(10.0_dp + ramp_up_factor &
                                              *(-15.0_dp+6.0_dp*ramp_up_factor))
                                  ! make transition smooth (quintic function)

   c_slide = c_slide * ramp_up_factor

end if

#endif

!  ------ Coefficients for the contribution of the basal water layer

#if (BASAL_HYDROLOGY==1)

#if (defined(C_HW_SLIDE))
  c_hw_slide = c_hw_slide
#else
  c_hw_slide = 0.0_dp
#endif

#if (defined(HW0_SLIDE))
  hw0_slide = hw0_slide
#else
  hw0_slide = 1.0_dp
#endif

  hw0_slide_inv = 1.0_dp/max(hw0_slide, eps_dp)

#else
  c_hw_slide    = 0.0_dp   ! dummy value
  hw0_slide     = 1.0_dp   ! dummy value
  hw0_slide_inv = 1.0_dp   ! dummy value
#endif

!-------- Computation of basal stresses --------

!  ------ Basal pressure p_b, basal water pressure p_b_w,
!         reduced pressure p_b_red

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

   if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
                     ! grounded ice, or floating ice at the grounding line

      p_b(j,i)     = rho*g*(h_c(j,i)+h_t(j,i))
      p_b_w(j,i)   = rho_sw*g*max((z_sl-zb(j,i)), 0.0_dp)
      p_b_red(j,i) = p_b(j,i)-p_b_w(j,i)
      p_b_red(j,i) = max(p_b_red(j,i), red_pres_limit_fact*p_b(j,i))   
                     ! in order to avoid very small values, which may lead to
                     ! huge sliding velocities

   else   ! maske(j,i) == 1, 2 or 3 away from the grounding line

      p_b(j,i)     = 0.0_dp
      p_b_w(j,i)   = 0.0_dp
      p_b_red(j,i) = 0.0_dp

   end if

end do
end do

!  ------ Absolute value of the basal shear stress, tau_b

do i=0, imax
do j=0, jmax
      tau_b(j,i) = p_b(j,i)*(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
end do
end do

!-------- Computation of d_help_b (defined on the grid points (i,j)) --------

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

   if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
                     ! grounded ice, or floating ice at the grounding line

!  ------ Abbreviations

#if (SLIDE_LAW==1)
      cvxy1 = c_slide(j,i) &
              * (tau_b(j,i)**(p_weert(j,i)-1)/p_b(j,i)**q_weert(j,i)) &
              * p_b(j,i)
      ctau1 = 1.0_dp/c_slide(j,i)**p_weert_inv(j,i) &
              * p_b(j,i)**(q_weert(j,i)*p_weert_inv(j,i))
              ! Basal sliding at pressure melting
#elif (SLIDE_LAW==2)
      cvxy1 = c_slide(j,i) &
              * (tau_b(j,i)**(p_weert(j,i)-1)/p_b_red(j,i)**q_weert(j,i)) &
              * p_b(j,i)
      ctau1 = 1.0_dp/c_slide(j,i)**p_weert_inv(j,i) &
              * p_b_red(j,i)**(q_weert(j,i)*p_weert_inv(j,i))
              ! Basal sliding at pressure melting
#endif

!  ------ d_help_b, c_drag

      if (n_cts(j,i) == -1) then   ! cold ice base

         if (sub_melt_flag(j,i)) then
            temp_diff = max((temp_c_m(0,j,i)-temp_c(0,j,i)), 0.0_dp)
            cvxy1a    = exp(-gamma_slide_inv(j,i)*temp_diff)  ! sub-melt sliding
            ctau1a    = 1.0_dp/cvxy1a**p_weert_inv(j,i)
         else
            cvxy1a    = 0.0_dp   ! no sub-melt sliding
            ctau1a    = 1.0_dp/eps**p_weert_inv(j,i)   ! dummy value
         end if

         d_help_b(j,i) = cvxy1*cvxy1a
         c_drag(j,i)   = ctau1*ctau1a

      else if (n_cts(j,i) == 0) then   ! temperate ice base

         d_help_b(j,i) = cvxy1   ! basal sliding
         c_drag(j,i)   = ctau1   ! (pressure-melting conditions)

      else   ! n_cts(j,i) == 1, temperate ice layer

         d_help_b(j,i) = cvxy1   ! basal sliding
         c_drag(j,i)   = ctau1   ! (pressure-melting conditions)

      end if

!    ---- Contribution of the basal water layer

#if (BASAL_HYDROLOGY==1)

      cvxy1b = 1.0_dp &
               + c_hw_slide*(1.0_dp-exp(-max(h_w(j,i),0.0_dp)*hw0_slide_inv))

      ctau1b = 1.0_dp/cvxy1b**p_weert_inv(j,i)

      d_help_b(j,i) = d_help_b(j,i) *cvxy1b
      c_drag(j,i)   = c_drag(j,i)   *ctau1b

#endif

   else   ! maske(j,i) == 1, 2 or 3 away from the grounding line

      d_help_b(j,i) = 0.0_dp
      c_drag(j,i)   = 0.0_dp

   end if

end do
end do

!-------- Computation of vx_b (defined at (i+1/2,j)) --------

do i=0, imax-1
do j=1, jmax-1
   vx_b(j,i) = -0.5_dp*(d_help_b(j,i)+d_help_b(j,i+1))*dzs_dxi(j,i)
end do
end do

!-------- Computation of vy_b (defined at (i,j+1/2)) --------

do i=1, imax-1
do j=0, jmax-1
   vy_b(j,i) = -0.5_dp*(d_help_b(j,i)+d_help_b(j+1,i))*dzs_deta(j,i)
end do
end do

!-------- Computation of vx_b_g and vy_b_g (defined at (i,j)) --------

do i=0, imax
do j=0, jmax
   vx_b_g(j,i) = -d_help_b(j,i)*dzs_dxi_g(j,i)
   vy_b_g(j,i) = -d_help_b(j,i)*dzs_deta_g(j,i)
end do
end do

!-------- Limitation of computed vx_b, vy_b, vx_b_g, vy_b_g to the interval
!         [-VH_MAX, VH_MAX] --------

vh_max = vh_max/year_sec

do i=0, imax
do j=0, jmax
   vx_b(j,i)   = max(vx_b(j,i),   -vh_max)
   vx_b(j,i)   = min(vx_b(j,i),    vh_max)
   vy_b(j,i)   = max(vy_b(j,i),   -vh_max)
   vy_b(j,i)   = min(vy_b(j,i),    vh_max)
   vx_b_g(j,i) = max(vx_b_g(j,i), -vh_max)
   vx_b_g(j,i) = min(vx_b_g(j,i),  vh_max)
   vy_b_g(j,i) = max(vy_b_g(j,i), -vh_max)
   vy_b_g(j,i) = min(vy_b_g(j,i),  vh_max)
end do
end do

end subroutine calc_vxy_b_sia

!-------------------------------------------------------------------------------
!> Computation of the shear stresses txz, tyz, the effective shear stress
!! sigma, the depth-averaged fluidity flui_ave_sia, the horizontal
!! velocity vx, vy and the horizontal volume flux qx, qy in the shallow ice
!! approximation.
!<------------------------------------------------------------------------------
subroutine calc_vxy_sia(dzeta_c, dzeta_t)

use ice_material_properties_m, only : ratefac_c, ratefac_t, ratefac_c_t, creep

implicit none

real(dp), intent(in) :: dzeta_c, dzeta_t

integer(i4b) :: i, j, kc, kt
real(dp) :: avxy3(0:kcmax), aqxy1(0:kcmax)
real(dp) :: ctxyz1(0:kcmax,0:jmax,0:imax), &
            ctxyz2(0:ktmax,0:jmax,0:imax)
real(dp) :: flui_t(0:ktmax), flui_c(0:kcmax)
real(dp) :: cflui0(0:ktmax), cflui1(0:kcmax)
real(dp) :: cvxy2(0:ktmax), cvxy3(0:kcmax)
real(dp) :: cqxy0(0:ktmax), cqxy1(0:kcmax)
real(dp) :: vh_max
real(dp) :: ratio_sl_threshold

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

do kc=0, kcmax
   if (flag_aa_nonzero) then
      avxy3(kc) = aa*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
      aqxy1(kc) = aa/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
   else
      avxy3(kc) = dzeta_c
      aqxy1(kc) = dzeta_c
   end if
end do

!-------- Computation of stresses --------

!  ------ Term abbreviations

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

   if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
                     ! grounded ice, or floating ice at the grounding line

      do kc=0, kcmax
         ctxyz1(kc,j,i) = rho*g*h_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
      end do

      if (n_cts(j,i) == 1) then   ! temperate layer

         do kt=0, ktmax
            ctxyz2(kt,j,i) = rho*g*h_t(j,i)*(1.0_dp-zeta_t(kt))
         end do

      else   ! cold base (-1), temperate base (0)

         do kt=0, ktmax
            ctxyz2(kt,j,i) = 0.0_dp
         end do

      end if

   else   ! maske(j,i) == 1, 2 or 3 away from the grounding line

      do kc=0, kcmax
         ctxyz1(kc,j,i) = 0.0_dp
      end do

      do kt=0, ktmax
         ctxyz2(kt,j,i) = 0.0_dp
      end do

   end if

end do
end do

!  ------ Shear stress txz (defined at (i+1/2,j,kc/t))

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

   do kc=0, kcmax
      txz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j,i+1)) &
                      *dzs_dxi(j,i)
   end do

   do kt=0, ktmax
      txz_t(kt,j,i) = txz_c(0,j,i) &
                      -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j,i+1)) &
                      *dzs_dxi(j,i)
   end do

end do
end do

!  ------ Shear stress tyz (defined at (i,j+1/2,kc/t))

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

   do kc=0, kcmax
      tyz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j+1,i)) &
                      *dzs_deta(j,i)
   end do

   do kt=0, ktmax
      tyz_t(kt,j,i) = tyz_c(0,j,i) &
                      -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j+1,i)) &
                      *dzs_deta(j,i)
   end do

end do
end do

!  ------ Effective shear stress sigma (defined at (i,j,kc/t))

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

   do kc=0, kcmax
      sigma_c(kc,j,i) = ctxyz1(kc,j,i) &
           *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
   end do

   do kt=0, ktmax
      sigma_t(kt,j,i) = sigma_c(0,j,i) &
         + ctxyz2(kt,j,i) &
           *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
   end do

end do
end do

!-------- Computation of the depth-averaged fluidity
!                 (defined on the grid points (i,j)) --------

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

   if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
                     ! grounded ice, or floating ice at the grounding line

!  ------ Fluidity, abbreviations

      do kt=0, ktmax
         flui_t(kt) = 2.0_dp &
                      *enh_t(kt,j,i) &
                      *ratefac_t(omega_t(kt,j,i)) &
                      *creep(sigma_t(kt,j,i))
         cflui0(kt) = h_t(j,i)*flui_t(kt)*dzeta_t
      end do

      do kc=0, kcmax
         flui_c(kc) = 2.0_dp &
                      *enh_c(kc,j,i) &
#if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
                      *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
#elif (CALCMOD==2 || CALCMOD==3)
                      *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), &
                                                   temp_c_m(kc,j,i)) &
#endif
                      *creep(sigma_c(kc,j,i))
         cflui1(kc) = aqxy1(kc)*h_c(j,i)*flui_c(kc)
      end do

!  ------ Depth average

      flui_ave_sia(j,i) = 0.0_dp

      if (n_cts(j,i) == 1) then

         do kt=0, ktmax-1
            flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui0(kt+1)+cflui0(kt))
         end do

      end if

      do kc=0, kcmax-1
         flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui1(kc+1)+cflui1(kc))
      end do

      flui_ave_sia(j,i) = flui_ave_sia(j,i)/(h_c(j,i)+h_t(j,i))

   else   ! maske(j,i) == 1, 2 or 3 away from the grounding line

      flui_ave_sia(j,i) = 0.0_dp

   end if

end do
end do

!-------- Computation of d_help_c/t
!         (defined on the grid points (i,j,kc/t)) --------

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

   if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
                     ! grounded ice, or floating ice at the grounding line

!  ------ Abbreviations

      do kt=0, ktmax
         cvxy2(kt) = 2.0_dp*h_t(j,i) &
                     *enh_t(kt,j,i) &
                     *ratefac_t(omega_t(kt,j,i)) &
                     *creep(sigma_t(kt,j,i)) &
                     *(ctxyz1(0,j,i)+ctxyz2(kt,j,i)) &
                     *dzeta_t
      end do

      do kc=0, kcmax
         cvxy3(kc) = 2.0_dp*avxy3(kc)*h_c(j,i) &
                     *enh_c(kc,j,i) &
#if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
                     *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
#elif (CALCMOD==2 || CALCMOD==3)
                     *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), &
                                                  temp_c_m(kc,j,i)) &
#endif
                     *creep(sigma_c(kc,j,i)) &
                     *ctxyz1(kc,j,i)
      end do

!  ------ d_help_c, d_help_t

      if (n_cts(j,i) == -1) then   ! cold ice base

         do kt=0, ktmax
            d_help_t(kt,j,i) = d_help_b(j,i)
         end do

         d_help_c(0,j,i) = d_help_t(ktmax,j,i)

         do kc=0, kcmax-1
            d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
                                +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
         end do

      else if (n_cts(j,i) == 0) then   ! temperate ice base

         do kt=0, ktmax
            d_help_t(kt,j,i) = d_help_b(j,i)
         end do

         d_help_c(0,j,i) = d_help_t(ktmax,j,i)

         do kc=0, kcmax-1
            d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
                                +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
         end do

      else   ! n_cts(j,i) == 1, temperate ice layer

         d_help_t(0,j,i) = d_help_b(j,i)

         do kt=0, ktmax-1
            d_help_t(kt+1,j,i) = d_help_t(kt,j,i) &
                                +0.5_dp*(cvxy2(kt+1)+cvxy2(kt))
         end do

         d_help_c(0,j,i) = d_help_t(ktmax,j,i)

         do kc=0, kcmax-1
            d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
                                +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
         end do

      end if

   else   ! maske(j,i) == 1, 2 or 3 away from the grounding line

      do kt=0, ktmax
         d_help_t(kt,j,i) = 0.0_dp
      end do

      do kc=0, kcmax
         d_help_c(kc,j,i) = 0.0_dp
      end do

   end if

end do
end do

!-------- Computation of vx_c/t (defined at (i+1/2,j,kc/t)) --------

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

   do kt=0, ktmax
      vx_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j,i+1)) &
                     *dzs_dxi(j,i)
   end do

   do kc=0, kcmax
      vx_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j,i+1)) &
                     *dzs_dxi(j,i)
   end do

end do
end do

!-------- Computation of vy_c/t (defined at (i,j+1/2,kc/t)) --------

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

   do kt=0, ktmax
      vy_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j+1,i)) &
                     *dzs_deta(j,i)
   end do

   do kc=0, kcmax
      vy_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j+1,i)) &
                     *dzs_deta(j,i)
   end do

end do
end do

!-------- Computation of the surface velocities vx_s_g and vy_s_g
!                                              (defined at (i,j)) --------

do i=0, imax
do j=0, jmax
   vx_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_dxi_g(j,i)
   vy_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_deta_g(j,i)
end do
end do

!-------- Limitation of computed vx_c/t, vy_c/t, vx_s_g, vy_s_g
!         to the interval [-VH_MAX, VH_MAX] --------

vh_max = vh_max/year_sec

do i=0, imax
do j=0, jmax
   vx_s_g(j,i) = max(vx_s_g(j,i), -vh_max)
   vx_s_g(j,i) = min(vx_s_g(j,i),  vh_max)
   vy_s_g(j,i) = max(vy_s_g(j,i), -vh_max)
   vy_s_g(j,i) = min(vy_s_g(j,i),  vh_max)
   do kt=0, ktmax
      vx_t(kt,j,i) = max(vx_t(kt,j,i), -vh_max)
      vx_t(kt,j,i) = min(vx_t(kt,j,i),  vh_max)
      vy_t(kt,j,i) = max(vy_t(kt,j,i), -vh_max)
      vy_t(kt,j,i) = min(vy_t(kt,j,i),  vh_max)
   end do
   do kc=0, kcmax
      vx_c(kc,j,i) = max(vx_c(kc,j,i), -vh_max)
      vx_c(kc,j,i) = min(vx_c(kc,j,i),  vh_max)
      vy_c(kc,j,i) = max(vy_c(kc,j,i), -vh_max)
      vy_c(kc,j,i) = min(vy_c(kc,j,i),  vh_max)
   end do
end do
end do

!-------- Computation of h_diff
!         (defined on the grid points (i,j)) --------

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

   if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
                     ! grounded ice, or floating ice at the grounding line

!  ------ Abbreviations

      do kt=0, ktmax
         cqxy0(kt) = h_t(j,i)*d_help_t(kt,j,i)*dzeta_t
      end do

      do kc=0, kcmax
         cqxy1(kc) = aqxy1(kc)*h_c(j,i)*d_help_c(kc,j,i)
      end do

!  ------ h_diff

      h_diff(j,i) = 0.0_dp

      if (n_cts(j,i) == 1) then

         do kt=0, ktmax-1
            h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy0(kt+1)+cqxy0(kt))
         end do

      end if

      do kc=0, kcmax-1
         h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy1(kc+1)+cqxy1(kc))
      end do

!  ------ Limitation of h_diff

      if (h_diff(j,i) < hd_min) h_diff(j,i) = 0.0_dp
      if (h_diff(j,i) > hd_max) h_diff(j,i) = hd_max

   else   ! maske(j,i) == 1, 2 or 3 away from the grounding line

      h_diff(j,i) = 0.0_dp

   end if

end do
end do

!-------- Computation of the horizontal volume flux
!                            and the depth-averaged velocity --------

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

   qx(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j,i+1))*dzs_dxi(j,i)

   if ( (maske(j,i)==0).or.(maske(j,i+1)==0) ) then   ! at least one neighbour
                                                      ! point is grounded ice

      vx_m(j,i)       = qx(j,i) &
                        / ( 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) )

      vx_m(j,i) = max(vx_m(j,i), -vh_max)   ! Limitation
      vx_m(j,i) = min(vx_m(j,i),  vh_max)   ! to the interval [-VH_MAX, VH_MAX]

      ratio_sl_x(j,i) = abs(vx_t(0,j,i)) / max(abs(vx_c(kcmax,j,i)), epsi)

   else 

      vx_m(j,i)       = 0.0_dp
      ratio_sl_x(j,i) = 0.0_dp

   end if

end do
end do

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

   qy(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j+1,i))*dzs_deta(j,i)

   if ( (maske(j,i)==0).or.(maske(j+1,i)==0) ) then   ! at least one neighbour
                                                      ! point is grounded ice

      vy_m(j,i)       = qy(j,i) &
                        / ( 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) )

      vy_m(j,i) = max(vy_m(j,i), -vh_max)   ! Limitation
      vy_m(j,i) = min(vy_m(j,i),  vh_max)   ! to the interval [-VH_MAX, VH_MAX]

      ratio_sl_y(j,i) = abs(vy_t(0,j,i)) / max(abs(vy_c(kcmax,j,i)), epsi)

   else 

      vy_m(j,i)       = 0.0_dp
      ratio_sl_y(j,i) = 0.0_dp

   end if

end do
end do

!-------- Detection of shelfy stream points --------

flag_shelfy_stream_x = .false.
flag_shelfy_stream_y = .false.
flag_shelfy_stream   = .false.

#if (DYNAMICS==0 || DYNAMICS==1)

ratio_sl_threshold = 1.11e+11_dp   ! dummy value

#elif (DYNAMICS==2)

#if ( defined(RATIO_SL_THRESH) )
ratio_sl_threshold = ratio_sl_thresh
#else
ratio_sl_threshold = 0.5_dp   ! default value
#endif

do i=0, imax-1
do j=0, jmax
   if (ratio_sl_x(j,i) > ratio_sl_threshold) flag_shelfy_stream_x(j,i) = .true.
end do
end do

do i=0, imax
do j=0, jmax-1
   if (ratio_sl_y(j,i) > ratio_sl_threshold) flag_shelfy_stream_y(j,i) = .true.
end do
end do

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

   if ( (maske(j,i) == 0) &   ! grounded ice
        .and. &
        (     flag_shelfy_stream_x(j,i-1)   &   ! at least
          .or.flag_shelfy_stream_x(j,i)     &   ! one neighbour
          .or.flag_shelfy_stream_y(j-1,i)   &   ! on the staggered grid
          .or.flag_shelfy_stream_y(j,i)   ) &   ! is a shelfy stream point
      ) then

      flag_shelfy_stream(j,i) = .true.

   end if

end do
end do

#else
stop ' >>> calc_vxy_sia: DYNAMICS must be 0, 1 or 2!'
#endif

!-------- Initialisation of the variable q_gl_g
!         (volume flux across the grounding line, to be
!         computed in the routine calc_vxy_ssa
!         if ice shelves are present)

q_gl_g = 0.0_dp

end subroutine calc_vxy_sia

!-------------------------------------------------------------------------------
!> Computation of the horizontal velocity vx, vy, the horizontal volume flux
!> qx, qy etc. for static ice.
!<------------------------------------------------------------------------------
subroutine calc_vxy_static()

implicit none

c_slide = 0.0_dp
p_weert = 0.0_dp
q_weert = 0.0_dp
p_b_w   = 0.0_dp

d_help_b = 0.0_dp
c_drag   = 0.0_dp

vx_b   = 0.0_dp
vy_b   = 0.0_dp
vx_b_g = 0.0_dp
vy_b_g = 0.0_dp

txz_c = 0.0_dp
txz_t = 0.0_dp

tyz_c = 0.0_dp
tyz_t = 0.0_dp

sigma_c = 0.0_dp
sigma_t = 0.0_dp

flui_ave_sia = 0.0_dp
de_ssa       = 0.0_dp
vis_int_g    = 0.0_dp

d_help_t = 0.0_dp
d_help_c = 0.0_dp

vx_c = 0.0_dp
vy_c = 0.0_dp

vx_t = 0.0_dp
vy_t = 0.0_dp

vx_s_g = 0.0_dp
vy_s_g = 0.0_dp

h_diff = 0.0_dp

qx = 0.0_dp
qy = 0.0_dp

vx_m = 0.0_dp
vy_m = 0.0_dp

ratio_sl_x = 0.0_dp
ratio_sl_y = 0.0_dp

flag_shelfy_stream_x = .false.
flag_shelfy_stream_y = .false.
flag_shelfy_stream   = .false.

q_gl_g = 0.0_dp

end subroutine calc_vxy_static

!-------------------------------------------------------------------------------
!> Computation of the horizontal velocity vx, vy, the horizontal volume flux
!! qx, qy and the flux across the grounding line q_gl_g in the shallow shelf
!! approximation (SSA) or the shelfy stream approximation (SStA).
!<------------------------------------------------------------------------------
subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)

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) :: z_sl, dxi, deta, dzeta_c, dzeta_t

integer(i4b) :: i, j, kc, kt, m
integer(i4b) :: iter_ssa
real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_sia, vy_m_sia
real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_prev, vy_m_prev
real(dp) :: rel_ssa
real(dp) :: dxi_inv, deta_inv
real(dp) :: vh_max
real(dp) :: ratio_sl_threshold, ratio_help
real(dp), dimension(0:JMAX,0:IMAX) :: weigh_ssta_sia_x, weigh_ssta_sia_y
real(dp) :: qx_gl_g, qy_gl_g

#if (MARGIN==3 || DYNAMICS==2)

!-------- Parameters for the relaxation scheme --------

#if (defined(N_ITER_SSA))
   iter_ssa = max(n_iter_ssa, 1)   ! number of iterations
#else
   iter_ssa = 2                    ! default value
#endif

#if (defined(RELAX_FACT_SSA))
   rel_ssa = relax_fact_ssa   ! relaxation factor
#else
   rel_ssa = 0.7_dp           ! default value
#endif

write(6,'(10x,a,i0,a,f6.3)') 'calc_vxy_ssa: iter_ssa = ', iter_ssa, &
                             ', rel_ssa =' , rel_ssa

!-------- 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) ) &   ! ocean point
      ) &
   flag_calving_front_1(j,i) = .true.

   if ( (maske(j,i)==2_i2b) &   ! ocean
        .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) &   ! ocean
        .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) &   ! ocean
        .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) &   ! ocean
        .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) &   ! ocean
        .and. (maske(j,i-1)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

end do

!-------- Save mean (depth-averaged) horizontal velocities from SIA --------

vx_m_sia = vx_m
vy_m_sia = vy_m

!-------- First iteration --------

m=1

!  ------ Depth-integrated viscosity vis_int_g

vis_int_g = 1.0e+15_dp*(h_c+h_t)
            ! viscosity of 10^15 Pa s assumed

!  ------ Horizontal velocity vx_m, vy_m

call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn, m, iter_ssa)

!-------- Further iterations --------

do m=2, iter_ssa

!  ------ Save velocities from previous iteration

   vx_m_prev = vx_m
   vy_m_prev = vy_m

!  ------ Depth-integrated viscosity vis_int_g

   call calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)

!  ------ Horizontal velocity vx_m, vy_m

   call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn, m, iter_ssa)

!  ------ Relaxation scheme

   vx_m = rel_ssa*vx_m + (1.0_dp-rel_ssa)*vx_m_prev
   vy_m = rel_ssa*vy_m + (1.0_dp-rel_ssa)*vy_m_prev

end do

!  ------ Depth-integrated viscosity vis_int_g

call calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)

!-------- Limitation of computed vx_m and vy_m to the interval
!                                            [-VH_MAX, VH_MAX] --------

vh_max = vh_max/year_sec

vx_m = max(vx_m, -vh_max)
vx_m = min(vx_m,  vh_max)
vy_m = max(vy_m, -vh_max)
vy_m = min(vy_m,  vh_max)

!-------- 3-D velocities, basal velocities and volume flux --------

#if (DYNAMICS==0 || DYNAMICS==1)

ratio_sl_threshold = 1.11e+11_dp   ! dummy value
ratio_help         = 0.0_dp

#elif (DYNAMICS==2)

#if ( defined(RATIO_SL_THRESH) )
ratio_sl_threshold = ratio_sl_thresh
#else
ratio_sl_threshold = 0.5_dp   ! default value
#endif

ratio_help = 1.0_dp/(1.0_dp-ratio_sl_threshold)

#else
stop ' >>> calc_vxy_ssa: DYNAMICS must be 0, 1 or 2!'
#endif

weigh_ssta_sia_x = 0.0_dp
weigh_ssta_sia_y = 0.0_dp

!  ------ x-component

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

   if (flag_shelfy_stream_x(j,i)) then
                               ! shelfy stream

      weigh_ssta_sia_x(j,i) = (ratio_sl_x(j,i)-ratio_sl_threshold)*ratio_help

      weigh_ssta_sia_x(j,i) = max(min(weigh_ssta_sia_x(j,i), 1.0_dp), 0.0_dp)
                              ! constrain to interval [0,1]

#if (SSTA_SIA_WEIGH_FCT==0)

      ! stick to the linear function set above

#elif (SSTA_SIA_WEIGH_FCT==1)

      weigh_ssta_sia_x(j,i) = weigh_ssta_sia_x(j,i)*weigh_ssta_sia_x(j,i) &
                                   *(3.0_dp-2.0_dp*weigh_ssta_sia_x(j,i))
                              ! make transition smooth (cubic function)

#elif (SSTA_SIA_WEIGH_FCT==2)

      weigh_ssta_sia_x(j,i) = weigh_ssta_sia_x(j,i)*weigh_ssta_sia_x(j,i) &
                                                   *weigh_ssta_sia_x(j,i) &
                                   *(10.0_dp + weigh_ssta_sia_x(j,i) &
                                       *(-15.0_dp+6.0_dp*weigh_ssta_sia_x(j,i)))
                              ! make transition even smoother (quintic function)

#else
      stop ' >>> calc_vxy_ssa: SSTA_SIA_WEIGH_FCT must be 0, 1 or 2!'
#endif

      do kt=0, ktmax
         vx_t(kt,j,i) = weigh_ssta_sia_x(j,i)*vx_m(j,i) &
                        + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_t(kt,j,i)
      end do

      do kc=0, kcmax
         vx_c(kc,j,i) = weigh_ssta_sia_x(j,i)*vx_m(j,i) &
                        + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_c(kc,j,i)
      end do

      vx_b(j,i) = vx_t(0,j,i)

      vx_m(j,i) = weigh_ssta_sia_x(j,i)*vx_m(j,i) &
                  + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_m_sia(j,i)

      qx(j,i)   = vx_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))

   else if ( (maske(j,i)==3_i2b).or.(maske(j,i+1)==3_i2b) ) then
                               ! at least one neighbour point is floating ice
      do kt=0, ktmax
         vx_t(kt,j,i) = vx_m(j,i)
      end do

      do kc=0, kcmax
         vx_c(kc,j,i) = vx_m(j,i)
      end do

      vx_b(j,i) = vx_m(j,i)

      qx(j,i)   = vx_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))

   end if

end do
end do

!  ------ y-component

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

   if (flag_shelfy_stream_y(j,i)) then
                               ! shelfy stream

      weigh_ssta_sia_y(j,i) = (ratio_sl_y(j,i)-ratio_sl_threshold)*ratio_help

      weigh_ssta_sia_y(j,i) = max(min(weigh_ssta_sia_y(j,i), 1.0_dp), 0.0_dp)
                              ! constrain to interval [0,1]

#if (SSTA_SIA_WEIGH_FCT==0)

      ! stick to the linear function set above

#elif (SSTA_SIA_WEIGH_FCT==1)

      weigh_ssta_sia_y(j,i) = weigh_ssta_sia_y(j,i)*weigh_ssta_sia_y(j,i) &
                                   *(3.0_dp-2.0_dp*weigh_ssta_sia_y(j,i))
                              ! make transition smooth (cubic function)

#elif (SSTA_SIA_WEIGH_FCT==2)

      weigh_ssta_sia_y(j,i) = weigh_ssta_sia_y(j,i)*weigh_ssta_sia_y(j,i) &
                                                   *weigh_ssta_sia_y(j,i) &
                                   *(10.0_dp + weigh_ssta_sia_y(j,i) &
                                       *(-15.0_dp+6.0_dp*weigh_ssta_sia_y(j,i)))
                              ! make transition even smoother (quintic function)

#else
      stop ' >>> calc_vxy_ssa: SSTA_SIA_WEIGH_FCT must be 0, 1 or 2!'
#endif

      do kt=0, ktmax
         vy_t(kt,j,i) = weigh_ssta_sia_y(j,i)*vy_m(j,i) &
                        + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_t(kt,j,i)
      end do

      do kc=0, kcmax
         vy_c(kc,j,i) = weigh_ssta_sia_y(j,i)*vy_m(j,i) &
                        + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_c(kc,j,i)
      end do

      vy_b(j,i) = vy_t(0,j,i)

      vy_m(j,i) = weigh_ssta_sia_y(j,i)*vy_m(j,i) &
                  + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_m_sia(j,i)

      qy(j,i)   = vy_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i))

   else if ( (maske(j,i)==3_i2b).or.(maske(j+1,i)==3_i2b) ) then
                               ! at least one neighbour point is floating ice
      do kt=0, ktmax
         vy_t(kt,j,i) = vy_m(j,i)
      end do

      do kc=0, kcmax
         vy_c(kc,j,i) = vy_m(j,i)
      end do

      vy_b(j,i) = vy_m(j,i)

      qy(j,i)   = vy_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i))

   end if

end do
end do

!-------- Surface and basal velocities vx_s_g vy_s_g, vx_b_g vy_b_g
!                                                (defined at (i,j)) --------

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

   if (flag_shelfy_stream(j,i)) then   ! shelfy stream

      vx_s_g(j,i) = 0.5_dp*(vx_c(kcmax,j,i-1)+vx_c(kcmax,j,i))
      vx_b_g(j,i) = 0.5_dp*(vx_b(      j,i-1)+vx_b(      j,i))

      vy_s_g(j,i) = 0.5_dp*(vy_c(kcmax,j-1,i)+vy_c(kcmax,j,i))
      vy_b_g(j,i) = 0.5_dp*(vy_b(      j-1,i)+vy_b(      j,i))

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

      vx_s_g(j,i) = 0.5_dp*(vx_m(j,i-1)+vx_m(j,i))
      vx_b_g(j,i) = vx_s_g(j,i)

      vy_s_g(j,i) = 0.5_dp*(vy_m(j-1,i)+vy_m(j,i))
      vy_b_g(j,i) = vy_s_g(j,i)

   end if

end do
end do

!-------- Computation of the flux across the grounding line q_gl_g

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

   if ( flag_grounding_line_1(j,i) ) then   ! grounding line

      qx_gl_g = 0.5_dp*(qx(j,i)+qx(j,i-1))
      qy_gl_g = 0.5_dp*(qy(j,i)+qy(j-1,i))

      q_gl_g(j,i) = sqrt(qx_gl_g*qx_gl_g+qy_gl_g*qy_gl_g)

   end if

end do
end do

#else

stop ' >>> calc_vxy_ssa: Only to be called for MARGIN==3 or DYNAMICS==2!'

#endif

end subroutine calc_vxy_ssa

!-------------------------------------------------------------------------------
!> Solution of the system of linear equations for the horizontal velocities
!! vx_m, vy_m in the shallow shelf approximation.
!<------------------------------------------------------------------------------
subroutine calc_vxy_ssa_matrix(z_sl,dxi, deta, ii, jj, nn, m, iter_ssa)

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)
integer(i4b), intent(in) :: m, iter_ssa
real(dp), intent(in) :: dxi, deta, z_sl

integer(i4b) :: i, j, k, n
integer(i4b) :: i1, j1
real(dp) :: inv_dxi, inv_deta, inv_dxi_deta, inv_dxi2, inv_deta2
real(dp) :: factor_rhs_1, factor_rhs_2, rhosw_rho_ratio
real(dp) :: H_mid, zl_mid
real(dp), dimension(0:JMAX,0:IMAX) :: vis_int_sgxy, beta_drag
character(len=256) :: ch_solver_set_option

#if (MARGIN==3 || DYNAMICS==2)

lis_integer :: ierr
lis_integer :: nc, nr
lis_integer :: iter
lis_matrix  :: lgs_a
lis_vector  :: lgs_b, lgs_x
lis_solver  :: solver

lis_integer, parameter                 :: nmax   =  2*(imax+1)*(jmax+1)
lis_integer, parameter                 :: n_sprs = 20*(imax+1)*(jmax+1)
lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
lis_scalar,  allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value

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

inv_dxi      = 1.0_dp/dxi
inv_deta     = 1.0_dp/deta
inv_dxi_deta = 1.0_dp/(dxi*deta)
inv_dxi2     = 1.0_dp/(dxi*dxi)
inv_deta2    = 1.0_dp/(deta*deta)

rhosw_rho_ratio = rho_sw/rho
factor_rhs_1 = 0.5_dp*rho*g
factor_rhs_2 = 0.5_dp*rho*g*(rho_sw-rho)/rho_sw

!-------- Depth-integrated viscosity on the staggered grid
!                                       [at (i+1/2,j+1/2)] --------

vis_int_sgxy = 0.0_dp   ! initialisation

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

   k=0

   if ((maske(j,i)==0_i2b).or.(maske(j,i)==3_i2b)) then
      k = k+1                              ! floating or grounded ice
      vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j,i)
   end if

   if ((maske(j,i+1)==0_i2b).or.(maske(j,i+1)==3_i2b)) then
      k = k+1                                  ! floating or grounded ice
      vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j,i+1)
   end if

   if ((maske(j+1,i)==0_i2b).or.(maske(j+1,i)==3_i2b)) then
      k = k+1                                  ! floating or grounded ice
      vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j+1,i)
   end if

   if ((maske(j+1,i+1)==0_i2b).or.(maske(j+1,i+1)==3_i2b)) then
      k = k+1                                      ! floating or grounded ice
      vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j+1,i+1)
   end if

   if (k>0) vis_int_sgxy(j,i) = vis_int_sgxy(j,i)/real(k,dp)

end do
end do

!-------- Basal drag parameter (for shelfy stream) --------

beta_drag = 0.0_dp   ! initialisation

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

   if (flag_shelfy_stream(j,i)) then

      if (m==1) then   ! first iteration

         beta_drag(j,i) = c_drag(j,i) &
                          / sqrt(vx_b_g(j,i)**2  &
                                +vy_b_g(j,i)**2) &
                                          **(1.0_dp-p_weert_inv(j,i))
                          ! computed with the SIA basal velocities
      else

         beta_drag(j,i) = c_drag(j,i) &
                          / sqrt((0.5_dp*(vx_m(j,i)+vx_m(j,i-1)))**2  &
                                +(0.5_dp*(vy_m(j,i)+vy_m(j-1,i)))**2) &
                                          **(1.0_dp-p_weert_inv(j,i))
                          ! computed with the SSA basal velocities
                          ! from the previous iteration
      end if

   end if

end do
end do

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

allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
allocate(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 n=1, nmax-1, 2

   i = ii((n+1)/2)
   j = jj((n+1)/2)

!  ------ Equations for vx_m (at (i+1/2,j))

   nr = n   ! row counter

   if ( (i /= imax).and.(j /= 0).and.(j /= jmax) ) then
      ! inner point on the staggered grid in x-direction

      h_mid  = 0.50_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))
      zl_mid = 0.50_dp*(zl(j,i)+zl(j,i+1))

      if ( &
           ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==3_i2b) ) &
           .or. &
           ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j,i+1) &
             .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
           .or. &
           ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j,i+1) &
             .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
         ) then
           ! both neighbours are floating ice,
           ! or one neighbour is floating ice and the other is grounded ice
           ! (grounding line)
           ! and floating conditions are satisfied;
           ! discretisation of the x-component of the PDE

         flag_shelfy_stream_x(j,i)=.false.
                                   ! make sure not to treat as shelfy stream

         nc = 2*nn(j-1,i)-1   ! smallest nc (column counter), for vx_m(j-1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j-1,i)     ! next nc (column counter), for vy_m(j-1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j-1,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i-1)-1   ! next nc (column counter), for vx_m(j,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j-1,i+1)   ! next nc (column counter), for vy_m(j-1,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)-1     ! next nc (column counter), for vx_m(j,i)
         if (nc /= nr) &      !                      (diagonal element)
            stop ' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!'
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -4.0_dp*inv_dxi2 &
                                 *(vis_int_g(j,i+1)+vis_int_g(j,i)) &
                          -inv_deta2 &
                                 *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)       ! next nc (column counter), for vy_m(j,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i+1)-1   ! next nc (column counter), for vx_m(j,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i+1)
         lgs_a_index(k) = nc

         nc = 2*nn(j,i+1)     ! next nc (column counter), for vy_m(j,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i)-1   ! largest nc (column counter), for vx_m(j+1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
         lgs_a_index(k) = nc

         lgs_b_value(nr) = factor_rhs_1 &
                           * (h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) &
                           * dzs_dxi(j,i)

         lgs_x_value(nr) = vx_m(j,i)
 
      else if (flag_shelfy_stream_x(j,i)) then
           ! shelfy stream (as determined by routine calc_vxy_sia)

         nc = 2*nn(j-1,i)-1   ! smallest nc (column counter), for vx_m(j-1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j-1,i)     ! next nc (column counter), for vy_m(j-1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j-1,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i-1)-1   ! next nc (column counter), for vx_m(j,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j-1,i+1)   ! next nc (column counter), for vy_m(j-1,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)-1     ! next nc (column counter), for vx_m(j,i)
         if (nc /= nr) &      !                      (diagonal element)
            stop ' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!'
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -4.0_dp*inv_dxi2 &
                                 *(vis_int_g(j,i+1)+vis_int_g(j,i)) &
                          -inv_deta2 &
                                 *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i)) &
                          -0.5_dp*(beta_drag(j,i+1)+beta_drag(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)       ! next nc (column counter), for vy_m(j,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i+1)-1   ! next nc (column counter), for vx_m(j,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i+1)
         lgs_a_index(k) = nc

         nc = 2*nn(j,i+1)     ! next nc (column counter), for vy_m(j,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i)-1   ! largest nc (column counter), for vx_m(j+1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
         lgs_a_index(k) = nc

         lgs_b_value(nr) = factor_rhs_1 &
                           * (h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) &
                           * dzs_dxi(j,i)

         lgs_x_value(nr) = vx_m(j,i)

      else if ( &
                ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j,i+1) &
                  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
                .or. &
                ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j,i+1) &
                  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
              ) then 
              ! one neighbour is floating ice and the other is grounded ice
              ! (grounding line),
              ! but floating conditions are not satisfied

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = vx_m(j,i)
         lgs_x_value(nr) = vx_m(j,i)

      else if ( &
                ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==1_i2b) ) &
                .or. &
                ( (maske(j,i)==1_i2b).and.(maske(j,i+1)==3_i2b) ) &
              ) then
              ! one neighbour is floating ice and the other is ice-free land;
              ! velocity assumed to be zero

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = 0.0_dp

         lgs_x_value(nr) = 0.0_dp

      else if ( &
                ( flag_calving_front_1(j,i).and.flag_calving_front_2(j,i+1) ) &
                .or. &
                ( flag_calving_front_2(j,i).and.flag_calving_front_1(j,i+1) ) &
              ) then
              ! one neighbour is floating ice and the other is ocean
              ! (calving front)

         if (maske(j,i)==3_i2b) then
            i1 = i     ! floating ice marker
         else   ! maske(j,i+1)==3_i2b
            i1 = i+1   ! floating ice marker
         end if

         if ( &
              ( (maske(j,i1-1)==3_i2b).or.(maske(j,i1+1)==3_i2b) ) &
              .or. &
              ( (maske(j,i1-1)==0_i2b).or.(maske(j,i1+1)==0_i2b) ) &
              .or. &
              ( (maske(j,i1-1)==1_i2b).or.(maske(j,i1+1)==1_i2b) ) &
            ) then
              ! discretisation of the x-component of the BC at the calving front

            nc = 2*nn(j-1,i1)  ! smallest nc (column counter), for vy_m(j-1,i1)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = -2.0_dp*inv_deta*vis_int_g(j,i1)
            lgs_a_index(k) = nc

            nc = 2*nn(j,i1-1)-1   ! next nc (column counter), for vx_m(j,i1-1)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = -4.0_dp*inv_dxi*vis_int_g(j,i1)
            lgs_a_index(k) = nc

            nc = 2*nn(j,i1)-1     ! next nc (column counter), for vx_m(j,i1)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = 4.0_dp*inv_dxi*vis_int_g(j,i1)
            lgs_a_index(k) = nc

            nc = 2*nn(j,i1)       ! largest nc (column counter), for vy_m(j,i1)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = 2.0_dp*inv_deta*vis_int_g(j,i1)
            lgs_a_index(k) = nc

            lgs_b_value(nr) = factor_rhs_2 &
                              *(h_c(j,i1)+h_t(j,i1))*(h_c(j,i1)+h_t(j,i1))

            lgs_x_value(nr) = vx_m(j,i)

         else   ! (maske(j,i1-1)==2_i2b).and.(maske(j,i1+1)==2_i2b);
                ! velocity assumed to be zero

            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = 1.0_dp   ! diagonal element only
            lgs_a_index(k) = nr

            lgs_b_value(nr) = 0.0_dp

            lgs_x_value(nr) = 0.0_dp

         end if

      else if ( (maske(j,i)==0_i2b).or.(maske(j,i+1)==0_i2b) ) then
           ! neither neighbour is floating ice, but at least one neighbour is
           ! grounded ice; velocity taken from the solution for grounded ice

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = vx_m(j,i)

         lgs_x_value(nr) = vx_m(j,i)

      else   ! neither neighbour is floating or grounded ice,
             ! velocity assumed to be zero

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = 0.0_dp

         lgs_x_value(nr) = 0.0_dp

      end if

   else   ! boundary condition, velocity assumed to be zero

      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
      lgs_a_value(k) = 1.0_dp   ! diagonal element only
      lgs_a_index(k) = nr

      lgs_b_value(nr) = 0.0_dp

      lgs_x_value(nr) = 0.0_dp

   end if

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

!  ------ Equations for vy_m (at (i,j+1/2))

   nr = n+1   ! row counter

   if ( (j /= jmax).and.(i /= 0).and.(i /= imax) ) then
      ! inner point on the staggered grid in y-direction

      h_mid  = 0.5_dp*(h_c(j+1,i)+h_t(j+1,i)+h_c(j,i)+h_t(j,i))
      zl_mid = 0.5_dp*(zl(j+1,i)+zl(j,i))
   
      if ( &
           ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==3_i2b) ) &
           .or. &
           ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j+1,i) &
             .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
           .or. &
           ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j+1,i) &
             .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
         ) then
           ! both neighbours are floating ice,
           ! or one neighbour is floating ice and the other is grounded ice
           ! (grounding line)
           ! and floating conditions are satisfied;
           ! discretisation of the x-component of the PDE

         flag_shelfy_stream_y(j,i)=.false.
                                   ! make sure not to treat as shelfy stream

         nc = 2*nn(j-1,i)     ! smallest nc (column counter), for vy_m(j-1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j,i-1)-1   ! next nc (column counter), for vx_m(j,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i-1))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i-1)     ! next nc (column counter), for vy_m(j,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)-1     ! next nc (column counter), for vx_m(j,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)       ! next nc (column counter), for vy_m(j,i)
         if (nc /= nr) &      !                      (diagonal element)
            stop ' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!'
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -4.0_dp*inv_deta2 &
                                 *(vis_int_g(j+1,i)+vis_int_g(j,i)) &
                          -inv_dxi2 &
                                 *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1))
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i-1)-1 ! next nc (column counter), for vx_m(j+1,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i+1)     ! next nc (column counter), for vy_m(j,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i)-1   ! next nc (column counter), for vx_m(j+1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i)     ! largest nc (column counter), for vy_m(j+1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j+1,i)
         lgs_a_index(k) = nc

         lgs_b_value(nr) = factor_rhs_1 &
                           * (h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) &
                           * dzs_deta(j,i)

         lgs_x_value(nr) = vy_m(j,i)
 
      else if (flag_shelfy_stream_y(j,i)) then
           ! shelfy stream (as determined by routine calc_vxy_sia)

         nc = 2*nn(j-1,i)     ! smallest nc (column counter), for vy_m(j-1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j,i-1)-1   ! next nc (column counter), for vx_m(j,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i-1))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i-1)     ! next nc (column counter), for vy_m(j,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)-1     ! next nc (column counter), for vx_m(j,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i)       ! next nc (column counter), for vy_m(j,i)
         if (nc /= nr) &      !                      (diagonal element)
            stop ' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!'
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -4.0_dp*inv_deta2 &
                                 *(vis_int_g(j+1,i)+vis_int_g(j,i)) &
                          -inv_dxi2 &
                                 *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1)) &
                          -0.5_dp*(beta_drag(j+1,i)+beta_drag(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i-1)-1 ! next nc (column counter), for vx_m(j+1,i-1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = -inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
         lgs_a_index(k) = nc

         nc = 2*nn(j,i+1)     ! next nc (column counter), for vy_m(j,i+1)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i)-1   ! next nc (column counter), for vx_m(j+1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = inv_dxi_deta &
                                 *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i))
         lgs_a_index(k) = nc

         nc = 2*nn(j+1,i)     ! largest nc (column counter), for vy_m(j+1,i)
         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j+1,i)
         lgs_a_index(k) = nc

         lgs_b_value(nr) = factor_rhs_1 &
                           * (h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) &
                           * dzs_deta(j,i)

         lgs_x_value(nr) = vy_m(j,i)

      else if ( &
                ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j+1,i) &
                  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
                .or. &
                ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j+1,i) &
                  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
              ) then
              ! one neighbour is floating ice and the other is grounded ice
              ! (grounding line),
              ! but floating conditions are not satisfied

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = vy_m(j,i)
         lgs_x_value(nr) = vy_m(j,i)

      else if ( &
                ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==1_i2b) ) &
                .or. &
                ( (maske(j,i)==1_i2b).and.(maske(j+1,i)==3_i2b) ) &
              ) then
           ! one neighbour is floating ice and the other is ice-free land;
           ! velocity assumed to be zero

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = 0.0_dp

         lgs_x_value(nr) = 0.0_dp

      else if ( &
                ( flag_calving_front_1(j,i).and.flag_calving_front_2(j+1,i) ) &
                .or. &
                ( flag_calving_front_2(j,i).and.flag_calving_front_1(j+1,i) ) &
              ) then
              ! one neighbour is floating ice and the other is ocean
              ! (calving front)

         if (maske(j,i)==3_i2b) then
            j1 = j     ! floating ice marker
         else   ! maske(j+1,i)==3_i2b
            j1 = j+1   ! floating ice marker
         end if

         if ( &
              ( (maske(j1-1,i)==3_i2b).or.(maske(j1+1,i)==3_i2b) ) &
              .or. &
              ( (maske(j1-1,i)==0_i2b).or.(maske(j1+1,i)==0_i2b) ) &
              .or. &
              ( (maske(j1-1,i)==1_i2b).or.(maske(j1+1,i)==1_i2b) ) &
            ) then
              ! discretisation of the y-component of the BC at the calving front

            nc = 2*nn(j1-1,i)  ! smallest nc (column counter), for vy_m(j1-1,i)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = -4.0_dp*inv_deta*vis_int_g(j1,i)
            lgs_a_index(k) = nc

            nc = 2*nn(j1,i-1)-1   ! next nc (column counter), for vx_m(j1,i-1)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = -2.0_dp*inv_dxi*vis_int_g(j1,i)
            lgs_a_index(k) = nc

            nc = 2*nn(j1,i)-1     ! next nc (column counter), for vx_m(j1,i)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = 2.0_dp*inv_dxi*vis_int_g(j1,i)
            lgs_a_index(k) = nc

            nc = 2*nn(j1,i)       ! largest nc (column counter), for vy_m(j1,i)
            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = 4.0_dp*inv_deta*vis_int_g(j1,i)
            lgs_a_index(k) = nc

            lgs_b_value(nr) = factor_rhs_2 &
                              *(h_c(j1,i)+h_t(j1,i))*(h_c(j1,i)+h_t(j1,i))

            lgs_x_value(nr) = vy_m(j,i)

         else   ! (maske(j1-1,i)==2_i2b).and.(maske(j1+1,i)==2_i2b);
                ! velocity assumed to be zero

            k  = k+1
            ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
            lgs_a_value(k) = 1.0_dp   ! diagonal element only
            lgs_a_index(k) = nr

            lgs_b_value(nr) = 0.0_dp

            lgs_x_value(nr) = 0.0_dp

         end if

      else if ( (maske(j,i)==0_i2b).or.(maske(j+1,i)==0_i2b) ) then
           ! neither neighbour is floating ice, but at least one neighbour is
           ! grounded ice; velocity taken from the solution for grounded ice

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = vy_m(j,i)

         lgs_x_value(nr) = vy_m(j,i)

      else   ! neither neighbour is floating or grounded ice,
             ! velocity assumed to be zero

         k  = k+1
         ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
         lgs_a_value(k) = 1.0_dp   ! diagonal element only
         lgs_a_index(k) = nr

         lgs_b_value(nr) = 0.0_dp
         lgs_x_value(nr) = 0.0_dp

      end if

   else   ! boundary condition, velocity assumed to be zero

      k  = k+1
      ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
      lgs_a_value(k) = 1.0_dp   ! diagonal element only
      lgs_a_index(k) = nr

      lgs_b_value(nr) = 0.0_dp
      lgs_x_value(nr) = 0.0_dp

   end if

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

end do

!-------- 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 of the system of linear equations with Lis --------

call lis_solver_create(solver, ierr)


#if (defined(LIS_OPTS))
    ch_solver_set_option = trim(lis_opts)
#else
    ch_solver_set_option = '-i bicgsafe -p jacobi '// &
                           '-maxiter 100 -tol 1.0e-4 -initx_zeros false'
#endif
                          
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)

if (iter_ssa == 1) then
   write(6,'(10x,a,i0)') 'calc_vxy_ssa_matrix: LIS iter = ', iter
else if (m == 1) then
   write(6,'(10x,a,i0)', advance='no') 'calc_vxy_ssa_matrix: LIS iter = ', iter
else if (m == iter_ssa) then
   write(6,'(a,i0)') ', ', iter
else
   write(6,'(a,i0)', advance='no') ', ', iter
end if

!!! call lis_solver_get_time(solver,solver_time,ierr)
!!! print *, 'calc_vxy_ssa_matrix: time (s) = ', solver_time

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 n=1, nmax-1, 2

   i = ii((n+1)/2)
   j = jj((n+1)/2)

   nr = n
   vx_m(j,i) = lgs_x_value(nr)

   nr = n+1
   vy_m(j,i) = lgs_x_value(nr)

end do

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

#else

stop ' >>> calc_vxy_ssa_matrix: Only to be called for MARGIN==3 or DYNAMICS==2!'

#endif

end subroutine calc_vxy_ssa_matrix

!-------------------------------------------------------------------------------
!> Computation of the depth-integrated viscosity vis_int_g in the
!! shallow shelf approximation.
!<------------------------------------------------------------------------------
subroutine calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)

use ice_material_properties_m, only : viscosity

implicit none

real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t

integer(i4b) :: i, j, kc, kt
real(dp) :: dxi_inv, deta_inv
real(dp) :: dvx_dxi, dvx_deta, dvy_dxi, dvy_deta
real(dp) :: aqxy1(0:kcmax)
real(dp) :: cvis0(0:ktmax), cvis1(0:kcmax)

#if (MARGIN==3 || DYNAMICS==2)

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

dxi_inv  = 1.0_dp/dxi
deta_inv = 1.0_dp/deta

do kc=0, kcmax
   if (flag_aa_nonzero) then
      aqxy1(kc) = aa/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
   else
      aqxy1(kc) = dzeta_c
   end if
end do

!-------- Computation of the depth-integrated viscosity --------

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

   if ((maske(j,i)==0_i2b).and.(.not.flag_shelfy_stream(j,i))) then
                                                   ! grounded ice, but
                                                   ! not shelfy stream
      de_ssa(j,i) = 0.0_dp   ! dummy value

      vis_int_g(j,i) = (h_c(j,i)+h_t(j,i)) / flui_ave_sia(j,i)

   else if ((maske(j,i)==1_i2b).or.(maske(j,i)==2_i2b)) then
                                                   ! ice-free land or ocean
      de_ssa(j,i) = 0.0_dp   ! dummy value

      vis_int_g(j,i) = 0.0_dp   ! dummy value

   else   ! (maske(j,i)==3_i2b).or.(flag_shelfy_stream(j,i)),
          ! floating ice or shelfy stream; 
          ! must not be at the margin of the computational domain

!  ------ Effective strain rate

      dvx_dxi  = (vx_m(j,i)-vx_m(j,i-1))*dxi_inv 
      dvy_deta = (vy_m(j,i)-vy_m(j-1,i))*deta_inv

      dvx_deta = 0.25_dp*deta_inv &
                 *(vx_m(j+1,i)+vx_m(j+1,i-1)-vx_m(j-1,i)-vx_m(j-1,i-1))
      dvy_dxi  = 0.25_dp*dxi_inv &
                 *(vy_m(j,i+1)+vy_m(j-1,i+1)-vy_m(j,i-1)-vy_m(j-1,i-1))

      de_ssa(j,i) = sqrt( dvx_dxi*dvx_dxi &
                          + dvy_deta*dvy_deta &
                          + dvx_dxi*dvy_deta &
                          + 0.25_dp*(dvx_deta+dvy_dxi)*(dvx_deta+dvy_dxi) )

!  ------ Term abbreviations

#if (DYNAMICS==2)
   if (.not.flag_shelfy_stream(j,i)) then
#endif

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
                        enh_c(kc,j,i), 0_i2b)
      end do
      ! Ice shelves (floating ice) are assumed to consist of cold ice only

#if (DYNAMICS==2)
   else   ! flag_shelfy_stream(j,i) == .true.

#if (CALCMOD==-1 || CALCMOD==0)

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
                        enh_c(kc,j,i), 0_i2b)
      end do

#elif (CALCMOD==1)

      do kt=0, ktmax
         cvis0(kt) = dzeta_t*h_t(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_t_m(kt,j,i), temp_t_m(kt,j,i), omega_t(kt,j,i), &
                        enh_t(kt,j,i), 1_i2b)
      end do

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
                        enh_c(kc,j,i), 0_i2b)
      end do

#elif (CALCMOD==2 || CALCMOD==3)

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), omega_c(kc,j,i), &
                        enh_c(kc,j,i), 2_i2b)
      end do

#else
   stop ' >>> calc_vis_ssa: CALCMOD must be -1, 0, 1, 2 or 3!'
#endif

   end if
#endif

!  ------ Depth-integrated viscosity

      vis_int_g(j,i) = 0.0_dp

#if (CALCMOD==1)
      do kt=0, ktmax-1
         vis_int_g(j,i) = vis_int_g(j,i)+0.5_dp*(cvis0(kt+1)+cvis0(kt))
      end do
#endif

      do kc=0, kcmax-1
         vis_int_g(j,i) = vis_int_g(j,i)+0.5_dp*(cvis1(kc+1)+cvis1(kc))
      end do

   end if
  
end do
end do

#else

stop ' >>> calc_vis_ssa: Only to be called for MARGIN==3 or DYNAMICS==2!'

#endif

end subroutine calc_vis_ssa

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

end module calc_vxy_m
!