calc_vz_ssa.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v z _ s s a
!
!> @file
!!
!! Computation of the vertical velocity vz 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 vertical velocity vz in the shallow shelf approximation.
!<------------------------------------------------------------------------------
subroutine calc_vz_ssa(z_sl, dxi, deta, dzeta_c)

use sico_types
use sico_variables

implicit none

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

integer(i4b) :: i, j, kt, kc
real(dp), dimension(0:JMAX,0:IMAX) :: vz_sl
real(dp), dimension(0:KCMAX) :: zeta_c_sgz, eaz_c_sgz, eaz_c_quotient_sgz
real(dp) :: dxi_inv, deta_inv
real(dp) :: dvx_dxi, dvy_deta

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

!-------- Abbreviations --------

do kc=0, kcmax-1
   zeta_c_sgz(kc) = (kc+0.5_dp)*dzeta_c
   eaz_c_sgz(kc)  = exp(deform*zeta_c_sgz(kc))
   eaz_c_quotient_sgz(kc) = (eaz_c_sgz(kc)-1.0_dp)/(ea-1.0_dp)
end do

!-------- Computation of vz --------

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

   if (maske(j,i)==3_i2b) then   ! floating ice

!  ------ Derivatives of the horizontal velocity

      dvx_dxi =  insq_g11_g(j,i)*insq_g22_g(j,i) &
                 *(vx_m(j,i)*sq_g22_sgx(j,i)-vx_m(j,i-1)*sq_g22_sgx(j,i-1)) &
                 *dxi_inv

      dvy_deta = insq_g11_g(j,i)*insq_g22_g(j,i) &
                 *(vy_m(j,i)*sq_g11_sgy(j,i)-vy_m(j-1,i)*sq_g11_sgy(j-1,i)) &
                 *deta_inv

!  ------ Basal velocity vz_b

      vz_b(j,i) = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1))*dzb_dxi_g(j,i) &
                  +0.5_dp*(vy_m(j,i)+vy_m(j-1,i))*dzb_deta_g(j,i) &
                  +dzb_dtau(j,i)-q_b_tot(j,i)
                  ! kinematic boundary condition at the ice base

!  ------ Velocity at sea level vz_sl

      vz_sl(j,i) = vz_b(j,i) - (z_sl-zb(j,i))*(dvx_dxi+dvy_deta)

!  ------ Surface velocity vz_s

      vz_s(j,i) = vz_sl(j,i) - (zs(j,i)-z_sl)*(dvx_dxi+dvy_deta)

!  ------ 3-D velocity vz_c and vz_t

      do kc=0, kcmax-1
         vz_c(kc,j,i) = vz_sl(j,i) &
                        -(zm(j,i)+eaz_c_quotient_sgz(kc)*h_c(j,i)-z_sl) &
                         *(dvx_dxi+dvy_deta)
      end do

      if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
                        ! cold ice base, temperate ice base

         do kt=0, ktmax-1
            vz_t(kt,j,i) = vz_b(j,i)
         end do

      else   ! n_cts(j,i) == 1, temperate ice layer

         do kt=0, ktmax-1
            vz_t(kt,j,i) = vz_sl(j,i) &
                           -(zb(j,i) &
                             +0.5_dp*(zeta_t(kt)+zeta_t(kt+1))*h_t(j,i) &
                             -z_sl) &
                            *(dvx_dxi+dvy_deta)
         end do

      end if

   end if

end do
end do

end subroutine calc_vz_ssa
!