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