calc_vxy_b_sia.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v x y _ b _s i a
!
!> @file
!!
!! Computation of the basal horizontal velocity vx_b, vy_b 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 basal horizontal velocity vx_b, vy_b in the shallow ice
!! approximation.
!<------------------------------------------------------------------------------
subroutine calc_vxy_b_sia(time, z_sl)

use sico_types
use sico_variables
use sico_vars

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, ctau1, ctau1a
real(dp) :: temp_diff
real(dp) :: vh_max
real(dp) :: year_sec_inv
logical, dimension(0:JMAX,0:IMAX) :: sub_melt_flag

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

year_sec_inv = 1.0_dp/year_sec

#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

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

   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
!