SICOPOLIS V5-dev  Revision 1173
calc_thk_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ t h k _ m
4 !
5 !> @file
6 !!
7 !! Computation of the ice thickness.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 Ralf Greve, Reinhard Calov, Tatsuru Sato
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 ice thickness.
34 !<------------------------------------------------------------------------------
35 module calc_thk_m
36 
37  use sico_types_m
39  use sico_vars_m
40 
41  implicit none
42 
43  real(dp), dimension(0:JMAX,0:IMAX), save :: h, h_neu
44  real(dp), dimension(0:JMAX,0:IMAX), save :: mb_source
45  logical, save :: flag_solver_explicit
46 
47  private
48  public :: calc_thk_init
50  public :: calc_thk_expl, calc_thk_impl
51  public :: calc_thk_mask_update
52 
53 contains
54 
55 !-------------------------------------------------------------------------------
56 !> Initialisations for the ice thickness computation.
57 !<------------------------------------------------------------------------------
58 subroutine calc_thk_init()
59 
60 implicit none
61 
62 #if (defined(CALCZS))
63  stop ' >>> calc_thk_init: Replace CALCZS by CALCTHK in header file!'
64 #elif (!defined(CALCTHK))
65  stop ' >>> calc_thk_init: Define CALCTHK in header file!'
66 #endif
67 
68 !-------- Computation/initialisation of the ice base topography
69 ! and its time derivative --------
70 
71 #if (MARGIN==1 || MARGIN==2) /* only grounded ice */
72 
73 zb_neu = zl_neu
75 
76 #elif (MARGIN==3) /* grounded and floating ice */
77 
78 where (maske <= 1_i2b) ! grounded ice or ice-free land
79  zb_neu = zl_neu
81 elsewhere ! (maske >= 2_i2b; ocean or floating ice)
82  zb_neu = zb ! initialisation,
83  dzb_dtau = 0.0_dp ! will be overwritten later
84 end where
85 
86 #else
87 stop ' >>> calc_thk_init: MARGIN must be either 1, 2 or 3!'
88 #endif
89 
90 !-------- Initialisation of the ice thickness
91 ! and surface topography --------
92 
93 h = h_c + h_t
94 
95 zs_neu = zs ! initialisation,
96 h_neu = h ! will be overwritten later
97 
98 !-------- Solver type --------
99 
100 #if (CALCTHK==1 || CALCTHK==4) /* explicit solver */
101 
102  flag_solver_explicit = .true.
103 
104 #elif (CALCTHK==2 || CALCTHK==3 \
105  || calcthk==5 || calcthk==6) /* implicit solver */
106 
107  flag_solver_explicit = .false.
108 
109 #else
110 stop ' >>> calc_thk_init: CALCTHK must be between 1 and 6!'
111 #endif
112 
113 !-------- Source term for the ice thickness equation --------
114 
115 #if (MB_ACCOUNT==0)
116 
118 
119 #elif (MB_ACCOUNT==1)
120 
121  if (flag_solver_explicit) then
122  mb_source = 0.0_dp ! source term will be applied separately
123  else
125  end if
126 
127 #else
128  stop ' >>> calc_thk_init: MB_ACCOUNT must be either 0 or 1!'
129 #endif
130 
131 end subroutine calc_thk_init
132 
133 !-------------------------------------------------------------------------------
134 !> Explicit solver for the diffusive SIA ice surface equation.
135 !<------------------------------------------------------------------------------
136 subroutine calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
138 implicit none
139 
140 real(dp), intent(in) :: time, dtime, dxi, deta
141 real(dp), intent(in) :: z_mar
142 
143 integer(i4b) :: i, j
144 real(dp) :: azs2, azs3
145 real(dp), dimension(0:JMAX,0:IMAX) :: czs2, czs3
146 
147 !-------- Abbreviations --------
148 
149 azs2 = dtime/(dxi*dxi)
150 azs3 = dtime/(deta*deta)
151 
152 czs2 = 0.0_dp
153 czs3 = 0.0_dp
154 
155 do i=0, imax-1
156 do j=0, jmax
157  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
158  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
159 end do
160 end do
161 
162 do i=0, imax
163 do j=0, jmax-1
164  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
165  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
166 end do
167 end do
168 
169 !-------- Solution of the explicit scheme --------
170 
171 do i=0, imax
172 do j=0, jmax
173 
174  if ( flag_inner_point(j,i) & ! inner point
175 #if (MARGIN==1 && MB_ACCOUNT==0)
176  .and.(maske(j,i) <= 1_i2b) & ! grounded ice or ice-free land
177 #endif
178  ) then
179 
180  zs_neu(j,i) = zs(j,i) &
181  + dtime*(mb_source(j,i)+dzb_dtau(j,i)) &
182  + ( czs2(j,i) *(zs(j,i+1)-zs(j,i) ) &
183  -czs2(j,i-1)*(zs(j,i) -zs(j,i-1)) &
184  +czs3(j,i) *(zs(j+1,i)-zs(j,i) ) &
185  -czs3(j-1,i)*(zs(j,i) -zs(j-1,i)) ) &
186  *insq_g11_g(j,i)*insq_g22_g(j,i)
187 
188  else
189  zs_neu(j,i) = zb_neu(j,i) ! zero-thickness boundary condition
190  end if
191 
192 end do
193 end do
194 
195 !-------- Ice thickness --------
196 
197 h_neu = zs_neu - zb_neu
198 
199 !-------- Applying the so-far neglected source term --------
200 
201 #if (MB_ACCOUNT==1)
202 call apply_mb_source(dtime, z_mar)
203 #endif
204 
205 !-------- Adjusting the ice thickness, if needed --------
206 
207 call thk_adjust(time, dtime)
208 
209 end subroutine calc_thk_sia_expl
210 
211 !-------------------------------------------------------------------------------
212 !> Over-implicit solver for the diffusive SIA ice surface equation.
213 !<------------------------------------------------------------------------------
214 subroutine calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, &
215  ii, jj, nn)
217 use sico_maths_m, only : sor_sprs
218 
219 implicit none
220 
221 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
222 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
223 real(dp), intent(in) :: time, dtime, dxi, deta
224 real(dp), intent(in) :: z_mar
225 real(dp), intent(in) :: mean_accum
226 
227 integer(i4b) :: i, j
228 integer(i4b) :: k, nnz
229 real(dp) :: azs2, azs3
230 real(dp), dimension(0:JMAX,0:IMAX) :: czs2, czs3
231 
232 #if (CALCTHK==2)
233 
234 integer(i4b) :: ierr
235 integer(i4b) :: iter
236 integer(i4b) :: nc, nr
237 integer(i4b), parameter :: nmax = (imax+1)*(jmax+1)
238 integer(i4b), parameter :: n_sprs = 10*(imax+1)*(jmax+1)
239 integer(i4b), allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
240 integer(i4b), allocatable, dimension(:) :: lgs_a_diag_index
241 integer(i4b), allocatable, dimension(:) :: lgs_a_index_trim
242 real(dp), allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
243 real(dp), allocatable, dimension(:) :: lgs_a_value_trim
244 real(dp) :: eps_sor
245 
246 #elif (CALCTHK==3)
247 
248 lis_integer :: ierr
249 lis_integer :: iter
250 lis_integer :: nc, nr
251 lis_integer, parameter :: nmax = (imax+1)*(jmax+1)
252 lis_integer, parameter :: n_sprs = 10*(imax+1)*(jmax+1)
253 lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
254 lis_integer, allocatable, dimension(:) :: lgs_a_diag_index
255 lis_matrix :: lgs_a
256 lis_vector :: lgs_b, lgs_x
257 lis_scalar, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
258 lis_solver :: solver
259 character(len=256) :: ch_solver_set_option
260 
261 #endif
262 
263 !-------- Abbreviations --------
264 
265 azs2 = dtime/(dxi*dxi)
266 azs3 = dtime/(deta*deta)
267 
268 czs2 = 0.0_dp
269 czs3 = 0.0_dp
270 
271 do i=0, imax-1
272 do j=0, jmax
273  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
274  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
275 end do
276 end do
277 
278 do i=0, imax
279 do j=0, jmax-1
280  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
281  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
282 end do
283 end do
284 
285 !-------- Assembly of the system of linear equations
286 ! (matrix storage: compressed sparse row CSR) --------
287 
288 #if (CALCTHK==2 || CALCTHK==3)
289 
290 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
291 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
292 
293 lgs_a_value = 0.0_dp
294 lgs_a_index = 0
295 lgs_a_ptr = 0
296 
297 lgs_b_value = 0.0_dp
298 lgs_x_value = 0.0_dp
299 
300 lgs_a_ptr(1) = 1
301 
302 k = 0
303 
304 do nr=1, nmax ! loop over rows
305 
306  i = ii(nr)
307  j = jj(nr)
308 
309  if ( flag_inner_point(j,i) & ! inner point
310 #if (MARGIN==1 && MB_ACCOUNT==0)
311  .and.(maske(j,i) <= 1_i2b) & ! grounded ice or ice-free land
312 #endif
313  ) then
314 
315  nc = nn(j-1,i) ! smallest nc (column counter), for zs(j-1,i)
316  k = k+1
317  ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
318  lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
319  *insq_g11_g(j,i)*insq_g22_g(j,i)
320  lgs_a_index(k) = nc
321 
322  nc = nn(j,i-1) ! next nc (column counter), for zs(j,i-1)
323  k = k+1
324  ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
325  lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
326  *insq_g11_g(j,i)*insq_g22_g(j,i)
327  lgs_a_index(k) = nc
328 
329  nc = nn(j,i) ! next nc (column counter), for zs(j,i)
330  ! if (nc /= nr) & (diagonal element)
331  ! stop ' >>> calc_thk_sia_impl: Check for diagonal element failed!'
332  k = k+1
333  ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
334  lgs_a_value(k) = 1.0_dp &
335  + (czs2(j,i)+czs2(j,i-1)+czs3(j,i)+czs3(j-1,i)) &
336  *ovi_weight &
337  *insq_g11_g(j,i)*insq_g22_g(j,i)
338  lgs_a_diag_index(nr) = k
339  lgs_a_index(k) = nc
340 
341  nc = nn(j,i+1) ! next nc (column counter), for zs(j,i+1)
342  k = k+1
343  ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
344  lgs_a_value(k) = -czs2(j,i)*ovi_weight &
345  *insq_g11_g(j,i)*insq_g22_g(j,i)
346  lgs_a_index(k) = nc
347 
348  nc = nn(j+1,i) ! largest nc (column counter), for zs(j+1,i)
349  k = k+1
350  ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
351  lgs_a_value(k) = -czs3(j,i)*ovi_weight &
352  *insq_g11_g(j,i)*insq_g22_g(j,i)
353  lgs_a_index(k) = nc
354 
355  lgs_b_value(nr) = zs(j,i) &
356  + dtime*(mb_source(j,i)+dzb_dtau(j,i)) &
357  + ( czs2(j,i)*(zs(j,i+1)-zs(j,i)) &
358  -czs2(j,i-1)*(zs(j,i)-zs(j,i-1)) &
359  +czs3(j,i)*(zs(j+1,i)-zs(j,i)) &
360  -czs3(j-1,i)*(zs(j,i)-zs(j-1,i)) ) &
361  *(1.0_dp-ovi_weight) &
362  *insq_g11_g(j,i)*insq_g22_g(j,i)
363  ! right-hand side
364 
365  else ! zero-thickness boundary condition
366 
367  k = k+1
368  ! if (k > n_sprs) stop ' >>> calc_thk_sia_impl: n_sprs too small!'
369  lgs_a_value(k) = 1.0_dp ! diagonal element only
370  lgs_a_diag_index(nr) = k
371  lgs_a_index(k) = nr
372  lgs_b_value(nr) = zb_neu(j,i)
373 
374  end if
375 
376  lgs_x_value(nr) = zs(j,i) ! old surface topography,
377  ! initial guess for solution vector
378 
379  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
380 
381 end do
382 
383 nnz = k ! number of non-zero elements of the matrix
384 
385 #else
386 stop ' >>> calc_thk_sia_impl: CALCTHK must be either 2 or 3!'
387 #endif
388 
389 !-------- Solution of the system of linear equations --------
390 
391 #if (CALCTHK==2)
392 
393 ! ------ Solution with the built-in SOR solver
394 
395 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
396 
397 do k=1, nnz ! relocate matrix to trimmed arrays
398  lgs_a_value_trim(k) = lgs_a_value(k)
399  lgs_a_index_trim(k) = lgs_a_index(k)
400 end do
401 
402 deallocate(lgs_a_value, lgs_a_index)
403 
404 eps_sor = 1.0e-05_dp*mean_accum*dtime ! convergence parameter
405 
406 call sor_sprs(lgs_a_value_trim, &
407  lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
408  lgs_b_value, &
409  nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
410 
411 do nr=1, nmax
412  i = ii(nr)
413  j = jj(nr)
414  zs_neu(j,i) = lgs_x_value(nr)
415 end do
416 
417 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
418 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
419 
420 #elif (CALCTHK==3)
421 
422 ! ------ Settings for Lis
423 
424 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
425 call lis_vector_create(lis_comm_world, lgs_b, ierr)
426 call lis_vector_create(lis_comm_world, lgs_x, ierr)
427 
428 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
429 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
430 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
431 
432 do nr=1, nmax
433 
434  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
435  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
436  lgs_a_value(nc), lgs_a, ierr)
437  end do
438 
439  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
440  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
441 
442 end do
443 
444 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
445 call lis_matrix_assemble(lgs_a, ierr)
446 
447 ! ------ Solution with Lis
448 
449 call lis_solver_create(solver, ierr)
450 
451 ch_solver_set_option = '-i bicg -p ilu '// &
452  '-maxiter 1000 -tol 1.0e-04 -initx_zeros false'
453 
454 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
455 call chkerr(ierr)
456 
457 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
458 call chkerr(ierr)
459 
460 call lis_solver_get_iter(solver, iter, ierr)
461 write(6,'(10x,a,i0)') 'calc_thk_sia_impl: iter = ', iter
462 
463 lgs_x_value = 0.0_dp
464 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
465 call lis_matrix_destroy(lgs_a, ierr)
466 call lis_vector_destroy(lgs_b, ierr)
467 call lis_vector_destroy(lgs_x, ierr)
468 call lis_solver_destroy(solver, ierr)
469 
470 do nr=1, nmax
471  i = ii(nr)
472  j = jj(nr)
473  zs_neu(j,i) = lgs_x_value(nr)
474 end do
475 
476 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
477 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
478 
479 #endif
480 
481 !-------- Ice thickness --------
482 
483 h_neu = zs_neu - zb_neu
484 
485 !-------- Applying the so-far neglected source term --------
486 
487 #if (MB_ACCOUNT==1)
488 call apply_mb_source(dtime, z_mar)
489 #endif
490 
491 !-------- Adjusting the ice thickness, if needed --------
492 
493 call thk_adjust(time, dtime)
494 
495 end subroutine calc_thk_sia_impl
496 
497 !-------------------------------------------------------------------------------
498 !> Explicit solver for the general ice thickness equation.
499 !<------------------------------------------------------------------------------
500 subroutine calc_thk_expl(time, dtime, dxi, deta, z_mar)
502 implicit none
503 
504 real(dp), intent(in) :: time, dtime, dxi, deta
505 real(dp), intent(in) :: z_mar
506 
507 integer(i4b) :: i, j
508 real(dp), dimension(0:JMAX,0:IMAX) :: dt_darea
509 real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
510 real(dp), dimension(0:JMAX,0:IMAX) :: upH_x_1, upH_x_2, upH_y_1, upH_y_2
511 real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_x_1, sq_g22_x_2
512 real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_y_1, sq_g11_y_2
513 
514 !-------- Abbreviations --------
515 
516 dt_darea = dtime/area
517 
518 do i=0, imax
519 do j=0, jmax
520 
521  if (flag_inner_point(j,i)) then
522 
523  vx_m_1(j,i) = vx_m(j,i-1)
524  vx_m_2(j,i) = vx_m(j,i)
525  vy_m_1(j,i) = vy_m(j-1,i)
526  vy_m_2(j,i) = vy_m(j,i)
527 
528  if (vx_m_1(j,i) >= 0.0_dp) then
529  uph_x_1(j,i) = h(j,i-1)
530  else
531  uph_x_1(j,i) = h(j,i)
532  end if
533 
534  if (vx_m_2(j,i) >= 0.0_dp) then
535  uph_x_2(j,i) = h(j,i)
536  else
537  uph_x_2(j,i) = h(j,i+1)
538  end if
539 
540  if (vy_m_1(j,i) >= 0.0_dp) then
541  uph_y_1(j,i) = h(j-1,i)
542  else
543  uph_y_1(j,i) = h(j,i)
544  end if
545 
546  if (vy_m_2(j,i) >= 0.0_dp) then
547  uph_y_2(j,i) = h(j,i)
548  else
549  uph_y_2(j,i) = h(j+1,i)
550  end if
551 
552  sq_g22_x_1(j,i) = sq_g22_sgx(j,i-1)
553  sq_g22_x_2(j,i) = sq_g22_sgx(j,i)
554  sq_g11_y_1(j,i) = sq_g11_sgy(j-1,i)
555  sq_g11_y_2(j,i) = sq_g11_sgy(j,i)
556 
557  else ! .not.(flag_inner_point(j,i))
558 
559  vx_m_1(j,i) = 0.0_dp
560  vx_m_2(j,i) = 0.0_dp
561  vy_m_1(j,i) = 0.0_dp
562  vy_m_2(j,i) = 0.0_dp
563 
564  uph_x_1(j,i) = 0.0_dp
565  uph_x_2(j,i) = 0.0_dp
566  uph_y_1(j,i) = 0.0_dp
567  uph_y_2(j,i) = 0.0_dp
568 
569  sq_g22_x_1(j,i) = 0.0_dp
570  sq_g22_x_2(j,i) = 0.0_dp
571  sq_g11_y_1(j,i) = 0.0_dp
572  sq_g11_y_2(j,i) = 0.0_dp
573 
574  end if
575 
576 end do
577 end do
578 
579 !-------- Solution of the explicit scheme --------
580 
581 where ( flag_inner_point & ! inner point
582 #if (MARGIN==1 && MB_ACCOUNT==0)
583  .and.(maske <= 1_i2b) & ! grounded ice or ice-free land
584 #endif
585  )
586 
587  h_neu = h + dtime*mb_source &
588  - dt_darea &
589  * ( ( vx_m_2*uph_x_2*sq_g22_x_2*deta &
590  -vx_m_1*uph_x_1*sq_g22_x_1*deta ) &
591  + ( vy_m_2*uph_y_2*sq_g11_y_2*dxi &
592  -vy_m_1*uph_y_1*sq_g11_y_1*dxi ) )
593 
594 elsewhere
595  h_neu = 0.0_dp ! zero-thickness boundary condition
596 end where
597 
598 !-------- Applying the so-far neglected source term --------
599 
600 #if (MB_ACCOUNT==1)
601 call apply_mb_source(dtime, z_mar)
602 #endif
603 
604 !-------- Adjusting the ice thickness, if needed --------
605 
606 call thk_adjust(time, dtime)
607 
608 end subroutine calc_thk_expl
609 
610 !-------------------------------------------------------------------------------
611 !> Over-implicit solver for the general ice thickness equation.
612 !<------------------------------------------------------------------------------
613 subroutine calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, &
614  ii, jj, nn)
616 use sico_maths_m, only : sor_sprs
617 
618 implicit none
619 
620 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
621 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
622 real(dp), intent(in) :: time, dtime, dxi, deta
623 real(dp), intent(in) :: z_mar
624 real(dp), intent(in) :: mean_accum
625 
626 integer(i4b) :: i, j
627 integer(i4b) :: k, nnz
628 real(dp), dimension(0:JMAX,0:IMAX) :: dt_darea
629 real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
630 real(dp), dimension(0:JMAX,0:IMAX) :: upH_x_1, upH_x_2, upH_y_1, upH_y_2
631 real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_x_1, sq_g22_x_2
632 real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_y_1, sq_g11_y_2
633 
634 #if (CALCTHK==5)
635 
636 integer(i4b) :: ierr
637 integer(i4b) :: iter
638 integer(i4b) :: nc, nr
639 integer(i4b), parameter :: nmax = (imax+1)*(jmax+1)
640 integer(i4b), parameter :: n_sprs = 10*(imax+1)*(jmax+1)
641 integer(i4b), allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
642 integer(i4b), allocatable, dimension(:) :: lgs_a_diag_index
643 integer(i4b), allocatable, dimension(:) :: lgs_a_index_trim
644 real(dp), allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
645 real(dp), allocatable, dimension(:) :: lgs_a_value_trim
646 real(dp) :: eps_sor
647 
648 #elif (CALCTHK==6)
649 
650 lis_integer :: ierr
651 lis_integer :: iter
652 lis_integer :: nc, nr
653 lis_integer, parameter :: nmax = (imax+1)*(jmax+1)
654 lis_integer, parameter :: n_sprs = 10*(imax+1)*(jmax+1)
655 lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
656 lis_integer, allocatable, dimension(:) :: lgs_a_diag_index
657 lis_matrix :: lgs_a
658 lis_vector :: lgs_b, lgs_x
659 lis_scalar, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
660 lis_solver :: solver
661 character(len=256) :: ch_solver_set_option
662 
663 #endif
664 
665 !-------- Abbreviations --------
666 
667 dt_darea = dtime/area
668 
669 do i=0, imax
670 do j=0, jmax
671 
672  if (flag_inner_point(j,i)) then
673 
674  vx_m_1(j,i) = vx_m(j,i-1)
675  vx_m_2(j,i) = vx_m(j,i)
676  vy_m_1(j,i) = vy_m(j-1,i)
677  vy_m_2(j,i) = vy_m(j,i)
678 
679  if (vx_m_1(j,i) >= 0.0_dp) then
680  uph_x_1(j,i) = h(j,i-1)
681  else
682  uph_x_1(j,i) = h(j,i)
683  end if
684 
685  if (vx_m_2(j,i) >= 0.0_dp) then
686  uph_x_2(j,i) = h(j,i)
687  else
688  uph_x_2(j,i) = h(j,i+1)
689  end if
690 
691  if (vy_m_1(j,i) >= 0.0_dp) then
692  uph_y_1(j,i) = h(j-1,i)
693  else
694  uph_y_1(j,i) = h(j,i)
695  end if
696 
697  if (vy_m_2(j,i) >= 0.0_dp) then
698  uph_y_2(j,i) = h(j,i)
699  else
700  uph_y_2(j,i) = h(j+1,i)
701  end if
702 
703  sq_g22_x_1(j,i) = sq_g22_sgx(j,i-1)
704  sq_g22_x_2(j,i) = sq_g22_sgx(j,i)
705  sq_g11_y_1(j,i) = sq_g11_sgy(j-1,i)
706  sq_g11_y_2(j,i) = sq_g11_sgy(j,i)
707 
708  else ! .not.(flag_inner_point(j,i))
709 
710  vx_m_1(j,i) = 0.0_dp
711  vx_m_2(j,i) = 0.0_dp
712  vy_m_1(j,i) = 0.0_dp
713  vy_m_2(j,i) = 0.0_dp
714 
715  uph_x_1(j,i) = 0.0_dp
716  uph_x_2(j,i) = 0.0_dp
717  uph_y_1(j,i) = 0.0_dp
718  uph_y_2(j,i) = 0.0_dp
719 
720  sq_g22_x_1(j,i) = 0.0_dp
721  sq_g22_x_2(j,i) = 0.0_dp
722  sq_g11_y_1(j,i) = 0.0_dp
723  sq_g11_y_2(j,i) = 0.0_dp
724 
725  end if
726 
727 end do
728 end do
729 
730 !-------- Assembly of the system of linear equations
731 ! (matrix storage: compressed sparse row CSR) --------
732 
733 #if (CALCTHK==5 || CALCTHK==6)
734 
735 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
736 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
737 
738 lgs_a_value = 0.0_dp
739 lgs_a_index = 0
740 lgs_a_ptr = 0
741 
742 lgs_b_value = 0.0_dp
743 lgs_x_value = 0.0_dp
744 
745 lgs_a_ptr(1) = 1
746 
747 k = 0
748 
749 do nr=1, nmax ! loop over rows
750 
751  i = ii(nr)
752  j = jj(nr)
753 
754  if ( flag_inner_point(j,i) & ! inner point
755 #if (MARGIN==1 && MB_ACCOUNT==0)
756  .and.(maske(j,i) <= 1_i2b) & ! grounded ice or ice-free land
757 #endif
758  ) then
759 
760  k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc ! for H(j-1,i)
761  if (vy_m_1(j,i) > 0.0_dp) &
762  lgs_a_value(k) = -dt_darea(j,i)*vy_m_1(j,i) &
763  *sq_g11_y_1(j,i)*dxi*ovi_weight
764 
765  k=k+1 ; nc=nn(j,i-1) ; lgs_a_index(k)=nc ! for H(j,i-1)
766  if (vx_m_1(j,i) > 0.0_dp) &
767  lgs_a_value(k) = -dt_darea(j,i)*vx_m_1(j,i) &
768  *sq_g22_x_1(j,i)*deta*ovi_weight
769 
770  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k ! for H(j,i)
771  lgs_a_value(k) = 1.0_dp ! (diagonal element)
772  if (vy_m_1(j,i) < 0.0_dp) &
773  lgs_a_value(k) = lgs_a_value(k) &
774  - dt_darea(j,i)*vy_m_1(j,i) &
775  *sq_g11_y_1(j,i)*dxi*ovi_weight
776  if (vx_m_1(j,i) < 0.0_dp) &
777  lgs_a_value(k) = lgs_a_value(k) &
778  - dt_darea(j,i)*vx_m_1(j,i) &
779  *sq_g22_x_1(j,i)*deta*ovi_weight
780  if (vx_m_2(j,i) > 0.0_dp) &
781  lgs_a_value(k) = lgs_a_value(k) &
782  + dt_darea(j,i)*vx_m_2(j,i) &
783  *sq_g22_x_2(j,i)*deta*ovi_weight
784  if (vy_m_2(j,i) > 0.0_dp) &
785  lgs_a_value(k) = lgs_a_value(k) &
786  + dt_darea(j,i)*vy_m_2(j,i) &
787  *sq_g11_y_2(j,i)*dxi*ovi_weight
788 
789  k=k+1 ; nc=nn(j,i+1) ; lgs_a_index(k)=nc ! for H(j,i+1)
790  if (vx_m_2(j,i) < 0.0_dp) &
791  lgs_a_value(k) = dt_darea(j,i)*vx_m_2(j,i) &
792  *sq_g22_x_2(j,i)*deta*ovi_weight
793 
794  k=k+1 ; nc=nn(j+1,i) ; lgs_a_index(k)=nc ! for H(j+1,i)
795  if (vy_m_2(j,i) < 0.0_dp) &
796  lgs_a_value(k) = dt_darea(j,i)*vy_m_2(j,i) &
797  *sq_g11_y_2(j,i)*dxi*ovi_weight
798 
799  lgs_b_value(nr) = h(j,i) &
800  +dtime*mb_source(j,i) &
801  -(1.0_dp-ovi_weight) &
802  * dt_darea(j,i) &
803  * ( ( vx_m_2(j,i)*uph_x_2(j,i) &
804  *sq_g22_x_2(j,i)*deta &
805  -vx_m_1(j,i)*uph_x_1(j,i) &
806  *sq_g22_x_1(j,i)*deta ) &
807  + ( vy_m_2(j,i)*uph_y_2(j,i) &
808  *sq_g11_y_2(j,i)*dxi &
809  -vy_m_1(j,i)*uph_y_1(j,i) &
810  *sq_g11_y_1(j,i)*dxi ) )
811  ! right-hand side
812 
813  else ! zero-thickness boundary condition
814 
815  k = k+1
816  lgs_a_value(k) = 1.0_dp ! diagonal element only
817  lgs_a_diag_index(nr) = k
818  lgs_a_index(k) = nr
819  lgs_b_value(nr) = 0.0_dp
820 
821  end if
822 
823  lgs_x_value(nr) = h(j,i) ! old ice thickness,
824  ! initial guess for solution vector
825 
826  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
827 
828 end do
829 
830 nnz = k ! number of non-zero elements of the matrix
831 
832 #else
833 stop ' >>> calc_thk_impl: CALCTHK must be either 5 or 6!'
834 #endif
835 
836 !-------- Solution of the system of linear equations --------
837 
838 #if (CALCTHK==5)
839 
840 ! ------ Solution with the built-in SOR solver
841 
842 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
843 
844 do k=1, nnz ! relocate matrix to trimmed arrays
845  lgs_a_value_trim(k) = lgs_a_value(k)
846  lgs_a_index_trim(k) = lgs_a_index(k)
847 end do
848 
849 deallocate(lgs_a_value, lgs_a_index)
850 
851 eps_sor = 1.0e-05_dp*mean_accum*dtime ! convergence parameter
852 
853 call sor_sprs(lgs_a_value_trim, &
854  lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
855  lgs_b_value, &
856  nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
857 
858 do nr=1, nmax
859  i = ii(nr)
860  j = jj(nr)
861  h_neu(j,i) = lgs_x_value(nr)
862 end do
863 
864 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
865 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
866 
867 #elif (CALCTHK==6)
868 
869 ! ------ Settings for Lis
870 
871 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
872 call lis_vector_create(lis_comm_world, lgs_b, ierr)
873 call lis_vector_create(lis_comm_world, lgs_x, ierr)
874 
875 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
876 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
877 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
878 
879 do nr=1, nmax
880 
881  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
882  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
883  lgs_a_value(nc), lgs_a, ierr)
884  end do
885 
886  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
887  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
888 
889 end do
890 
891 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
892 call lis_matrix_assemble(lgs_a, ierr)
893 
894 ! ------ Solution with Lis
895 
896 call lis_solver_create(solver, ierr)
897 
898 ch_solver_set_option = '-i bicg -p ilu '// &
899  '-maxiter 1000 -tol 1.0e-04 -initx_zeros false'
900 
901 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
902 call chkerr(ierr)
903 
904 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
905 call chkerr(ierr)
906 
907 call lis_solver_get_iter(solver, iter, ierr)
908 write(6,'(10x,a,i0)') 'calc_thk_impl: iter = ', iter
909 
910 lgs_x_value = 0.0_dp
911 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
912 call lis_matrix_destroy(lgs_a, ierr)
913 call lis_vector_destroy(lgs_b, ierr)
914 call lis_vector_destroy(lgs_x, ierr)
915 call lis_solver_destroy(solver, ierr)
916 
917 do nr=1, nmax
918  i = ii(nr)
919  j = jj(nr)
920  h_neu(j,i) = lgs_x_value(nr)
921 end do
922 
923 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
924 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
925 
926 #endif
927 
928 !-------- Applying the so-far neglected source term --------
929 
930 #if (MB_ACCOUNT==1)
931 call apply_mb_source(dtime, z_mar)
932 #endif
933 
934 !-------- Adjusting the ice thickness, if needed --------
935 
936 call thk_adjust(time, dtime)
937 
938 end subroutine calc_thk_impl
939 
940 !-------------------------------------------------------------------------------
941 !> Ice thickness evolution due to the source term (surface mass balance,
942 !! basal mass balance and calving).
943 !! Determination of the several components of the mass balance.
944 !<------------------------------------------------------------------------------
945 subroutine apply_mb_source(dtime, z_mar)
947  implicit none
948 
949  real(dp), intent(in) :: dtime
950  real(dp), intent(in) :: z_mar
951 
952  real(dp), dimension(0:JMAX,0:IMAX) :: H_smb, H_adv
953  logical, dimension(0:JMAX,0:IMAX) :: flag_ocean_point
954 
955  if (flag_solver_explicit) then
957  ! define source term for the ice thickness equation
958  else
959  where (flag_inner_point) h_neu = h_neu - dtime*mb_source
960  ! subtract source term from the ice thickness solution
961  end if
962 
963 ! Store thickness as local variable for updating
964  h_smb = h_neu
965 
966 ! ### At the end, forbid appearance of ice past two outer grid points
967 ! ### (this relates to the one outer points in calc_top and, therefore, applies
968 ! ### for non-cyclic RB for domain boundary of topography only)
969  h_smb(0,:) = 0.0_dp; h_smb(1,:) = 0.0_dp
970  h_smb(jmax-1,:) = 0.0_dp; h_smb(jmax,:) = 0.0_dp
971  h_smb(:,0) = 0.0_dp; h_smb(:,1) = 0.0_dp
972  h_smb(:,imax-1) = 0.0_dp; h_smb(:,imax) = 0.0_dp
973 
974 ! Now, for inner points *two* steps from the boundaries
975 ! ### Apply mass balance source term ###
976  h_smb(2:jmax-2,2:imax-2) &
977  = max( h_neu(2:jmax-2,2:imax-2) + dtime*mb_source(2:jmax-2,2:imax-2), &
978  0.0_dp)
979 
980 ! Note: Still the ice mass which went past the continental shelf and to the second
981 ! outer points will by counted later in this scheme by using H_adv and H_smb. This
982 ! is why the next loop runs from 1 to IMAX and from 1 to JMAX.
983 
984 #if (MARGIN==1)
985  ! Ocean kill: if ocean point and ice exists, remove the ice
986  where (maske == 2_i2b .and. h_smb > eps_dp) h_smb = 0.0_dp
987 #else /* margin = 2 or 3 */
988  ! If outside continental shelf, remove the ice
989  where (zl0 < z_mar) h_smb = 0.0_dp !!! Is this OK? ... !!!
990 #endif
991 
992  ! Forbid negative thickness
993  where (h_smb < 0.0_dp) h_smb = 0.0_dp
994 
995  ! Set mass balance and melt quatities zero
996  mb_source_apl = 0.0_dp
997  as_perp_apl = 0.0_dp
998  runoff_top = 0.0_dp
999  disc_top = 0.0_dp
1000  runoff_bot = 0.0_dp
1001  mask_run = 0_i2b
1002 
1003  ! For inner points *one* step from the boundaries. All mass which went here by
1004  ! ice flow will be counted.
1005 
1006  h_adv = h_neu
1007 
1008  where (flag_inner_point)
1009 
1010  where (maske == 2_i2b)
1011  flag_ocean_point = .true.
1012  elsewhere
1013  flag_ocean_point = .false.
1014  end where
1015 
1016  where (h_smb >= eps_dp) ! ## CASE 1: glaciated point ##
1017 
1018  ! Actual mass balance, based on temporal change in ice sheet
1019  mb_source_apl = (h_smb-h_adv)/dtime
1020 
1021  ! Store melt quantities here for later accounting
1022  runoff_top = runoff
1025 
1026  ! Actual SMB
1028 
1029  ! Store melting mask value. This became an ice point. Regarded
1030  ! as land always and the melt will be attributed to land.
1031  mask_run = 1_i2b ! grounded ice
1032 
1033  elsewhere (h_smb < eps_dp .and. h_adv >= eps_dp)
1034  ! ## CASE 2: ice disappeared ##
1035 
1036  ! Actual mass balance, based on temporal change in ice sheet
1037  mb_source_apl = (h_smb-h_adv)/dtime
1038 
1039  ! Assume all hidden melt comes from the surface
1040  ! (since no basal calculations were made at these points)
1041  where (runoff+calv_grounded > 0.0_dp)
1042  ! share by the respective melting continents by their fractions
1045  elsewhere (runoff > 0.0_dp .and. calv_grounded <= 0.0_dp)
1046  ! only ordinary surface melt applied
1048  elsewhere (calv_grounded > 0.0_dp .and. runoff <= 0.0_dp)
1049  ! only surface discharge applied
1051  elsewhere
1052  ! safety measure (hidden would be small here anyhow)
1053  runoff_top = -0.5_dp*mb_source_apl
1054  disc_top = -0.5_dp*mb_source_apl
1055  end where
1056 
1057  ! Actual SMB
1059 
1060  ! Store melting mask value. This grid point is ice free now, therefore it
1061  ! can be either ocean or land (using actual maske value).
1062  where (flag_ocean_point) ! Ocean point (hidden melt)
1063  mask_run = -2_i2b
1064  elsewhere ! Land (hidden melt)
1065  mask_run = -1_i2b
1066  end where
1067 
1068  end where
1069 
1070  end where
1071 
1072 ! Update the ice thickness
1073  h_neu = h_smb
1074 
1075 end subroutine apply_mb_source
1076 
1077 !-------------------------------------------------------------------------------
1078 !> Adjustment of the newly computed ice thickness distribution.
1079 !<------------------------------------------------------------------------------
1080  subroutine thk_adjust(time, dtime)
1081 
1082  implicit none
1083 
1084  real(dp), intent(in) :: time, dtime
1085 
1086  real(dp), dimension(0:JMAX,0:IMAX) :: H_neu_tmp
1087  real(dp) :: dtime_inv
1088 
1089 !-------- Saving computed H_neu before any adjustments --------
1090 
1091  h_neu_tmp = h_neu
1092 
1093 !-------- Correct negative thickness values --------
1094 
1095  where (h_neu < 0.0_dp) h_neu = 0.0_dp
1096 
1097 !-------- Further adjustments --------
1098 
1099 #if (THK_EVOL==0)
1100 
1101 ! ------ No evolution of the ice thickness
1102 
1103  h_neu = h ! newly computed ice thickness is discarded
1104 
1105 #elif (THK_EVOL==1)
1106 
1107 ! ------ No adjustment, ice thickness evolves freely;
1108 ! thus nothing to be done
1109 
1110 #elif (THK_EVOL==2)
1111 
1112 ! ------ Adjustment due to prescribed target topography
1113 
1114  if (time >= time_target_topo_final) then
1115 
1116  h_neu = h_target
1117 
1118  else if (time >= time_target_topo_init) then
1119 
1120  target_topo_tau = (time_target_topo_final-time)/3.0_dp
1121 
1122  h_neu = ( target_topo_tau*h_neu + dtime*h_target ) &
1123  / ( target_topo_tau + dtime )
1124 
1125  end if
1126 
1127 #elif (THK_EVOL==3)
1128 
1129 ! ------ Adjustment due to prescribed target topography with relaxation time
1130 
1131  ! relaxation time target_topo_tau defined in routine sico_init
1132 
1133  h_neu = ( target_topo_tau*h_neu + dtime*h_target ) &
1134  / ( target_topo_tau + dtime )
1135 
1136 #elif (THK_EVOL==4)
1137 
1138 ! ------ Maximum ice extent constrained by prescribed mask
1139 
1140  where (maske_maxextent == 0_i2b) & ! not allowed to glaciate
1141  h_neu = 0.0_dp
1142 
1143 #else
1144  stop ' >>> thk_adjust: THK_EVOL must be between 0 and 4!'
1145 #endif
1146 
1147 !-------- Computation of the mass balance adjustment --------
1148 
1149  dtime_inv = 1.0_dp/dtime
1150 
1151  smb_corr = (h_neu-h_neu_tmp)*dtime_inv
1152 
1153 end subroutine thk_adjust
1154 
1155 !-------------------------------------------------------------------------------
1156 !> Update of the ice-land-ocean mask etc.
1157 !<------------------------------------------------------------------------------
1158 subroutine calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, &
1159  n_calc_thk_mask_update_aux)
1161 use topograd_m
1162 
1163 implicit none
1164 
1165 real(dp), intent(in) :: time
1166 real(dp), intent(in) :: dtime, dxi, deta
1167 real(dp), intent(in) :: z_sl, z_mar
1168 integer(i2b), intent(in) :: n_calc_thk_mask_update_aux
1169 
1170 integer(i4b) :: i, j
1171 real(dp) :: year_sec_inv
1172 real(dp), dimension(0:JMAX,0:IMAX) :: H_neu_tmp
1173 real(dp) :: dtime_inv
1174 
1175 !-------- Term abbreviations --------
1176 
1177 year_sec_inv = 1.0_dp/year_sec
1178 
1179 !-------- Saving computed H_neu before any modifications --------
1180 
1181 h_neu_tmp = h_neu
1182 
1183 !-------- Update of the mask --------
1184 
1185 dtime_inv = 1.0_dp/dtime
1186 
1187 if (n_calc_thk_mask_update_aux == 1_i2b) then
1188  call calc_thk_mask_update_aux1(time, dtime, dtime_inv, z_sl, z_mar)
1189 else if (n_calc_thk_mask_update_aux == 2_i2b) then
1190  call calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)
1191 else if (n_calc_thk_mask_update_aux == 3_i2b) then
1192  call calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)
1193 else
1194  write(6, fmt='(a)') ' >>> calc_thk_mask_update:'
1195  write(6, fmt='(a)') ' n_calc_thk_mask_update_aux has no valid value!'
1196  stop
1197 end if
1198 
1199 !-------- Time derivatives --------
1200 
1201 dzs_dtau = (zs_neu-zs)*dtime_inv
1202 dzb_dtau = (zb_neu-zb)*dtime_inv
1205 
1206 #if (THK_EVOL==2)
1207 if ( abs((time-time_target_topo_final)*year_sec_inv) < eps ) then
1208  dzs_dtau = 0.0_dp ! Introduced
1209  dzb_dtau = 0.0_dp ! by
1210  dzm_dtau = 0.0_dp ! Tatsuru Sato
1211  dh_c_dtau = 0.0_dp ! for
1212  dh_t_dtau = 0.0_dp ! stability reasons
1213 end if
1214 #endif
1215 
1216 !-------- New gradients --------
1217 
1218 #if (TOPOGRAD==0)
1219 call topograd_1(dxi, deta, 2)
1220 #elif (TOPOGRAD==1)
1221 call topograd_2(dxi, deta, 2)
1222 #endif
1223 
1224 !-------- Compute the volume flux across the CTS, am_perp --------
1225 
1226 #if (CALCMOD==1)
1227 
1228 do i=1, imax-1
1229 do j=1, jmax-1
1230 
1231  if ( ((maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b)) &
1232  .and.(n_cts(j,i) == 1_i2b)) then
1233 
1234  am_perp_st(j,i) = &
1235  0.5_dp*(vx_c(0,j,i)+vx_c(0,j,i-1))*dzm_dxi_g(j,i) &
1236  + 0.5_dp*(vy_c(0,j,i)+vy_c(0,j-1,i))*dzm_deta_g(j,i) &
1237  - 0.5_dp*(vz_c(0,j,i)+vz_t(ktmax-1,j,i))
1238  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
1239 
1240  else
1241  am_perp_st(j,i) = 0.0_dp
1242  am_perp(j,i) = 0.0_dp
1243  end if
1244 
1245 end do
1246 end do
1247 
1248 #endif
1249 
1250 !-------- Compute the applied surface mass balance as_perp_apl --------
1251 
1252 #if (MB_ACCOUNT==0)
1253 
1254 where ( (maske==0_i2b).or.(maske_old==0_i2b) &
1255  .or.(maske==3_i2b).or.(maske_old==3_i2b) )
1256  ! grounded or floating ice before or after the time step
1257 
1259 
1260 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1261 
1262  as_perp_apl = 0.0_dp
1263 
1264 end where
1265 
1266 #endif
1267 
1268 end subroutine calc_thk_mask_update
1269 
1270 !-------------------------------------------------------------------------------
1271 !> Update of the ice-land-ocean mask for SIA-only dynamics of ice sheets
1272 !! without ice shelves.
1273 !<------------------------------------------------------------------------------
1274 subroutine calc_thk_mask_update_aux1(time, dtime, dtime_inv, z_sl, z_mar)
1276 implicit none
1277 
1278 real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1279 
1280 integer(i4b) :: i, j, kc, kt
1281 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1282 logical :: flag_calving_event
1283 
1284 real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1285 
1286 real(dp), parameter :: H_isol_max = 1.0e+03_dp
1287  ! Thickness limit of isolated glaciated points
1288 
1289 !-------- Term abbreviations --------
1290 
1291 rhosw_rho_ratio = rho_sw/rho
1292 rho_rhosw_ratio = rho/rho_sw
1293 
1294 !-------- Update of the mask --------
1295 
1296 zs_neu = zb_neu + h_neu ! ice surface topography
1297 h_sea_neu = z_sl - zl_neu ! sea depth
1298 h_balance = h_sea_neu*rhosw_rho_ratio
1299 
1300 #if (THK_EVOL>=1)
1301 
1302 do i=0, imax
1303 do j=0, jmax
1304 
1305  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
1306 
1307  if (h_neu(j,i) >= eps_dp) then
1308  maske(j,i) = 0_i2b ! grounded ice
1309  else
1310  maske(j,i) = 1_i2b ! ice-free land
1311  end if
1312 
1313 #if (MARGIN==2)
1314 
1315  else ! (maske(j,i) == 2_i2b); sea
1316 
1317  if (h_neu(j,i) >= eps_dp) then
1318 
1319  if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then
1320 
1321 #if (MARINE_ICE_FORMATION==1)
1322  maske(j,i) = 2_i2b ! floating ice cut off -> sea
1323 #elif (MARINE_ICE_FORMATION==2)
1324  maske(j,i) = 0_i2b ! "underwater ice"
1325 #endif
1326  else
1327  maske(j,i) = 0_i2b ! grounded ice
1328  end if
1329 
1330 #if (MARINE_ICE_CALVING==2 \
1331  || marine_ice_calving==4 \
1332  || marine_ice_calving==6)
1333  if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1334 #elif (MARINE_ICE_CALVING==3 \
1335  || marine_ice_calving==5 \
1336  || marine_ice_calving==7)
1337  if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1338 #endif
1339 
1340  else
1341  maske(j,i) = 2_i2b ! sea
1342  end if
1343 
1344 #endif
1345 
1346  end if
1347 
1348 end do
1349 end do
1350 
1351 #endif
1352 
1353 ! ------ Adjustment due to prescribed target topography
1354 
1355 #if (THK_EVOL==2)
1357 #endif
1358 
1359 !-------- Detection of the grounding line and the calving front
1360 ! (not relevant for SIA-only dynamics) --------
1361 
1362 flag_grounding_line_1 = .false.
1363 flag_grounding_line_2 = .false.
1364 
1365 flag_calving_front_1 = .false.
1366 flag_calving_front_2 = .false.
1367 
1368 !-------- Correction of zs_neu and zb_neu for ice-free land and sea --------
1369 
1370 do i=0, imax
1371 do j=0, jmax
1372 
1373  if (maske(j,i) == 1_i2b) then ! ice-free land
1374 
1375  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1376  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1377 
1378  else if (maske(j,i) == 2_i2b) then ! sea
1379 
1380  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1381  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1382 
1383  else if (maske(j,i) == 3_i2b) then ! floating ice
1384 
1385  write(6, fmt='(a)') ' >>> calc_thk_mask_update_aux1:'
1386  write(6, fmt='(a)') ' maske(j,i)==3 not allowed for'
1387  write(6, fmt='(a)') ' SIA-only dynamics!'
1388  stop
1389 
1390  end if
1391 
1392 end do
1393 end do
1394 
1395 !-------- Computation of further quantities --------
1396 
1397 do i=0, imax
1398 do j=0, jmax
1399 
1400  if (maske(j,i) == 0_i2b) then ! grounded ice
1401 
1402 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
1403 
1404  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
1405  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
1406  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
1407  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
1408  ) then ! all nearest neighbours ice-free
1409  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
1410  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1411  end if
1412 
1413 ! ------ Further stuff
1414 
1415  if (n_cts(j,i) == 1) then
1416  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
1417  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
1418  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1419  else
1420  h_t_neu(j,i) = 0.0_dp
1421  zm_neu(j,i) = zb_neu(j,i)
1422  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1423  end if
1424 
1425  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
1426 
1427  zs_neu(j,i) = zb_neu(j,i)
1428  zm_neu(j,i) = zb_neu(j,i)
1429  h_c_neu(j,i) = 0.0_dp
1430  h_t_neu(j,i) = 0.0_dp
1431  h_neu(j,i) = 0.0_dp
1432 
1433  end if
1434 
1435 end do
1436 end do
1437 
1438 end subroutine calc_thk_mask_update_aux1
1439 
1440 !-------------------------------------------------------------------------------
1441 !> Update of the ice-land-ocean mask for hybrid SIA/SStA dynamics of ice sheets
1442 !! without ice shelves.
1443 !<------------------------------------------------------------------------------
1444 subroutine calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)
1445 
1446 implicit none
1447 
1448 real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1449 
1450 integer(i4b) :: i, j, kc, kt
1451 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1452 logical :: flag_calving_event
1453 
1454 real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1455 
1456 real(dp), parameter :: H_isol_max = 1.0e+03_dp
1457  ! Thickness limit of isolated glaciated points
1458 
1459 !-------- Term abbreviations --------
1460 
1461 rhosw_rho_ratio = rho_sw/rho
1462 rho_rhosw_ratio = rho/rho_sw
1463 
1464 !-------- Update of the mask --------
1465 
1466 zs_neu = zb_neu + h_neu ! ice surface topography
1467 h_sea_neu = z_sl - zl_neu ! sea depth
1468 h_balance = h_sea_neu*rhosw_rho_ratio
1469 
1470 #if (THK_EVOL>=1)
1471 
1472 do i=0, imax
1473 do j=0, jmax
1474 
1475  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
1476 
1477  if (h_neu(j,i) >= eps_dp) then
1478  maske(j,i) = 0_i2b ! grounded ice
1479  else
1480  maske(j,i) = 1_i2b ! ice-free land
1481  end if
1482 
1483 #if (MARGIN==2)
1484 
1485  else ! (maske(j,i) == 2_i2b); sea
1486 
1487  if (h_neu(j,i) >= eps_dp) then
1488 
1489  if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then
1490 
1491 #if (MARINE_ICE_FORMATION==1)
1492  maske(j,i) = 2_i2b ! floating ice cut off -> sea
1493 #elif (MARINE_ICE_FORMATION==2)
1494  maske(j,i) = 0_i2b ! "underwater ice"
1495 #endif
1496  else
1497  maske(j,i) = 0_i2b ! grounded ice
1498  end if
1499 
1500 #if (MARINE_ICE_CALVING==2 \
1501  || marine_ice_calving==4 \
1502  || marine_ice_calving==6)
1503  if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1504 #elif (MARINE_ICE_CALVING==3 \
1505  || marine_ice_calving==5 \
1506  || marine_ice_calving==7)
1507  if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1508 #endif
1509 
1510  else
1511  maske(j,i) = 2_i2b ! sea
1512  end if
1513 
1514 #endif
1515 
1516  end if
1517 
1518 end do
1519 end do
1520 
1521 #endif
1522 
1523 ! ------ Adjustment due to prescribed target topography
1524 
1525 #if (THK_EVOL==2)
1527 #endif
1528 
1529 !-------- Detection of the grounding line and the calving front
1530 ! (not relevant for hybrid SIA/SStA dynamics
1531 ! of ice sheets without ice shelves) --------
1532 
1533 flag_grounding_line_1 = .false.
1534 flag_grounding_line_2 = .false.
1535 
1536 flag_calving_front_1 = .false.
1537 flag_calving_front_2 = .false.
1538 
1539 !-------- Correction of zs_neu and zb_neu for ice-free land and sea --------
1540 
1541 do i=0, imax
1542 do j=0, jmax
1543 
1544  if (maske(j,i) == 1_i2b) then ! ice-free land
1545 
1546  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1547  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1548 
1549  else if (maske(j,i) == 2_i2b) then ! sea
1550 
1551  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1552  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1553 
1554  else if (maske(j,i) == 3_i2b) then ! floating ice
1555 
1556  write(6, fmt='(a)') ' >>> calc_thk_mask_update_aux2:'
1557  write(6, fmt='(a)') ' maske(j,i)==3 not allowed for'
1558  write(6, fmt='(a)') ' hybrid SIA/SStA dynamics of ice sheets'
1559  write(6, fmt='(a)') ' without ice shelves!'
1560  stop
1561 
1562  end if
1563 
1564 end do
1565 end do
1566 
1567 !-------- Computation of further quantities --------
1568 
1569 do i=0, imax
1570 do j=0, jmax
1571 
1572  if (maske(j,i) == 0_i2b) then ! grounded ice
1573 
1574 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
1575 
1576  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
1577  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
1578  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
1579  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
1580  ) then ! all nearest neighbours ice-free
1581  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
1582  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1583  end if
1584 
1585 ! ------ Further stuff
1586 
1587  if (n_cts(j,i) == 1) then
1588  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
1589  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
1590  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1591  else
1592  h_t_neu(j,i) = 0.0_dp
1593  zm_neu(j,i) = zb_neu(j,i)
1594  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1595  end if
1596 
1597  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
1598 
1599  zs_neu(j,i) = zb_neu(j,i)
1600  zm_neu(j,i) = zb_neu(j,i)
1601  h_c_neu(j,i) = 0.0_dp
1602  h_t_neu(j,i) = 0.0_dp
1603  h_neu(j,i) = 0.0_dp
1604 
1605  end if
1606 
1607 end do
1608 end do
1609 
1610 end subroutine calc_thk_mask_update_aux2
1611 
1612 !-------------------------------------------------------------------------------
1613 !> Update of the ice-land-ocean mask for coupled SIA/SSA or
1614 !! SIA/SStA/SSA dynamics of ice sheets with ice shelves.
1615 !<------------------------------------------------------------------------------
1616 subroutine calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)
1617 
1618 implicit none
1619 
1620 real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1621 
1622 integer(i4b) :: i, j, kc, kt
1623 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1624 logical :: flag_calving_event
1625 
1626 real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1627 
1628 real(dp), parameter :: H_isol_max = 1.0e+03_dp
1629  ! Thickness limit of isolated glaciated points
1630 
1631 !-------- Term abbreviations --------
1632 
1633 rhosw_rho_ratio = rho_sw/rho
1634 rho_rhosw_ratio = rho/rho_sw
1635 
1636 !-------- Update of the mask --------
1637 
1638 zs_neu = zb_neu + h_neu ! ice surface topography
1639 h_sea_neu = z_sl - zl_neu ! sea depth
1640 h_balance = h_sea_neu*rhosw_rho_ratio
1641 
1642 #if (THK_EVOL>=1)
1643 
1644 do i=1, imax-1
1645 do j=1, jmax-1
1646 
1647  ! grounding_line migration check
1648 
1649  if ( ( maske(j,i) <= 1_i2b &
1650  .and. (maske(j,i+1)>1_i2b.or.maske(j,i-1)>1_i2b &
1651  .or.maske(j+1,i)>1_i2b.or.maske(j-1,i)>1_i2b) ) &
1652  .or. &
1653  ( maske(j,i)>=2_i2b &
1654  .and. (maske(j,i+1)==0_i2b.or.maske(j,i-1)==0_i2b &
1655  .or.maske(j+1,i)==0_i2b.or.maske(j-1,i)==0_i2b) ) ) then
1656 
1657  if (h_neu(j,i) >= eps_dp) then
1658 
1659  if (h_neu(j,i)<h_balance(j,i).and.zl_neu(j,i)<z_sl) then
1660  maske(j,i) = 3_i2b
1661  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_neu(j,i)
1662  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
1663  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1664  else
1665  maske(j,i) = 0_i2b
1666  zb_neu(j,i) = zl_neu(j,i)
1667  dzb_dtau(j,i) = dzl_dtau(j,i)
1668  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1669  end if
1670 
1671  else ! if (H_neu(j,i) <= eps_dp) then
1672 
1673  if (zl_neu(j,i)>z_sl) then
1674  maske(j,i) = 1_i2b
1675  zb_neu(j,i) = zl_neu(j,i)
1676  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
1677  zs_neu(j,i) = zb_neu(j,i)
1678  else
1679  maske(j,i) = 2_i2b
1680  zb_neu(j,i) = z_sl
1681  dzb_dtau(j,i) = 0.0_dp
1682  zs_neu(j,i) = z_sl
1683  end if
1684 
1685  end if
1686 
1687  else if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
1688 
1689  if (h_neu(j,i) >= eps_dp) then ! can change maske 0 or 1
1690  maske(j,i) = 0_i2b
1691  zb_neu(j,i) = zl_neu(j,i)
1692  dzb_dtau(j,i) = dzl_dtau(j,i)
1693  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1694  else
1695  maske(j,i) = 1_i2b
1696  zb_neu(j,i) = zl_neu(j,i)
1697  dzb_dtau(j,i) = dzl_dtau(j,i)
1698  zs_neu(j,i) = zl_neu(j,i)
1699  end if
1700 
1701  else ! if (maske(j,i)==2.or.maske(j,i)==3) then
1702 
1703  if (h_neu(j,i) > eps_dp) then
1704 
1705  if (h_neu(j,i)<h_balance(j,i).and.zl_neu(j,i)<z_sl) then
1706  maske(j,i) = 3_i2b
1707  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_neu(j,i)
1708  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
1709  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1710  else
1711  maske(j,i) = 0_i2b
1712  zb_neu(j,i) = zl_neu(j,i)
1713  dzb_dtau(j,i) = dzl_dtau(j,i)
1714  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1715  end if
1716 
1717  else ! if (H_neu(j,i) <= eps_dp) then
1718 
1719  if (zl_neu(j,i)>z_sl) then
1720  maske(j,i) = 1_i2b
1721  zb_neu(j,i) = zl_neu(j,i)
1722  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
1723  zs_neu(j,i) = zb_neu(j,i)
1724  else
1725  maske(j,i) = 2_i2b
1726  zb_neu(j,i) = z_sl
1727  dzb_dtau(j,i) = 0.0_dp
1728  zs_neu(j,i) = z_sl
1729  end if
1730 
1731  end if
1732 
1733  end if
1734 
1735 end do
1736 end do
1737 
1738 #if (ICE_SHELF_CALVING==2)
1739 
1740 do
1741 
1742  flag_calving_front_1 = .false.
1743  flag_calving_event = .false.
1744 
1745  do i=1, imax-1
1746  do j=1, jmax-1
1747 
1748  if ( (maske(j,i)==3_i2b) & ! floating ice
1749  .and. &
1750  ( (maske(j,i+1)==2_i2b) & ! with
1751  .or.(maske(j,i-1)==2_i2b) & ! one
1752  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1753  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
1754  ) &
1755  flag_calving_front_1(j,i) = .true. ! preliminary detection
1756  ! of the calving front
1757 
1758  if ( flag_calving_front_1(j,i).and.(h_neu(j,i) < h_calv) ) then
1759  flag_calving_event = .true. ! calving event,
1760  maske(j,i) = 2_i2b ! floating ice point changes to sea point
1761  end if
1762 
1763  end do
1764  end do
1765 
1766  if (.not.flag_calving_event) exit
1767 
1768 end do
1769 
1770 #elif (ICE_SHELF_CALVING==9)
1771 
1772 #if (defined(MISMIPP))
1773 
1774 do i=0, imax
1775 do j=0, jmax
1776 
1777  if ((maske(j,i)==3_i2b).and.(xi(i) > lx)) then
1778  maske(j,i) = 2_i2b ! floating ice point changes to sea point
1779  end if
1780 
1781 end do
1782 end do
1783 
1784 #else
1785 write(6, fmt='(a)') ' >>> calc_thk_mask_update_aux3:'
1786 write(6, fmt='(a)') ' Option ICE_SHELF_CALVING==9'
1787 write(6, fmt='(a)') ' only defined for MISMIP+!'
1788 stop
1789 #endif
1790 
1791 #endif
1792 
1793 #endif
1794 
1795 ! ------ Adjustment due to prescribed target topography
1796 
1797 #if (THK_EVOL==2)
1799 #endif
1800 
1801 !-------- Detection of the grounding line and the calving front --------
1802 
1803 flag_grounding_line_1 = .false.
1804 flag_grounding_line_2 = .false.
1805 
1806 flag_calving_front_1 = .false.
1807 flag_calving_front_2 = .false.
1808 
1809 do i=1, imax-1
1810 do j=1, jmax-1
1811 
1812  if ( (maske(j,i)==0_i2b) & ! grounded ice
1813  .and. &
1814  ( (maske(j,i+1)==3_i2b) & ! with
1815  .or.(maske(j,i-1)==3_i2b) & ! one
1816  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1817  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1818  ) &
1819  flag_grounding_line_1(j,i) = .true.
1820 
1821  if ( (maske(j,i)==3_i2b) & ! floating ice
1822  .and. &
1823  ( (maske(j,i+1)==0_i2b) & ! with
1824  .or.(maske(j,i-1)==0_i2b) & ! one
1825  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1826  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1827  ) &
1828  flag_grounding_line_2(j,i) = .true.
1829 
1830  if ( (maske(j,i)==3_i2b) & ! floating ice
1831  .and. &
1832  ( (maske(j,i+1)==2_i2b) & ! with
1833  .or.(maske(j,i-1)==2_i2b) & ! one
1834  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1835  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
1836  ) &
1837  flag_calving_front_1(j,i) = .true.
1838 
1839  if ( (maske(j,i)==2_i2b) & ! sea
1840  .and. &
1841  ( (maske(j,i+1)==3_i2b) & ! with
1842  .or.(maske(j,i-1)==3_i2b) & ! one
1843  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1844  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1845  ) &
1846  flag_calving_front_2(j,i) = .true.
1847 
1848 end do
1849 end do
1850 
1851 do i=1, imax-1
1852 
1853  j=0
1854  if ( (maske(j,i)==2_i2b) & ! sea
1855  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1856  ) & ! floating ice point
1857  flag_calving_front_2(j,i) = .true.
1858 
1859  j=jmax
1860  if ( (maske(j,i)==2_i2b) & ! sea
1861  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1862  ) & ! floating ice point
1863  flag_calving_front_2(j,i) = .true.
1864 
1865 end do
1866 
1867 do j=1, jmax-1
1868 
1869  i=0
1870  if ( (maske(j,i)==2_i2b) & ! sea
1871  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1872  ) & ! floating ice point
1873  flag_calving_front_2(j,i) = .true.
1874 
1875  i=imax
1876  if ( (maske(j,i)==2_i2b) & ! sea
1877  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1878  ) & ! floating ice point
1879  flag_calving_front_2(j,i) = .true.
1880 
1881 end do
1882 
1883 !-------- Correction of zs_neu and zb_neu
1884 ! for ice-free land, sea and floating ice --------
1885 
1886 do i=0, imax
1887 do j=0, jmax
1888 
1889  if (maske(j,i) == 1_i2b) then ! ice-free land
1890 
1891  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1892  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1893 
1894  else if (maske(j,i) == 2_i2b) then ! sea
1895 
1896  zs_neu(j,i) = z_sl ! both zs_neu(j,i) and zb_neu(j,i)
1897  zb_neu(j,i) = z_sl ! set to the current sea level stand
1898  h_neu(j,i) = 0.0_dp
1899 
1900  else if (maske(j,i) == 3_i2b) then ! floating ice
1901 
1902  h_neu(j,i) = zs_neu(j,i)-zb_neu(j,i) ! ice thickness
1903 
1904  zb_neu(j,i) = z_sl - rho_rhosw_ratio*h_neu(j,i) ! floating condition
1905  zs_neu(j,i) = zb_neu(j,i) + h_neu(j,i)
1906 
1907  end if
1908 
1909 end do
1910 end do
1911 
1912 !-------- Computation of further quantities --------
1913 
1914 do i=0, imax
1915 do j=0, jmax
1916 
1917  if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
1918  ! grounded or floating ice
1919 
1920 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
1921 
1922  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
1923  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
1924  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
1925  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
1926  ) then ! all nearest neighbours ice-free
1927  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
1928  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1929  end if
1930 
1931 ! ------ Further stuff
1932 
1933  if (n_cts(j,i) == 1) then
1934  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
1935  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
1936  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1937  else
1938  h_t_neu(j,i) = 0.0_dp
1939  zm_neu(j,i) = zb_neu(j,i)
1940  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1941  end if
1942 
1943  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
1944 
1945  zs_neu(j,i) = zb_neu(j,i)
1946  zm_neu(j,i) = zb_neu(j,i)
1947  h_c_neu(j,i) = 0.0_dp
1948  h_t_neu(j,i) = 0.0_dp
1949  h_neu(j,i) = 0.0_dp
1950 
1951  end if
1952 
1953 end do
1954 end do
1955 
1956 end subroutine calc_thk_mask_update_aux3
1957 
1958 !-------------------------------------------------------------------------------
1959 
1960 end module calc_thk_m
1961 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
Computation of the ice thickness.
Definition: calc_thk_m.F90:35
subroutine, public calc_thk_init()
Initialisations for the ice thickness computation.
Definition: calc_thk_m.F90:59
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, lgs_b_value, nnz, nmax, omega, eps_sor, lgs_x_value, ierr)
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:jmax, 0:imax) insq_g22_sgy
insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) vx_m
vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) mb_source_apl
mb_source_apl(j,i): Applied mass balance source (SMB, BMB, calving)
real(dp), parameter eps
eps: Small number
integer(i2b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
subroutine calc_thk_mask_update_aux1(time, dtime, dtime_inv, z_sl, z_mar)
Update of the ice-land-ocean mask for SIA-only dynamics of ice sheets without ice shelves...
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:jmax, 0:imax) insq_g11_sgx
insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
Several mathematical tools used by SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
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) sq_g22_sgx
sq_g22_sgx(j,i): Square root of g22, at (i+1/2,j)
subroutine, public calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
Over-implicit solver for the general ice thickness equation.
Definition: calc_thk_m.F90:615
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
logical, save flag_solver_explicit
Definition: calc_thk_m.F90:45
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
subroutine, public calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
Explicit solver for the diffusive SIA ice surface equation.
Definition: calc_thk_m.F90:137
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_m.F90:56
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
subroutine, public calc_thk_expl(time, dtime, dxi, deta, z_mar)
Explicit solver for the general ice thickness equation.
Definition: calc_thk_m.F90:501
real(dp), dimension(0:jmax, 0:imax) vy_m
vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax), save h
Definition: calc_thk_m.F90:43
real(dp), dimension(0:jmax, 0:imax) runoff_bot
runoff_bot(j,i): Runoff rate (accounting at ice bottom)
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
real(dp), dimension(0:jmax, 0:imax) h_diff
h_diff(j,i): Diffusivity of the SIA evolution equation of the ice surface
real(dp), dimension(0:jmax, 0:imax) runoff_top
runoff_top(j,i): Runoff rate (accounting at ice surface)
subroutine apply_mb_source(dtime, z_mar)
Ice thickness evolution due to the source term (surface mass balance, basal mass balance and calving)...
Definition: calc_thk_m.F90:946
subroutine, public calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, n_calc_thk_mask_update_aux)
Update of the ice-land-ocean mask etc.
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
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), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
Definition: topograd_m.F90:201
real(dp), dimension(0:jmax, 0:imax), save mb_source
Definition: calc_thk_m.F90:44
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:jmax, 0:imax) sq_g11_sgy
sq_g11_sgy(j,i): Square root of g11, at (i,j+1/2)
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
real(dp), dimension(0:jmax, 0:imax) zs_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax), save h_neu
Definition: calc_thk_m.F90:43
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) disc_top
disc_top(j,i): Ice discharge rate (accounting at ice surface)
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
Calculation of topography gradients on the staggered grid and on the grid points (including length re...
Definition: topograd_m.F90:39
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
real(dp) rho_sw
RHO_SW: Density of sea water.
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
real(dp), dimension(0:jmax, 0:imax) zb_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
integer(i2b), dimension(0:jmax, 0:imax) mask_run
mask_run(j,i): mask indicating melt type. 2: visible (ocean, for later developments), 1: visible (grounded ice), -1: hidden on land, -2: hidden in ocean
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:jmax, 0:imax) h_target
H_target(j,i): Target topography (ice thickness)
subroutine, public calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
Over-implicit solver for the diffusive SIA ice surface equation.
Definition: calc_thk_m.F90:216
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax) am_perp_st
am_perp_st(j,i): Steady-state part of am_perp (without contribution of dzm_dtau)
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske_target
maske_target(j,i): Target topography (ice-land-ocean mask)
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) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)