topograd_2.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t o p o g r a d _ 2
!
!> @file
!!
!! Calculation of topography gradients on the staggered grid and on the grid
!! points (the latter by fourth-order discretization).
!! Length rescaling with the corresponding components of the metric tensor is
!! already included.
!!
!! @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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Calculation of topography gradients on the staggered grid and on the grid
!! points (the latter by fourth-order discretization).
!! Length rescaling with the corresponding components of the metric tensor is
!! already included.
!<------------------------------------------------------------------------------
subroutine topograd_2(dxi, deta, n_switch)

use sico_types
use sico_variables
use sico_vars

implicit none

integer(i4b), intent(in) :: n_switch
real(dp),     intent(in) :: dxi, deta

integer(i4b) :: i, j
real(dp)     :: dxi_inv, deta_inv, dxi12_inv, deta12_inv
real(dp)     :: zs_help(0:jmax,0:imax), zm_help(0:jmax,0:imax), &
                zb_help(0:jmax,0:imax)

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

dxi12_inv  = 1.0_dp/(12.0_dp*dxi)
deta12_inv = 1.0_dp/(12.0_dp*deta)

!-------- Distinguish between old and new topography data --------

if (n_switch == 1) then

   do i=0, imax
   do j=0, jmax
      zs_help(j,i) = zs(j,i)
      zm_help(j,i) = zm(j,i)
      zb_help(j,i) = zb(j,i)
   end do
   end do

else if (n_switch == 2) then

   do i=0, imax
   do j=0, jmax
      zs_help(j,i) = zs_neu(j,i)
      zm_help(j,i) = zm_neu(j,i)
      zb_help(j,i) = zb_neu(j,i)
   end do
   end do

else
   stop ' topograd_2: Wrong value for n_switch!'
end if

!-------- Topography gradients on the staggered grid --------

!  ------ x-derivatives

do i=0, imax-1
do j=0, jmax
   dzs_dxi(j,i)  = (zs_help(j,i+1)-zs_help(j,i))*dxi_inv &
                   *insq_g11_sgx(j,i)
   dzm_dxi(j,i)  = (zm_help(j,i+1)-zm_help(j,i))*dxi_inv &
                   *insq_g11_sgx(j,i)
   dzb_dxi(j,i)  = (zb_help(j,i+1)-zb_help(j,i))*dxi_inv &
                   *insq_g11_sgx(j,i)
   dh_c_dxi(j,i) = dzs_dxi(j,i)-dzm_dxi(j,i)
   dh_t_dxi(j,i) = dzm_dxi(j,i)-dzb_dxi(j,i)
end do
end do

!  ------ y-derivatives

do i=0, imax
do j=0, jmax-1
   dzs_deta(j,i)  = (zs_help(j+1,i)-zs_help(j,i))*deta_inv &
                   *insq_g22_sgy(j,i)
   dzm_deta(j,i)  = (zm_help(j+1,i)-zm_help(j,i))*deta_inv &
                   *insq_g22_sgy(j,i)
   dzb_deta(j,i)  = (zb_help(j+1,i)-zb_help(j,i))*deta_inv &
                   *insq_g22_sgy(j,i)
   dh_c_deta(j,i) = dzs_deta(j,i)-dzm_deta(j,i)
   dh_t_deta(j,i) = dzm_deta(j,i)-dzb_deta(j,i)
end do
end do

!-------- Topography gradients on the grid points --------

!  ------ x-derivatives

do i=2, imax-2
do j=0, jmax
   dzs_dxi_g(j,i) &
      = (        -zs_help(j,i+2) + 8.0_dp*zs_help(j,i+1) &
          -8.0_dp*zs_help(j,i-1) +        zs_help(j,i-2) ) &
        *dxi12_inv &
        *insq_g11_g(j,i)
   dzm_dxi_g(j,i) &
      = (      -zm_help(j,i+2)   + 8.0_dp*zm_help(j,i+1) &
          -8.0_dp*zm_help(j,i-1) +        zm_help(j,i-2) ) &
        *dxi12_inv &
        *insq_g11_g(j,i)
   dzb_dxi_g(j,i) &
      = (     -zb_help(j,i+2)    + 8.0_dp*zb_help(j,i+1) &
          -8.0_dp*zb_help(j,i-1) +        zb_help(j,i-2) ) &
        *dxi12_inv &
        *insq_g11_g(j,i)
   dh_c_dxi_g(j,i) = dzs_dxi_g(j,i)-dzm_dxi_g(j,i)
   dh_t_dxi_g(j,i) = dzm_dxi_g(j,i)-dzb_dxi_g(j,i)
end do
end do

do j=0, jmax

   dzs_dxi_g(j,0)       = (zs_help(j,1)-zs_help(j,0))*dxi_inv &
                          *insq_g11_g(j,0)
   dzm_dxi_g(j,0)       = (zm_help(j,1)-zm_help(j,0))*dxi_inv &
                          *insq_g11_g(j,0)
   dzb_dxi_g(j,0)       = (zb_help(j,1)-zb_help(j,0))*dxi_inv &
                          *insq_g11_g(j,0)
   dh_c_dxi_g(j,0)      = dzs_dxi_g(j,0)-dzm_dxi_g(j,0)
   dh_t_dxi_g(j,0)      = dzm_dxi_g(j,0)-dzb_dxi_g(j,0)

   dzs_dxi_g(j,1)       = (zs_help(j,2)-zs_help(j,0)) &
                          *0.5_dp*dxi_inv &
                          *insq_g11_g(j,1)
   dzm_dxi_g(j,1)       = (zm_help(j,2)-zm_help(j,0)) &
                          *0.5_dp*dxi_inv &
                          *insq_g11_g(j,1)
   dzb_dxi_g(j,1)       = (zb_help(j,2)-zb_help(j,0)) &
                          *0.5_dp*dxi_inv &
                          *insq_g11_g(j,1)
   dh_c_dxi_g(j,1)      = dzs_dxi_g(j,1)-dzm_dxi_g(j,1)
   dh_t_dxi_g(j,1)      = dzm_dxi_g(j,1)-dzb_dxi_g(j,1)

   dzs_dxi_g(j,imax-1)  = (zs_help(j,imax)-zs_help(j,imax-2)) &
                          *0.5_dp*dxi_inv &
                          *insq_g11_g(j,imax-1)
   dzm_dxi_g(j,imax-1)  = (zm_help(j,imax)-zm_help(j,imax-2)) &
                          *0.5_dp*dxi_inv &
                          *insq_g11_g(j,imax-1)
   dzb_dxi_g(j,imax-1)  = (zb_help(j,imax)-zb_help(j,imax-2)) &
                          *0.5_dp*dxi_inv &
                          *insq_g11_g(j,imax-1)
   dh_c_dxi_g(j,imax-1) = dzs_dxi_g(j,imax-1) &
                          -dzm_dxi_g(j,imax-1)
   dh_t_dxi_g(j,imax-1) = dzm_dxi_g(j,imax-1) &
                          -dzb_dxi_g(j,imax-1)

   dzs_dxi_g(j,imax)    = (zs_help(j,imax)-zs_help(j,imax-1)) &
                          *dxi_inv &
                          *insq_g11_g(j,imax)
   dzm_dxi_g(j,imax)    = (zm_help(j,imax)-zm_help(j,imax-1)) &
                          *dxi_inv &
                          *insq_g11_g(j,imax)
   dzb_dxi_g(j,imax)    = (zb_help(j,imax)-zb_help(j,imax-1)) &
                          *dxi_inv &
                          *insq_g11_g(j,imax)
   dh_c_dxi_g(j,imax)   = dzs_dxi_g(j,imax)-dzm_dxi_g(j,imax)
   dh_t_dxi_g(j,imax)   = dzm_dxi_g(j,imax)-dzb_dxi_g(j,imax)

end do

!  ------ y-derivatives

do i=0, imax
do j=2, jmax-2
   dzs_deta_g(j,i) &
      = (     -zs_help(j+2,i)    + 8.0_dp*zs_help(j+1,i) &
          -8.0_dp*zs_help(j-1,i) +        zs_help(j-2,i) ) &
        *deta12_inv &
        *insq_g22_g(j,i)
   dzm_deta_g(j,i) &
      = (     -zm_help(j+2,i)    + 8.0_dp*zm_help(j+1,i) &
          -8.0_dp*zm_help(j-1,i) +        zm_help(j-2,i) ) &
        *deta12_inv &
        *insq_g22_g(j,i)
   dzb_deta_g(j,i) &
      = (     -zb_help(j+2,i)    + 8.0_dp*zb_help(j+1,i) &
          -8.0_dp*zb_help(j-1,i) +        zb_help(j-2,i) ) &
        *deta12_inv &
        *insq_g22_g(j,i)
   dh_c_deta_g(j,i) = dzs_deta_g(j,i)-dzm_deta_g(j,i)
   dh_t_deta_g(j,i) = dzm_deta_g(j,i)-dzb_deta_g(j,i)
end do
end do

do i=0, imax

   dzs_deta_g(0,i)       = (zs_help(1,i)-zs_help(0,i))*deta_inv &
                           *insq_g22_g(0,i)
   dzm_deta_g(0,i)       = (zm_help(1,i)-zm_help(0,i))*deta_inv &
                           *insq_g22_g(0,i)
   dzb_deta_g(0,i)       = (zb_help(1,i)-zb_help(0,i))*deta_inv &
                           *insq_g22_g(0,i)
   dh_c_deta_g(0,i)      = dzs_deta_g(0,i)-dzm_deta_g(0,i)
   dh_t_deta_g(0,i)      = dzm_deta_g(0,i)-dzb_deta_g(0,i)

   dzs_deta_g(1,i)       = (zs_help(2,i)-zs_help(0,i)) &
                           *0.5_dp*deta_inv &
                           *insq_g22_g(1,i)
   dzm_deta_g(1,i)       = (zm_help(2,i)-zm_help(0,i)) &
                           *0.5_dp*deta_inv &
                           *insq_g22_g(1,i)
   dzb_deta_g(1,i)       = (zb_help(2,i)-zb_help(0,i)) &
                           *0.5_dp*deta_inv &
                           *insq_g22_g(1,i)
   dh_c_deta_g(1,i)      = dzs_deta_g(1,i)-dzm_deta_g(1,i)
   dh_t_deta_g(1,i)      = dzm_deta_g(1,i)-dzb_deta_g(1,i)

   dzs_deta_g(jmax-1,i)  = (zs_help(jmax,i)-zs_help(jmax-2,i)) &
                           *0.5_dp*deta_inv &
                           *insq_g22_g(jmax-1,i)
   dzm_deta_g(jmax-1,i)  = (zm_help(jmax,i)-zm_help(jmax-2,i)) &
                           *0.5_dp*deta_inv &
                           *insq_g22_g(jmax-1,i)
   dzb_deta_g(jmax-1,i)  = (zb_help(jmax,i)-zb_help(jmax-2,i)) &
                           *0.5_dp*deta_inv &
                           *insq_g22_g(jmax-1,i)
   dh_c_deta_g(jmax-1,i) = dzs_deta_g(jmax-1,i) &
                           -dzm_deta_g(jmax-1,i)
   dh_t_deta_g(jmax-1,i) = dzm_deta_g(jmax-1,i) &
                           -dzb_deta_g(jmax-1,i)

   dzs_deta_g(jmax,i)    = (zs_help(jmax,i)-zs_help(jmax-1,i)) &
                           *deta_inv &
                           *insq_g22_g(jmax,i)
   dzm_deta_g(jmax,i)    = (zm_help(jmax,i)-zm_help(jmax-1,i)) &
                           *deta_inv &
                           *insq_g22_g(jmax,i)
   dzb_deta_g(jmax,i)    = (zb_help(jmax,i)-zb_help(jmax-1,i)) &
                           *deta_inv &
                           *insq_g22_g(jmax,i)
   dh_c_deta_g(jmax,i)   = dzs_deta_g(jmax,i)-dzm_deta_g(jmax,i)
   dh_t_deta_g(jmax,i)   = dzm_deta_g(jmax,i)-dzb_deta_g(jmax,i)

end do

end subroutine topograd_2
!