SICOPOLIS V5-dev  Revision 1420
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-2019 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  use error_m
42 
43  implicit none
44 
45  private
46  public :: calc_gia
47 
48 contains
49 
50 !-------------------------------------------------------------------------------
51 !> Main subroutine of calc_gia_m:
52 !! Computation of the glacial isostatic adjustment of the lithosphere surface.
53 !<------------------------------------------------------------------------------
54 subroutine calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)
55 
56 implicit none
57 
58 integer(i4b), intent(in) :: itercount, iter_wss
59 real(dp), intent(in) :: time
60 real(dp), intent(in) :: z_sl
61 real(dp), intent(in) :: dtime, dxi, deta
62 
63 integer(i4b) :: i, j
64 real(dp) :: tldt_inv(0:jmax,0:imax)
65 real(dp) :: dtime_inv
66 real(dp) :: rho_rhoa_ratio, rhosw_rhoa_ratio
67 real(dp) :: time_dimless
68 real(dp) :: target_topo_tau_factor, target_topo_tau
69 
70 !-------- Term abbreviations --------
71 
72 do i=0, imax
73 do j=0, jmax
74  tldt_inv(j,i) = 1.0_dp/(time_lag_asth(j,i)+dtime)
75 end do
76 end do
77 
78 rho_rhoa_ratio = rho/rho_a
79 rhosw_rhoa_ratio = rho_sw/rho_a
80 
81 dtime_inv = 1.0_dp/dtime
82 
83 !-------- Computation of zl_neu and its time derivative --------
84 
85 #if (REBOUND==0)
86 
87 do i=0, imax
88 do j=0, jmax
89  zl_neu(j,i) = zl(j,i)
90  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
91 end do
92 end do
93 
94 #elif (REBOUND==1)
95 
96 do i=0, imax
97 do j=0, jmax
98 
99  if (maske(j,i) <= 1_i1b) then ! grounded ice or ice-free land
100 
101  zl_neu(j,i) = tldt_inv(j,i)*( time_lag_asth(j,i)*zl(j,i) &
102  + dtime*(zl0(j,i) &
103  -frac_llra*rho_rhoa_ratio*(h_c(j,i)+h_t(j,i))) )
104 
105  else ! (maske(j,i) >= 2_i1b)
106 
107  zl_neu(j,i) = tldt_inv(j,i)*( time_lag_asth(j,i)*zl(j,i) &
108  + dtime*(zl0(j,i) &
109  -frac_llra*rhosw_rhoa_ratio*z_sl) )
110 
111  ! Water load relative to the present sea-level stand (0 m)
112  ! -> can be positive or negative
113 
114  end if
115 
116  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
117 
118 end do
119 end do
120 
121 #elif (REBOUND==2)
122 
123 if ((mod(itercount, iter_wss)==0).or.(itercount==1)) then
124  write(6, fmt='(10x,a)') 'Computation of wss'
125  call calc_elra(z_sl, dxi, deta) ! Update of the steady-state displacement
126  ! of the lithosphere (wss)
127 end if
128 
129 do i=0, imax
130 do j=0, jmax
131  zl_neu(j,i) = tldt_inv(j,i) &
132  * ( time_lag_asth(j,i)*zl(j,i) + dtime*(zl0(j,i)-wss(j,i)) )
133  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
134 end do
135 end do
136 
137 #endif
138 
139 ! ------ Adjustment due to prescribed target topography
140 
141 #if (THK_EVOL==2)
142 
143 if (time >= time_target_topo_final) then
144 
145  zl_neu = zl_target
146  dzl_dtau = 0.0_dp
147 
148 else if (time > time_target_topo_init) then
149 
150  time_dimless = (time -time_target_topo_init) &
152 
153  time_dimless = max(min(time_dimless, 1.0_dp), 0.0_dp)
154  ! constrain to interval [0,1]
155 
156  target_topo_tau_factor = time_dimless*time_dimless*time_dimless &
157  *(10.0_dp + time_dimless &
158  *(-15.0_dp+6.0_dp*time_dimless))
159  ! make transition smooth (quintic function)
160 
161  target_topo_tau_factor = 1.0_dp-target_topo_tau_factor
162 
163  target_topo_tau = target_topo_tau_0 * target_topo_tau_factor
164 
165  zl_neu = ( target_topo_tau*zl_neu + dtime*zl_target ) &
166  / ( target_topo_tau + dtime )
167 
168  dzl_dtau = (zl_neu-zl)*dtime_inv
169 
170 end if
171 
172 #elif (THK_EVOL==3)
173 
174 target_topo_tau = target_topo_tau_0
175 
176 zl_neu = ( target_topo_tau*zl_neu + dtime*zl_target ) &
177  / ( target_topo_tau + dtime )
178 
179 dzl_dtau = (zl_neu-zl)*dtime_inv
180 
181 #endif
182 
183 !======== Generation of an isostatically relaxed
184 ! lithosphere surface topography. Not to be used regularly! ========
185 
186 #if (defined(EXEC_MAKE_ZL0))
187 
188 call make_zl0()
189 
190 errormsg = ' >>> calc_gia: Routine make_zl0 successfully completed,' &
191  // end_of_line &
192  //' isostatically relaxed lithosphere' &
193  // end_of_line &
194  //' topography written on file' &
195  // end_of_line &
196  //' (in directory specified by OUTPATH).' &
197  // end_of_line &
198  //' Execution of SICOPOLIS stopped.'
199 call error(errormsg) ! actually not an error,
200  ! just a regular stop with an info message
201 
202 #endif
203 
204 end subroutine calc_gia
205 
206 !-------------------------------------------------------------------------------
207 !> Computation of the isostatic steady-state displacement of the lithosphere
208 !! (wss) for the ELRA model.
209 !<------------------------------------------------------------------------------
210 subroutine calc_elra(z_sl, dxi, deta)
211 
212 #if defined(ALLOW_OPENAD) /* OpenAD */
213  use ctrl_m, only: myfloor
214 #endif /* OpenAD */
215 
216 implicit none
217 
218 real(dp), intent(in) :: z_sl, dxi, deta
219 
220 integer(i4b) :: i, j, ir, jr, il, jl, n
221 integer(i4b) :: ir_max, jr_max, min_imax_jmax
222 integer(i4b) :: il_begin, il_end, jl_begin, jl_end
223 real(dp) :: rho_g, rho_sw_g, rho_a_g_inv
224 real(dp) :: dxi_inv, deta_inv
225 real(dp) :: kei_r_incr_inv
226 real(dp), dimension(0:JMAX,0:IMAX) :: l_r, l_r_inv, fac_wss
227 real(dp) :: ra_max
228 
229 #if !defined(ALLOW_OPENAD) /* Normal */
230 real(dp), allocatable, dimension(:,:) :: f_0
231 #else /* OpenAD */
232 real(dp), dimension(-JMAX:2*JMAX,-IMAX:2*IMAX) :: f_0
233 #endif /* Normal vs. OpenAD */
234 
235 real(dp), parameter :: r_infl = 8.0_dp
236  ! Radius of non-locality influence of the elastic
237  ! lithosphere plate (in units of l_r)
238 
239 !-------- Initialisations --------
240 
241 rho_g = rho*g
242 rho_sw_g = rho_sw*g
243 rho_a_g_inv = 1.0_dp/(rho_a*g)
244 
245 dxi_inv = 1.0_dp/dxi
246 deta_inv = 1.0_dp/deta
247 
248 kei_r_incr_inv = 1.0_dp/kei_r_incr
249 
250 do i=0, imax
251 do j=0, jmax
252  l_r(j,i) = (flex_rig_lith(j,i)*rho_a_g_inv)**0.25_dp
253  l_r_inv(j,i) = 1.0_dp/l_r(j,i)
254  fac_wss(j,i) = l_r(j,i)*l_r(j,i)/(2.0_dp*pi*flex_rig_lith(j,i))
255 end do
256 end do
257 
258 ra_max = r_infl*maxval(l_r) ! Radius of non-locality influence (in m)
259 
260 #if !defined(ALLOW_OPENAD) /* Normal */
261 
262 ir_max = floor(ra_max*dxi_inv)
263 jr_max = floor(ra_max*deta_inv)
264  ! Radius of non-locality influence (in grid points)
265 
266 #else /* OpenAD */
267 
268 call myfloor(ra_max*dxi_inv, ir_max)
269 call myfloor(ra_max*deta_inv, jr_max)
270  ! Radius of non-locality influence (in grid points)
271 
272 if (ir_max.gt.imax) then
273  write(*,*) "ERROR: unhandled case ir_max greater than IMAX for f_0 dimensions"
274 end if
275 if (jr_max.gt.jmax) then
276  write(*,*) "ERROR: unhandled case jr_max greater than JMAX for f_0 dimensions"
277 end if
278 
279 #endif /* Normal vs. OpenAD */
280 
281 min_imax_jmax = min(imax, jmax)
282 ir_max = min(ir_max, min_imax_jmax)
283 jr_max = min(jr_max, min_imax_jmax)
284  ! ir_max and jr_max constrained by size of domain
285 
286 il_begin = -ir_max
287 il_end = ir_max+imax
288 jl_begin = -jr_max
289 jl_end = jr_max+jmax
290 
291 !-------- Ice/water load --------
292 
293 #if !defined(ALLOW_OPENAD) /* Normal */
294 allocate(f_0(jl_begin:jl_end, il_begin:il_end))
295 #endif /* Normal */
296 
297 do il=il_begin, il_end
298 do jl=jl_begin, jl_end
299 
300  i = min(max(il, 0), imax)
301  j = min(max(jl, 0), jmax)
302 
303  if (maske(j,i)==0_i1b) then
304  f_0(jl,il) = rho_g * area(j,i) * (h_c(j,i) + h_t(j,i))
305  else if (maske(j,i)==1_i1b) then
306  f_0(jl,il) = 0.0_dp
307  else ! (maske(j,i)>=2_i1b)
308  f_0(jl,il) = rho_sw_g * area(j,i) * z_sl
309  ! Water load relative to the present sea-level stand (0 m)
310  ! -> can be positive or negative
311  end if
312 
313 end do
314 end do
315 
316 !-------- Computation of the steady-state displacement of the lithosphere
317 ! (wss, positive downward) --------
318 
319 do i=0, imax
320 do j=0, jmax
321 
322  wss(j,i) = 0.0_dp
323 
324  do ir=-ir_max, ir_max
325  do jr=-jr_max, jr_max
326 
327 #if (GRID==0)
328  n = nint( dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
329 #elif (GRID==1)
330  n = nint( sq_g11_g(j,i)*dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
331  ! The correction with g11 only is OK because g11=g22 for this case
332 #elif (GRID==2)
333  n = nint( dist_dxdy(jr,ir)*l_r_inv(j,i)*kei_r_incr_inv )
334  ! The distortion due to g11 and g22 is already accounted for in
335  ! the variable dist_dxdy(jr,ir)
336 #endif
337 
338  wss(j,i) = wss(j,i) - fac_wss(j,i)*f_0(jr+j,ir+i)*kei(n)
339 
340  end do
341  end do
342 
343 end do
344 end do
345 
346 #if !defined(ALLOW_OPENAD) /* Normal */
347 deallocate (f_0)
348 #endif /* Normal */
349 
350 end subroutine calc_elra
351 
352 !-------------------------------------------------------------------------------
353 !> Generation of an isostatically relaxed lithosphere surface topography
354 !! for either the rigid lithosphere, the local lithosphere or the elastic
355 !! lithosphere model (depending on the setting of the parameter REBOUND).
356 !! This routine is not to be used regularly, and it is only executed if the
357 !! parameter EXEC_MAKE_ZL0 is defined in the header file.
358 !<------------------------------------------------------------------------------
359 subroutine make_zl0()
360 
361 #if defined(ALLOW_OPENAD) /* OpenAD */
362  use ctrl_m, only: myceiling
363 #endif /* OpenAD */
364 
365 implicit none
366 
367 integer(i4b) :: i, j, m, n, i_f, j_f, n_filter
368 integer(i4b) :: ios
369 real(dp), dimension(0:JMAX,0:IMAX) :: zl0_raw, zl0_smoothed
370 real(dp) :: dx
371 real(dp) :: rho_ratio
372 real(dp) :: filter_width, sigma_filter
373 real(dp) :: dist, weigh, sum_weigh
374 character(len= 8) :: ch_resolution
375 character(len= 64) :: ch_model
376 character(len=256) :: filename, filename_with_path
377 
378 !-------- Checking of the grid --------
379 
380 #if (GRID==0 || GRID==1)
381 
382 dx = real(DX, dp) ! horizontal grid spacing (in km)
383 
384 #else
385 
386 dx = 0.0_dp ! dummy value
387 
388 errormsg = ' >>> make_zl0: Routine works only for GRID 0 or 1!'
389 call error(errormsg)
390 
391 #endif
392 
393 !-------- Computation of the raw (unsmoothed) topography --------
394 
395 rho_ratio = rho/rho_a
396 
397 #if (REBOUND==0)
398 
399 zl0_raw = zl ! rigid lithosphere
400 
401 #elif (REBOUND==1)
402 
403 do i=0, imax
404 do j=0, jmax
405 
406  if (maske(j,i) == 0_i1b) then
407  zl0_raw(j,i) = zl(j,i) + rho_ratio*(h_c(j,i)+h_t(j,i))
408  else
409  zl0_raw(j,i) = zl(j,i)
410  end if ! local lithosphere
411 
412 end do
413 end do
414 
415 #elif (REBOUND==2)
416 
417 zl0_raw = zl + wss ! elastic lithosphere
418 
419 #else
420 
421 errormsg = ' >>> make_zl0: REBOUND must be either 0, 1 or 2!'
422 call error(errormsg)
423 
424 #endif
425 
426 !-------- Smoothing of the topography (Gaussian filter) --------
427 
428 filter_width = dx ! half span of filtered area, in km
429  ! (chosen as one grid spacing)
430 
431 sigma_filter = filter_width/dx ! half span of filtered area,
432  ! in grid points
433 
434 #if !defined(ALLOW_OPENAD) /* Normal */
435 n_filter = ceiling(2.0_dp*sigma_filter)
436 #else /* OpenAD */
437 call myceiling(2.0_dp*sigma_filter, n_filter)
438 #endif /* Normal vs. OpenAD */
439 
440 n_filter = max(n_filter, 5)
441 
442 do i=0, imax
443 do j=0, jmax
444 
445  sum_weigh = 0.0_dp
446 
447  do m=-n_filter, n_filter
448  do n=-n_filter, n_filter
449 
450  i_f = i+m
451  j_f = j+n
452 
453  if (i_f < 0) i_f = 0
454  if (i_f > imax) i_f = imax
455 
456  if (j_f < 0) j_f = 0
457  if (j_f > jmax) j_f = jmax
458 
459  dist = sqrt(real(m,dp)**2+real(n,dp)**2)
460  weigh = exp(-(dist/sigma_filter)**2)
461  sum_weigh = sum_weigh + weigh
462 
463  zl0_smoothed(j,i) = zl0_smoothed(j,i) + weigh*zl0_raw(j_f,i_f)
464 
465  end do
466  end do
467 
468  zl0_smoothed(j,i) = zl0_smoothed(j,i)/sum_weigh
469 
470 end do
471 end do
472 
473 !-------- Writing on file --------
474 
475 if ( abs(dx-nint(dx)) < eps ) then
476  write(ch_resolution, fmt='(i8)') nint(dx)
477 else
478  write(ch_resolution, fmt='(f8.2)') dx
479 end if
480 
481 #if !defined(ALLOW_OPENAD) /* Normal */
482 ch_resolution = adjustl(ch_resolution)
483 #endif /* Normal */
484 
485 #if (REBOUND==0)
486 ch_model = 'rigid_lithosphere'
487 #elif (REBOUND==1)
488 ch_model = 'local_lithosphere'
489 #elif (REBOUND==2)
490 ch_model = 'elastic_lithosphere'
491 #endif
492 
493 filename = trim(ch_domain_short)//'_'//trim(ch_resolution)
494 filename = trim(filename)//'_zl0_'//trim(ch_model)//'.dat'
495 filename_with_path = trim(outpath)//'/'//trim(filename)
496 
497 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='replace')
498 
499 if (ios /= 0) then
500  errormsg = ' >>> make_zl0: Error when opening the new zl0 file!'
501  call error(errormsg)
502 end if
503 
504 write(23,'(a)') &
505 '%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
506 write(23,'(a)') '% '//trim(ch_domain_long)//':'
507 write(23,'(a)') &
508 '% Relaxed lithosphere surface topography without ice load, in m AMSL.'
509 write(23,'(a)') &
510 '% Horizontal resolution '//trim(ch_resolution)//' km.'
511 write(23,'(a,i3,a,i3,a,i3,a,i3,a)') &
512 '% ', jmax+1, ' records [j = ', &
513 jmax, ' (-1) 0] with ', imax+1, &
514 ' values [i = 0 (1) ', imax, '] in each record.'
515 write(23,'(a)') &
516 '%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
517 
518 do j=jmax, 0, -1
519  do i=0, imax-1
520  zl0_smoothed(j,i) = max(zl0_smoothed(j,i), no_value_neg_2)
521  write(23,'(f8.1)', advance='no') zl0_smoothed(j,i)
522  end do
523  i=imax
524  zl0_smoothed(j,i) = max(zl0_smoothed(j,i), no_value_neg_2)
525  write(23,'(f8.1)') zl0_smoothed(j,i)
526 end do
527 
528 close(23, status='keep')
529 
530 end subroutine make_zl0
531 
532 !-------------------------------------------------------------------------------
533 
534 end module calc_gia_m
535 !
Definition: ctrl_m.F90:9
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
real(dp), parameter eps
eps: Small number
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), 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...
real(dp) target_topo_tau_0
target_topo_tau_0: Relaxation time for target-topography adjustment
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:55
Writing of error messages and stopping execution.
Definition: error_m.F90:35
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
character, parameter end_of_line
end_of_line: End-of-line string
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
integer(i1b), 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(-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.
character(len=256) errormsg
errormsg: Error-message string
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)