calc_elra.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ e l r a
!
!> @file
!!
!! Computation of the isostatic steady-state displacement of the lithosphere
!! (wss) for the ELRA model.
!!
!! @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 isostatic steady-state displacement of the lithosphere
!! (wss) for the ELRA model.
!<------------------------------------------------------------------------------
subroutine calc_elra(z_sl, dxi, deta)

use sico_types
use sico_variables
use sico_vars

implicit none

real(dp), intent(in) :: z_sl, dxi, deta

integer(i4b) :: i, j, ir, jr, il, jl, n
integer(i4b) :: ir_max, jr_max, min_imax_jmax
integer(i4b) :: il_begin, il_end, jl_begin, jl_end
real(dp)                              :: rho_g, rho_sw_g, rho_a_g_inv
real(dp)                              :: dxi_inv, deta_inv
real(dp)                              :: kei_r_incr_inv
real(dp), dimension(0:JMAX,0:IMAX)    :: l_r, l_r_inv, fac_wss
real(dp)                              :: ra_max
real(dp), allocatable, dimension(:,:) :: f_0

real(dp), parameter :: r_infl = 8.0_dp
                       ! Radius of non-locality influence of the elastic
                       ! lithosphere plate (in units of l_r)

!-------- Initialisations --------

rho_g       = rho*g
rho_sw_g    = rho_sw*g
rho_a_g_inv = 1.0_dp/(rho_a*g)

dxi_inv   = 1.0_dp/dxi
deta_inv  = 1.0_dp/deta

kei_r_incr_inv = 1.0_dp/kei_r_incr

do i=0, imax
do j=0, jmax
   l_r(j,i)     = (flex_rig_lith(j,i)*rho_a_g_inv)**0.25_dp
   l_r_inv(j,i) = 1.0_dp/l_r(j,i)
   fac_wss(j,i) = l_r(j,i)*l_r(j,i)/(2.0_dp*pi*flex_rig_lith(j,i))
end do
end do

ra_max  = r_infl*maxval(l_r)   ! Radius of non-locality influence (in m)

ir_max = floor(ra_max*dxi_inv)
jr_max = floor(ra_max*deta_inv)
         ! Radius of non-locality influence (in grid points)

min_imax_jmax = min(imax, jmax)
ir_max        = min(ir_max, min_imax_jmax)
jr_max        = min(jr_max, min_imax_jmax)
                ! ir_max and jr_max constrained by size of domain

il_begin = -ir_max
il_end   = ir_max+imax
jl_begin = -jr_max
jl_end   = jr_max+jmax

!-------- Ice/water load --------

allocate(f_0(jl_begin:jl_end, il_begin:il_end))

do il=il_begin, il_end
do jl=jl_begin, jl_end

   i = min(max(il, 0), imax)
   j = min(max(jl, 0), jmax)

   if (maske(j,i)==0) then
      f_0(jl,il) = rho_g * area(j,i) * (h_c(j,i) + h_t(j,i))
   else if (maske(j,i)==1) then
      f_0(jl,il) = 0.0_dp
   else   ! (maske(j,i)>=2)
      f_0(jl,il) = rho_sw_g * area(j,i) * z_sl
                   ! Water load relative to the present sea-level stand (0 m)
                   ! -> can be positive or negative
   end if

end do
end do

!-------- Computation of the steady-state displacement of the lithosphere
!         (wss, positive downward) --------

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

   wss(j,i) = 0.0_dp

   do ir=-ir_max, ir_max
   do jr=-jr_max, jr_max

#if GRID==0
      n = nint( dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
#elif GRID==1
      n = nint( sq_g11_g(j,i)*dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
              ! The correction with g11 only is OK because g11=g22 for this case
#elif GRID==2
      n = nint( dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
              ! The distortion due to g11 and g22 is already accounted for in
              ! the variable dist_dxdy(jr,ir)
#endif

      wss(j,i) = wss(j,i) - fac_wss(j,i)*f_0(jr+j,ir+i)*kei(n)

   end do
   end do

end do
end do

deallocate (f_0)

end subroutine calc_elra
!