topograd_1.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t o p o g r a d _ 1
!
!> @file
!!
!! Calculation of topography gradients on the staggered grid and on the grid
!! points (the latter by second-order discretization).
!! Length rescaling with the corresponding components of the metric tensor is
!! already included.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 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 second-order discretization).
!! Length rescaling with the corresponding components of the metric tensor is
!! already included.
!<------------------------------------------------------------------------------
subroutine topograd_1(dxi, deta, n_switch)

use sico_types
use sico_variables

implicit none

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

integer(i4b) :: i, j
real(dp)     :: dxi_inv, deta_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

!-------- 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_1: 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=1, imax-1
do j=0, jmax
   dzs_dxi_g(j,i)  = (zs_help(j,i+1)-zs_help(j,i-1))*0.5_dp*dxi_inv &
                     *insq_g11_g(j,i)
   dzm_dxi_g(j,i)  = (zm_help(j,i+1)-zm_help(j,i-1))*0.5_dp*dxi_inv &
                     *insq_g11_g(j,i)
   dzb_dxi_g(j,i)  = (zb_help(j,i+1)-zb_help(j,i-1))*0.5_dp*dxi_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,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=1, jmax-1
   dzs_deta_g(j,i)  = (zs_help(j+1,i)-zs_help(j-1,i)) &
                      *0.5_dp*deta_inv &
                      *insq_g22_g(j,i)
   dzm_deta_g(j,i)  = (zm_help(j+1,i)-zm_help(j-1,i)) &
                      *0.5_dp*deta_inv &
                      *insq_g22_g(j,i)
   dzb_deta_g(j,i)  = (zb_help(j+1,i)-zb_help(j-1,i)) &
                      *0.5_dp*deta_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(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_1
!