calc_vz_sia.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ v z _ s i a
!
!> @file
!!
!! Computation of the vertical velocity vz 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 vertical velocity vz in the shallow ice approximation.
!<------------------------------------------------------------------------------
subroutine calc_vz_sia(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)     :: avz3(0:kcmax)
real(dp)     :: cvz00, cvz0(0:ktmax), cvz1(0:ktmax), cvz2(0:ktmax), &
                cvz3(0:kcmax), cvz4(0:kcmax), cvz5(0:kcmax)
real(dp)     :: dxi_inv, deta_inv

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

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

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

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

   if (maske(j,i) == 0) then   ! glaciated land

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

      cvz00 = 0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1))*dzb_dxi_g(j,i) &
              +0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i))*dzb_deta_g(j,i) &
              +dzb_dtau(j,i)-q_b_tot(j,i)

      kt=0
      cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
              *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
                 -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
                +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
                 -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
              *dzeta_t
      cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
              *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
               -vx_t(kt,j,i)  -vx_t(kt,j,i-1))*0.5_dp
      cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
              *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
               -vy_t(kt,j,i)  -vy_t(kt,j-1,i))*0.5_dp

      do kt=1, ktmax-1
         cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
              *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
                 -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
                +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
                 -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
              *dzeta_t
         cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
              *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
               -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.25_dp
         cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
              *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
               -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.25_dp
      end do

      kt=ktmax
      cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
              *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
                 -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
                +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
                 -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
              *dzeta_t
      cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
              *(vx_t(kt,j,i)  +vx_t(kt,j,i-1) &
               -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.5_dp
      cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
              *(vy_t(kt,j,i)  +vy_t(kt,j-1,i) &
               -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.5_dp

      kc=0
      cvz3(kc) = avz3(kc)*h_c(j,i) &
                 *insq_g11_g(j,i)*insq_g22_g(j,i) &
              *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
                 -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
                +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
                 -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
      cvz4(kc) = (dzm_dxi_g(j,i) &
              +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
              *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
               -vx_c(kc,j,i)  -vx_c(kc,j,i-1))*0.5_dp
      cvz5(kc) = (dzm_deta_g(j,i) &
              +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
              *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
               -vy_c(kc,j,i)  -vy_c(kc,j-1,i))*0.5_dp

      do kc=1, kcmax-1
         cvz3(kc) = avz3(kc)*h_c(j,i) &
                    *insq_g11_g(j,i)*insq_g22_g(j,i) &
              *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
                 -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
                +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
                 -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
         cvz4(kc) = (dzm_dxi_g(j,i) &
              +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
              *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
               -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.25_dp
         cvz5(kc) = (dzm_deta_g(j,i) &
              +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
              *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
               -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.25_dp
      end do

      kc=kcmax
      cvz3(kc) = avz3(kc)*h_c(j,i) &
                 *insq_g11_g(j,i)*insq_g22_g(j,i) &
              *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
                 -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
                +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
                 -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
      cvz4(kc) = (dzm_dxi_g(j,i) &
              +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
              *(vx_c(kc,j,i)  +vx_c(kc,j,i-1) &
               -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.5_dp
      cvz5(kc) = (dzm_deta_g(j,i) &
              +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
              *(vy_c(kc,j,i)  +vy_c(kc,j-1,i) &
               -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.5_dp

!-------- Computation of vz_b --------

      vz_b(j,i) = cvz00

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

      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

         vz_m(j,i) = vz_b(j,i)

         vz_c(0,j,i) = vz_m(j,i) &
                       +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))

         do kc=1, kcmax-1
            vz_c(kc,j,i) = vz_c(kc-1,j,i) &
                           -cvz3(kc)+cvz4(kc)+cvz5(kc)
         end do

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

         vz_t(0,j,i) = vz_b(j,i) &
                       +0.5_dp*(-cvz0(0)+cvz1(0)+cvz2(0))

         do kt=1, ktmax-1
            vz_t(kt,j,i) = vz_t(kt-1,j,i) &
                           -cvz0(kt)+cvz1(kt)+cvz2(kt)
         end do

         vz_m(j,i) = vz_t(ktmax-1,j,i) &
                     +0.5_dp*(-cvz0(ktmax)+cvz1(ktmax)+cvz2(ktmax))

         vz_c(0,j,i) = vz_m(j,i) &
                       +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))

         do kc=1, kcmax-1
            vz_c(kc,j,i) = vz_c(kc-1,j,i) &
                           -cvz3(kc)+cvz4(kc)+cvz5(kc)
         end do

      end if

!-------- Computation of vz_s --------

      kc=kcmax

      vz_s(j,i) = vz_c(kc-1,j,i) &
                  +0.5_dp*(-cvz3(kc)+cvz4(kc)+cvz5(kc))

   else   ! maske(j,i) == (1 or 2)

      vz_b(j,i) = 0.0_dp

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

      vz_m(j,i) = 0.0_dp

      do kc=0, kcmax-1
         vz_c(kc,j,i) = 0.0_dp
      end do

      vz_s(j,i) = 0.0_dp

   end if

end do
end do

end subroutine calc_vz_sia
!