calc_vis_ssa.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v i s _ s s a
!
!> @file
!!
!! Computation of the depth-integrated viscosity vis_int_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 depth-integrated viscosity vis_int_g in the
!! shallow shelf approximation.
!<------------------------------------------------------------------------------
subroutine calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)

use sico_types
use sico_variables
use sico_vars

implicit none

real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t

integer(i4b) :: i, j, kc, kt
real(dp) :: dxi_inv, deta_inv
real(dp) :: dvx_dxi, dvx_deta, dvy_dxi, dvy_deta
real(dp) :: aqxy1(0:kcmax)
real(dp) :: cvis0(0:ktmax), cvis1(0:kcmax)

!-------- Term abbreviations --------

dxi_inv  = 1.0_dp/dxi
deta_inv = 1.0_dp/deta

do kc=0, kcmax
   if (flag_aa_nonzero) then
      aqxy1(kc) = aa/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
   else
      aqxy1(kc) = dzeta_c
   end if
end do

!-------- Computation of the depth-integrated viscosity --------

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

   if ((maske(j,i)==0_i2b).and.(.not.flag_shelfy_stream(j,i))) then
                                                   ! grounded ice, but
                                                   ! not shelfy stream
      de_ssa(j,i) = 0.0_dp   ! dummy value

      vis_int_g(j,i) = (h_c(j,i)+h_t(j,i)) / flui_ave_sia(j,i)

   else if ((maske(j,i)==1_i2b).or.(maske(j,i)==2_i2b)) then
                                                   ! ice-free land or ocean
      de_ssa(j,i) = 0.0_dp   ! dummy value

      vis_int_g(j,i) = 0.0_dp   ! dummy value

   else   ! (maske(j,i)==3_i2b).or.(flag_shelfy_stream(j,i)),
          ! floating ice or shelfy stream; 
          ! must not be at the margin of the computational domain

!  ------ Effective strain rate

      dvx_dxi  = (vx_m(j,i)-vx_m(j,i-1))*dxi_inv 
      dvy_deta = (vy_m(j,i)-vy_m(j-1,i))*deta_inv

      dvx_deta = 0.25_dp*deta_inv &
                 *(vx_m(j+1,i)+vx_m(j+1,i-1)-vx_m(j-1,i)-vx_m(j-1,i-1))
      dvy_dxi  = 0.25_dp*dxi_inv &
                 *(vy_m(j,i+1)+vy_m(j-1,i+1)-vy_m(j,i-1)-vy_m(j-1,i-1))

      de_ssa(j,i) = sqrt( dvx_dxi*dvx_dxi &
                          + dvy_deta*dvy_deta &
                          + dvx_dxi*dvy_deta &
                          + 0.25_dp*(dvx_deta+dvy_dxi)*(dvx_deta+dvy_dxi) )

!  ------ Term abbreviations

#if (DYNAMICS==2)
   if (.not.flag_shelfy_stream(j,i)) then
#endif

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
                        enh_c(kc,j,i), 0_i2b)
      end do
      ! Ice shelves (floating ice) are assumed to consist of cold ice only

#if (DYNAMICS==2)
   else   ! flag_shelfy_stream(j,i) == .true.

#if (CALCMOD==-1 || CALCMOD==0)

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
                        enh_c(kc,j,i), 0_i2b)
      end do

#elif (CALCMOD==1)

      do kt=0, ktmax
         cvis0(kt) = dzeta_t*h_t(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_t_m(kt,j,i), temp_t_m(kt,j,i), omega_t(kt,j,i), &
                        enh_t(kt,j,i), 1_i2b)
      end do

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
                        enh_c(kc,j,i), 0_i2b)
      end do

#elif (CALCMOD==2 || CALCMOD==3)

      do kc=0, kcmax
         cvis1(kc) = aqxy1(kc)*h_c(j,i) &
                      *viscosity(de_ssa(j,i), &
                        temp_c(kc,j,i), temp_c_m(kc,j,i), omega_c(kc,j,i), &
                        enh_c(kc,j,i), 2_i2b)
      end do

#else
   stop ' calc_vis_ssa: CALCMOD must be either -1, 0, 1, 2 or 3!'
#endif

   end if
#endif

!  ------ Depth-integrated viscosity

      vis_int_g(j,i) = 0.0_dp

#if (CALCMOD==1)
      do kt=0, ktmax-1
         vis_int_g(j,i) = vis_int_g(j,i)+0.5_dp*(cvis0(kt+1)+cvis0(kt))
      end do
#endif

      do kc=0, kcmax-1
         vis_int_g(j,i) = vis_int_g(j,i)+0.5_dp*(cvis1(kc+1)+cvis1(kc))
      end do

   end if
  
end do
end do

end subroutine calc_vis_ssa
!