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