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