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