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