SICOPOLIS V5-dev  Revision 1105
calc_gia_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ g i a _ m
4 !
5 !> @file
6 !!
7 !! Computation of the glacial isostatic adjustment of the lithosphere surface.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 Ralf Greve, Sascha Knell
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of the glacial isostatic adjustment of the lithosphere surface.
34 !<------------------------------------------------------------------------------
35 
36 module calc_gia_m
37 
38  use sico_types_m
40  use sico_vars_m
41 
42  implicit none
43 
44  private
45  public :: calc_gia
46 
47 contains
48 
49 !-------------------------------------------------------------------------------
50 !> Main subroutine of calc_gia_m:
51 !! Computation of the glacial isostatic adjustment of the lithosphere surface.
52 !<------------------------------------------------------------------------------
53 subroutine calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)
54 
55 implicit none
56 
57 integer(i4b), intent(in) :: itercount, iter_wss
58 real(dp), intent(in) :: time
59 real(dp), intent(in) :: z_sl
60 real(dp), intent(in) :: dtime, dxi, deta
61 
62 integer(i4b) :: i, j
63 real(dp) :: tldt_inv(0:jmax,0:imax)
64 real(dp) :: dtime_inv
65 real(dp) :: rho_rhoa_ratio, rhosw_rhoa_ratio
66 
67 !-------- Term abbreviations --------
68 
69 do i=0, imax
70 do j=0, jmax
71  tldt_inv(j,i) = 1.0_dp/(time_lag_asth(j,i)+dtime)
72 end do
73 end do
74 
75 rho_rhoa_ratio = rho/rho_a
76 rhosw_rhoa_ratio = rho_sw/rho_a
77 
78 dtime_inv = 1.0_dp/dtime
79 
80 !-------- Computation of zl_neu and its time derivative --------
81 
82 #if (REBOUND==0)
83 
84 do i=0, imax
85 do j=0, jmax
86  zl_neu(j,i) = zl(j,i)
87  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
88 end do
89 end do
90 
91 #elif (REBOUND==1)
92 
93 do i=0, imax
94 do j=0, jmax
95 
96  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
97 
98  zl_neu(j,i) = tldt_inv(j,i)*( time_lag_asth(j,i)*zl(j,i) &
99  + dtime*(zl0(j,i) &
100  -frac_llra*rho_rhoa_ratio*(h_c(j,i)+h_t(j,i))) )
101 
102  else ! (maske(j,i) >= 2_i2b)
103 
104  zl_neu(j,i) = tldt_inv(j,i)*( time_lag_asth(j,i)*zl(j,i) &
105  + dtime*(zl0(j,i) &
106  -frac_llra*rhosw_rhoa_ratio*z_sl) )
107 
108  ! Water load relative to the present sea-level stand (0 m)
109  ! -> can be positive or negative
110 
111  end if
112 
113  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
114 
115 end do
116 end do
117 
118 #elif (REBOUND==2)
119 
120 if ((mod(itercount, iter_wss)==0).or.(itercount==1)) then
121  write(6, fmt='(10x,a)') 'Computation of wss'
122  call calc_elra(z_sl, dxi, deta) ! Update of the steady-state displacement
123  ! of the lithosphere (wss)
124 end if
125 
126 do i=0, imax
127 do j=0, jmax
128  zl_neu(j,i) = tldt_inv(j,i) &
129  * ( time_lag_asth(j,i)*zl(j,i) + dtime*(zl0(j,i)-wss(j,i)) )
130  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
131 end do
132 end do
133 
134 #endif
135 
136 ! ------ Adjustment due to prescribed target topography
137 
138 #if (THK_EVOL==2)
139 
140 if (time >= time_target_topo_final) then
141 
142  zl_neu = zl_target
143  dzl_dtau = 0.0_dp ! previously this was (zl_neu-zl)*dtime_inv
144 
145 else if (time >= time_target_topo_init) then
146 
148 
149  zl_neu = ( target_topo_tau*zl_neu + dtime*zl_target ) &
150  / ( target_topo_tau + dtime )
151 
152  dzl_dtau = (zl_neu-zl)*dtime_inv
153 
154 end if
155 
156 #elif (THK_EVOL==3)
157 
158 ! relaxation time target_topo_tau defined in routine sico_init
159 
160 zl_neu = ( target_topo_tau*zl_neu + dtime*zl_target ) &
161  / ( target_topo_tau + dtime )
162 
163 dzl_dtau = (zl_neu-zl)*dtime_inv
164 
165 #endif
166 
167 !======== Generation of an isostatically relaxed
168 ! lithosphere surface topography. Not to be used regularly! ========
169 
170 #if (defined(EXEC_MAKE_ZL0))
171 
172 call make_zl0()
173 
174 write(6, fmt='(a)') ' '
175 write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
176 write(6, fmt='(a)') ' >>> calc_gia: Routine make_zl0 successfully completed, '
177 write(6, fmt='(a)') ' isostatically relaxed lithosphere '
178 write(6, fmt='(a)') ' topography written on file '
179 write(6, fmt='(a)') ' (in directory specified by OUTPATH). '
180 write(6, fmt='(a)') ' Execution of SICOPOLIS stopped. '
181 write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
182 write(6, fmt='(a)') ' '
183 
184 stop
185 
186 #endif
187 
188 end subroutine calc_gia
189 
190 !-------------------------------------------------------------------------------
191 !> Computation of the isostatic steady-state displacement of the lithosphere
192 !! (wss) for the ELRA model.
193 !<------------------------------------------------------------------------------
194 subroutine calc_elra(z_sl, dxi, deta)
195 
196 implicit none
197 
198 real(dp), intent(in) :: z_sl, dxi, deta
199 
200 integer(i4b) :: i, j, ir, jr, il, jl, n
201 integer(i4b) :: ir_max, jr_max, min_imax_jmax
202 integer(i4b) :: il_begin, il_end, jl_begin, jl_end
203 real(dp) :: rho_g, rho_sw_g, rho_a_g_inv
204 real(dp) :: dxi_inv, deta_inv
205 real(dp) :: kei_r_incr_inv
206 real(dp), dimension(0:JMAX,0:IMAX) :: l_r, l_r_inv, fac_wss
207 real(dp) :: ra_max
208 real(dp), allocatable, dimension(:,:) :: f_0
209 
210 real(dp), parameter :: r_infl = 8.0_dp
211  ! Radius of non-locality influence of the elastic
212  ! lithosphere plate (in units of l_r)
213 
214 !-------- Initialisations --------
215 
216 rho_g = rho*g
217 rho_sw_g = rho_sw*g
218 rho_a_g_inv = 1.0_dp/(rho_a*g)
219 
220 dxi_inv = 1.0_dp/dxi
221 deta_inv = 1.0_dp/deta
222 
223 kei_r_incr_inv = 1.0_dp/kei_r_incr
224 
225 do i=0, imax
226 do j=0, jmax
227  l_r(j,i) = (flex_rig_lith(j,i)*rho_a_g_inv)**0.25_dp
228  l_r_inv(j,i) = 1.0_dp/l_r(j,i)
229  fac_wss(j,i) = l_r(j,i)*l_r(j,i)/(2.0_dp*pi*flex_rig_lith(j,i))
230 end do
231 end do
232 
233 ra_max = r_infl*maxval(l_r) ! Radius of non-locality influence (in m)
234 
235 ir_max = floor(ra_max*dxi_inv)
236 jr_max = floor(ra_max*deta_inv)
237  ! Radius of non-locality influence (in grid points)
238 
239 min_imax_jmax = min(imax, jmax)
240 ir_max = min(ir_max, min_imax_jmax)
241 jr_max = min(jr_max, min_imax_jmax)
242  ! ir_max and jr_max constrained by size of domain
243 
244 il_begin = -ir_max
245 il_end = ir_max+imax
246 jl_begin = -jr_max
247 jl_end = jr_max+jmax
248 
249 !-------- Ice/water load --------
250 
251 allocate(f_0(jl_begin:jl_end, il_begin:il_end))
252 
253 do il=il_begin, il_end
254 do jl=jl_begin, jl_end
255 
256  i = min(max(il, 0), imax)
257  j = min(max(jl, 0), jmax)
258 
259  if (maske(j,i)==0) then
260  f_0(jl,il) = rho_g * area(j,i) * (h_c(j,i) + h_t(j,i))
261  else if (maske(j,i)==1) then
262  f_0(jl,il) = 0.0_dp
263  else ! (maske(j,i)>=2)
264  f_0(jl,il) = rho_sw_g * area(j,i) * z_sl
265  ! Water load relative to the present sea-level stand (0 m)
266  ! -> can be positive or negative
267  end if
268 
269 end do
270 end do
271 
272 !-------- Computation of the steady-state displacement of the lithosphere
273 ! (wss, positive downward) --------
274 
275 do i=0, imax
276 do j=0, jmax
277 
278  wss(j,i) = 0.0_dp
279 
280  do ir=-ir_max, ir_max
281  do jr=-jr_max, jr_max
282 
283 #if (GRID==0)
284  n = nint( dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
285 #elif (GRID==1)
286  n = nint( sq_g11_g(j,i)*dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
287  ! The correction with g11 only is OK because g11=g22 for this case
288 #elif (GRID==2)
289  n = nint( dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
290  ! The distortion due to g11 and g22 is already accounted for in
291  ! the variable dist_dxdy(jr,ir)
292 #endif
293 
294  wss(j,i) = wss(j,i) - fac_wss(j,i)*f_0(jr+j,ir+i)*kei(n)
295 
296  end do
297  end do
298 
299 end do
300 end do
301 
302 deallocate (f_0)
303 
304 end subroutine calc_elra
305 
306 !-------------------------------------------------------------------------------
307 !> Generation of an isostatically relaxed lithosphere surface topography
308 !! for either the rigid lithosphere, the local lithosphere or the elastic
309 !! lithosphere model (depending on the setting of the parameter REBOUND).
310 !! This routine is not to be used regularly, and it is only executed if the
311 !! parameter EXEC_MAKE_ZL0 is defined in the header file.
312 !<------------------------------------------------------------------------------
313 subroutine make_zl0()
314 
315 implicit none
316 
317 integer(i4b) :: i, j, m, n, i_f, j_f, n_filter
318 integer(i4b) :: ios
319 real(dp), dimension(0:JMAX,0:IMAX) :: zl0_raw, zl0_smoothed
320 real(dp) :: dx
321 real(dp) :: rho_ratio
322 real(dp) :: filter_width, sigma_filter
323 real(dp) :: dist, weigh, sum_weigh
324 character(len= 8) :: ch_resolution
325 character(len= 64) :: ch_model
326 character(len=256) :: filename, filename_with_path
327 
328 !-------- Checking of the grid --------
329 
330 #if (GRID==0 || GRID==1)
331 
332 dx = real(DX, dp) ! horizontal grid spacing (in km)
333 
334 #else
335 
336 dx = 0.0_dp ! dummy value
337 stop ' >>> make_zl0: Routine works only for GRID 0 or 1!'
338 
339 #endif
340 
341 !-------- Computation of the raw (unsmoothed) topography --------
342 
343 rho_ratio = rho/rho_a
344 
345 #if (REBOUND==0)
346 
347 zl0_raw = zl ! rigid lithosphere
348 
349 #elif (REBOUND==1)
350 
351 do i=0, imax
352 do j=0, jmax
353 
354  if (maske(j,i) == 0) then
355  zl0_raw(j,i) = zl(j,i) + rho_ratio*(h_c(j,i)+h_t(j,i))
356  else
357  zl0_raw(j,i) = zl(j,i)
358  end if ! local lithosphere
359 
360 end do
361 end do
362 
363 #elif (REBOUND==2)
364 
365 zl0_raw = zl + wss ! elastic lithosphere
366 
367 #else
368 
369 stop ' >>> make_zl0: REBOUND must be either 0, 1 or 2!'
370 
371 #endif
372 
373 !-------- Smoothing of the topography (Gaussian filter) --------
374 
375 filter_width = dx ! half span of filtered area, in km
376  ! (chosen as one grid spacing)
377 
378 sigma_filter = filter_width/dx ! half span of filtered area,
379  ! in grid points
380 
381 n_filter = ceiling(2.0_dp*sigma_filter)
382 n_filter = max(n_filter, 5)
383 
384 do i=0, imax
385 do j=0, jmax
386 
387  sum_weigh = 0.0_dp
388 
389  do m=-n_filter, n_filter
390  do n=-n_filter, n_filter
391 
392  i_f = i+m
393  j_f = j+n
394 
395  if (i_f < 0) i_f = 0
396  if (i_f > imax) i_f = imax
397 
398  if (j_f < 0) j_f = 0
399  if (j_f > jmax) j_f = jmax
400 
401  dist = sqrt(real(m,dp)**2+real(n,dp)**2)
402  weigh = exp(-(dist/sigma_filter)**2)
403  sum_weigh = sum_weigh + weigh
404 
405  zl0_smoothed(j,i) = zl0_smoothed(j,i) + weigh*zl0_raw(j_f,i_f)
406 
407  end do
408  end do
409 
410  zl0_smoothed(j,i) = zl0_smoothed(j,i)/sum_weigh
411 
412 end do
413 end do
414 
415 !-------- Writing on file --------
416 
417 if ( abs(dx-nint(dx)) < eps ) then
418  write(ch_resolution, fmt='(i8)') nint(dx)
419 else
420  write(ch_resolution, fmt='(f8.2)') dx
421 end if
422 
423 ch_resolution = adjustl(ch_resolution)
424 
425 #if (REBOUND==0)
426 ch_model = 'rigid_lithosphere'
427 #elif (REBOUND==1)
428 ch_model = 'local_lithosphere'
429 #elif (REBOUND==2)
430 ch_model = 'elastic_lithosphere'
431 #endif
432 
433 filename = trim(ch_domain_short)//'_'//trim(ch_resolution)
434 filename = trim(filename)//'_zl0_'//trim(ch_model)//'.dat'
435 filename_with_path = trim(outpath)//'/'//trim(filename)
436 
437 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='replace')
438 
439 if (ios /= 0) stop ' >>> make_zl0: Error when opening the new zl0 file!'
440 
441 write(23,'(a)') &
442 '%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
443 write(23,'(a)') '% '//trim(ch_domain_long)//':'
444 write(23,'(a)') &
445 '% Relaxed lithosphere surface topography without ice load, in m AMSL.'
446 write(23,'(a)') &
447 '% Horizontal resolution '//trim(ch_resolution)//' km.'
448 write(23,'(a,i3,a,i3,a,i3,a,i3,a)') &
449 '% ', jmax+1, ' records [j = ', &
450 jmax, ' (-1) 0] with ', imax+1, &
451 ' values [i = 0 (1) ', imax, '] in each record.'
452 write(23,'(a)') &
453 '%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
454 
455 do j=jmax, 0, -1
456  do i=0, imax-1
457  zl0_smoothed(j,i) = max(zl0_smoothed(j,i), no_value_neg_2)
458  write(23,'(f8.1)', advance='no') zl0_smoothed(j,i)
459  end do
460  i=imax
461  zl0_smoothed(j,i) = max(zl0_smoothed(j,i), no_value_neg_2)
462  write(23,'(f8.1)') zl0_smoothed(j,i)
463 end do
464 
465 close(23, status='keep')
466 
467 end subroutine make_zl0
468 
469 !-------------------------------------------------------------------------------
470 
471 end module calc_gia_m
472 !
real(dp), parameter eps
eps: Small number
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) zl_target
zl_target(j,i): Target topography (lithosphere surface)
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
real(dp) g
G: Acceleration due to gravity.
real(dp), dimension(0:jmax, 0:imax) zl_neu
(.)_neu: New value of quantity (.) computed during an integration step
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
real(dp), parameter no_value_neg_2
no_value_neg_2: Negative no-value parameter
real(dp), dimension(0:jmax, 0:imax) wss
wss(j,i): Isostatic steady-state displacement of the lithosphere
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
subroutine, public calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)
Main subroutine of calc_gia_m: Computation of the glacial isostatic adjustment of the lithosphere sur...
Definition: calc_gia_m.F90:54
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
real(dp), dimension(-jmax:jmax,-imax:imax) dist_dxdy
dist_dxdy(jr,ir): Distance between grid points with delta_i=ir, delta_j=jr
real(dp), parameter pi
pi: Constant pi
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp) rho_sw
RHO_SW: Density of sea water.
real(dp) rho
RHO: Density of ice.
real(dp) rho_a
RHO_A: Density of the asthenosphere.
integer, parameter dp
Double-precision reals.
Computation of the glacial isostatic adjustment of the lithosphere surface.
Definition: calc_gia_m.F90:36
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) flex_rig_lith
flex_rig_lith(j,i): Flexural rigidity of the lithosphere
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)