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
real(dp)     :: xi0, eta0
real(dp)     :: h_ice, freeboard_ratio
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//'/asf/'//zs_present_file, &
     recl=16384, status='old')

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

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

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

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

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

open(24, iostat=ios, &
     file=inpath//'/asf/'//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

freeboard_ratio = (rho_sw-rho)/rho_sw

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

   if (maske(j,i) <= 1 ) then

      zb(j,i) = zl(j,i)   ! ensure consistency

   else if (maske(j,i) == 2 ) then

#if ( MARGIN==1 || MARGIN==2 )
      zs(j,i) = zl(j,i)   ! ensure
      zb(j,i) = zl(j,i)   ! consistency
#elif ( MARGIN==3 )
      zs(j,i) = 0.0_dp    ! present-day
      zb(j,i) = 0.0_dp    ! sea level
#endif

   else if (maske(j,i) == 3 ) then

#if ( MARGIN==1 || ( MARGIN==2 && MARINE_ICE_FORMATION==1 ) )
      maske(j,i) = 2     ! floating ice cut off
      zs(j,i) = zl(j,i)
      zb(j,i) = zl(j,i)
#elif ( MARGIN==2 && MARINE_ICE_FORMATION==2 )
      maske(j,i) = 0     ! floating ice becomes "underwater ice"
      h_ice   = zs(j,i)-zb(j,i)   ! ice thickness
      zs(j,i) = zl(j,i)+h_ice
      zb(j,i) = zl(j,i)
#elif ( MARGIN==3 )
      h_ice = zs(j,i)-zb(j,i)   ! ice thickness
      zs(j,i) = freeboard_ratio*h_ice   ! ensure properly
      zb(j,i) = zs(j,i)-h_ice           ! floating ice
#endif

   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) = 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
!