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, n
integer(i4b) :: ios
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(23, iostat=ios, &
     file=inpath//'/nhem/'//zl0_file, &
     recl=16384, status='old')

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

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

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

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

do j=jmax, 0, -1
   read(23, fmt=*) (zl0(j,i), i=0,imax)
end do

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

do j=jmax, 0, -1
   read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
end do

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
      maske(j,i) = 1
      zs(j,i) = zl0(j,i)
      zb(j,i) = zl0(j,i)
      zl(j,i) = zl0(j,i)
   else   ! (maske(j,i) >= 2 )
      maske(j,i) = 2
#if ( MARGIN==1 || MARGIN==2 )
      zs(j,i) = zl0(j,i)
      zb(j,i) = zl0(j,i)
#elif ( MARGIN==3 )
      zs(j,i) = 0.0_dp   ! present-day
      zb(j,i) = 0.0_dp   ! sea level
#endif
      zl(j,i) = zl0(j,i)
   end if

   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
!