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