calc_vxy_sia.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v x y _ s i a
!
!> @file
!!
!! 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.
!!
!! @section Copyright
!!
!! Copyright 2009-2016 Ralf Greve
!!
!! @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 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 sico_types
use sico_variables
use sico_vars
use ice_material_quantities, only : ratefac_c, ratefac_t, ratefac_c_t

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