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