topography1.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t o p o g r a p h y 1
!
!> @file
!!
!! Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For present-day initial topography.
!!
!! @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 present-day initial topography.
!<------------------------------------------------------------------------------
subroutine topography1(dxi, deta)

use sico_types
use sico_variables

implicit none

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

integer(i4b) :: i, j, n
integer(i4b) :: ios, n_dummy
real(dp)     :: d_dummy
real(dp)     :: xi0, eta0
character    :: ch_dummy

character(len=  8) :: ch_imax
character(len=128) :: fmt4

write(ch_imax, fmt='(i8)') imax
write(fmt4,    fmt='(a)')  '('//trim(adjustl(ch_imax))//'(i1),i1)'

!-------- Read topography --------

open(21, iostat=ios, &
     file=inpath//'/nmars/'//zs_present_file, &
     recl=8192, status='old')

if (ios.ne.0) stop ' topography1: Error when opening the zs file!'

open(22, iostat=ios, &
     file=inpath//'/nmars/'//zl_present_file, &
     recl=8192, status='old')

if (ios.ne.0) stop ' topography1: Error when opening the zl file!'

open(23, iostat=ios, &
     file=inpath//'/nmars/'//zl0_file, &
     recl=8192, status='old')

if (ios.ne.0) stop ' topography1: Error when opening the zl0 file!'

open(24, iostat=ios, &
     file=inpath//'/nmars/'//mask_present_file, &
     recl=1024, status='old')

if (ios.ne.0) stop ' topography1: Error when opening the mask file!'

do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
do n=1, 6; read(24, fmt='(a)') ch_dummy; end do

do j=jmax, 0, -1
   read(21, fmt=*)          (zs(j,i),    i=0,imax)
   read(22, fmt=*)          (zl(j,i),    i=0,imax)
   read(23, fmt=*)          (zl0(j,i),   i=0,imax)
   read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
end do

close(21, status='keep')
close(22, status='keep')
close(23, status='keep')
close(24, status='keep')

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

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

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

   if (maske(j,i) <= 1 ) then
      zb(j,i) = zl(j,i)
   else   ! (maske(j,i)>=2)
      stop ' topography1: maske(j,i)>=2 not allowed for initial topography!'
   end if

   zs(j,i)  = zs(j,i) *1000.0_dp
   zb(j,i)  = zb(j,i) *1000.0_dp    ! km --> m
   zl(j,i)  = zl(j,i) *1000.0_dp
   zl0(j,i) = zl0(j,i)*1000.0_dp

   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) = zs(j,i)-zm(j,i)
   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 topography1
!