7 !! Computation of the isostatic steady-state displacement of the lithosphere
8 !! (wss) for the ELRA model.
12 !! Copyright 2009-2013 Ralf Greve, Sascha Knell
16 !! This file is part of SICOPOLIS.
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 !> Computation of the isostatic steady-state displacement of the lithosphere
35 !! (wss) for the ELRA model.
36 !<------------------------------------------------------------------------------
44 real(dp),
intent(in) :: z_sl, dxi, deta
46 integer(i4b) :: i, j, ir, jr, il, jl, n
47 integer(i4b) :: ir_max, jr_max, min_imax_jmax
48 integer(i4b) :: il_begin, il_end, jl_begin, jl_end
49 real(dp) :: l_r, l_r_inv, rho_g, rho_sw_g, fac_wss
51 real(dp) :: dxi_inv, deta_inv
52 real(dp) :: kei_r_incr_inv
53 real(dp),
allocatable,
dimension(:,:) :: f_0
55 real(dp),
parameter :: r_infl = 8.0_dp
61 l_r = (flex_rig_l/(rho_a*g))**0.25_dp
69 fac_wss = l_r**2.0_dp/(2.0_dp*pi*flex_rig_l)
72 deta_inv = 1.0_dp/deta
74 kei_r_incr_inv = 1.0_dp/kei_r_incr
76 ir_max = floor(ra_max*dxi_inv)
77 jr_max = floor(ra_max*deta_inv)
80 min_imax_jmax = min(imax, jmax)
81 ir_max = min(ir_max, min_imax_jmax)
82 jr_max = min(jr_max, min_imax_jmax)
92 allocate(f_0(jl_begin:jl_end, il_begin:il_end))
94 do il=il_begin, il_end
95 do jl=jl_begin, jl_end
97 i = min(max(il, 0), imax)
98 j = min(max(jl, 0), jmax)
100 if (maske(j,i)==0)
then
101 f_0(jl,il) = rho_g * area(j,i) * (h_c(j,i) + h_t(j,i))
102 else if (maske(j,i)==1)
then
105 f_0(jl,il) = rho_sw_g * area(j,i) * z_sl
121 do ir=-ir_max, ir_max
122 do jr=-jr_max, jr_max
125 n = nint( dist_dxdy(jr,ir)*l_r_inv*kei_r_incr_inv )
127 n = nint( sq_g11_g(j,i)*dist_dxdy(jr,ir)*l_r_inv*kei_r_incr_inv )
130 n = nint( dist_dxdy(jr,ir)*l_r_inv*kei_r_incr_inv )
135 wss(j,i) = wss(j,i) - fac_wss*f_0(jr+j,ir+i)*kei(n)