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-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 depth-integrated viscosity vis_int_g in the
!! shallow shelf approximation.
!<------------------------------------------------------------------------------
subroutine calc_vis_ssa(dxi, deta, dzeta_c)

use sico_types
use sico_variables

implicit none

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

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

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

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

do kc=0, kcmax
   aqxy1(kc) = deform/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
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_sf(j,i)) then   ! grounded ice

      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, flag_sf(j,i)==.true., floating ice; 
          ! 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 abbreviation

      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), .true.)

      end do

!  ------ Depth-integrated viscosity

      vis_int_g(j,i) = 0.0_dp

      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
!