make_zl0.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  m a k e _ z l 0
!
!> @file
!!
!! Generation of an isostatically relaxed lithosphere surface topography
!! for either the rigid lithosphere, the local lithosphere or the elastic
!! lithosphere model (depending on the setting of the parameter REBOUND).
!! This routine is not to be used regularly, and it is only executed if the
!! parameter EXEC_MAKE_ZL0 is defined in the header file.
!!
!! @section Copyright
!!
!! Copyright 2014-2016 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Generation of an isostatically relaxed lithosphere surface topography
!! for either the rigid lithosphere, the local lithosphere or the elastic
!! lithosphere model (depending on the setting of the parameter REBOUND).
!! This routine is not to be used regularly, and it is only executed if the
!! parameter EXEC_MAKE_ZL0 is defined in the header file.
!<------------------------------------------------------------------------------
subroutine make_zl0()

use sico_types
use sico_variables
use sico_vars

implicit none

integer(i4b)                       :: i, j, m, n, i_f, j_f, n_filter
integer(i4b)                       :: ios
real(dp), dimension(0:JMAX,0:IMAX) :: zl0_raw, zl0_smoothed
real(dp)                           :: rho_ratio
real(dp)                           :: filter_width, sigma_filter
real(dp)                           :: dist, weigh, sum_weigh
character(len=  8)                 :: ch_resolution
character(len= 64)                 :: ch_model
character(len=256)                 :: ch_filename

!-------- Checking of the grid --------

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

continue

#else

stop ' make_zl0: Routine works only for GRID 0 or 1!'

#endif

!-------- Computation of the raw (unsmoothed) topography --------

rho_ratio = rho/rho_a

#if (REBOUND==0)

zl0_raw = zl   ! rigid lithosphere

#elif (REBOUND==1)

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

   if (maske(j,i) == 0) then
      zl0_raw(j,i) = zl(j,i) + rho_ratio*(h_c(j,i)+h_t(j,i))
   else
      zl0_raw(j,i) = zl(j,i)
   end if            ! local lithosphere
   
end do
end do

#elif (REBOUND==2)

zl0_raw = zl + wss   ! elastic lithosphere

#else

stop ' make_zl0: REBOUND must be either 0, 1 or 2!'

#endif

!-------- Smoothing of the topography (Gaussian filter) --------

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

   filter_width = real(DX, dp)    ! half span of filtered area, in km
                                  ! (chosen as one grid spacing)

   sigma_filter = filter_width/real(DX, dp)   ! half span of filtered area,
                                              ! in grid points

   n_filter     = ceiling(2.0_dp*sigma_filter)
   n_filter     = max(n_filter, 5)

   sum_weigh = 0.0_dp

   do m=-n_filter, n_filter
   do n=-n_filter, n_filter

      i_f = i+m
      j_f = j+n

      if (i_f <    0) i_f =    0
      if (i_f > imax) i_f = imax

      if (j_f <    0) j_f =    0
      if (j_f > jmax) j_f = jmax

      dist      = sqrt(real(m,dp)**2+real(n,dp)**2)
      weigh     = exp(-(dist/sigma_filter)**2)
      sum_weigh = sum_weigh + weigh

      zl0_smoothed(j,i) = zl0_smoothed(j,i) + weigh*zl0_raw(j_f,i_f)

   end do
   end do

   zl0_smoothed(j,i) = zl0_smoothed(j,i)/sum_weigh

end do
end do

!-------- Writing on file --------

if ( abs(real(dx,dp)-nint(real(dx,dp))) < eps ) then
   write(ch_resolution, fmt='(i8)') nint(real(dx,dp))
else
   write(ch_resolution, fmt='(f8.2)') real(dx,dp)
end if

ch_resolution = adjustl(ch_resolution)

#if (REBOUND==0)
ch_model = 'rigid_lithosphere'
#elif (REBOUND==1)
ch_model = 'local_lithosphere'
#elif (REBOUND==2)
ch_model = 'elastic_lithosphere'
#endif

ch_filename = trim(ch_domain_short)//'_'//trim(ch_resolution)
ch_filename = trim(ch_filename)//'_zl0_'//trim(ch_model)//'.dat'

open(23, iostat=ios, &
         file=outpath//'/'//trim(ch_filename), &
         recl=8192, status='replace')
if (ios /= 0) stop ' make_zl0: Error when opening the new zl0 file!'

write(23,'(a)') &
'%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
write(23,'(a)') '% '//trim(ch_domain_long)//':'
write(23,'(a)') &
'% Relaxed lithosphere surface topography without ice load, in m AMSL.'
write(23,'(a)') &
'% Horizontal resolution '//trim(ch_resolution)//' km.'
write(23,'(a,i3,a,i3,a,i3,a,i3,a)') &
'% ', jmax+1, ' records [j = ', &
jmax, ' (-1) 0] with ', imax+1, &
' values [i = 0 (1) ', imax, '] in each record.'
write(23,'(a)') &
'%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'

do j=jmax, 0, -1
   do i=0, imax-1
      if ( zl0_smoothed(j,i) < -9999.0_dp ) zl0_smoothed(j,i) = -9999.0_dp
      write(23,'(f8.1)', advance='no') zl0_smoothed(j,i)
   end do
   i=imax
   if ( zl0_smoothed(j,i) < -9999.0_dp ) zl0_smoothed(j,i) = -9999.0_dp
   write(23,'(f8.1)') zl0_smoothed(j,i)
end do

close(23, status='keep')

end subroutine make_zl0
!