SICOPOLIS V5-dev  Revision 1279
calc_vxy_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ v x y _ m
4 !
5 !> @file
6 !!
7 !! Computation of the horizontal velocity vx, vy.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 Ralf Greve, Tatsuru Sato, Thomas Goelles, Jorge Bernales,
12 !! Fuyuki Saito
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of the horizontal velocity vx, vy.
35 !<------------------------------------------------------------------------------
36 module calc_vxy_m
37 
38  use sico_types_m
40  use sico_vars_m
41  use error_m
42 
43  implicit none
44 
45  private
47 
48 contains
49 
50 !-------------------------------------------------------------------------------
51 !> Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice
52 !! approximation.
53 !<------------------------------------------------------------------------------
54 subroutine calc_vxy_b_sia(time, z_sl)
55 
56 implicit none
57 
58 real(dp), intent(in) :: time, z_sl
59 
60 integer(i4b) :: i, j
61 integer(i4b) :: r_smw_1, s_smw_1, r_smw_2, s_smw_2
62 real(dp), dimension(0:JMAX,0:IMAX) :: gamma_slide_inv
63 real(dp) :: gamma_slide_inv_1, gamma_slide_inv_2
64 real(dp) :: smw_coeff_1, smw_coeff_2
65 real(dp), dimension(0:JMAX,0:IMAX) :: p_b, p_b_red, tau_b
66 real(dp) :: cvxy1, cvxy1a, cvxy1b, ctau1, ctau1a, ctau1b
67 real(dp) :: temp_diff
68 real(dp) :: c_Hw_slide, Hw0_slide, Hw0_slide_inv
69 real(dp) :: vh_max
70 real(dp) :: year_sec_inv, time_in_years
71 real(dp) :: ramp_up_factor
72 logical, dimension(0:JMAX,0:IMAX) :: sub_melt_flag
73 
74 year_sec_inv = 1.0_dp/year_sec
75 time_in_years = time*year_sec_inv
76 
77 !-------- Sliding-law coefficients --------
78 
79 #if (defined(ANT) && (!defined(SEDI_SLIDE) || SEDI_SLIDE==1))
80  ! Antarctica without soft sediment,
81  ! one sliding law for the entire modelling domain
82 
83 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
84 
85 do i=0, imax
86 do j=0, jmax
87  c_slide(j,i) = c_slide*year_sec_inv
88  p_weert(j,i) = p_weert
89  q_weert(j,i) = q_weert
90  gamma_slide_inv(j,i) = gamma_slide_inv_1
91  sub_melt_flag(j,i) = (gamma_slide >= eps)
92 end do
93 end do
94 
95 #elif (defined(ANT) && SEDI_SLIDE==2)
96  ! Antarctica with soft sediment,
97  ! different sliding laws for hard rock and soft sediment
98 
99 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
100 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)
101 
102 do i=0, imax
103 do j=0, jmax
104 
105  if (maske_sedi(j,i) /= 7) then ! no soft sediment
106 
107 #if (defined(TRANS_LENGTH_SL))
108  c_slide(j,i) = ( (1.0_dp-r_mask_sedi(j,i))*c_slide &
109  +r_mask_sedi(j,i) *c_slide_sedi ) &
110  *year_sec_inv
111 #else
112  c_slide(j,i) = c_slide*year_sec_inv
113 #endif
114 
115  p_weert(j,i) = p_weert
116  q_weert(j,i) = q_weert
117  gamma_slide_inv(j,i) = gamma_slide_inv_1
118  sub_melt_flag(j,i) = (gamma_slide >= eps)
119 
120  else ! maske_sedi(j,i) == 7, soft sediment
121 
122 #if (defined(TRANS_LENGTH_SL))
123  c_slide(j,i) = ( (1.0_dp-r_mask_sedi(j,i))*c_slide &
124  +r_mask_sedi(j,i) *c_slide_sedi ) &
125  *year_sec_inv
126 #else
127  c_slide(j,i) = c_slide_sedi*year_sec_inv
128 #endif
129 
130  p_weert(j,i) = p_weert_sedi
131  q_weert(j,i) = q_weert_sedi
132  gamma_slide_inv(j,i) = gamma_slide_inv_2
133  sub_melt_flag(j,i) = (gamma_slide_sedi >= eps)
134 
135  end if
136 
137 end do
138 end do
139 
140 #elif (defined(GRL) && ICE_STREAM==1 \
141  || defined(asf) && ice_stream==1)
142  ! Greenland or Austfonna without ice streams,
143  ! one sliding law for the entire modelling domain,
144  ! with surface meltwater effect
145 
146 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
147 
148 r_smw_1 = r_smw
149 s_smw_1 = s_smw
150 smw_coeff_1 = smw_coeff*year_sec**s_smw_1
151 
152 do i=0, imax
153 do j=0, jmax
154  c_slide(j,i) = c_slide*year_sec_inv &
155  *(1.0_dp+smw_coeff_1*runoff(j,i)**s_smw_1 &
156  /max((h_c(j,i)+h_t(j,i))**r_smw_1, eps))
157  p_weert(j,i) = p_weert
158  q_weert(j,i) = q_weert
159  gamma_slide_inv(j,i) = gamma_slide_inv_1
160  sub_melt_flag(j,i) = (gamma_slide >= eps)
161 end do
162 end do
163 
164 #elif (defined(GRL) && ICE_STREAM==2 \
165  || defined(asf) && ice_stream==2)
166  ! Greenland or Austfonna with ice streams,
167  ! different sliding laws for hard rock and ice streams
168 
169 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
170 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)
171 
172 r_smw_1 = r_smw
173 s_smw_1 = s_smw
174 smw_coeff_1 = smw_coeff *year_sec**s_smw_1
175 r_smw_2 = r_smw_sedi
176 s_smw_2 = s_smw_sedi
177 smw_coeff_2 = smw_coeff_sedi*year_sec**s_smw_2
178 
179 do i=0, imax
180 do j=0, jmax
181 
182  if (maske_sedi(j,i) /= 7) then ! no ice stream
183  c_slide(j,i) = c_slide*year_sec_inv &
184  *(1.0_dp+smw_coeff_1*runoff(j,i)**s_smw_1 &
185  /max((h_c(j,i)+h_t(j,i))**r_smw_1, eps))
186  p_weert(j,i) = p_weert
187  q_weert(j,i) = q_weert
188  gamma_slide_inv(j,i) = gamma_slide_inv_1
189  sub_melt_flag(j,i) = (gamma_slide >= eps)
190  else ! maske_sedi(j,i) == 7, ice stream
191  c_slide(j,i) = c_slide_sedi*year_sec_inv &
192  *(1.0_dp+smw_coeff_2*runoff(j,i)**s_smw_2 &
193  /max((h_c(j,i)+h_t(j,i))**r_smw_2, eps))
194  p_weert(j,i) = p_weert_sedi
195  q_weert(j,i) = q_weert_sedi
196  gamma_slide_inv(j,i) = gamma_slide_inv_2
197  sub_melt_flag(j,i) = (gamma_slide_sedi >= eps)
198  end if
199 
200 end do
201 end do
202 
203 #elif (defined(XYZ) && defined(HEINO))
204  ! ISMIP HEINO,
205  ! different sliding laws for hard rock and soft sediment
206 
207 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
208 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)
209 
210 do i=0, imax
211 do j=0, jmax
212 
213  if (maske_sedi(j,i) /= 7) then ! no soft sediment
214  c_slide(j,i) = c_slide*year_sec_inv
215  p_weert(j,i) = p_weert
216  q_weert(j,i) = q_weert
217  gamma_slide_inv(j,i) = gamma_slide_inv_1
218  sub_melt_flag(j,i) = (gamma_slide >= eps)
219  else ! maske_sedi(j,i) == 7, soft sediment
220  c_slide(j,i) = c_slide_sedi*year_sec_inv
221  p_weert(j,i) = p_weert_sedi
222  q_weert(j,i) = q_weert_sedi
223  gamma_slide_inv(j,i) = gamma_slide_inv_2
224  sub_melt_flag(j,i) = (gamma_slide_sedi >= eps)
225  end if
226 
227 end do
228 end do
229 
230 #elif (defined(SCAND) \
231  || defined(nhem) \
232  || defined(tibet) \
233  || defined(nmars) \
234  || defined(smars) \
235  || defined(emtp2sge) \
236  || defined(xyz))
237  ! one sliding law for the entire modelling domain
238 
239 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
240 
241 do i=0, imax
242 do j=0, jmax
243  c_slide(j,i) = c_slide*year_sec_inv
244  p_weert(j,i) = p_weert
245  q_weert(j,i) = q_weert
246  gamma_slide_inv(j,i) = gamma_slide_inv_1
247  sub_melt_flag(j,i) = (gamma_slide >= eps)
248 end do
249 end do
250 
251 #else
252 
253 errormsg = ' >>> calc_vxy_b_sia: No valid domain specified!'
254 call error(errormsg)
255 
256 #endif
257 
258 do i=0, imax
259 do j=0, jmax
260  p_weert_inv(j,i) = 1.0_dp/max(real(p_weert(j,i),dp), eps)
261 end do
262 end do
263 
264 ! ------ Ramping up basal sliding
265 
266 ramp_up_factor = 1.0_dp
267 
268 #if (defined(TIME_RAMP_UP_SLIDE))
269 
270 if ( (time_ramp_up_slide > eps_dp) &
271  .and. &
272  ((time_in_years-(time_init0)) < (time_ramp_up_slide)) ) then
273 
274  ramp_up_factor = (time_in_years-(time_init0))/(time_ramp_up_slide)
275 
276  ramp_up_factor = max(min(ramp_up_factor, 1.0_dp), 0.0_dp)
277  ! constrain to interval [0,1]
278 
279  ramp_up_factor = ramp_up_factor*ramp_up_factor*ramp_up_factor &
280  *(10.0_dp + ramp_up_factor &
281  *(-15.0_dp+6.0_dp*ramp_up_factor))
282  ! make transition smooth (quintic function)
283 
284  c_slide = c_slide * ramp_up_factor
285 
286 end if
287 
288 #endif
289 
290 ! ------ Coefficients for the contribution of the basal water layer
291 
292 #if (BASAL_HYDROLOGY==1)
293 
294 #if (defined(C_HW_SLIDE))
295  c_hw_slide = c_hw_slide
296 #else
297  c_hw_slide = 0.0_dp
298 #endif
299 
300 #if (defined(HW0_SLIDE))
301  hw0_slide = hw0_slide
302 #else
303  hw0_slide = 1.0_dp
304 #endif
305 
306  hw0_slide_inv = 1.0_dp/max(hw0_slide, eps_dp)
307 
308 #else
309  c_hw_slide = 0.0_dp ! dummy value
310  hw0_slide = 1.0_dp ! dummy value
311  hw0_slide_inv = 1.0_dp ! dummy value
312 #endif
313 
314 !-------- Computation of basal stresses --------
315 
316 ! ------ Basal pressure p_b, basal water pressure p_b_w,
317 ! reduced pressure p_b_red
318 
319 do i=0, imax
320 do j=0, jmax
321 
322  if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
323  ! grounded ice, or floating ice at the grounding line
324 
325  p_b(j,i) = max(rho*g*(h_c(j,i)+h_t(j,i)), 0.0_dp)
326  p_b_w(j,i) = rho_sw*g*max((z_sl-zb(j,i)), 0.0_dp)
327  p_b_red(j,i) = p_b(j,i)-p_b_w(j,i)
328  p_b_red(j,i) = max(p_b_red(j,i), red_pres_limit_fact*p_b(j,i))
329  ! in order to avoid very small values, which may lead to
330  ! huge sliding velocities
331 
332  else ! maske(j,i) == 1, 2 or 3 away from the grounding line
333 
334  p_b(j,i) = 0.0_dp
335  p_b_w(j,i) = 0.0_dp
336  p_b_red(j,i) = 0.0_dp
337 
338  end if
339 
340 end do
341 end do
342 
343 ! ------ Absolute value of the basal shear stress, tau_b
344 
345 do i=0, imax
346 do j=0, jmax
347 
348 #ifndef ALLOW_OPENAD /* Normal */
349 
350  tau_b(j,i) = p_b(j,i)*sqrt(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)
351 
352 #else /* OpenAD: guarding against non-differentiable sqrt(0) */
353 
354  if ((dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2) > 0) then
355  tau_b(j,i) = p_b(j,i)*sqrt(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)
356  else
357  tau_b(j,i) = 0.0_dp
358  end if
359 
360 #endif /* Normal vs. OpenAD */
361 
362 end do
363 end do
364 
365 !-------- Computation of d_help_b (defined on the grid points (i,j)) --------
366 
367 do i=0, imax
368 do j=0, jmax
369 
370  if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
371  ! grounded ice, or floating ice at the grounding line
372 
373 ! ------ Abbreviations
374 
375 #if (SLIDE_LAW==1)
376  cvxy1 = c_slide(j,i) &
377  * ( (tau_b(j,i)+eps_dp)**(p_weert(j,i)-1) &
378  /(p_b(j,i)+eps_dp)**q_weert(j,i) ) &
379  * p_b(j,i)
380  ctau1 = 1.0_dp/(c_slide(j,i)+eps_dp)**p_weert_inv(j,i) &
381  * (p_b(j,i)+eps_dp)**(q_weert(j,i)*p_weert_inv(j,i))
382  ! Basal sliding at pressure melting
383 #elif (SLIDE_LAW==2)
384  cvxy1 = c_slide(j,i) &
385  * ( (tau_b(j,i)+eps_dp)**(p_weert(j,i)-1) &
386  /(p_b_red(j,i)+eps_dp)**q_weert(j,i) ) &
387  * p_b(j,i)
388  ctau1 = 1.0_dp/(c_slide(j,i)+eps_dp)**p_weert_inv(j,i) &
389  * (p_b_red(j,i)+eps_dp)**(q_weert(j,i)*p_weert_inv(j,i))
390  ! Basal sliding at pressure melting
391 #endif
392 
393 ! ------ d_help_b, c_drag
394 
395  if (n_cts(j,i) == -1) then ! cold ice base
396 
397  if (sub_melt_flag(j,i)) then
398  temp_diff = max((temp_c_m(0,j,i)-temp_c(0,j,i)), 0.0_dp)
399  cvxy1a = exp(-gamma_slide_inv(j,i)*temp_diff) ! sub-melt sliding
400  ctau1a = 1.0_dp/cvxy1a**p_weert_inv(j,i)
401  else
402  cvxy1a = 0.0_dp ! no sub-melt sliding
403  ctau1a = 1.0_dp/eps**p_weert_inv(j,i) ! dummy value
404  end if
405 
406  d_help_b(j,i) = cvxy1*cvxy1a
407  c_drag(j,i) = ctau1*ctau1a
408 
409  else if (n_cts(j,i) == 0) then ! temperate ice base
410 
411  d_help_b(j,i) = cvxy1 ! basal sliding
412  c_drag(j,i) = ctau1 ! (pressure-melting conditions)
413 
414  else ! n_cts(j,i) == 1, temperate ice layer
415 
416  d_help_b(j,i) = cvxy1 ! basal sliding
417  c_drag(j,i) = ctau1 ! (pressure-melting conditions)
418 
419  end if
420 
421 ! ---- Contribution of the basal water layer
422 
423 #if (BASAL_HYDROLOGY==1)
424 
425  cvxy1b = (1.0_dp-c_hw_slide) &
426  + c_hw_slide*(1.0_dp-exp(-max(h_w(j,i),0.0_dp)*hw0_slide_inv))
427 
428  ctau1b = 1.0_dp/cvxy1b**p_weert_inv(j,i)
429 
430  d_help_b(j,i) = d_help_b(j,i) *cvxy1b
431  c_drag(j,i) = c_drag(j,i) *ctau1b
432 
433 #endif
434 
435  else ! maske(j,i) == 1, 2 or 3 away from the grounding line
436 
437  d_help_b(j,i) = 0.0_dp
438  c_drag(j,i) = 0.0_dp
439 
440  end if
441 
442 end do
443 end do
444 
445 !-------- Computation of vx_b (defined at (i+1/2,j)) --------
446 
447 do i=0, imax-1
448 do j=1, jmax-1
449  vx_b(j,i) = -0.5_dp*(d_help_b(j,i)+d_help_b(j,i+1))*dzs_dxi(j,i)
450 end do
451 end do
452 
453 !-------- Computation of vy_b (defined at (i,j+1/2)) --------
454 
455 do i=1, imax-1
456 do j=0, jmax-1
457  vy_b(j,i) = -0.5_dp*(d_help_b(j,i)+d_help_b(j+1,i))*dzs_deta(j,i)
458 end do
459 end do
460 
461 !-------- Computation of vx_b_g and vy_b_g (defined at (i,j)) --------
462 
463 do i=0, imax
464 do j=0, jmax
465  vx_b_g(j,i) = -d_help_b(j,i)*dzs_dxi_g(j,i)
466  vy_b_g(j,i) = -d_help_b(j,i)*dzs_deta_g(j,i)
467 end do
468 end do
469 
470 !-------- Limitation of computed vx_b, vy_b, vx_b_g, vy_b_g to the interval
471 ! [-VH_MAX, VH_MAX] --------
472 
473 vh_max = vh_max/year_sec
474 
475 do i=0, imax
476 do j=0, jmax
477  vx_b(j,i) = max(vx_b(j,i), -vh_max)
478  vx_b(j,i) = min(vx_b(j,i), vh_max)
479  vy_b(j,i) = max(vy_b(j,i), -vh_max)
480  vy_b(j,i) = min(vy_b(j,i), vh_max)
481  vx_b_g(j,i) = max(vx_b_g(j,i), -vh_max)
482  vx_b_g(j,i) = min(vx_b_g(j,i), vh_max)
483  vy_b_g(j,i) = max(vy_b_g(j,i), -vh_max)
484  vy_b_g(j,i) = min(vy_b_g(j,i), vh_max)
485 end do
486 end do
487 
488 end subroutine calc_vxy_b_sia
489 
490 !-------------------------------------------------------------------------------
491 !> Computation of the shear stresses txz, tyz, the effective shear stress
492 !! sigma, the depth-averaged fluidity flui_ave_sia, the horizontal
493 !! velocity vx, vy and the horizontal volume flux qx, qy in the shallow ice
494 !! approximation.
495 !<------------------------------------------------------------------------------
496 subroutine calc_vxy_sia(dzeta_c, dzeta_t)
499 
500 implicit none
501 
502 real(dp), intent(in) :: dzeta_c, dzeta_t
503 
504 integer(i4b) :: i, j, kc, kt
505 real(dp) :: avxy3(0:kcmax), aqxy1(0:kcmax)
506 real(dp) :: ctxyz1(0:kcmax,0:jmax,0:imax), &
507  ctxyz2(0:ktmax,0:jmax,0:imax)
508 real(dp) :: flui_t(0:ktmax), flui_c(0:kcmax)
509 real(dp) :: cflui0(0:ktmax), cflui1(0:kcmax)
510 real(dp) :: cvxy2(0:ktmax), cvxy3(0:kcmax)
511 real(dp) :: cqxy0(0:ktmax), cqxy1(0:kcmax)
512 real(dp) :: vh_max
513 real(dp) :: ratio_sl_threshold
514 
515 !-------- Term abbreviations --------
516 
517 do kc=0, kcmax
518  if (flag_aa_nonzero) then
519  avxy3(kc) = aa*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
520  aqxy1(kc) = aa/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
521  else
522  avxy3(kc) = dzeta_c
523  aqxy1(kc) = dzeta_c
524  end if
525 end do
526 
527 !-------- Computation of stresses --------
528 
529 ! ------ Term abbreviations
530 
531 do i=0, imax
532 do j=0, jmax
533 
534  if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
535  ! grounded ice, or floating ice at the grounding line
536 
537  do kc=0, kcmax
538  ctxyz1(kc,j,i) = rho*g*h_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
539  end do
540 
541  if (n_cts(j,i) == 1) then ! temperate layer
542 
543  do kt=0, ktmax
544  ctxyz2(kt,j,i) = rho*g*h_t(j,i)*(1.0_dp-zeta_t(kt))
545  end do
546 
547  else ! cold base (-1), temperate base (0)
548 
549  do kt=0, ktmax
550  ctxyz2(kt,j,i) = 0.0_dp
551  end do
552 
553  end if
554 
555  else ! maske(j,i) == 1, 2 or 3 away from the grounding line
556 
557  do kc=0, kcmax
558  ctxyz1(kc,j,i) = 0.0_dp
559  end do
560 
561  do kt=0, ktmax
562  ctxyz2(kt,j,i) = 0.0_dp
563  end do
564 
565  end if
566 
567 end do
568 end do
569 
570 ! ------ Shear stress txz (defined at (i+1/2,j,kc/t))
571 
572 do i=0, imax-1
573 do j=0, jmax
574 
575  do kc=0, kcmax
576  txz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j,i+1)) &
577  *dzs_dxi(j,i)
578  end do
579 
580  do kt=0, ktmax
581  txz_t(kt,j,i) = txz_c(0,j,i) &
582  -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j,i+1)) &
583  *dzs_dxi(j,i)
584  end do
585 
586 end do
587 end do
588 
589 ! ------ Shear stress tyz (defined at (i,j+1/2,kc/t))
590 
591 do i=0, imax
592 do j=0, jmax-1
593 
594  do kc=0, kcmax
595  tyz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j+1,i)) &
596  *dzs_deta(j,i)
597  end do
598 
599  do kt=0, ktmax
600  tyz_t(kt,j,i) = tyz_c(0,j,i) &
601  -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j+1,i)) &
602  *dzs_deta(j,i)
603  end do
604 
605 end do
606 end do
607 
608 ! ------ Effective shear stress sigma (defined at (i,j,kc/t))
609 
610 do i=0, imax
611 do j=0, jmax
612 
613  do kc=0, kcmax
614 
615 #ifndef ALLOW_OPENAD /* Normal */
616 
617  sigma_c(kc,j,i) = ctxyz1(kc,j,i) &
618  *sqrt(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)
619 
620 #else /* OpenAD: guarding against non-differentiable sqrt(0) */
621 
622  if ( (dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2) > 0 ) then
623  sigma_c(kc,j,i) = ctxyz1(kc,j,i) &
624  *sqrt(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)
625  else
626  sigma_c(kc,j,i) = 0.0_dp
627  end if
628 
629 #endif /* Normal vs. OpenAD */
630 
631  end do
632 
633  do kt=0, ktmax
634 
635 #ifndef ALLOW_OPENAD /* Normal */
636 
637  sigma_t(kt,j,i) = sigma_c(0,j,i) &
638  + ctxyz2(kt,j,i) &
639  *sqrt(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)
640 
641 #else /* OpenAD: guarding against non-differentiable sqrt(0) */
642 
643  if ( (dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2) > 0 ) then
644  sigma_t(kt,j,i) = sigma_c(0,j,i) &
645  + ctxyz2(kt,j,i) &
646  *sqrt(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)
647  else
648  sigma_t(kt,j,i) = sigma_c(0,j,i)
649  end if
650 
651 #endif /* Normal vs. OpenAD */
652 
653  end do
654 
655 end do
656 end do
657 
658 !-------- Computation of the depth-averaged fluidity
659 ! (defined on the grid points (i,j)) --------
660 
661 do i=0, imax
662 do j=0, jmax
663 
664  if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
665  ! grounded ice, or floating ice at the grounding line
666 
667 ! ------ Fluidity, abbreviations
668 
669  do kt=0, ktmax
670  flui_t(kt) = 2.0_dp &
671  *enh_t(kt,j,i) &
672  *ratefac_t(omega_t(kt,j,i)) &
673  *creep(sigma_t(kt,j,i))
674  cflui0(kt) = h_t(j,i)*flui_t(kt)*dzeta_t
675  end do
676 
677  do kc=0, kcmax
678  flui_c(kc) = 2.0_dp &
679  *enh_c(kc,j,i) &
680 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
681  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
682 #elif (CALCMOD==2 || CALCMOD==3)
683  *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), &
684  temp_c_m(kc,j,i)) &
685 #endif
686  *creep(sigma_c(kc,j,i))
687  cflui1(kc) = aqxy1(kc)*h_c(j,i)*flui_c(kc)
688  end do
689 
690 ! ------ Depth average
691 
692  flui_ave_sia(j,i) = 0.0_dp
693 
694  if (n_cts(j,i) == 1) then
695 
696  do kt=0, ktmax-1
697  flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui0(kt+1)+cflui0(kt))
698  end do
699 
700  end if
701 
702  do kc=0, kcmax-1
703  flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui1(kc+1)+cflui1(kc))
704  end do
705 
706  flui_ave_sia(j,i) = flui_ave_sia(j,i)/(h_c(j,i)+h_t(j,i))
707 
708  else ! maske(j,i) == 1, 2 or 3 away from the grounding line
709 
710  flui_ave_sia(j,i) = 0.0_dp
711 
712  end if
713 
714 end do
715 end do
716 
717 !-------- Computation of d_help_c/t
718 ! (defined on the grid points (i,j,kc/t)) --------
719 
720 do i=0, imax
721 do j=0, jmax
722 
723  if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
724  ! grounded ice, or floating ice at the grounding line
725 
726 ! ------ Abbreviations
727 
728  do kt=0, ktmax
729  cvxy2(kt) = 2.0_dp*h_t(j,i) &
730  *enh_t(kt,j,i) &
731  *ratefac_t(omega_t(kt,j,i)) &
732  *creep(sigma_t(kt,j,i)) &
733  *(ctxyz1(0,j,i)+ctxyz2(kt,j,i)) &
734  *dzeta_t
735  end do
736 
737  do kc=0, kcmax
738  cvxy3(kc) = 2.0_dp*avxy3(kc)*h_c(j,i) &
739  *enh_c(kc,j,i) &
740 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
741  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
742 #elif (CALCMOD==2 || CALCMOD==3)
743  *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), &
744  temp_c_m(kc,j,i)) &
745 #endif
746  *creep(sigma_c(kc,j,i)) &
747  *ctxyz1(kc,j,i)
748  end do
749 
750 ! ------ d_help_c, d_help_t
751 
752  if (n_cts(j,i) == -1) then ! cold ice base
753 
754  do kt=0, ktmax
755  d_help_t(kt,j,i) = d_help_b(j,i)
756  end do
757 
758  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
759 
760  do kc=0, kcmax-1
761  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
762  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
763  end do
764 
765  else if (n_cts(j,i) == 0) then ! temperate ice base
766 
767  do kt=0, ktmax
768  d_help_t(kt,j,i) = d_help_b(j,i)
769  end do
770 
771  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
772 
773  do kc=0, kcmax-1
774  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
775  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
776  end do
777 
778  else ! n_cts(j,i) == 1, temperate ice layer
779 
780  d_help_t(0,j,i) = d_help_b(j,i)
781 
782  do kt=0, ktmax-1
783  d_help_t(kt+1,j,i) = d_help_t(kt,j,i) &
784  +0.5_dp*(cvxy2(kt+1)+cvxy2(kt))
785  end do
786 
787  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
788 
789  do kc=0, kcmax-1
790  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
791  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
792  end do
793 
794  end if
795 
796  else ! maske(j,i) == 1, 2 or 3 away from the grounding line
797 
798  do kt=0, ktmax
799  d_help_t(kt,j,i) = 0.0_dp
800  end do
801 
802  do kc=0, kcmax
803  d_help_c(kc,j,i) = 0.0_dp
804  end do
805 
806  end if
807 
808 end do
809 end do
810 
811 !-------- Computation of vx_c/t (defined at (i+1/2,j,kc/t)) --------
812 
813 do i=0, imax-1
814 do j=1, jmax-1
815 
816  do kt=0, ktmax
817  vx_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j,i+1)) &
818  *dzs_dxi(j,i)
819  end do
820 
821  do kc=0, kcmax
822  vx_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j,i+1)) &
823  *dzs_dxi(j,i)
824  end do
825 
826 end do
827 end do
828 
829 !-------- Computation of vy_c/t (defined at (i,j+1/2,kc/t)) --------
830 
831 do i=1, imax-1
832 do j=0, jmax-1
833 
834  do kt=0, ktmax
835  vy_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j+1,i)) &
836  *dzs_deta(j,i)
837  end do
838 
839  do kc=0, kcmax
840  vy_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j+1,i)) &
841  *dzs_deta(j,i)
842  end do
843 
844 end do
845 end do
846 
847 !-------- Computation of the surface velocities vx_s_g and vy_s_g
848 ! (defined at (i,j)) --------
849 
850 do i=0, imax
851 do j=0, jmax
852  vx_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_dxi_g(j,i)
853  vy_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_deta_g(j,i)
854 end do
855 end do
856 
857 !-------- Limitation of computed vx_c/t, vy_c/t, vx_s_g, vy_s_g
858 ! to the interval [-VH_MAX, VH_MAX] --------
859 
860 vh_max = vh_max/year_sec
861 
862 do i=0, imax
863 do j=0, jmax
864  vx_s_g(j,i) = max(vx_s_g(j,i), -vh_max)
865  vx_s_g(j,i) = min(vx_s_g(j,i), vh_max)
866  vy_s_g(j,i) = max(vy_s_g(j,i), -vh_max)
867  vy_s_g(j,i) = min(vy_s_g(j,i), vh_max)
868  do kt=0, ktmax
869  vx_t(kt,j,i) = max(vx_t(kt,j,i), -vh_max)
870  vx_t(kt,j,i) = min(vx_t(kt,j,i), vh_max)
871  vy_t(kt,j,i) = max(vy_t(kt,j,i), -vh_max)
872  vy_t(kt,j,i) = min(vy_t(kt,j,i), vh_max)
873  end do
874  do kc=0, kcmax
875  vx_c(kc,j,i) = max(vx_c(kc,j,i), -vh_max)
876  vx_c(kc,j,i) = min(vx_c(kc,j,i), vh_max)
877  vy_c(kc,j,i) = max(vy_c(kc,j,i), -vh_max)
878  vy_c(kc,j,i) = min(vy_c(kc,j,i), vh_max)
879  end do
880 end do
881 end do
882 
883 !-------- Computation of h_diff
884 ! (defined on the grid points (i,j)) --------
885 
886 do i=0, imax
887 do j=0, jmax
888 
889  if ((maske(j,i) == 0).or.flag_grounding_line_2(j,i)) then
890  ! grounded ice, or floating ice at the grounding line
891 
892 ! ------ Abbreviations
893 
894  do kt=0, ktmax
895  cqxy0(kt) = h_t(j,i)*d_help_t(kt,j,i)*dzeta_t
896  end do
897 
898  do kc=0, kcmax
899  cqxy1(kc) = aqxy1(kc)*h_c(j,i)*d_help_c(kc,j,i)
900  end do
901 
902 ! ------ h_diff
903 
904  h_diff(j,i) = 0.0_dp
905 
906  if (n_cts(j,i) == 1) then
907 
908  do kt=0, ktmax-1
909  h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy0(kt+1)+cqxy0(kt))
910  end do
911 
912  end if
913 
914  do kc=0, kcmax-1
915  h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy1(kc+1)+cqxy1(kc))
916  end do
917 
918 ! ------ Limitation of h_diff
919 
920  if (h_diff(j,i) < hd_min) h_diff(j,i) = 0.0_dp
921  if (h_diff(j,i) > hd_max) h_diff(j,i) = hd_max
922 
923  else ! maske(j,i) == 1, 2 or 3 away from the grounding line
924 
925  h_diff(j,i) = 0.0_dp
926 
927  end if
928 
929 end do
930 end do
931 
932 !-------- Computation of the horizontal volume flux
933 ! and the depth-averaged velocity --------
934 
935 do i=0, imax-1
936 do j=0, jmax
937 
938  qx(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j,i+1))*dzs_dxi(j,i)
939 
940  if ( (maske(j,i)==0).or.(maske(j,i+1)==0) ) then ! at least one neighbour
941  ! point is grounded ice
942 
943  vx_m(j,i) = qx(j,i) &
944  / ( 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) )
945 
946  vx_m(j,i) = max(vx_m(j,i), -vh_max) ! Limitation
947  vx_m(j,i) = min(vx_m(j,i), vh_max) ! to the interval [-VH_MAX, VH_MAX]
948 
949  ratio_sl_x(j,i) = abs(vx_t(0,j,i)) / max(abs(vx_c(kcmax,j,i)), eps_dp)
950 
951  else
952 
953  vx_m(j,i) = 0.0_dp
954  ratio_sl_x(j,i) = 0.0_dp
955 
956  end if
957 
958 end do
959 end do
960 
961 do i=0, imax
962 do j=0, jmax-1
963 
964  qy(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j+1,i))*dzs_deta(j,i)
965 
966  if ( (maske(j,i)==0).or.(maske(j+1,i)==0) ) then ! at least one neighbour
967  ! point is grounded ice
968 
969  vy_m(j,i) = qy(j,i) &
970  / ( 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) )
971 
972  vy_m(j,i) = max(vy_m(j,i), -vh_max) ! Limitation
973  vy_m(j,i) = min(vy_m(j,i), vh_max) ! to the interval [-VH_MAX, VH_MAX]
974 
975  ratio_sl_y(j,i) = abs(vy_t(0,j,i)) / max(abs(vy_c(kcmax,j,i)), eps_dp)
976 
977  else
978 
979  vy_m(j,i) = 0.0_dp
980  ratio_sl_y(j,i) = 0.0_dp
981 
982  end if
983 
984 end do
985 end do
986 
987 !-------- Detection of shelfy stream points --------
988 
989 flag_shelfy_stream_x = .false.
990 flag_shelfy_stream_y = .false.
991 flag_shelfy_stream = .false.
992 
993 #if (DYNAMICS==0 || DYNAMICS==1)
994 
995 ratio_sl_threshold = 1.11e+11_dp ! dummy value
996 
997 #elif (DYNAMICS==2)
998 
999 #if ( defined(RATIO_SL_THRESH) )
1000 ratio_sl_threshold = ratio_sl_thresh
1001 #else
1002 ratio_sl_threshold = 0.5_dp ! default value
1003 #endif
1004 
1005 do i=0, imax-1
1006 do j=0, jmax
1007  if (ratio_sl_x(j,i) > ratio_sl_threshold) flag_shelfy_stream_x(j,i) = .true.
1008 end do
1009 end do
1010 
1011 do i=0, imax
1012 do j=0, jmax-1
1013  if (ratio_sl_y(j,i) > ratio_sl_threshold) flag_shelfy_stream_y(j,i) = .true.
1014 end do
1015 end do
1016 
1017 do i=1, imax-1
1018 do j=1, jmax-1
1019 
1020  if ( (maske(j,i) == 0) & ! grounded ice
1021  .and. &
1022  ( flag_shelfy_stream_x(j,i-1) & ! at least
1023  .or.flag_shelfy_stream_x(j,i) & ! one neighbour
1024  .or.flag_shelfy_stream_y(j-1,i) & ! on the staggered grid
1025  .or.flag_shelfy_stream_y(j,i) ) & ! is a shelfy stream point
1026  ) then
1027 
1028  flag_shelfy_stream(j,i) = .true.
1029 
1030  end if
1031 
1032 end do
1033 end do
1034 
1035 #else
1036 
1037 errormsg = ' >>> calc_vxy_sia: DYNAMICS must be 0, 1 or 2!'
1038 call error(errormsg)
1039 
1040 #endif
1041 
1042 !-------- Initialisation of the variable q_gl_g
1043 ! (volume flux across the grounding line, to be
1044 ! computed in the routine calc_vxy_ssa
1045 ! if ice shelves are present)
1046 
1047 q_gl_g = 0.0_dp
1048 
1049 end subroutine calc_vxy_sia
1050 
1051 !-------------------------------------------------------------------------------
1052 !> Computation of the horizontal velocity vx, vy, the horizontal volume flux
1053 !> qx, qy etc. for static ice.
1054 !<------------------------------------------------------------------------------
1055 subroutine calc_vxy_static()
1057 implicit none
1058 
1059 c_slide = 0.0_dp
1060 p_weert = 0.0_dp
1061 q_weert = 0.0_dp
1062 p_b_w = 0.0_dp
1063 
1064 d_help_b = 0.0_dp
1065 c_drag = 0.0_dp
1066 
1067 vx_b = 0.0_dp
1068 vy_b = 0.0_dp
1069 vx_b_g = 0.0_dp
1070 vy_b_g = 0.0_dp
1071 
1072 txz_c = 0.0_dp
1073 txz_t = 0.0_dp
1074 
1075 tyz_c = 0.0_dp
1076 tyz_t = 0.0_dp
1077 
1078 sigma_c = 0.0_dp
1079 sigma_t = 0.0_dp
1080 
1081 flui_ave_sia = 0.0_dp
1082 de_ssa = 0.0_dp
1083 vis_int_g = 0.0_dp
1084 
1085 d_help_t = 0.0_dp
1086 d_help_c = 0.0_dp
1087 
1088 vx_c = 0.0_dp
1089 vy_c = 0.0_dp
1090 
1091 vx_t = 0.0_dp
1092 vy_t = 0.0_dp
1093 
1094 vx_s_g = 0.0_dp
1095 vy_s_g = 0.0_dp
1096 
1097 h_diff = 0.0_dp
1098 
1099 qx = 0.0_dp
1100 qy = 0.0_dp
1101 
1102 vx_m = 0.0_dp
1103 vy_m = 0.0_dp
1104 
1105 ratio_sl_x = 0.0_dp
1106 ratio_sl_y = 0.0_dp
1107 
1108 flag_shelfy_stream_x = .false.
1109 flag_shelfy_stream_y = .false.
1110 flag_shelfy_stream = .false.
1111 
1112 q_gl_g = 0.0_dp
1113 
1114 end subroutine calc_vxy_static
1115 
1116 !-------------------------------------------------------------------------------
1117 !> Computation of the horizontal velocity vx, vy, the horizontal volume flux
1118 !! qx, qy and the flux across the grounding line q_gl_g in the shallow shelf
1119 !! approximation (SSA) or the shelfy stream approximation (SStA).
1120 !<------------------------------------------------------------------------------
1121 subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1123 implicit none
1124 
1125 #ifndef ALLOW_OPENAD /* Normal */
1126 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
1127 #else /* OpenAD */
1128 integer(i4b), intent(inout) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
1129 #endif /* Normal vs. OpenAD */
1130 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
1131 real(dp), intent(in) :: z_sl, dxi, deta, dzeta_c, dzeta_t
1132 
1133 integer(i4b) :: i, j, kc, kt, m
1134 integer(i4b) :: iter_ssa
1135 real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_sia, vy_m_sia
1136 real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_prev, vy_m_prev
1137 real(dp) :: rel_ssa
1138 real(dp) :: dxi_inv, deta_inv
1139 real(dp) :: vh_max
1140 real(dp) :: ratio_sl_threshold, ratio_help
1141 real(dp), dimension(0:JMAX,0:IMAX) :: weigh_ssta_sia_x, weigh_ssta_sia_y
1142 real(dp) :: qx_gl_g, qy_gl_g
1143 
1144 #if (MARGIN==3 || DYNAMICS==2)
1145 
1146 !-------- Parameters for the relaxation scheme --------
1147 
1148 #if (defined(N_ITER_SSA))
1149  iter_ssa = max(n_iter_ssa, 1) ! number of iterations
1150 #else
1151  iter_ssa = 2 ! default value
1152 #endif
1153 
1154 #if (defined(RELAX_FACT_SSA))
1155  rel_ssa = relax_fact_ssa ! relaxation factor
1156 #else
1157  rel_ssa = 0.7_dp ! default value
1158 #endif
1159 
1160 write(6,'(10x,a,i0,a,f6.3)') 'calc_vxy_ssa: iter_ssa = ', iter_ssa, &
1161  ', rel_ssa =' , rel_ssa
1162 
1163 !-------- Detection of the grounding line and the calving front --------
1164 
1165 flag_grounding_line_1 = .false.
1166 flag_grounding_line_2 = .false.
1167 
1168 flag_calving_front_1 = .false.
1169 flag_calving_front_2 = .false.
1170 
1171 do i=1, imax-1
1172 do j=1, jmax-1
1173 
1174  if ( (maske(j,i)==0_i2b) & ! grounded ice
1175  .and. &
1176  ( (maske(j,i+1)==3_i2b) & ! with
1177  .or.(maske(j,i-1)==3_i2b) & ! one
1178  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1179  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1180  ) &
1181  flag_grounding_line_1(j,i) = .true.
1182 
1183  if ( (maske(j,i)==3_i2b) & ! floating ice
1184  .and. &
1185  ( (maske(j,i+1)==0_i2b) & ! with
1186  .or.(maske(j,i-1)==0_i2b) & ! one
1187  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1188  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1189  ) &
1190  flag_grounding_line_2(j,i) = .true.
1191 
1192  if ( (maske(j,i)==3_i2b) & ! floating ice
1193  .and. &
1194  ( (maske(j,i+1)==2_i2b) & ! with
1195  .or.(maske(j,i-1)==2_i2b) & ! one
1196  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1197  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1198  ) &
1199  flag_calving_front_1(j,i) = .true.
1200 
1201  if ( (maske(j,i)==2_i2b) & ! ocean
1202  .and. &
1203  ( (maske(j,i+1)==3_i2b) & ! with
1204  .or.(maske(j,i-1)==3_i2b) & ! one
1205  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1206  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1207  ) &
1208  flag_calving_front_2(j,i) = .true.
1209 
1210 end do
1211 end do
1212 
1213 do i=1, imax-1
1214 
1215  j=0
1216  if ( (maske(j,i)==2_i2b) & ! ocean
1217  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1218  ) & ! floating ice point
1219  flag_calving_front_2(j,i) = .true.
1220 
1221  j=jmax
1222  if ( (maske(j,i)==2_i2b) & ! ocean
1223  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1224  ) & ! floating ice point
1225  flag_calving_front_2(j,i) = .true.
1226 
1227 end do
1228 
1229 do j=1, jmax-1
1230 
1231  i=0
1232  if ( (maske(j,i)==2_i2b) & ! ocean
1233  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1234  ) & ! floating ice point
1235  flag_calving_front_2(j,i) = .true.
1236 
1237  i=imax
1238  if ( (maske(j,i)==2_i2b) & ! ocean
1239  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1240  ) & ! floating ice point
1241  flag_calving_front_2(j,i) = .true.
1242 
1243 end do
1244 
1245 !-------- Save mean (depth-averaged) horizontal velocities from SIA --------
1246 
1247 vx_m_sia = vx_m
1248 vy_m_sia = vy_m
1249 
1250 !-------- First iteration --------
1251 
1252 m=1
1253 
1254 ! ------ Depth-integrated viscosity vis_int_g
1255 
1256 vis_int_g = 1.0e+15_dp*(h_c+h_t)
1257  ! viscosity of 10^15 Pa s assumed
1258 
1259 ! ------ Horizontal velocity vx_m, vy_m
1260 
1261 call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn, m, iter_ssa)
1262 
1263 !-------- Further iterations --------
1264 
1265 do m=2, iter_ssa
1266 
1267 ! ------ Save velocities from previous iteration
1268 
1269  vx_m_prev = vx_m
1270  vy_m_prev = vy_m
1271 
1272 ! ------ Depth-integrated viscosity vis_int_g
1273 
1274  call calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)
1275 
1276 ! ------ Horizontal velocity vx_m, vy_m
1277 
1278  call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn, m, iter_ssa)
1279 
1280 ! ------ Relaxation scheme
1281 
1282  vx_m = rel_ssa*vx_m + (1.0_dp-rel_ssa)*vx_m_prev
1283  vy_m = rel_ssa*vy_m + (1.0_dp-rel_ssa)*vy_m_prev
1284 
1285 end do
1286 
1287 ! ------ Depth-integrated viscosity vis_int_g
1288 
1289 call calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)
1290 
1291 !-------- Limitation of computed vx_m and vy_m to the interval
1292 ! [-VH_MAX, VH_MAX] --------
1293 
1294 vh_max = vh_max/year_sec
1295 
1296 vx_m = max(vx_m, -vh_max)
1297 vx_m = min(vx_m, vh_max)
1298 vy_m = max(vy_m, -vh_max)
1299 vy_m = min(vy_m, vh_max)
1300 
1301 !-------- 3-D velocities, basal velocities and volume flux --------
1302 
1303 #if (DYNAMICS==0 || DYNAMICS==1)
1304 
1305 ratio_sl_threshold = 1.11e+11_dp ! dummy value
1306 ratio_help = 0.0_dp
1307 
1308 #elif (DYNAMICS==2)
1309 
1310 #if ( defined(RATIO_SL_THRESH) )
1311 ratio_sl_threshold = ratio_sl_thresh
1312 #else
1313 ratio_sl_threshold = 0.5_dp ! default value
1314 #endif
1315 
1316 ratio_help = 1.0_dp/(1.0_dp-ratio_sl_threshold)
1317 
1318 #else
1319 
1320 errormsg = ' >>> calc_vxy_ssa: DYNAMICS must be 0, 1 or 2!'
1321 call error(errormsg)
1322 
1323 #endif
1324 
1325 weigh_ssta_sia_x = 0.0_dp
1326 weigh_ssta_sia_y = 0.0_dp
1327 
1328 ! ------ x-component
1329 
1330 do i=0, imax-1
1331 do j=0, jmax
1332 
1333  if (flag_shelfy_stream_x(j,i)) then
1334  ! shelfy stream
1335 
1336  weigh_ssta_sia_x(j,i) = (ratio_sl_x(j,i)-ratio_sl_threshold)*ratio_help
1337 
1338  weigh_ssta_sia_x(j,i) = max(min(weigh_ssta_sia_x(j,i), 1.0_dp), 0.0_dp)
1339  ! constrain to interval [0,1]
1340 
1341 #if (SSTA_SIA_WEIGH_FCT==0)
1342 
1343  ! stick to the linear function set above
1344 
1345 #elif (SSTA_SIA_WEIGH_FCT==1)
1346 
1347  weigh_ssta_sia_x(j,i) = weigh_ssta_sia_x(j,i)*weigh_ssta_sia_x(j,i) &
1348  *(3.0_dp-2.0_dp*weigh_ssta_sia_x(j,i))
1349  ! make transition smooth (cubic function)
1350 
1351 #elif (SSTA_SIA_WEIGH_FCT==2)
1352 
1353  weigh_ssta_sia_x(j,i) = weigh_ssta_sia_x(j,i)*weigh_ssta_sia_x(j,i) &
1354  *weigh_ssta_sia_x(j,i) &
1355  *(10.0_dp + weigh_ssta_sia_x(j,i) &
1356  *(-15.0_dp+6.0_dp*weigh_ssta_sia_x(j,i)))
1357  ! make transition even smoother (quintic function)
1358 
1359 #else
1360  errormsg = ' >>> calc_vxy_ssa: SSTA_SIA_WEIGH_FCT must be 0, 1 or 2!'
1361  call error(errormsg)
1362 #endif
1363 
1364  do kt=0, ktmax
1365  vx_t(kt,j,i) = weigh_ssta_sia_x(j,i)*vx_m(j,i) &
1366  + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_t(kt,j,i)
1367  end do
1368 
1369  do kc=0, kcmax
1370  vx_c(kc,j,i) = weigh_ssta_sia_x(j,i)*vx_m(j,i) &
1371  + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_c(kc,j,i)
1372  end do
1373 
1374  vx_b(j,i) = vx_t(0,j,i)
1375 
1376  vx_m(j,i) = weigh_ssta_sia_x(j,i)*vx_m(j,i) &
1377  + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_m_sia(j,i)
1378 
1379  qx(j,i) = vx_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))
1380 
1381  else if ( (maske(j,i)==3_i2b).or.(maske(j,i+1)==3_i2b) ) then
1382  ! at least one neighbour point is floating ice
1383  do kt=0, ktmax
1384  vx_t(kt,j,i) = vx_m(j,i)
1385  end do
1386 
1387  do kc=0, kcmax
1388  vx_c(kc,j,i) = vx_m(j,i)
1389  end do
1390 
1391  vx_b(j,i) = vx_m(j,i)
1392 
1393  qx(j,i) = vx_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))
1394 
1395  end if
1396 
1397 end do
1398 end do
1399 
1400 ! ------ y-component
1401 
1402 do i=0, imax
1403 do j=0, jmax-1
1404 
1405  if (flag_shelfy_stream_y(j,i)) then
1406  ! shelfy stream
1407 
1408  weigh_ssta_sia_y(j,i) = (ratio_sl_y(j,i)-ratio_sl_threshold)*ratio_help
1409 
1410  weigh_ssta_sia_y(j,i) = max(min(weigh_ssta_sia_y(j,i), 1.0_dp), 0.0_dp)
1411  ! constrain to interval [0,1]
1412 
1413 #if (SSTA_SIA_WEIGH_FCT==0)
1414 
1415  ! stick to the linear function set above
1416 
1417 #elif (SSTA_SIA_WEIGH_FCT==1)
1418 
1419  weigh_ssta_sia_y(j,i) = weigh_ssta_sia_y(j,i)*weigh_ssta_sia_y(j,i) &
1420  *(3.0_dp-2.0_dp*weigh_ssta_sia_y(j,i))
1421  ! make transition smooth (cubic function)
1422 
1423 #elif (SSTA_SIA_WEIGH_FCT==2)
1424 
1425  weigh_ssta_sia_y(j,i) = weigh_ssta_sia_y(j,i)*weigh_ssta_sia_y(j,i) &
1426  *weigh_ssta_sia_y(j,i) &
1427  *(10.0_dp + weigh_ssta_sia_y(j,i) &
1428  *(-15.0_dp+6.0_dp*weigh_ssta_sia_y(j,i)))
1429  ! make transition even smoother (quintic function)
1430 
1431 #else
1432  errormsg = ' >>> calc_vxy_ssa: SSTA_SIA_WEIGH_FCT must be 0, 1 or 2!'
1433  call error(errormsg)
1434 #endif
1435 
1436  do kt=0, ktmax
1437  vy_t(kt,j,i) = weigh_ssta_sia_y(j,i)*vy_m(j,i) &
1438  + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_t(kt,j,i)
1439  end do
1440 
1441  do kc=0, kcmax
1442  vy_c(kc,j,i) = weigh_ssta_sia_y(j,i)*vy_m(j,i) &
1443  + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_c(kc,j,i)
1444  end do
1445 
1446  vy_b(j,i) = vy_t(0,j,i)
1447 
1448  vy_m(j,i) = weigh_ssta_sia_y(j,i)*vy_m(j,i) &
1449  + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_m_sia(j,i)
1450 
1451  qy(j,i) = vy_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i))
1452 
1453  else if ( (maske(j,i)==3_i2b).or.(maske(j+1,i)==3_i2b) ) then
1454  ! at least one neighbour point is floating ice
1455  do kt=0, ktmax
1456  vy_t(kt,j,i) = vy_m(j,i)
1457  end do
1458 
1459  do kc=0, kcmax
1460  vy_c(kc,j,i) = vy_m(j,i)
1461  end do
1462 
1463  vy_b(j,i) = vy_m(j,i)
1464 
1465  qy(j,i) = vy_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i))
1466 
1467  end if
1468 
1469 end do
1470 end do
1471 
1472 !-------- Surface and basal velocities vx_s_g vy_s_g, vx_b_g vy_b_g
1473 ! (defined at (i,j)) --------
1474 
1475 do i=1, imax-1
1476 do j=1, jmax-1
1477 
1478  if (flag_shelfy_stream(j,i)) then ! shelfy stream
1479 
1480  vx_s_g(j,i) = 0.5_dp*(vx_c(kcmax,j,i-1)+vx_c(kcmax,j,i))
1481  vx_b_g(j,i) = 0.5_dp*(vx_b( j,i-1)+vx_b( j,i))
1482 
1483  vy_s_g(j,i) = 0.5_dp*(vy_c(kcmax,j-1,i)+vy_c(kcmax,j,i))
1484  vy_b_g(j,i) = 0.5_dp*(vy_b( j-1,i)+vy_b( j,i))
1485 
1486  else if (maske(j,i)==3_i2b) then ! floating ice
1487 
1488  vx_s_g(j,i) = 0.5_dp*(vx_m(j,i-1)+vx_m(j,i))
1489  vx_b_g(j,i) = vx_s_g(j,i)
1490 
1491  vy_s_g(j,i) = 0.5_dp*(vy_m(j-1,i)+vy_m(j,i))
1492  vy_b_g(j,i) = vy_s_g(j,i)
1493 
1494  end if
1495 
1496 end do
1497 end do
1498 
1499 !-------- Computation of the flux across the grounding line q_gl_g
1500 
1501 do i=1, imax-1
1502 do j=1, jmax-1
1503 
1504  if ( flag_grounding_line_1(j,i) ) then ! grounding line
1505 
1506  qx_gl_g = 0.5_dp*(qx(j,i)+qx(j,i-1))
1507  qy_gl_g = 0.5_dp*(qy(j,i)+qy(j-1,i))
1508 
1509 #ifndef ALLOW_OPENAD /* Normal */
1510 
1511  q_gl_g(j,i) = sqrt(qx_gl_g*qx_gl_g+qy_gl_g*qy_gl_g)
1512 
1513 #else /* OpenAD: guarding against non-differentiable sqrt(0) */
1514 
1515  if ( (qx_gl_g*qx_gl_g+qy_gl_g*qy_gl_g) > 0 ) then
1516  q_gl_g(j,i) = sqrt(qx_gl_g*qx_gl_g+qy_gl_g*qy_gl_g)
1517  else
1518  q_gl_g(j,i) = 0.0_dp
1519  end if
1520 
1521 #endif /* Normal vs. OpenAD */
1522 
1523  end if
1524 
1525 end do
1526 end do
1527 
1528 #else
1529 
1530 errormsg = ' >>> calc_vxy_ssa: Only to be called for MARGIN==3 or DYNAMICS==2!'
1531 call error(errormsg)
1532 
1533 #endif
1534 
1535 end subroutine calc_vxy_ssa
1536 
1537 !-------------------------------------------------------------------------------
1538 !> Solution of the system of linear equations for the horizontal velocities
1539 !! vx_m, vy_m in the shallow shelf approximation.
1540 !<------------------------------------------------------------------------------
1541 subroutine calc_vxy_ssa_matrix(z_sl,dxi, deta, ii, jj, nn, m, iter_ssa)
1542 
1543 #if (MARGIN==3 || DYNAMICS==2)
1544 #ifdef ALLOW_OPENAD /* Normal */
1545 use sico_maths_m
1546 #endif /* Normal */
1547 #endif
1548 
1549 implicit none
1550 
1551 
1552 #ifndef ALLOW_OPENAD /* Normal */
1553 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
1554 #else /* OpenAD */
1555 integer(i4b), intent(inout) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
1556 #endif /* Normal vs. OpenAD */
1557 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
1558 integer(i4b), intent(in) :: m, iter_ssa
1559 real(dp), intent(in) :: dxi, deta, z_sl
1560 
1561 integer(i4b) :: i, j, k, n
1562 integer(i4b) :: i1, j1
1563 real(dp) :: inv_dxi, inv_deta, inv_dxi_deta, inv_dxi2, inv_deta2
1564 real(dp) :: factor_rhs_1, factor_rhs_2, rhosw_rho_ratio
1565 real(dp) :: H_mid, zl_mid
1566 real(dp), dimension(0:JMAX,0:IMAX) :: vis_int_sgxy, beta_drag
1567 character(len=256) :: ch_solver_set_option
1568 
1569 #if (MARGIN==3 || DYNAMICS==2)
1570 
1571 #ifndef ALLOW_OPENAD /* Normal */
1572 
1573 lis_integer :: ierr
1574 lis_integer :: nc, nr
1575 lis_integer :: iter
1576 lis_matrix :: lgs_a
1577 lis_vector :: lgs_b, lgs_x
1578 lis_solver :: solver
1579 
1580 lis_integer, parameter :: nmax = 2*(imax+1)*(jmax+1)
1581 lis_integer, parameter :: n_sprs = 20*(imax+1)*(jmax+1)
1582 lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
1583 lis_scalar, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
1584 
1585 #else /* OpenAD */
1586 
1587 integer(i4b) :: ierr
1588 integer(i4b) :: nc, nr
1589 integer(i4b) :: iter
1590 
1591 integer(i4b), parameter :: nmax = 2*(imax+1)*(jmax+1)
1592 integer(i4b), parameter :: n_sprs = 20*(imax+1)*(jmax+1)
1593 integer(i4b), dimension(nmax+1) :: lgs_a_ptr
1594 real(dp), dimension(n_sprs) :: lgs_a_value
1595 real(dp), dimension(n_sprs) :: lgs_a_index
1596 integer(i4b), dimension(n_sprs) :: lgs_a_index_pass
1597 real(dp), dimension(nmax) :: lgs_b_value, lgs_x_value
1598 
1599 #endif /* Normal vs. OpenAD */
1600 
1601 !-------- Abbreviations --------
1602 
1603 inv_dxi = 1.0_dp/dxi
1604 inv_deta = 1.0_dp/deta
1605 inv_dxi_deta = 1.0_dp/(dxi*deta)
1606 inv_dxi2 = 1.0_dp/(dxi*dxi)
1607 inv_deta2 = 1.0_dp/(deta*deta)
1608 
1609 rhosw_rho_ratio = rho_sw/rho
1610 factor_rhs_1 = 0.5_dp*rho*g
1611 factor_rhs_2 = 0.5_dp*rho*g*(rho_sw-rho)/rho_sw
1612 
1613 !-------- Depth-integrated viscosity on the staggered grid
1614 ! [at (i+1/2,j+1/2)] --------
1615 
1616 vis_int_sgxy = 0.0_dp ! initialisation
1617 
1618 do i=0, imax-1
1619 do j=0, jmax-1
1620 
1621  k=0
1622 
1623  if ((maske(j,i)==0_i2b).or.(maske(j,i)==3_i2b)) then
1624  k = k+1 ! floating or grounded ice
1625  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j,i)
1626  end if
1627 
1628  if ((maske(j,i+1)==0_i2b).or.(maske(j,i+1)==3_i2b)) then
1629  k = k+1 ! floating or grounded ice
1630  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j,i+1)
1631  end if
1632 
1633  if ((maske(j+1,i)==0_i2b).or.(maske(j+1,i)==3_i2b)) then
1634  k = k+1 ! floating or grounded ice
1635  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j+1,i)
1636  end if
1637 
1638  if ((maske(j+1,i+1)==0_i2b).or.(maske(j+1,i+1)==3_i2b)) then
1639  k = k+1 ! floating or grounded ice
1640  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j+1,i+1)
1641  end if
1642 
1643  if (k>0) vis_int_sgxy(j,i) = vis_int_sgxy(j,i)/real(k,dp)
1644 
1645 end do
1646 end do
1647 
1648 !-------- Basal drag parameter (for shelfy stream) --------
1649 
1650 beta_drag = 0.0_dp ! initialisation
1651 
1652 do i=1, imax-1
1653 do j=1, jmax-1
1654 
1655  if (flag_shelfy_stream(j,i)) then
1656 
1657  if (m==1) then ! first iteration
1658 
1659  beta_drag(j,i) = c_drag(j,i) &
1660  / sqrt( ( vx_b_g(j,i)**2 &
1661  + vy_b_g(j,i)**2 ) &
1662  + eps_dp**2 ) &
1663  **(1.0_dp-p_weert_inv(j,i))
1664  ! computed with the SIA basal velocities
1665  else
1666 
1667  beta_drag(j,i) = c_drag(j,i) &
1668  / sqrt( ( (0.5_dp*(vx_m(j,i)+vx_m(j,i-1)))**2 &
1669  + (0.5_dp*(vy_m(j,i)+vy_m(j-1,i)))**2 ) &
1670  + eps_dp**2 ) &
1671  **(1.0_dp-p_weert_inv(j,i))
1672  ! computed with the SSA basal velocities
1673  ! from the previous iteration
1674  end if
1675 
1676  end if
1677 
1678 end do
1679 end do
1680 
1681 !-------- Assembly of the system of linear equations
1682 ! (matrix storage: compressed sparse row CSR) --------
1683 
1684 #ifndef ALLOW_OPENAD /* Normal */
1685 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
1686 allocate(lgs_b_value(nmax), lgs_x_value(nmax))
1687 #endif /* Normal */
1688 
1689 lgs_a_value = 0.0_dp
1690 #ifndef ALLOW_OPENAD /* Normal */
1691 lgs_a_index = 0
1692 #else /* OpenAD */
1693 lgs_a_index = 0.0_dp
1694 #endif /* Normal vs. OpenAD */
1695 lgs_a_ptr = 0
1696 
1697 lgs_b_value = 0.0_dp
1698 lgs_x_value = 0.0_dp
1699 
1700 lgs_a_ptr(1) = 1
1701 
1702 k = 0
1703 
1704 do n=1, nmax-1, 2
1705 
1706  i = ii((n+1)/2)
1707  j = jj((n+1)/2)
1708 
1709 ! ------ Equations for vx_m (at (i+1/2,j))
1710 
1711  nr = n ! row counter
1712 
1713  if ( (i /= imax).and.(j /= 0).and.(j /= jmax) ) then
1714  ! inner point on the staggered grid in x-direction
1715 
1716  h_mid = 0.50_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))
1717  zl_mid = 0.50_dp*(zl(j,i)+zl(j,i+1))
1718 
1719  if ( &
1720  ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==3_i2b) ) &
1721  .or. &
1722  ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j,i+1) &
1723  .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1724  .or. &
1725  ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j,i+1) &
1726  .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1727  ) then
1728  ! both neighbours are floating ice,
1729  ! or one neighbour is floating ice and the other is grounded ice
1730  ! (grounding line)
1731  ! and floating conditions are satisfied;
1732  ! discretisation of the x-component of the PDE
1733 
1734  flag_shelfy_stream_x(j,i)=.false.
1735  ! make sure not to treat as shelfy stream
1736 
1737  nc = 2*nn(j-1,i)-1 ! smallest nc (column counter), for vx_m(j-1,i)
1738  k = k+1
1739  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1740  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
1741  lgs_a_index(k) = nc
1742 
1743  nc = 2*nn(j-1,i) ! next nc (column counter), for vy_m(j-1,i)
1744  k = k+1
1745  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1746  lgs_a_value(k) = inv_dxi_deta &
1747  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j-1,i))
1748  lgs_a_index(k) = nc
1749 
1750  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
1751  k = k+1
1752  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1753  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i)
1754  lgs_a_index(k) = nc
1755 
1756  nc = 2*nn(j-1,i+1) ! next nc (column counter), for vy_m(j-1,i+1)
1757  k = k+1
1758  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1759  lgs_a_value(k) = -inv_dxi_deta &
1760  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
1761  lgs_a_index(k) = nc
1762 
1763  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
1764  if (nc /= nr) then ! (diagonal element)
1765  errormsg = ' >>> calc_vxy_ssa_matrix: ' &
1766  //'Check for diagonal element failed!'
1767  call error(errormsg)
1768  end if
1769  k = k+1
1770  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1771  lgs_a_value(k) = -4.0_dp*inv_dxi2 &
1772  *(vis_int_g(j,i+1)+vis_int_g(j,i)) &
1773  -inv_deta2 &
1774  *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i))
1775  lgs_a_index(k) = nc
1776 
1777  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
1778  k = k+1
1779  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1780  lgs_a_value(k) = -inv_dxi_deta &
1781  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
1782  lgs_a_index(k) = nc
1783 
1784  nc = 2*nn(j,i+1)-1 ! next nc (column counter), for vx_m(j,i+1)
1785  k = k+1
1786  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1787  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i+1)
1788  lgs_a_index(k) = nc
1789 
1790  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
1791  k = k+1
1792  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1793  lgs_a_value(k) = inv_dxi_deta &
1794  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j,i))
1795  lgs_a_index(k) = nc
1796 
1797  nc = 2*nn(j+1,i)-1 ! largest nc (column counter), for vx_m(j+1,i)
1798  k = k+1
1799  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1800  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
1801  lgs_a_index(k) = nc
1802 
1803  lgs_b_value(nr) = factor_rhs_1 &
1804  * (h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) &
1805  * dzs_dxi(j,i)
1806 
1807  lgs_x_value(nr) = vx_m(j,i)
1808 
1809  else if (flag_shelfy_stream_x(j,i)) then
1810  ! shelfy stream (as determined by routine calc_vxy_sia)
1811 
1812  nc = 2*nn(j-1,i)-1 ! smallest nc (column counter), for vx_m(j-1,i)
1813  k = k+1
1814  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1815  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
1816  lgs_a_index(k) = nc
1817 
1818  nc = 2*nn(j-1,i) ! next nc (column counter), for vy_m(j-1,i)
1819  k = k+1
1820  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1821  lgs_a_value(k) = inv_dxi_deta &
1822  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j-1,i))
1823  lgs_a_index(k) = nc
1824 
1825  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
1826  k = k+1
1827  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1828  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i)
1829  lgs_a_index(k) = nc
1830 
1831  nc = 2*nn(j-1,i+1) ! next nc (column counter), for vy_m(j-1,i+1)
1832  k = k+1
1833  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1834  lgs_a_value(k) = -inv_dxi_deta &
1835  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
1836  lgs_a_index(k) = nc
1837 
1838  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
1839  if (nc /= nr) then ! (diagonal element)
1840  errormsg = ' >>> calc_vxy_ssa_matrix: ' &
1841  //'Check for diagonal element failed!'
1842  call error(errormsg)
1843  end if
1844  k = k+1
1845  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1846  lgs_a_value(k) = -4.0_dp*inv_dxi2 &
1847  *(vis_int_g(j,i+1)+vis_int_g(j,i)) &
1848  -inv_deta2 &
1849  *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i)) &
1850  -0.5_dp*(beta_drag(j,i+1)+beta_drag(j,i))
1851  lgs_a_index(k) = nc
1852 
1853  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
1854  k = k+1
1855  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1856  lgs_a_value(k) = -inv_dxi_deta &
1857  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
1858  lgs_a_index(k) = nc
1859 
1860  nc = 2*nn(j,i+1)-1 ! next nc (column counter), for vx_m(j,i+1)
1861  k = k+1
1862  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1863  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i+1)
1864  lgs_a_index(k) = nc
1865 
1866  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
1867  k = k+1
1868  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1869  lgs_a_value(k) = inv_dxi_deta &
1870  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j,i))
1871  lgs_a_index(k) = nc
1872 
1873  nc = 2*nn(j+1,i)-1 ! largest nc (column counter), for vx_m(j+1,i)
1874  k = k+1
1875  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1876  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
1877  lgs_a_index(k) = nc
1878 
1879  lgs_b_value(nr) = factor_rhs_1 &
1880  * (h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) &
1881  * dzs_dxi(j,i)
1882 
1883  lgs_x_value(nr) = vx_m(j,i)
1884 
1885  else if ( &
1886  ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j,i+1) &
1887  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1888  .or. &
1889  ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j,i+1) &
1890  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1891  ) then
1892  ! one neighbour is floating ice and the other is grounded ice
1893  ! (grounding line),
1894  ! but floating conditions are not satisfied
1895 
1896  k = k+1
1897  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1898  lgs_a_value(k) = 1.0_dp ! diagonal element only
1899  lgs_a_index(k) = nr
1900 
1901  lgs_b_value(nr) = vx_m(j,i)
1902  lgs_x_value(nr) = vx_m(j,i)
1903 
1904  else if ( &
1905  ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==1_i2b) ) &
1906  .or. &
1907  ( (maske(j,i)==1_i2b).and.(maske(j,i+1)==3_i2b) ) &
1908  ) then
1909  ! one neighbour is floating ice and the other is ice-free land;
1910  ! velocity assumed to be zero
1911 
1912  k = k+1
1913  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1914  lgs_a_value(k) = 1.0_dp ! diagonal element only
1915  lgs_a_index(k) = nr
1916 
1917  lgs_b_value(nr) = 0.0_dp
1918 
1919  lgs_x_value(nr) = 0.0_dp
1920 
1921  else if ( &
1922  ( flag_calving_front_1(j,i).and.flag_calving_front_2(j,i+1) ) &
1923  .or. &
1924  ( flag_calving_front_2(j,i).and.flag_calving_front_1(j,i+1) ) &
1925  ) then
1926  ! one neighbour is floating ice and the other is ocean
1927  ! (calving front)
1928 
1929  if (maske(j,i)==3_i2b) then
1930  i1 = i ! floating ice marker
1931  else ! maske(j,i+1)==3_i2b
1932  i1 = i+1 ! floating ice marker
1933  end if
1934 
1935  if ( &
1936  ( (maske(j,i1-1)==3_i2b).or.(maske(j,i1+1)==3_i2b) ) &
1937  .or. &
1938  ( (maske(j,i1-1)==0_i2b).or.(maske(j,i1+1)==0_i2b) ) &
1939  .or. &
1940  ( (maske(j,i1-1)==1_i2b).or.(maske(j,i1+1)==1_i2b) ) &
1941  ) then
1942  ! discretisation of the x-component of the BC at the calving front
1943 
1944  nc = 2*nn(j-1,i1) ! smallest nc (column counter), for vy_m(j-1,i1)
1945  k = k+1
1946  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1947  lgs_a_value(k) = -2.0_dp*inv_deta*vis_int_g(j,i1)
1948  lgs_a_index(k) = nc
1949 
1950  nc = 2*nn(j,i1-1)-1 ! next nc (column counter), for vx_m(j,i1-1)
1951  k = k+1
1952  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1953  lgs_a_value(k) = -4.0_dp*inv_dxi*vis_int_g(j,i1)
1954  lgs_a_index(k) = nc
1955 
1956  nc = 2*nn(j,i1)-1 ! next nc (column counter), for vx_m(j,i1)
1957  k = k+1
1958  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1959  lgs_a_value(k) = 4.0_dp*inv_dxi*vis_int_g(j,i1)
1960  lgs_a_index(k) = nc
1961 
1962  nc = 2*nn(j,i1) ! largest nc (column counter), for vy_m(j,i1)
1963  k = k+1
1964  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1965  lgs_a_value(k) = 2.0_dp*inv_deta*vis_int_g(j,i1)
1966  lgs_a_index(k) = nc
1967 
1968  lgs_b_value(nr) = factor_rhs_2 &
1969  *(h_c(j,i1)+h_t(j,i1))*(h_c(j,i1)+h_t(j,i1))
1970 
1971  lgs_x_value(nr) = vx_m(j,i)
1972 
1973  else ! (maske(j,i1-1)==2_i2b).and.(maske(j,i1+1)==2_i2b);
1974  ! velocity assumed to be zero
1975 
1976  k = k+1
1977  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1978  lgs_a_value(k) = 1.0_dp ! diagonal element only
1979  lgs_a_index(k) = nr
1980 
1981  lgs_b_value(nr) = 0.0_dp
1982 
1983  lgs_x_value(nr) = 0.0_dp
1984 
1985  end if
1986 
1987  else if ( (maske(j,i)==0_i2b).or.(maske(j,i+1)==0_i2b) ) then
1988  ! neither neighbour is floating ice, but at least one neighbour is
1989  ! grounded ice; velocity taken from the solution for grounded ice
1990 
1991  k = k+1
1992  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
1993  lgs_a_value(k) = 1.0_dp ! diagonal element only
1994  lgs_a_index(k) = nr
1995 
1996  lgs_b_value(nr) = vx_m(j,i)
1997 
1998  lgs_x_value(nr) = vx_m(j,i)
1999 
2000  else ! neither neighbour is floating or grounded ice,
2001  ! velocity assumed to be zero
2002 
2003  k = k+1
2004  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2005  lgs_a_value(k) = 1.0_dp ! diagonal element only
2006  lgs_a_index(k) = nr
2007 
2008  lgs_b_value(nr) = 0.0_dp
2009 
2010  lgs_x_value(nr) = 0.0_dp
2011 
2012  end if
2013 
2014  else ! boundary condition, velocity assumed to be zero
2015 
2016  k = k+1
2017  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2018  lgs_a_value(k) = 1.0_dp ! diagonal element only
2019  lgs_a_index(k) = nr
2020 
2021  lgs_b_value(nr) = 0.0_dp
2022 
2023  lgs_x_value(nr) = 0.0_dp
2024 
2025  end if
2026 
2027  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
2028 
2029 ! ------ Equations for vy_m (at (i,j+1/2))
2030 
2031  nr = n+1 ! row counter
2032 
2033  if ( (j /= jmax).and.(i /= 0).and.(i /= imax) ) then
2034  ! inner point on the staggered grid in y-direction
2035 
2036  h_mid = 0.5_dp*(h_c(j+1,i)+h_t(j+1,i)+h_c(j,i)+h_t(j,i))
2037  zl_mid = 0.5_dp*(zl(j+1,i)+zl(j,i))
2038 
2039  if ( &
2040  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==3_i2b) ) &
2041  .or. &
2042  ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j+1,i) &
2043  .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
2044  .or. &
2045  ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j+1,i) &
2046  .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
2047  ) then
2048  ! both neighbours are floating ice,
2049  ! or one neighbour is floating ice and the other is grounded ice
2050  ! (grounding line)
2051  ! and floating conditions are satisfied;
2052  ! discretisation of the x-component of the PDE
2053 
2054  flag_shelfy_stream_y(j,i)=.false.
2055  ! make sure not to treat as shelfy stream
2056 
2057  nc = 2*nn(j-1,i) ! smallest nc (column counter), for vy_m(j-1,i)
2058  k = k+1
2059  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2060  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j,i)
2061  lgs_a_index(k) = nc
2062 
2063  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
2064  k = k+1
2065  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2066  lgs_a_value(k) = inv_dxi_deta &
2067  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i-1))
2068  lgs_a_index(k) = nc
2069 
2070  nc = 2*nn(j,i-1) ! next nc (column counter), for vy_m(j,i-1)
2071  k = k+1
2072  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2073  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
2074  lgs_a_index(k) = nc
2075 
2076  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
2077  k = k+1
2078  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2079  lgs_a_value(k) = -inv_dxi_deta &
2080  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
2081  lgs_a_index(k) = nc
2082 
2083  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
2084  if (nc /= nr) then ! (diagonal element)
2085  errormsg = ' >>> calc_vxy_ssa_matrix: ' &
2086  //'Check for diagonal element failed!'
2087  call error(errormsg)
2088  end if
2089  k = k+1
2090  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2091  lgs_a_value(k) = -4.0_dp*inv_deta2 &
2092  *(vis_int_g(j+1,i)+vis_int_g(j,i)) &
2093  -inv_dxi2 &
2094  *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1))
2095  lgs_a_index(k) = nc
2096 
2097  nc = 2*nn(j+1,i-1)-1 ! next nc (column counter), for vx_m(j+1,i-1)
2098  k = k+1
2099  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2100  lgs_a_value(k) = -inv_dxi_deta &
2101  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
2102  lgs_a_index(k) = nc
2103 
2104  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
2105  k = k+1
2106  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2107  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
2108  lgs_a_index(k) = nc
2109 
2110  nc = 2*nn(j+1,i)-1 ! next nc (column counter), for vx_m(j+1,i)
2111  k = k+1
2112  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2113  lgs_a_value(k) = inv_dxi_deta &
2114  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i))
2115  lgs_a_index(k) = nc
2116 
2117  nc = 2*nn(j+1,i) ! largest nc (column counter), for vy_m(j+1,i)
2118  k = k+1
2119  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2120  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j+1,i)
2121  lgs_a_index(k) = nc
2122 
2123  lgs_b_value(nr) = factor_rhs_1 &
2124  * (h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) &
2125  * dzs_deta(j,i)
2126 
2127  lgs_x_value(nr) = vy_m(j,i)
2128 
2129  else if (flag_shelfy_stream_y(j,i)) then
2130  ! shelfy stream (as determined by routine calc_vxy_sia)
2131 
2132  nc = 2*nn(j-1,i) ! smallest nc (column counter), for vy_m(j-1,i)
2133  k = k+1
2134  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2135  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j,i)
2136  lgs_a_index(k) = nc
2137 
2138  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
2139  k = k+1
2140  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2141  lgs_a_value(k) = inv_dxi_deta &
2142  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i-1))
2143  lgs_a_index(k) = nc
2144 
2145  nc = 2*nn(j,i-1) ! next nc (column counter), for vy_m(j,i-1)
2146  k = k+1
2147  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2148  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
2149  lgs_a_index(k) = nc
2150 
2151  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
2152  k = k+1
2153  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2154  lgs_a_value(k) = -inv_dxi_deta &
2155  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
2156  lgs_a_index(k) = nc
2157 
2158  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
2159  if (nc /= nr) then ! (diagonal element)
2160  errormsg = ' >>> calc_vxy_ssa_matrix: ' &
2161  //'Check for diagonal element failed!'
2162  call error(errormsg)
2163  end if
2164  k = k+1
2165  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2166  lgs_a_value(k) = -4.0_dp*inv_deta2 &
2167  *(vis_int_g(j+1,i)+vis_int_g(j,i)) &
2168  -inv_dxi2 &
2169  *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1)) &
2170  -0.5_dp*(beta_drag(j+1,i)+beta_drag(j,i))
2171  lgs_a_index(k) = nc
2172 
2173  nc = 2*nn(j+1,i-1)-1 ! next nc (column counter), for vx_m(j+1,i-1)
2174  k = k+1
2175  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2176  lgs_a_value(k) = -inv_dxi_deta &
2177  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
2178  lgs_a_index(k) = nc
2179 
2180  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
2181  k = k+1
2182  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2183  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
2184  lgs_a_index(k) = nc
2185 
2186  nc = 2*nn(j+1,i)-1 ! next nc (column counter), for vx_m(j+1,i)
2187  k = k+1
2188  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2189  lgs_a_value(k) = inv_dxi_deta &
2190  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i))
2191  lgs_a_index(k) = nc
2192 
2193  nc = 2*nn(j+1,i) ! largest nc (column counter), for vy_m(j+1,i)
2194  k = k+1
2195  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2196  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j+1,i)
2197  lgs_a_index(k) = nc
2198 
2199  lgs_b_value(nr) = factor_rhs_1 &
2200  * (h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) &
2201  * dzs_deta(j,i)
2202 
2203  lgs_x_value(nr) = vy_m(j,i)
2204 
2205  else if ( &
2206  ( flag_grounding_line_1(j,i).and.flag_grounding_line_2(j+1,i) &
2207  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
2208  .or. &
2209  ( flag_grounding_line_2(j,i).and.flag_grounding_line_1(j+1,i) &
2210  .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
2211  ) then
2212  ! one neighbour is floating ice and the other is grounded ice
2213  ! (grounding line),
2214  ! but floating conditions are not satisfied
2215 
2216  k = k+1
2217  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2218  lgs_a_value(k) = 1.0_dp ! diagonal element only
2219  lgs_a_index(k) = nr
2220 
2221  lgs_b_value(nr) = vy_m(j,i)
2222  lgs_x_value(nr) = vy_m(j,i)
2223 
2224  else if ( &
2225  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==1_i2b) ) &
2226  .or. &
2227  ( (maske(j,i)==1_i2b).and.(maske(j+1,i)==3_i2b) ) &
2228  ) then
2229  ! one neighbour is floating ice and the other is ice-free land;
2230  ! velocity assumed to be zero
2231 
2232  k = k+1
2233  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2234  lgs_a_value(k) = 1.0_dp ! diagonal element only
2235  lgs_a_index(k) = nr
2236 
2237  lgs_b_value(nr) = 0.0_dp
2238 
2239  lgs_x_value(nr) = 0.0_dp
2240 
2241  else if ( &
2242  ( flag_calving_front_1(j,i).and.flag_calving_front_2(j+1,i) ) &
2243  .or. &
2244  ( flag_calving_front_2(j,i).and.flag_calving_front_1(j+1,i) ) &
2245  ) then
2246  ! one neighbour is floating ice and the other is ocean
2247  ! (calving front)
2248 
2249  if (maske(j,i)==3_i2b) then
2250  j1 = j ! floating ice marker
2251  else ! maske(j+1,i)==3_i2b
2252  j1 = j+1 ! floating ice marker
2253  end if
2254 
2255  if ( &
2256  ( (maske(j1-1,i)==3_i2b).or.(maske(j1+1,i)==3_i2b) ) &
2257  .or. &
2258  ( (maske(j1-1,i)==0_i2b).or.(maske(j1+1,i)==0_i2b) ) &
2259  .or. &
2260  ( (maske(j1-1,i)==1_i2b).or.(maske(j1+1,i)==1_i2b) ) &
2261  ) then
2262  ! discretisation of the y-component of the BC at the calving front
2263 
2264  nc = 2*nn(j1-1,i) ! smallest nc (column counter), for vy_m(j1-1,i)
2265  k = k+1
2266  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2267  lgs_a_value(k) = -4.0_dp*inv_deta*vis_int_g(j1,i)
2268  lgs_a_index(k) = nc
2269 
2270  nc = 2*nn(j1,i-1)-1 ! next nc (column counter), for vx_m(j1,i-1)
2271  k = k+1
2272  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2273  lgs_a_value(k) = -2.0_dp*inv_dxi*vis_int_g(j1,i)
2274  lgs_a_index(k) = nc
2275 
2276  nc = 2*nn(j1,i)-1 ! next nc (column counter), for vx_m(j1,i)
2277  k = k+1
2278  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2279  lgs_a_value(k) = 2.0_dp*inv_dxi*vis_int_g(j1,i)
2280  lgs_a_index(k) = nc
2281 
2282  nc = 2*nn(j1,i) ! largest nc (column counter), for vy_m(j1,i)
2283  k = k+1
2284  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2285  lgs_a_value(k) = 4.0_dp*inv_deta*vis_int_g(j1,i)
2286  lgs_a_index(k) = nc
2287 
2288  lgs_b_value(nr) = factor_rhs_2 &
2289  *(h_c(j1,i)+h_t(j1,i))*(h_c(j1,i)+h_t(j1,i))
2290 
2291  lgs_x_value(nr) = vy_m(j,i)
2292 
2293  else ! (maske(j1-1,i)==2_i2b).and.(maske(j1+1,i)==2_i2b);
2294  ! velocity assumed to be zero
2295 
2296  k = k+1
2297  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2298  lgs_a_value(k) = 1.0_dp ! diagonal element only
2299  lgs_a_index(k) = nr
2300 
2301  lgs_b_value(nr) = 0.0_dp
2302 
2303  lgs_x_value(nr) = 0.0_dp
2304 
2305  end if
2306 
2307  else if ( (maske(j,i)==0_i2b).or.(maske(j+1,i)==0_i2b) ) then
2308  ! neither neighbour is floating ice, but at least one neighbour is
2309  ! grounded ice; velocity taken from the solution for grounded ice
2310 
2311  k = k+1
2312  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2313  lgs_a_value(k) = 1.0_dp ! diagonal element only
2314  lgs_a_index(k) = nr
2315 
2316  lgs_b_value(nr) = vy_m(j,i)
2317 
2318  lgs_x_value(nr) = vy_m(j,i)
2319 
2320  else ! neither neighbour is floating or grounded ice,
2321  ! velocity assumed to be zero
2322 
2323  k = k+1
2324  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2325  lgs_a_value(k) = 1.0_dp ! diagonal element only
2326  lgs_a_index(k) = nr
2327 
2328  lgs_b_value(nr) = 0.0_dp
2329  lgs_x_value(nr) = 0.0_dp
2330 
2331  end if
2332 
2333  else ! boundary condition, velocity assumed to be zero
2334 
2335  k = k+1
2336  ! if (k > n_sprs) stop ' >>> calc_vxy_ssa_matrix: n_sprs too small!'
2337  lgs_a_value(k) = 1.0_dp ! diagonal element only
2338  lgs_a_index(k) = nr
2339 
2340  lgs_b_value(nr) = 0.0_dp
2341  lgs_x_value(nr) = 0.0_dp
2342 
2343  end if
2344 
2345  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
2346 
2347 end do
2348 
2349 !-------- Settings for Lis --------
2350 
2351 #ifndef ALLOW_OPENAD /* Normal */
2352 
2353 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
2354 call lis_vector_create(lis_comm_world, lgs_b, ierr)
2355 call lis_vector_create(lis_comm_world, lgs_x, ierr)
2356 
2357 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
2358 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
2359 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
2360 
2361 do nr=1, nmax
2362 
2363  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
2364  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
2365  lgs_a_value(nc), lgs_a, ierr)
2366  end do
2367 
2368  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
2369  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
2370 
2371 end do
2372 
2373 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
2374 call lis_matrix_assemble(lgs_a, ierr)
2375 
2376 !-------- Solution of the system of linear equations with Lis --------
2377 
2378 call lis_solver_create(solver, ierr)
2379 
2380 #if (defined(LIS_OPTS))
2381  ch_solver_set_option = trim(lis_opts)
2382 #else
2383  ch_solver_set_option = '-i bicgsafe -p jacobi '// &
2384  '-maxiter 100 -tol 1.0e-4 -initx_zeros false'
2385 #endif
2386 
2387 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
2388 call chkerr(ierr)
2389 
2390 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
2391 call chkerr(ierr)
2392 
2393 call lis_solver_get_iter(solver, iter, ierr)
2394 
2395 if (iter_ssa == 1) then
2396  write(6,'(10x,a,i0)') 'calc_vxy_ssa_matrix: LIS iter = ', iter
2397 else if (m == 1) then
2398  write(6,'(10x,a,i0)', advance='no') 'calc_vxy_ssa_matrix: LIS iter = ', iter
2399 else if (m == iter_ssa) then
2400  write(6,'(a,i0)') ', ', iter
2401 else
2402  write(6,'(a,i0)', advance='no') ', ', iter
2403 end if
2404 
2405 !!! call lis_solver_get_time(solver,solver_time,ierr)
2406 !!! print *, 'calc_vxy_ssa_matrix: time (s) = ', solver_time
2407 
2408 lgs_x_value = 0.0_dp
2409 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
2410 call lis_matrix_destroy(lgs_a, ierr)
2411 call lis_vector_destroy(lgs_b, ierr)
2412 call lis_vector_destroy(lgs_x, ierr)
2413 call lis_solver_destroy(solver, ierr)
2414 
2415 #else /* OpenAD */
2416 
2417 lgs_a_index_pass = int(lgs_a_index)
2418 call sico_lis_solver(nmax, n_sprs, &
2419  lgs_a_ptr, lgs_a_index_pass, &
2420  lgs_a_value, lgs_b_value, lgs_x_value)
2421 
2422 #endif /* Normal vs. OpenAD */
2423 
2424 do n=1, nmax-1, 2
2425 
2426  i = ii((n+1)/2)
2427  j = jj((n+1)/2)
2428 
2429  nr = n
2430  vx_m(j,i) = lgs_x_value(nr)
2431 
2432  nr = n+1
2433  vy_m(j,i) = lgs_x_value(nr)
2434 
2435 end do
2436 
2437 #ifndef ALLOW_OPENAD /* Normal */
2438 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
2439 deallocate(lgs_b_value, lgs_x_value)
2440 #endif /* Normal */
2441 
2442 #else
2443 
2444 errormsg = ' >>> calc_vxy_ssa_matrix: ' &
2445  //'Only to be called for MARGIN==3 or DYNAMICS==2!'
2446 call error(errormsg)
2447 
2448 #endif
2449 
2450 end subroutine calc_vxy_ssa_matrix
2451 
2452 !-------------------------------------------------------------------------------
2453 !> Computation of the depth-integrated viscosity vis_int_g in the
2454 !! shallow shelf approximation.
2455 !<------------------------------------------------------------------------------
2456 subroutine calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)
2457 
2459 
2460 implicit none
2461 
2462 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t
2463 
2464 integer(i4b) :: i, j, kc, kt
2465 real(dp) :: dxi_inv, deta_inv
2466 real(dp) :: dvx_dxi, dvx_deta, dvy_dxi, dvy_deta
2467 real(dp) :: aqxy1(0:kcmax)
2468 real(dp) :: cvis0(0:ktmax), cvis1(0:kcmax)
2469 
2470 #if (MARGIN==3 || DYNAMICS==2)
2471 
2472 !-------- Term abbreviations --------
2473 
2474 dxi_inv = 1.0_dp/dxi
2475 deta_inv = 1.0_dp/deta
2476 
2477 do kc=0, kcmax
2478  if (flag_aa_nonzero) then
2479  aqxy1(kc) = aa/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
2480  else
2481  aqxy1(kc) = dzeta_c
2482  end if
2483 end do
2484 
2485 !-------- Computation of the depth-integrated viscosity --------
2486 
2487 do i=0, imax
2488 do j=0, jmax
2489 
2490  if ((maske(j,i)==0_i2b).and.(.not.flag_shelfy_stream(j,i))) then
2491  ! grounded ice, but
2492  ! not shelfy stream
2493  de_ssa(j,i) = 0.0_dp ! dummy value
2494 
2495  vis_int_g(j,i) = (h_c(j,i)+h_t(j,i)) / flui_ave_sia(j,i)
2496 
2497  else if ((maske(j,i)==1_i2b).or.(maske(j,i)==2_i2b)) then
2498  ! ice-free land or ocean
2499  de_ssa(j,i) = 0.0_dp ! dummy value
2500 
2501  vis_int_g(j,i) = 0.0_dp ! dummy value
2502 
2503  else ! (maske(j,i)==3_i2b).or.(flag_shelfy_stream(j,i)),
2504  ! floating ice or shelfy stream;
2505  ! must not be at the margin of the computational domain
2506 
2507 ! ------ Effective strain rate
2508 
2509  dvx_dxi = (vx_m(j,i)-vx_m(j,i-1))*dxi_inv
2510  dvy_deta = (vy_m(j,i)-vy_m(j-1,i))*deta_inv
2511 
2512  dvx_deta = 0.25_dp*deta_inv &
2513  *(vx_m(j+1,i)+vx_m(j+1,i-1)-vx_m(j-1,i)-vx_m(j-1,i-1))
2514  dvy_dxi = 0.25_dp*dxi_inv &
2515  *(vy_m(j,i+1)+vy_m(j-1,i+1)-vy_m(j,i-1)-vy_m(j-1,i-1))
2516 
2517 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
2518 
2519  de_ssa(j,i) = sqrt( dvx_dxi*dvx_dxi &
2520  + dvy_deta*dvy_deta &
2521  + dvx_dxi*dvy_deta &
2522  + 0.25_dp*(dvx_deta+dvy_dxi)*(dvx_deta+dvy_dxi) )
2523 
2524 #else /* OpenAD: guarding against non-differentiable sqrt(0) */
2525 
2526  if ( ( dvx_dxi*dvx_dxi &
2527  + dvy_deta*dvy_deta &
2528  + dvx_dxi*dvy_deta &
2529  + 0.25_dp*(dvx_deta+dvy_dxi)*(dvx_deta+dvy_dxi) ) > 0 ) then
2530  de_ssa(j,i) = sqrt( dvx_dxi*dvx_dxi &
2531  + dvy_deta*dvy_deta &
2532  + dvx_dxi*dvy_deta &
2533  + 0.25_dp*(dvx_deta+dvy_dxi)*(dvx_deta+dvy_dxi) )
2534  else
2535  de_ssa(j,i) = 0.0_dp
2536  end if
2537 
2538 #endif /* Normal vs. OpenAD */
2539 
2540 ! ------ Term abbreviations
2541 
2542 #if (DYNAMICS==2)
2543  if (.not.flag_shelfy_stream(j,i)) then
2544 #endif
2545 
2546  do kc=0, kcmax
2547  cvis1(kc) = aqxy1(kc)*h_c(j,i) &
2548  *viscosity(de_ssa(j,i), &
2549  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2550  enh_c(kc,j,i), 0_i2b)
2551  end do
2552  ! Ice shelves (floating ice) are assumed to consist of cold ice only
2553 
2554 #if (DYNAMICS==2)
2555  else ! flag_shelfy_stream(j,i) == .true.
2556 
2557 #if (CALCMOD==-1 || CALCMOD==0)
2558 
2559  do kc=0, kcmax
2560  cvis1(kc) = aqxy1(kc)*h_c(j,i) &
2561  *viscosity(de_ssa(j,i), &
2562  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2563  enh_c(kc,j,i), 0_i2b)
2564  end do
2565 
2566 #elif (CALCMOD==1)
2567 
2568  do kt=0, ktmax
2569  cvis0(kt) = dzeta_t*h_t(j,i) &
2570  *viscosity(de_ssa(j,i), &
2571  temp_t_m(kt,j,i), temp_t_m(kt,j,i), omega_t(kt,j,i), &
2572  enh_t(kt,j,i), 1_i2b)
2573  end do
2574 
2575  do kc=0, kcmax
2576  cvis1(kc) = aqxy1(kc)*h_c(j,i) &
2577  *viscosity(de_ssa(j,i), &
2578  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2579  enh_c(kc,j,i), 0_i2b)
2580  end do
2581 
2582 #elif (CALCMOD==2 || CALCMOD==3)
2583 
2584  do kc=0, kcmax
2585  cvis1(kc) = aqxy1(kc)*h_c(j,i) &
2586  *viscosity(de_ssa(j,i), &
2587  temp_c(kc,j,i), temp_c_m(kc,j,i), omega_c(kc,j,i), &
2588  enh_c(kc,j,i), 2_i2b)
2589  end do
2590 
2591 #else
2592  errormsg = ' >>> calc_vis_ssa: CALCMOD must be -1, 0, 1, 2 or 3!'
2593  call error(errormsg)
2594 #endif
2595 
2596  end if
2597 
2598 #endif
2599 
2600 ! ------ Depth-integrated viscosity
2601 
2602  vis_int_g(j,i) = 0.0_dp
2603 
2604 #if (CALCMOD==1)
2605  do kt=0, ktmax-1
2606  vis_int_g(j,i) = vis_int_g(j,i)+0.5_dp*(cvis0(kt+1)+cvis0(kt))
2607  end do
2608 #endif
2609 
2610  do kc=0, kcmax-1
2611  vis_int_g(j,i) = vis_int_g(j,i)+0.5_dp*(cvis1(kc+1)+cvis1(kc))
2612  end do
2613 
2614  end if
2615 
2616 end do
2617 end do
2618 
2619 #else
2620 
2621 errormsg = ' >>> calc_vis_ssa: Only to be called for MARGIN==3 or DYNAMICS==2!'
2622 call error(errormsg)
2623 
2624 #endif
2625 
2626 end subroutine calc_vis_ssa
2627 
2628 !-------------------------------------------------------------------------------
2629 
2630 end module calc_vxy_m
2631 !
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)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
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:kcmax, 0:jmax, 0:imax) d_help_c
d_help_c(kc,j,i): Auxiliary quantity for the computation of vx, vy und zs
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
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)
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
real(dp), dimension(0:jmax, 0:imax) vx_s_g
vx_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
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...
Several mathematical tools used by SICOPOLIS.
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) p_b_w
p_b_w(j,i): Basal water pressure
real(dp), dimension(0:jmax, 0:imax) vy_b
vy_b(j,i): Velocity in y-direction at the ice base, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_x
flag_shelfy_stream_x(j,i): Shelfy stream flag in x-direction, at (i+1/2,j). .true.: shelfy stream point .false.: otherwise
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) txz_c
txz_c(kc,j,i): Shear stress txz in the upper (kc) ice domain (at (i+1/2,j,kc))
integer(i2b), dimension(0:jmax, 0:imax) maske_sedi
maske_sedi(j,i): Sediment mask. 1: hard rock, 7: soft sediment, 2: ocean.
Definition: sico_vars_m.F90:46
real(dp), dimension(0:jmax, 0:imax) dzs_dxi
dzs_dxi(j,i): Derivative of zs by xi (at i+1/2,j)
real(dp) function, public creep(sigma_val)
Creep response function for ice.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
real(dp) g
G: Acceleration due to gravity.
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
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)
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp), dimension(0:jmax, 0:imax) vis_int_g
vis_int_g(j,i): Depth-integrated viscosity of the SSA, at (i,j)
real(dp), dimension(0:jmax, 0:imax) ratio_sl_y
ratio_sl_y(j,i): Ratio of basal to surface velocity (slip ratio) in y-direction, at (i...
Computation of the horizontal velocity vx, vy.
Definition: calc_vxy_m.F90:36
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
real(dp), dimension(0:jmax, 0:imax) d_help_b
d_help_b(j,i): Auxiliary quantity for the computation of vx_b and vy_b
real(dp) function, public ratefac_c(temp_val, temp_m_val)
Rate factor for cold ice: Linear interpolation of tabulated values in RF(.).
real(dp) ea
ea: Abbreviation for exp(aa)
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
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:ktmax, 0:jmax, 0:imax) d_help_t
d_help_t(kt,j,i): Auxiliary quantity for the computation of vx, vy und zs
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) flui_ave_sia
flui_ave_sia(j,i): Depth-averaged fluidity of the SIA
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) txz_t
txz_t(kt,j,i): Shear stress txz in the lower (kt) ice domain (at (i+1/2,j,kt))
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
real(dp) function, public ratefac_t(omega_val)
Rate factor for temperate ice.
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) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_m
temp_c_m(kc,j,i): Melting temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) qx
qx(j,i): Volume flux in x-direction (depth-integrated vx, at (i+1/2,j))
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
real(dp), dimension(0:jmax, 0:imax) dzs_deta
dzs_deta(j,i): Derivative of zs by eta (at i,j+1/2)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) tyz_c
tyz_c(kc,j,i): Shear stress tyz in the upper (kc) ice domain (at (i,j+1/2,kc))
real(dp), dimension(0:jmax, 0:imax) q_gl_g
q_gl_g(j,i): Volume flux across the grounding line, at (i,j)
real(dp), dimension(0:jmax, 0:imax) vy_s_g
vy_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_m.F90:497
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
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) ratio_sl_x
ratio_sl_x(j,i): Ratio of basal to surface velocity (slip ratio) in x-direction, at (i+1/2...
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
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), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
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) function, public viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, i_flag_cold_temp)
Ice viscosity as a function of the effective strain rate and the temperature (in cold ice) or the wat...
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_y
flag_shelfy_stream_y(j,i): Shelfy stream flag in y-direction, at (i,j+1/2). .true.: shelfy stream point .false.: otherwise
real(dp) function, public ratefac_c_t(temp_val, omega_val, temp_m_val)
Rate factor for cold and temperate ice: Combination of ratefac_c and ratefac_t (only for the enthalpy...
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) vx_b
vx_b(j,i): Velocity in x-direction at the ice base, at (i+1/2,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) sigma_t
sigma_t(kt,j,i): Effective stress in the lower (kt) ice domain
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
integer, parameter dp
Double-precision reals.
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) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
subroutine, public calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
Definition: calc_vxy_m.F90:55
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) tyz_t
tyz_t(kt,j,i): Shear stress tyz in the lower (kt) ice domain (at (i,j+1/2,kt))
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
real(dp), dimension(0:jmax, 0:imax) qy
qy(j,i): Volume flux in y-direction (depth-integrated vy, at (i,j+1/2))