calc_vxy_ssa.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v x y _ s s a
!
!> @file
!!
!! 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.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> 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.
!<------------------------------------------------------------------------------
subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)

use sico_types
use sico_variables

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

integer(i4b) :: i, j, kc, kt, m
real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_prev, vy_m_prev
real(dp) :: dxi_inv, deta_inv
real(dp) :: vh_max
real(dp) :: qx_gl_g, qy_gl_g

integer(i4b), parameter :: iter_ssa = 2       ! number of iterations
real(dp),     parameter :: rel      = 0.7_dp   ! relaxation factor

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

flag_grounding_line = .false.
flag_calving_front  = .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) ) &   ! ice shelf point
      ) &
   flag_grounding_line(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(j,i) = .true.

end do
end do

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

m=1

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

do i=0, imax
do j=0, jmax
   vis_int_g(j,i) = 1.0e+15_dp*(h_c(j,i)+h_t(j,i))
                    ! viscosity of 10^15 Pa s assumed
end do
end do

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

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

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

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

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

!  ------ Relaxation scheme

   do i=0, imax
   do j=0, jmax
      vx_m(j,i) = rel*vx_m(j,i) + (1.0_dp-rel)*vx_m_prev(j,i)
      vy_m(j,i) = rel*vy_m(j,i) + (1.0_dp-rel)*vy_m_prev(j,i)
   end do
   end do

end do

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

vh_max = vh_max/year_sec

do i=0, imax
do j=0, jmax
   vx_m(j,i) = max(vx_m(j,i), -vh_max)
   vx_m(j,i) = min(vx_m(j,i),  vh_max)
   vy_m(j,i) = max(vy_m(j,i), -vh_max)
   vy_m(j,i) = min(vy_m(j,i),  vh_max)
end do
end do

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

!  ------ x-component

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

   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 ( (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 (maske(j,i)==3_i2b) then

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

end subroutine calc_vxy_ssa
!