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