calc_vxy_ssa_matrix.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v x y _ s s a _ m a t r i x
!
!> @file
!!
!! Solution of the system of linear equations for the horizontal velocities
!! vx_m, vy_m in the shallow shelf approximation.
!!
!! @section Copyright
!!
!! Copyright 2009-2016 Ralf Greve, 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> 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)

use sico_types
use sico_variables
use sico_vars

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
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

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)

ch_solver_set_option = '-i bicgstabl -p ilu '// &
                       '-maxiter 1000 -tol 1.0e-11 -initx_zeros false'
call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)

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

! call lis_solver_get_iters(solver, iter, ierr)
! write(6,'(11x,a,i5)') 'calc_vxy_ssa_matrix: 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 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)

end subroutine calc_vxy_ssa_matrix
!