calc_gia.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ g i a
!
!> @file
!!
!! Computation of the glacial isostatic adjustment of the lithosphere surface.
!!
!! @section Copyright
!!
!! Copyright 2009-2016 Ralf Greve, Sascha Knell
!!
!! @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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Computation of the glacial isostatic adjustment of the lithosphere surface.
!<------------------------------------------------------------------------------
subroutine calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)

use sico_types
use sico_variables
use sico_vars

implicit none

integer(i4b), intent(in) :: itercount, iter_wss
real(dp),     intent(in) :: time
real(dp),     intent(in) :: z_sl
real(dp),     intent(in) :: dtime, dxi, deta

integer(i4b) :: i, j
real(dp) :: tldt_inv(0:jmax,0:imax)
real(dp) :: dtime_inv
real(dp) :: rho_rhoa_ratio, rhosw_rhoa_ratio
real(dp) :: target_topo_time_lag

!-------- Term abbreviations --------

do i=0, imax
do j=0, jmax
   tldt_inv(j,i) = 1.0_dp/(time_lag_asth(j,i)+dtime)
end do
end do

rho_rhoa_ratio   = rho/rho_a
rhosw_rhoa_ratio = rho_sw/rho_a

dtime_inv = 1.0_dp/dtime

!-------- Computation of zl_neu and its time derivative --------

#if (REBOUND==0)

do i=0, imax
do j=0, jmax
   zl_neu(j,i)   = zl(j,i)
   dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
end do
end do

#elif (REBOUND==1)

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

   if (maske(j,i) <= 1_i2b) then   ! grounded ice or ice-free land

      zl_neu(j,i) = tldt_inv(j,i)*( time_lag_asth(j,i)*zl(j,i) &
                    + dtime*(zl0(j,i) &
                             -frac_llra*rho_rhoa_ratio*(h_c(j,i)+h_t(j,i))) )

   else   ! (maske(j,i) >= 2_i2b)

      zl_neu(j,i) = tldt_inv(j,i)*( time_lag_asth(j,i)*zl(j,i) &
                    + dtime*(zl0(j,i) &
                             -frac_llra*rhosw_rhoa_ratio*z_sl) )

                   ! Water load relative to the present sea-level stand (0 m)
                   ! -> can be positive or negative

   end if

   dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv

end do
end do

#elif (REBOUND==2)

if ((mod(itercount, iter_wss)==0).or.(itercount==1)) then
   write(6,*) ' -------- Computation of wss'
   call calc_elra(z_sl, dxi, deta)   ! Update of the steady-state displacement
                                     ! of the lithosphere (wss)
end if

do i=0, imax
do j=0, jmax
   zl_neu(j,i) = tldt_inv(j,i) &
                 * ( time_lag_asth(j,i)*zl(j,i) + dtime*(zl0(j,i)-wss(j,i)) )
   dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
end do
end do

#endif

!  ------ Adjustment due to prescribed target topography

#if (ZS_EVOL==2)

if (time >= time_target_topo_final) then

   zl_neu   = zl_target
   dzl_dtau = 0.0_dp   ! previously this was (zl_neu-zl)*dtime_inv

else if (time >= time_target_topo_init) then

   target_topo_time_lag = (time_target_topo_final-time)/3.0_dp

   zl_neu =   ( target_topo_time_lag*zl_neu + dtime*zl_target ) &
            / ( target_topo_time_lag        + dtime )

   dzl_dtau = (zl_neu-zl)*dtime_inv

end if

#endif

!======== Generation of an isostatically relaxed
!         lithosphere surface topography. Not to be used regularly! ========

#if (defined(EXEC_MAKE_ZL0))

call make_zl0()

write(6, fmt='(a)') ' '
write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
write(6, fmt='(a)') ' calc_gia: Routine make_zl0 successfully completed,      '
write(6, fmt='(a)') '           isostatically relaxed lithosphere topography  '
write(6, fmt='(a)') '           written on file                               '
write(6, fmt='(a)') '           (in directory specified by OUTPATH).          '
write(6, fmt='(a)') '           Execution of SICOPOLIS stopped.               '
write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
write(6, fmt='(a)') ' '

stop

#endif

end subroutine calc_gia
!