SICOPOLIS V5-dev  Revision 1279
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 logical :: flag_calving_event
1702 
1703 real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1704 
1705 real(dp), parameter :: H_isol_max = 1.0e+03_dp
1706  ! Thickness limit of isolated glaciated points
1707 
1708 !-------- Term abbreviations --------
1709 
1710 rhosw_rho_ratio = rho_sw/rho
1711 rho_rhosw_ratio = rho/rho_sw
1712 
1713 !-------- Update of the mask --------
1714 
1715 zs_neu = zb_neu + h_neu ! ice surface topography
1716 h_sea_neu = z_sl - zl_neu ! sea depth
1717 h_balance = h_sea_neu*rhosw_rho_ratio
1718 
1719 #if (THK_EVOL>=1)
1720 
1721 do i=0, imax
1722 do j=0, jmax
1723 
1724  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
1725 
1726  if (h_neu(j,i) >= eps_dp) then
1727  maske(j,i) = 0_i2b ! grounded ice
1728  else
1729  maske(j,i) = 1_i2b ! ice-free land
1730  end if
1731 
1732 #if (MARGIN==2)
1733 
1734  else ! (maske(j,i) == 2_i2b); sea
1735 
1736  if (h_neu(j,i) >= eps_dp) then
1737 
1738  if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then
1739 
1740 #if (MARINE_ICE_FORMATION==1)
1741  maske(j,i) = 2_i2b ! floating ice cut off -> sea
1742 #elif (MARINE_ICE_FORMATION==2)
1743  maske(j,i) = 0_i2b ! "underwater ice"
1744 #endif
1745  else
1746  maske(j,i) = 0_i2b ! grounded ice
1747  end if
1748 
1749 #if (MARINE_ICE_CALVING==2 \
1750  || marine_ice_calving==4 \
1751  || marine_ice_calving==6)
1752  if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1753 #elif (MARINE_ICE_CALVING==3 \
1754  || marine_ice_calving==5 \
1755  || marine_ice_calving==7)
1756  if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1757 #endif
1758 
1759  else
1760  maske(j,i) = 2_i2b ! sea
1761  end if
1762 
1763 #endif
1764 
1765  end if
1766 
1767 end do
1768 end do
1769 
1770 #endif
1771 
1772 ! ------ Adjustment due to prescribed target topography
1773 
1774 #if (THK_EVOL==2)
1776 #endif
1777 
1778 !-------- Detection of the grounding line and the calving front
1779 ! (not relevant for SIA-only dynamics) --------
1780 
1781 flag_grounding_line_1 = .false.
1782 flag_grounding_line_2 = .false.
1783 
1784 flag_calving_front_1 = .false.
1785 flag_calving_front_2 = .false.
1786 
1787 !-------- Correction of zs_neu and zb_neu for ice-free land and sea --------
1788 
1789 do i=0, imax
1790 do j=0, jmax
1791 
1792  if (maske(j,i) == 1_i2b) then ! ice-free land
1793 
1794  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1795  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1796 
1797  else if (maske(j,i) == 2_i2b) then ! sea
1798 
1799  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1800  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1801 
1802  else if (maske(j,i) == 3_i2b) then ! floating ice
1803 
1804  errormsg = ' >>> calc_thk_mask_update_aux1:' &
1805  // end_of_line &
1806  //' maske(j,i)==3 not allowed for' &
1807  // end_of_line &
1808  //' SIA-only dynamics!'
1809  call error(errormsg)
1810 
1811  end if
1812 
1813 end do
1814 end do
1815 
1816 !-------- Computation of further quantities --------
1817 
1818 do i=0, imax
1819 do j=0, jmax
1820 
1821  if (maske(j,i) == 0_i2b) then ! grounded ice
1822 
1823 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
1824 
1825  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
1826  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
1827  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
1828  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
1829  ) then ! all nearest neighbours ice-free
1830  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
1831  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
1832  end if
1833 
1834 ! ------ Further stuff
1835 
1836  if (n_cts(j,i) == 1) then
1837  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
1838  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
1839  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1840  else
1841  h_t_neu(j,i) = 0.0_dp
1842  zm_neu(j,i) = zb_neu(j,i)
1843  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1844  end if
1845 
1846  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
1847 
1848  zs_neu(j,i) = zb_neu(j,i)
1849  zm_neu(j,i) = zb_neu(j,i)
1850  h_c_neu(j,i) = 0.0_dp
1851  h_t_neu(j,i) = 0.0_dp
1852  h_neu(j,i) = 0.0_dp
1853 
1854  end if
1855 
1856 end do
1857 end do
1858 
1859 end subroutine calc_thk_mask_update_aux1
1860 
1861 !-------------------------------------------------------------------------------
1862 !> Update of the ice-land-ocean mask for hybrid SIA/SStA dynamics of ice sheets
1863 !! without ice shelves.
1864 !<------------------------------------------------------------------------------
1865 subroutine calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)
1866 
1867 implicit none
1868 
1869 real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1870 
1871 integer(i4b) :: i, j, kc, kt
1872 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1873 logical :: flag_calving_event
1874 
1875 real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1876 
1877 real(dp), parameter :: H_isol_max = 1.0e+03_dp
1878  ! Thickness limit of isolated glaciated points
1879 
1880 !-------- Term abbreviations --------
1881 
1882 rhosw_rho_ratio = rho_sw/rho
1883 rho_rhosw_ratio = rho/rho_sw
1884 
1885 !-------- Update of the mask --------
1886 
1887 zs_neu = zb_neu + h_neu ! ice surface topography
1888 h_sea_neu = z_sl - zl_neu ! sea depth
1889 h_balance = h_sea_neu*rhosw_rho_ratio
1890 
1891 #if (THK_EVOL>=1)
1892 
1893 do i=0, imax
1894 do j=0, jmax
1895 
1896  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
1897 
1898  if (h_neu(j,i) >= eps_dp) then
1899  maske(j,i) = 0_i2b ! grounded ice
1900  else
1901  maske(j,i) = 1_i2b ! ice-free land
1902  end if
1903 
1904 #if (MARGIN==2)
1905 
1906  else ! (maske(j,i) == 2_i2b); sea
1907 
1908  if (h_neu(j,i) >= eps_dp) then
1909 
1910  if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then
1911 
1912 #if (MARINE_ICE_FORMATION==1)
1913  maske(j,i) = 2_i2b ! floating ice cut off -> sea
1914 #elif (MARINE_ICE_FORMATION==2)
1915  maske(j,i) = 0_i2b ! "underwater ice"
1916 #endif
1917  else
1918  maske(j,i) = 0_i2b ! grounded ice
1919  end if
1920 
1921 #if (MARINE_ICE_CALVING==2 \
1922  || marine_ice_calving==4 \
1923  || marine_ice_calving==6)
1924  if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1925 #elif (MARINE_ICE_CALVING==3 \
1926  || marine_ice_calving==5 \
1927  || marine_ice_calving==7)
1928  if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
1929 #endif
1930 
1931  else
1932  maske(j,i) = 2_i2b ! sea
1933  end if
1934 
1935 #endif
1936 
1937  end if
1938 
1939 end do
1940 end do
1941 
1942 #endif
1943 
1944 ! ------ Adjustment due to prescribed target topography
1945 
1946 #if (THK_EVOL==2)
1948 #endif
1949 
1950 !-------- Detection of the grounding line and the calving front
1951 ! (not relevant for hybrid SIA/SStA dynamics
1952 ! of ice sheets without ice shelves) --------
1953 
1954 flag_grounding_line_1 = .false.
1955 flag_grounding_line_2 = .false.
1956 
1957 flag_calving_front_1 = .false.
1958 flag_calving_front_2 = .false.
1959 
1960 !-------- Correction of zs_neu and zb_neu for ice-free land and sea --------
1961 
1962 do i=0, imax
1963 do j=0, jmax
1964 
1965  if (maske(j,i) == 1_i2b) then ! ice-free land
1966 
1967  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1968  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1969 
1970  else if (maske(j,i) == 2_i2b) then ! sea
1971 
1972  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1973  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
1974 
1975  else if (maske(j,i) == 3_i2b) then ! floating ice
1976 
1977  errormsg = ' >>> calc_thk_mask_update_aux2:' &
1978  // end_of_line &
1979  //' maske(j,i)==3 not allowed for' &
1980  // end_of_line &
1981  //' hybrid SIA/SStA dynamics of ice sheets' &
1982  // end_of_line &
1983  //' without ice shelves!'
1984  call error(errormsg)
1985 
1986  end if
1987 
1988 end do
1989 end do
1990 
1991 !-------- Computation of further quantities --------
1992 
1993 do i=0, imax
1994 do j=0, jmax
1995 
1996  if (maske(j,i) == 0_i2b) then ! grounded ice
1997 
1998 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
1999 
2000  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
2001  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
2002  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
2003  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
2004  ) then ! all nearest neighbours ice-free
2005  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
2006  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2007  end if
2008 
2009 ! ------ Further stuff
2010 
2011  if (n_cts(j,i) == 1) then
2012  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
2013  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
2014  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
2015  else
2016  h_t_neu(j,i) = 0.0_dp
2017  zm_neu(j,i) = zb_neu(j,i)
2018  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
2019  end if
2020 
2021  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
2022 
2023  zs_neu(j,i) = zb_neu(j,i)
2024  zm_neu(j,i) = zb_neu(j,i)
2025  h_c_neu(j,i) = 0.0_dp
2026  h_t_neu(j,i) = 0.0_dp
2027  h_neu(j,i) = 0.0_dp
2028 
2029  end if
2030 
2031 end do
2032 end do
2033 
2034 end subroutine calc_thk_mask_update_aux2
2035 
2036 !-------------------------------------------------------------------------------
2037 !> Update of the ice-land-ocean mask for coupled SIA/SSA or
2038 !! SIA/SStA/SSA dynamics of ice sheets with ice shelves.
2039 !<------------------------------------------------------------------------------
2040 subroutine calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)
2041 
2042 implicit none
2043 
2044 real(dp), intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
2045 
2046 integer(i4b) :: i, j, kc, kt
2047 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
2048 logical :: flag_calving_event
2049 
2050 real(dp), dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
2051 
2052 real(dp), parameter :: H_isol_max = 1.0e+03_dp
2053  ! Thickness limit of isolated glaciated points
2054 
2055 !-------- Term abbreviations --------
2056 
2057 rhosw_rho_ratio = rho_sw/rho
2058 rho_rhosw_ratio = rho/rho_sw
2059 
2060 !-------- Update of the mask --------
2061 
2062 zs_neu = zb_neu + h_neu ! ice surface topography
2063 h_sea_neu = z_sl - zl_neu ! sea depth
2064 h_balance = h_sea_neu*rhosw_rho_ratio
2065 
2066 #if (THK_EVOL>=1)
2067 
2068 do i=1, imax-1
2069 do j=1, jmax-1
2070 
2071  ! grounding_line migration check
2072 
2073  if ( ( maske(j,i) <= 1_i2b &
2074  .and. (maske(j,i+1)>1_i2b.or.maske(j,i-1)>1_i2b &
2075  .or.maske(j+1,i)>1_i2b.or.maske(j-1,i)>1_i2b) ) &
2076  .or. &
2077  ( maske(j,i)>=2_i2b &
2078  .and. (maske(j,i+1)==0_i2b.or.maske(j,i-1)==0_i2b &
2079  .or.maske(j+1,i)==0_i2b.or.maske(j-1,i)==0_i2b) ) ) then
2080 
2081  if (h_neu(j,i) >= eps_dp) then
2082 
2083  if (h_neu(j,i)<h_balance(j,i).and.zl_neu(j,i)<z_sl) then
2084  maske(j,i) = 3_i2b
2085  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_neu(j,i)
2086  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
2087  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2088  else
2089  maske(j,i) = 0_i2b
2090  zb_neu(j,i) = zl_neu(j,i)
2091  dzb_dtau(j,i) = dzl_dtau(j,i)
2092  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2093  end if
2094 
2095  else ! if (H_neu(j,i) <= eps_dp) then
2096 
2097  if (zl_neu(j,i)>z_sl) then
2098  maske(j,i) = 1_i2b
2099  zb_neu(j,i) = zl_neu(j,i)
2100  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
2101  zs_neu(j,i) = zb_neu(j,i)
2102  else
2103  maske(j,i) = 2_i2b
2104  zb_neu(j,i) = z_sl
2105  dzb_dtau(j,i) = 0.0_dp
2106  zs_neu(j,i) = z_sl
2107  end if
2108 
2109  end if
2110 
2111  else if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
2112 
2113  if (h_neu(j,i) >= eps_dp) then ! can change maske 0 or 1
2114  maske(j,i) = 0_i2b
2115  zb_neu(j,i) = zl_neu(j,i)
2116  dzb_dtau(j,i) = dzl_dtau(j,i)
2117  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2118  else
2119  maske(j,i) = 1_i2b
2120  zb_neu(j,i) = zl_neu(j,i)
2121  dzb_dtau(j,i) = dzl_dtau(j,i)
2122  zs_neu(j,i) = zl_neu(j,i)
2123  end if
2124 
2125  else ! if (maske(j,i)==2.or.maske(j,i)==3) then
2126 
2127  if (h_neu(j,i) > eps_dp) then
2128 
2129  if (h_neu(j,i)<h_balance(j,i).and.zl_neu(j,i)<z_sl) then
2130  maske(j,i) = 3_i2b
2131  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_neu(j,i)
2132  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
2133  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2134  else
2135  maske(j,i) = 0_i2b
2136  zb_neu(j,i) = zl_neu(j,i)
2137  dzb_dtau(j,i) = dzl_dtau(j,i)
2138  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2139  end if
2140 
2141  else ! if (H_neu(j,i) <= eps_dp) then
2142 
2143  if (zl_neu(j,i)>z_sl) then
2144  maske(j,i) = 1_i2b
2145  zb_neu(j,i) = zl_neu(j,i)
2146  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
2147  zs_neu(j,i) = zb_neu(j,i)
2148  else
2149  maske(j,i) = 2_i2b
2150  zb_neu(j,i) = z_sl
2151  dzb_dtau(j,i) = 0.0_dp
2152  zs_neu(j,i) = z_sl
2153  end if
2154 
2155  end if
2156 
2157  end if
2158 
2159 end do
2160 end do
2161 
2162 #if (ICE_SHELF_CALVING==2)
2163 
2164 #ifndef ALLOW_OPENAD /* Normal */
2165 
2166 do
2167 
2168  flag_calving_front_1 = .false.
2169  flag_calving_event = .false.
2170 
2171  do i=1, imax-1
2172  do j=1, jmax-1
2173 
2174  if ( (maske(j,i)==3_i2b) & ! floating ice
2175  .and. &
2176  ( (maske(j,i+1)==2_i2b) & ! with
2177  .or.(maske(j,i-1)==2_i2b) & ! one
2178  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
2179  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
2180  ) &
2181  flag_calving_front_1(j,i) = .true. ! preliminary detection
2182  ! of the calving front
2183 
2184  if ( flag_calving_front_1(j,i).and.(h_neu(j,i) < h_calv) ) then
2185  flag_calving_event = .true. ! calving event,
2186  maske(j,i) = 2_i2b ! floating ice point changes to sea point
2187  end if
2188 
2189  end do
2190  end do
2191 
2192  if (.not.flag_calving_event) exit
2193 
2194 end do
2195 
2196 #else /* OpenAD */
2197 
2198 ! note: this block needs LIS to be tested
2199 ! just setting to .true. to get the loop started
2200  flag_calving_event = .true.
2201 
2202 do while (flag_calving_event)
2203 
2204  flag_calving_front_1 = .false.
2205  flag_calving_event = .false.
2206 
2207  do i=1, imax-1
2208  do j=1, jmax-1
2209 
2210  if ( (maske(j,i)==3_i2b) & ! floating ice
2211  .and. &
2212  ( (maske(j,i+1)==2_i2b) & ! with
2213  .or.(maske(j,i-1)==2_i2b) & ! one
2214  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
2215  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
2216  ) &
2217  flag_calving_front_1(j,i) = .true. ! preliminary detection
2218  ! of the calving front
2219 
2220  if ( (flag_calving_front_1(j,i)).and.(h_neu(j,i) < h_calv) ) then
2221  flag_calving_event = .true. ! calving event,
2222  maske(j,i) = 2_i2b ! floating ice point changes to sea point
2223  end if
2224 
2225  end do
2226  end do
2227 
2228 end do
2229 
2230 #endif /* Normal vs. OpenAD */
2231 
2232 #elif (ICE_SHELF_CALVING==3)
2233 
2234 do i=0, imax
2235 do j=0, jmax
2236 
2237  if (maske(j,i)==3_i2b) maske(j,i) = 2_i2b
2238  ! float-kill: all floating ice points changed to sea points
2239 
2240 end do
2241 end do
2242 
2243 #elif (ICE_SHELF_CALVING==9)
2244 
2245 #if (defined(MISMIPP))
2246 
2247 do i=0, imax
2248 do j=0, jmax
2249 
2250  if ((maske(j,i)==3_i2b).and.(xi(i) > lx)) then
2251  maske(j,i) = 2_i2b ! floating ice point changes to sea point
2252  end if
2253 
2254 end do
2255 end do
2256 
2257 #else
2258 
2259 errormsg = ' >>> calc_thk_mask_update_aux3:' &
2260  // end_of_line &
2261  //' Option ICE_SHELF_CALVING==9' &
2262  // end_of_line &
2263  //' only defined for MISMIP+!'
2264 call error(errormsg)
2265 
2266 #endif
2267 
2268 #endif
2269 
2270 #endif
2271 
2272 ! ------ Adjustment due to prescribed target topography
2273 
2274 #if (THK_EVOL==2)
2276 #endif
2277 
2278 !-------- Detection of the grounding line and the calving front --------
2279 
2280 flag_grounding_line_1 = .false.
2281 flag_grounding_line_2 = .false.
2282 
2283 flag_calving_front_1 = .false.
2284 flag_calving_front_2 = .false.
2285 
2286 do i=1, imax-1
2287 do j=1, jmax-1
2288 
2289  if ( (maske(j,i)==0_i2b) & ! grounded ice
2290  .and. &
2291  ( (maske(j,i+1)==3_i2b) & ! with
2292  .or.(maske(j,i-1)==3_i2b) & ! one
2293  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
2294  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
2295  ) &
2296  flag_grounding_line_1(j,i) = .true.
2297 
2298  if ( (maske(j,i)==3_i2b) & ! floating ice
2299  .and. &
2300  ( (maske(j,i+1)==0_i2b) & ! with
2301  .or.(maske(j,i-1)==0_i2b) & ! one
2302  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
2303  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
2304  ) &
2305  flag_grounding_line_2(j,i) = .true.
2306 
2307  if ( (maske(j,i)==3_i2b) & ! floating ice
2308  .and. &
2309  ( (maske(j,i+1)==2_i2b) & ! with
2310  .or.(maske(j,i-1)==2_i2b) & ! one
2311  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
2312  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
2313  ) &
2314  flag_calving_front_1(j,i) = .true.
2315 
2316  if ( (maske(j,i)==2_i2b) & ! sea
2317  .and. &
2318  ( (maske(j,i+1)==3_i2b) & ! with
2319  .or.(maske(j,i-1)==3_i2b) & ! one
2320  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
2321  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
2322  ) &
2323  flag_calving_front_2(j,i) = .true.
2324 
2325 end do
2326 end do
2327 
2328 do i=1, imax-1
2329 
2330  j=0
2331  if ( (maske(j,i)==2_i2b) & ! sea
2332  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
2333  ) & ! floating ice point
2334  flag_calving_front_2(j,i) = .true.
2335 
2336  j=jmax
2337  if ( (maske(j,i)==2_i2b) & ! sea
2338  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
2339  ) & ! floating ice point
2340  flag_calving_front_2(j,i) = .true.
2341 
2342 end do
2343 
2344 do j=1, jmax-1
2345 
2346  i=0
2347  if ( (maske(j,i)==2_i2b) & ! sea
2348  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
2349  ) & ! floating ice point
2350  flag_calving_front_2(j,i) = .true.
2351 
2352  i=imax
2353  if ( (maske(j,i)==2_i2b) & ! sea
2354  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
2355  ) & ! floating ice point
2356  flag_calving_front_2(j,i) = .true.
2357 
2358 end do
2359 
2360 !-------- Correction of zs_neu and zb_neu
2361 ! for ice-free land, sea and floating ice --------
2362 
2363 do i=0, imax
2364 do j=0, jmax
2365 
2366  if (maske(j,i) == 1_i2b) then ! ice-free land
2367 
2368  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
2369  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
2370 
2371  else if (maske(j,i) == 2_i2b) then ! sea
2372 
2373  zs_neu(j,i) = z_sl ! both zs_neu(j,i) and zb_neu(j,i)
2374  zb_neu(j,i) = z_sl ! set to the current sea level stand
2375  h_neu(j,i) = 0.0_dp
2376 
2377  else if (maske(j,i) == 3_i2b) then ! floating ice
2378 
2379  h_neu(j,i) = zs_neu(j,i)-zb_neu(j,i) ! ice thickness
2380 
2381  zb_neu(j,i) = z_sl - rho_rhosw_ratio*h_neu(j,i) ! floating condition
2382  zs_neu(j,i) = zb_neu(j,i) + h_neu(j,i)
2383 
2384  end if
2385 
2386 end do
2387 end do
2388 
2389 !-------- Computation of further quantities --------
2390 
2391 do i=0, imax
2392 do j=0, jmax
2393 
2394  if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
2395  ! grounded or floating ice
2396 
2397 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
2398 
2399  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
2400  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
2401  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
2402  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
2403  ) then ! all nearest neighbours ice-free
2404  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
2405  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
2406  end if
2407 
2408 ! ------ Further stuff
2409 
2410  if (n_cts(j,i) == 1) then
2411  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
2412  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
2413  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
2414  else
2415  h_t_neu(j,i) = 0.0_dp
2416  zm_neu(j,i) = zb_neu(j,i)
2417  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
2418  end if
2419 
2420  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
2421 
2422  zs_neu(j,i) = zb_neu(j,i)
2423  zm_neu(j,i) = zb_neu(j,i)
2424  h_c_neu(j,i) = 0.0_dp
2425  h_t_neu(j,i) = 0.0_dp
2426  h_neu(j,i) = 0.0_dp
2427 
2428  end if
2429 
2430 end do
2431 end do
2432 
2433 end subroutine calc_thk_mask_update_aux3
2434 
2435 !-------------------------------------------------------------------------------
2436 !> Enforce connectivity of the ocean.
2437 !<------------------------------------------------------------------------------
2438 subroutine ocean_connect()
2439 
2440 implicit none
2441 
2442 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal only for now */
2443 
2444 integer(i4b) :: i, j
2445 integer(i1b), dimension(0:JMAX,0:IMAX) :: mask_connect
2446 integer(i1b), dimension(0:JMAX,0:IMAX) :: mask_connect_save, mask_connect_diff
2447 logical :: flag_change
2448 
2449 !-------- Determine connected area allowed to be ocean --------
2450 
2451 mask_connect = 0_i1b
2452 
2453 mask_connect(0:1 , : ) = 1_i1b
2454 mask_connect(jmax-1:jmax , : ) = 1_i1b
2455 mask_connect(: , 0:1 ) = 1_i1b
2456 mask_connect(: , imax-1:imax) = 1_i1b
2457 
2458 flag_change = .true.
2459 
2460 do while (flag_change)
2461 
2462  mask_connect_save = mask_connect
2463 
2464  do i=1, imax-1
2465  do j=1, jmax-1
2466 
2467  if (mask_connect_save(j,i) == 1_i1b) then
2468  if (maske(j ,i+1) >= 2_i2b) mask_connect(j ,i+1) = 1_i1b
2469  if (maske(j ,i-1) >= 2_i2b) mask_connect(j ,i-1) = 1_i1b
2470  if (maske(j+1,i ) >= 2_i2b) mask_connect(j+1,i ) = 1_i1b
2471  if (maske(j-1,i ) >= 2_i2b) mask_connect(j-1,i ) = 1_i1b
2472  if (maske(j+1,i+1) >= 2_i2b) mask_connect(j+1,i+1) = 1_i1b
2473  if (maske(j+1,i-1) >= 2_i2b) mask_connect(j+1,i-1) = 1_i1b
2474  if (maske(j-1,i+1) >= 2_i2b) mask_connect(j-1,i+1) = 1_i1b
2475  if (maske(j-1,i-1) >= 2_i2b) mask_connect(j-1,i-1) = 1_i1b
2476  end if
2477 
2478  end do
2479  end do
2480 
2481  mask_connect_diff = abs(mask_connect-mask_connect_save)
2482 
2483  if (maxval(mask_connect_diff) > 0_i1b) then
2484  flag_change = .true.
2485  else
2486  flag_change = .false.
2487  end if
2488 
2489 end do
2490 
2491 !-------- Reset disconnected "ocean islands" to ice-free land --------
2492 
2493 where ((maske == 2_i2b).and.(mask_connect == 0_i1b)) maske = 1_i2b
2494 
2495 #endif /* Normal only for now */
2496 
2497 end subroutine ocean_connect
2498 
2499 !-------------------------------------------------------------------------------
2500 
2501 end module calc_thk_m
2502 !
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)