topography2.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t o p o g r a p h y 2
!
!> @file
!!
!! Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For ice-free initial topography with relaxed lithosphere.
!!
!! @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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For ice-free initial topography with relaxed lithosphere.
!<------------------------------------------------------------------------------
subroutine topography2(dxi, deta)

use sico_types
use sico_variables

implicit none

real(dp), intent(out) :: dxi, deta

integer(i4b) :: i, j
integer(i4b) :: ios
real(dp)     :: xi0, eta0

!-------- Set topography --------

zl0 = 0.0_dp

maske = 1

maske(0,:)    = 2
maske(jmax,:) = 2

maske(:,0)    = 2
maske(:,imax) = 2

!-------- Further stuff --------

#if (GRID==0 || GRID==1)

dxi  = dx *1000.0_dp   ! km -> m
deta = dx *1000.0_dp   ! km -> m

xi0  = x0 *1000.0_dp   ! km -> m
eta0 = y0 *1000.0_dp   ! km -> m

#elif GRID==2

stop ' topography2: GRID==2 not allowed for this application!'

#endif

do i=0, imax
do j=0, jmax

   zs(j,i) = zl0(j,i)
   zb(j,i) = zl0(j,i)
   zl(j,i) = zl0(j,i)

   xi(i)  = xi0  + real(i,dp)*dxi
   eta(j) = eta0 + real(j,dp)*deta

   call geo_coord(phi(j,i), lambda(j,i), xi(i), eta(j))

   zm(j,i) = zb(j,i)
   n_cts(j,i) = -1

   h_c(j,i) = 0.0_dp
   h_t(j,i) = 0.0_dp

   dzs_dtau(j,i)  = 0.0_dp
   dzm_dtau(j,i)  = 0.0_dp
   dzb_dtau(j,i)  = 0.0_dp
   dzl_dtau(j,i)  = 0.0_dp
   dh_c_dtau(j,i) = 0.0_dp
   dh_t_dtau(j,i) = 0.0_dp

end do
end do

!-------- Metric tensor, gradients of the topography --------

call metric()

#if TOPOGRAD==0
call topograd_1(dxi, deta, 1)
#elif TOPOGRAD==1
call topograd_2(dxi, deta, 1)
#endif

!-------- Corresponding area of grid points --------

do i=0, imax
do j=0, jmax
   area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
end do
end do

end subroutine topography2
!