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-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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> 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, dzeta_t, ii, jj, nn)

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)
real(dp), intent(in) :: z_sl, dxi, deta, dzeta_c, dzeta_t

integer(i4b) :: i, j, kc, kt, m
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) :: dxi_inv, deta_inv
real(dp) :: vh_max
real(dp) :: ratio_sl_threshold, ratio_help, weigh_ssa_sia
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_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

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, m)

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

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

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

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

#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 either 0, 1 or 2!'
#endif

!  ------ x-component

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

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

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

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

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

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

      vx_m(j,i) = weigh_ssa_sia*vx_m(j,i) &
                  + (1.0_dp-weigh_ssa_sia)*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_ssa_sia = (ratio_sl_y(j,i)-ratio_sl_threshold)*ratio_help

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

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

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

      vy_m(j,i) = weigh_ssa_sia*vy_m(j,i) &
                  + (1.0_dp-weigh_ssa_sia)*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

end subroutine calc_vxy_ssa
!