SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_vxy_ssa_matrix.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ v x y _ s s a _ m a t r i x
4 !
5 !> @file
6 !!
7 !! Solution of the system of linear equations for the horizontal velocities
8 !! vx_m, vy_m in the shallow shelf approximation.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2013 Ralf Greve, Tatsuru Sato
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 !> Solution of the system of linear equations for the horizontal velocities
35 !! vx_m, vy_m in the shallow shelf approximation.
36 !<------------------------------------------------------------------------------
37 subroutine calc_vxy_ssa_matrix(z_sl,dxi, deta, ii, jj, nn)
38 
39 use sico_types
40 use sico_variables
41 
42 implicit none
43 
44 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
45 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
46 real(dp), intent(in) :: dxi, deta, z_sl
47 
48 integer(i4b) :: i, j, k, n
49 integer(i4b) :: i1, j1
50 real(dp) :: inv_dxi, inv_deta, inv_dxi_deta, inv_dxi2, inv_deta2
51 real(dp) :: factor_rhs_1, factor_rhs_2, tau_b, rhosw_rho_ratio
52 real(dp) :: h_mid, zl_mid
53 real(dp), dimension(0:JMAX,0:IMAX) :: vis_int_sgxy
54 character(len=256) :: ch_solver_set_option
55 
56 lis_integer :: ierr
57 lis_integer :: nc, nr
58 ! LIS_INTEGER :: iter
59 lis_matrix :: lgs_a
60 lis_vector :: lgs_b, lgs_x
61 lis_solver :: solver
62 
63 lis_integer, parameter :: nmax = 2*(imax+1)*(jmax+1)
64 lis_integer, parameter :: n_sprs = 20*(imax+1)*(jmax+1)
65 lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
66 lis_scalar, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
67 
68 !-------- Abbreviations --------
69 
70 inv_dxi = 1.0_dp/dxi
71 inv_deta = 1.0_dp/deta
72 inv_dxi_deta = 1.0_dp/(dxi*deta)
73 inv_dxi2 = 1.0_dp/(dxi*dxi)
74 inv_deta2 = 1.0_dp/(deta*deta)
75 
76 rhosw_rho_ratio = rho_sw/rho
77 factor_rhs_1 = 0.5_dp*rho*g
78 factor_rhs_2 = 0.5_dp*rho*g*(rho_sw-rho)/rho_sw
79 
80 !-------- Depth-integrated viscosity on the staggered grid
81 ! [at (i+1/2,j+1/2)] --------
82 
83 vis_int_sgxy = 0.0_dp ! initialisation
84 
85 do i=0, imax-1
86 do j=0, jmax-1
87 
88  k=0
89 
90  if ((maske(j,i)==0_i2b).or.(maske(j,i)==3_i2b)) then
91  k = k+1 ! floating or grounded ice
92  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j,i)
93  end if
94 
95  if ((maske(j,i+1)==0_i2b).or.(maske(j,i+1)==3_i2b)) then
96  k = k+1 ! floating or grounded ice
97  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j,i+1)
98  end if
99 
100  if ((maske(j+1,i)==0_i2b).or.(maske(j+1,i)==3_i2b)) then
101  k = k+1 ! floating or grounded ice
102  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j+1,i)
103  end if
104 
105  if ((maske(j+1,i+1)==0_i2b).or.(maske(j+1,i+1)==3_i2b)) then
106  k = k+1 ! floating or grounded ice
107  vis_int_sgxy(j,i) = vis_int_sgxy(j,i) + vis_int_g(j+1,i+1)
108  end if
109 
110  if (k>0) vis_int_sgxy(j,i) = vis_int_sgxy(j,i)/real(k,dp)
111 
112 end do
113 end do
114 
115 !-------- Assembly of the system of linear equations
116 ! (matrix storage: compressed sparse row CSR) --------
117 
118 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
119 allocate(lgs_b_value(nmax), lgs_x_value(nmax))
120 
121 lgs_a_value = 0.0_dp
122 lgs_a_index = 0
123 lgs_a_ptr = 0
124 
125 lgs_b_value = 0.0_dp
126 lgs_x_value = 0.0_dp
127 
128 lgs_a_ptr(1) = 1
129 
130 k = 0
131 
132 do n=1, nmax-1, 2
133 
134  i = ii((n+1)/2)
135  j = jj((n+1)/2)
136 
137 ! ------ Equations for vx_m
138 
139  nr = n ! row counter
140 
141  if ( (i /= imax).and.(j /= 0).and.(j /= jmax) ) then
142 
143  h_mid=0.50_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))
144  zl_mid=0.50_dp*(zl(j,i)+zl(j,i+1))
145 
146  ! inner point on the staggered grid in x-direction
147  if ( &
148  ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==3_i2b) ) &
149  .or. &
150  ((maske(j,i)==3_i2b).and.(maske(j,i+1)==0_i2b) &
151  .and.(h_mid<(z_sl-zl_mid)*rhosw_rho_ratio)) &
152  .or. &
153  ((maske(j,i)==0_i2b).and.(maske(j,i+1)==3_i2b) &
154  .and.(h_mid<(z_sl-zl_mid)*rhosw_rho_ratio)) &
155  ) then
156  ! both neighbours are floating ice,
157  ! or one neighbour is floating ice and the other is grounded ice
158  ! and floating conditions are satisfied as a shelf;
159  ! discretisation of the x-component of the PDE
160 
161  nc = 2*nn(j-1,i)-1 ! smallest nc (column counter), for vx_m(j-1,i)
162  k = k+1
163  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
164  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
165  lgs_a_index(k) = nc
166 
167  nc = 2*nn(j-1,i) ! next nc (column counter), for vy_m(j-1,i)
168  k = k+1
169  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
170  lgs_a_value(k) = inv_dxi_deta &
171  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j-1,i))
172  lgs_a_index(k) = nc
173 
174  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
175  k = k+1
176  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
177  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i)
178  lgs_a_index(k) = nc
179 
180  nc = 2*nn(j-1,i+1) ! next nc (column counter), for vy_m(j-1,i+1)
181  k = k+1
182  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
183  lgs_a_value(k) = -inv_dxi_deta &
184  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
185  lgs_a_index(k) = nc
186 
187  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
188  if (nc /= nr) & ! (diagonal element)
189  stop ' calc_vxy_ssa_matrix: Check for diagonal element failed!'
190  k = k+1
191  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
192  lgs_a_value(k) = -4.0_dp*inv_dxi2 &
193  *(vis_int_g(j,i+1)+vis_int_g(j,i)) &
194  -inv_deta2 &
195  *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i))
196  lgs_a_index(k) = nc
197 
198  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
199  k = k+1
200  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
201  lgs_a_value(k) = -inv_dxi_deta &
202  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
203  lgs_a_index(k) = nc
204 
205  nc = 2*nn(j,i+1)-1 ! next nc (column counter), for vx_m(j,i+1)
206  k = k+1
207  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
208  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i+1)
209  lgs_a_index(k) = nc
210 
211  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
212  k = k+1
213  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
214  lgs_a_value(k) = inv_dxi_deta &
215  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j,i))
216  lgs_a_index(k) = nc
217 
218  nc = 2*nn(j+1,i)-1 ! largest nc (column counter), for vx_m(j+1,i)
219  k = k+1
220  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
221  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
222  lgs_a_index(k) = nc
223 
224  lgs_b_value(nr) = factor_rhs_1 &
225  * (h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) &
226  * dzs_dxi(j,i)
227 
228  lgs_x_value(nr) = vx_m(j,i)
229 
230 ! if floating conditions are not satisfied as a shelf
231 ! use ice sheet value or shelfy stream
232  else if ( ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==0_i2b) &
233  .and.(h_mid>(z_sl-zl(j,i)+zl(j,i+1))*rhosw_rho_ratio)) &
234  .or. &
235  ( (maske(j,i)==0_i2b).and.(maske(j,i+1)==3_i2b) &
236  .and.(h_mid>(z_sl-zl_mid)*rhosw_rho_ratio) ) &
237  ) then
238 
239 #if ( !defined(SHS) || SHS==0 )
240  k = k+1
241  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
242  lgs_a_value(k) = 1.0_dp ! diagonal element only
243  lgs_a_index(k) = nr
244 
245  lgs_b_value(nr) = vx_m(j,i)
246  lgs_x_value(nr) = vx_m(j,i)
247 
248 #elif SHS==1
249  ! definition of basal conditions
250  tau_b=((year_sec/c_slide_ssa)**(1.0_dp/p_ssa)) &
251  *(rho*g*(h_c(j,i)+h_t(j,i))**(q_ssa/p_ssa)) &
252  /((vx_m(j,i)*vx_m(j,i)+vy_m(j,i)*vy_m(j,i))**(0.50_dp*(1.0_dp-1.0_dp/p_ssa)))
253 
254 ! define flag shelfy flow
255  flag_sf(j,i)=.true.
256  flag_sf(j,i-1)=.true.
257 
258  nc = 2*nn(j-1,i)-1 ! smallest nc (column counter), for vx_m(j-1,i)
259  k = k+1
260  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
261  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
262  lgs_a_index(k) = nc
263 
264  nc = 2*nn(j-1,i) ! next nc (column counter), for vy_m(j-1,i)
265  k = k+1
266  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
267  lgs_a_value(k) = inv_dxi_deta &
268  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j-1,i))
269  lgs_a_index(k) = nc
270 
271  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
272  k = k+1
273  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
274  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i)
275  lgs_a_index(k) = nc
276 
277  nc = 2*nn(j-1,i+1) ! next nc (column counter), for vy_m(j-1,i+1)
278  k = k+1
279  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
280  lgs_a_value(k) = -inv_dxi_deta &
281  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
282  lgs_a_index(k) = nc
283 
284  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
285  if (nc /= nr) & ! (diagonal element)
286  stop ' calc_vxy_ssa_matrix: Check for diagonal element failed!'
287  k = k+1
288  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
289  lgs_a_value(k) = -4.0_dp*inv_dxi2 &
290  *(vis_int_g(j,i+1)+vis_int_g(j,i)) &
291  -inv_deta2 &
292  *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i))-tau_b
293  lgs_a_index(k) = nc
294 
295  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
296  k = k+1
297  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
298  lgs_a_value(k) = -inv_dxi_deta &
299  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
300  lgs_a_index(k) = nc
301 
302  nc = 2*nn(j,i+1)-1 ! next nc (column counter), for vx_m(j,i+1)
303  k = k+1
304  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
305  lgs_a_value(k) = 4.0_dp*inv_dxi2*vis_int_g(j,i+1)
306  lgs_a_index(k) = nc
307 
308  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
309  k = k+1
310  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
311  lgs_a_value(k) = inv_dxi_deta &
312  *(2.0_dp*vis_int_g(j,i+1)+vis_int_sgxy(j,i))
313  lgs_a_index(k) = nc
314 
315  nc = 2*nn(j+1,i)-1 ! largest nc (column counter), for vx_m(j+1,i)
316  k = k+1
317  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
318  lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
319  lgs_a_index(k) = nc
320 
321  lgs_b_value(nr) = factor_rhs_1 &
322  * (h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) &
323  * dzs_dxi(j,i)
324 
325  lgs_x_value(nr) = vx_m(j,i)
326 
327 #endif
328 
329  else if ( &
330  ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==1_i2b) ) &
331  .or. &
332  ( (maske(j,i)==1_i2b).and.(maske(j,i+1)==3_i2b) ) &
333  ) then
334  ! one neighbour is floating ice and the other is ice-free land;
335  ! velocity assumed to be zero
336 
337 
338  k = k+1
339  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
340  lgs_a_value(k) = 1.0_dp ! diagonal element only
341  lgs_a_index(k) = nr
342 
343  lgs_b_value(nr) = 0.0_dp
344 
345  lgs_x_value(nr) = 0.0_dp
346 
347  else if ( &
348  ( (maske(j,i)==3_i2b).and.(maske(j,i+1)==2_i2b) ) &
349  .or. &
350  ( (maske(j,i)==2_i2b).and.(maske(j,i+1)==3_i2b) ) &
351  ) then
352  ! one neighbour is floating ice and the other is ocean
353 
354  if (maske(j,i)==3_i2b) then
355  i1 = i ! floating ice marker
356  else ! maske(j,i+1)==3_i2b
357  i1 = i+1 ! floating ice marker
358  end if
359 
360  if ( &
361  ( (maske(j,i1-1)==3_i2b).or.(maske(j,i1+1)==3_i2b) ) &
362  .or. &
363  ( (maske(j,i1-1)==0_i2b).or.(maske(j,i1+1)==0_i2b) ) &
364  .or. &
365  ( (maske(j,i1-1)==1_i2b).or.(maske(j,i1+1)==1_i2b) ) &
366  ) then
367  ! discretisation of the x-component of the BC at the calving front
368 
369  nc = 2*nn(j-1,i1) ! smallest nc (column counter), for vy_m(j-1,i1)
370  k = k+1
371  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
372  lgs_a_value(k) = -2.0_dp*inv_deta*vis_int_g(j,i1)
373  lgs_a_index(k) = nc
374 
375  nc = 2*nn(j,i1-1)-1 ! next nc (column counter), for vx_m(j,i1-1)
376  k = k+1
377  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
378  lgs_a_value(k) = -4.0_dp*inv_dxi*vis_int_g(j,i1)
379  lgs_a_index(k) = nc
380 
381  nc = 2*nn(j,i1)-1 ! next nc (column counter), for vx_m(j,i1)
382  k = k+1
383  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
384  lgs_a_value(k) = 4.0_dp*inv_dxi*vis_int_g(j,i1)
385  lgs_a_index(k) = nc
386 
387  nc = 2*nn(j,i1) ! largest nc (column counter), for vy_m(j,i1)
388  k = k+1
389  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
390  lgs_a_value(k) = 2.0_dp*inv_deta*vis_int_g(j,i1)
391  lgs_a_index(k) = nc
392 
393  lgs_b_value(nr) = factor_rhs_2 &
394  *(h_c(j,i1)+h_t(j,i1))*(h_c(j,i1)+h_t(j,i1))
395 
396  lgs_x_value(nr) = vx_m(j,i)
397 
398  else ! (maske(j,i1-1)==2_i2b).and.(maske(j,i1+1)==2_i2b);
399  ! velocity assumed to be zero
400 
401  k = k+1
402  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
403  lgs_a_value(k) = 1.0_dp ! diagonal element only
404  lgs_a_index(k) = nr
405 
406  lgs_b_value(nr) = 0.0_dp
407 
408  lgs_x_value(nr) = 0.0_dp
409 
410  end if
411 
412  else if ( (maske(j,i)==0_i2b).or.(maske(j,i+1)==0_i2b) ) then
413  ! neither neighbour is floating ice, but at least one neighbour is
414  ! grounded ice; velocity taken from the solution for grounded ice
415 
416  k = k+1
417  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
418  lgs_a_value(k) = 1.0_dp ! diagonal element only
419  lgs_a_index(k) = nr
420 
421  lgs_b_value(nr) = vx_m(j,i)
422 
423  lgs_x_value(nr) = vx_m(j,i)
424 
425  else ! neither neighbour is floating or grounded ice,
426  ! velocity assumed to be zero
427 
428  k = k+1
429  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
430  lgs_a_value(k) = 1.0_dp ! diagonal element only
431  lgs_a_index(k) = nr
432 
433  lgs_b_value(nr) = 0.0_dp
434 
435  lgs_x_value(nr) = 0.0_dp
436 
437  end if
438 
439  else ! boundary condition, velocity assumed to be zero
440 
441  k = k+1
442  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
443  lgs_a_value(k) = 1.0_dp ! diagonal element only
444  lgs_a_index(k) = nr
445 
446  lgs_b_value(nr) = 0.0_dp
447 
448  lgs_x_value(nr) = 0.0_dp
449 
450  end if
451 
452  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
453 
454 ! ------ Equations for vy_m
455 
456  nr = n+1 ! row counter
457 
458  if ( (j /= jmax).and.(i /= 0).and.(i /= imax) ) then
459  ! inner point on the staggered grid in y-direction
460  h_mid=0.50_dp*(h_c(j+1,i)+h_t(j+1,i)+h_c(j,i)+h_t(j,i))
461  zl_mid=0.50_dp*(zl(j+1,i)+zl(j,i))
462 
463  if ( &
464  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==3_i2b) ) &
465  .or. &
466  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==0_i2b) &
467  .and.(h_mid<(z_sl-zl_mid)*rhosw_rho_ratio)) &
468  .or. &
469  ( (maske(j,i)==0_i2b).and.(maske(j+1,i)==3_i2b) &
470  .and.(h_mid<(z_sl-zl_mid)*rhosw_rho_ratio)) &
471  ) then
472  ! both neighbours are floating ice,
473  ! or one neighbour is floating ice and the other is grounded ice;
474  ! discretisation of the y-component of the PDE
475 
476  nc = 2*nn(j-1,i) ! smallest nc (column counter), for vy_m(j-1,i)
477  k = k+1
478  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
479  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j,i)
480  lgs_a_index(k) = nc
481 
482  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
483  k = k+1
484  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
485  lgs_a_value(k) = inv_dxi_deta &
486  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i-1))
487  lgs_a_index(k) = nc
488 
489  nc = 2*nn(j,i-1) ! next nc (column counter), for vy_m(j,i-1)
490  k = k+1
491  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
492  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
493  lgs_a_index(k) = nc
494 
495  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
496  k = k+1
497  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
498  lgs_a_value(k) = -inv_dxi_deta &
499  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
500  lgs_a_index(k) = nc
501 
502  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
503  if (nc /= nr) & ! (diagonal element)
504  stop ' calc_vxy_ssa_matrix: Check for diagonal element failed!'
505  k = k+1
506  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
507  lgs_a_value(k) = -4.0_dp*inv_deta2 &
508  *(vis_int_g(j+1,i)+vis_int_g(j,i)) &
509  -inv_dxi2 &
510  *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1))
511  lgs_a_index(k) = nc
512 
513  nc = 2*nn(j+1,i-1)-1 ! next nc (column counter), for vx_m(j+1,i-1)
514  k = k+1
515  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
516  lgs_a_value(k) = -inv_dxi_deta &
517  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
518  lgs_a_index(k) = nc
519 
520  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
521  k = k+1
522  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
523  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
524  lgs_a_index(k) = nc
525 
526  nc = 2*nn(j+1,i)-1 ! next nc (column counter), for vx_m(j+1,i)
527  k = k+1
528  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
529  lgs_a_value(k) = inv_dxi_deta &
530  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i))
531  lgs_a_index(k) = nc
532 
533  nc = 2*nn(j+1,i) ! largest nc (column counter), for vy_m(j+1,i)
534  k = k+1
535  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
536  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j+1,i)
537  lgs_a_index(k) = nc
538 
539  lgs_b_value(nr) = factor_rhs_1 &
540  * (h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) &
541  * dzs_deta(j,i)
542 
543  lgs_x_value(nr) = vy_m(j,i)
544 
545 ! floating conditions satisfied as ice sheet
546  else if ( &
547  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==0_i2b) &
548  .and.(h_mid>=(z_sl-zl_mid)*rhosw_rho_ratio)) &
549  .or. &
550  ( (maske(j,i)==0_i2b).and.(maske(j+1,i)==3_i2b) &
551  .and.(h_mid>=(z_sl-zl_mid)*rhosw_rho_ratio)) &
552  ) then
553 
554 #if ( !defined(SHS) || SHS==0 )
555  k = k+1
556  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
557  lgs_a_value(k) = 1.0_dp ! diagonal element only
558  lgs_a_index(k) = nr
559 
560  lgs_b_value(nr) = vy_m(j,i)
561  lgs_x_value(nr) = vy_m(j,i)
562 
563 #elif SHS==1
564  ! definition of basal conditions
565  tau_b=((year_sec/c_slide_ssa)**(1.0_dp/p_ssa)) &
566  *(rho*g*(h_c(j,i)+h_t(j,i))**(q_ssa/p_ssa)) &
567  /((vx_m(j,i)*vx_m(j,i)+vy_m(j,i)*vy_m(j,i)) &
568  **(0.50_dp*(1.0_dp-1.0_dp/p_ssa)))
569 
570  ! flag shelfy flow
571  flag_sf(j,i)=.true.
572  flag_sf(j-1,i)=.true.
573 
574  nc = 2*nn(j-1,i) ! smallest nc (column counter), for vy_m(j-1,i)
575  k = k+1
576  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
577  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j,i)
578  lgs_a_index(k) = nc
579 
580  nc = 2*nn(j,i-1)-1 ! next nc (column counter), for vx_m(j,i-1)
581  k = k+1
582  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
583  lgs_a_value(k) = inv_dxi_deta &
584  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i-1))
585  lgs_a_index(k) = nc
586 
587  nc = 2*nn(j,i-1) ! next nc (column counter), for vy_m(j,i-1)
588  k = k+1
589  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
590  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
591  lgs_a_index(k) = nc
592 
593  nc = 2*nn(j,i)-1 ! next nc (column counter), for vx_m(j,i)
594  k = k+1
595  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
596  lgs_a_value(k) = -inv_dxi_deta &
597  *(2.0_dp*vis_int_g(j,i)+vis_int_sgxy(j,i))
598  lgs_a_index(k) = nc
599 
600  nc = 2*nn(j,i) ! next nc (column counter), for vy_m(j,i)
601  if (nc /= nr) & ! (diagonal element)
602  stop ' calc_vxy_ssa_matrix: Check for diagonal element failed!'
603  k = k+1
604  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
605  lgs_a_value(k) = -4.0_dp*inv_deta2 &
606  *(vis_int_g(j+1,i)+vis_int_g(j,i)) &
607  -inv_dxi2 &
608  *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1))-tau_b
609  lgs_a_index(k) = nc
610 
611  nc = 2*nn(j+1,i-1)-1 ! next nc (column counter), for vx_m(j+1,i-1)
612  k = k+1
613  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
614  lgs_a_value(k) = -inv_dxi_deta &
615  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
616  lgs_a_index(k) = nc
617 
618  nc = 2*nn(j,i+1) ! next nc (column counter), for vy_m(j,i+1)
619  k = k+1
620  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
621  lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
622  lgs_a_index(k) = nc
623 
624  nc = 2*nn(j+1,i)-1 ! next nc (column counter), for vx_m(j+1,i)
625  k = k+1
626  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
627  lgs_a_value(k) = inv_dxi_deta &
628  *(2.0_dp*vis_int_g(j+1,i)+vis_int_sgxy(j,i))
629  lgs_a_index(k) = nc
630 
631  nc = 2*nn(j+1,i) ! largest nc (column counter), for vy_m(j+1,i)
632  k = k+1
633  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
634  lgs_a_value(k) = 4.0_dp*inv_deta2*vis_int_g(j+1,i)
635  lgs_a_index(k) = nc
636 
637  lgs_b_value(nr) = factor_rhs_1 &
638  * (h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) &
639  * dzs_deta(j,i)
640 
641  lgs_x_value(nr) = vy_m(j,i)
642 
643 #endif
644 
645  else if ( &
646  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==1_i2b) ) &
647  .or. &
648  ( (maske(j,i)==1_i2b).and.(maske(j+1,i)==3_i2b) ) &
649  ) then
650  ! one neighbour is floating ice and the other is ice-free land;
651  ! velocity assumed to be zero
652 
653  k = k+1
654  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
655  lgs_a_value(k) = 1.0_dp ! diagonal element only
656  lgs_a_index(k) = nr
657 
658  lgs_b_value(nr) = 0.0_dp
659 
660  lgs_x_value(nr) = 0.0_dp
661 
662  else if ( &
663  ( (maske(j,i)==3_i2b).and.(maske(j+1,i)==2_i2b) ) &
664  .or. &
665  ( (maske(j,i)==2_i2b).and.(maske(j+1,i)==3_i2b) ) &
666  ) then
667  ! one neighbour is floating ice and the other is ocean
668 
669  if (maske(j,i)==3_i2b) then
670  j1 = j ! floating ice marker
671  else ! maske(j+1,i)==3_i2b
672  j1 = j+1 ! floating ice marker
673  end if
674 
675  if ( &
676  ( (maske(j1-1,i)==3_i2b).or.(maske(j1+1,i)==3_i2b) ) &
677  .or. &
678  ( (maske(j1-1,i)==0_i2b).or.(maske(j1+1,i)==0_i2b) ) &
679  .or. &
680  ( (maske(j1-1,i)==1_i2b).or.(maske(j1+1,i)==1_i2b) ) &
681  ) then
682  ! discretisation of the y-component of the BC at the calving front
683 
684  nc = 2*nn(j1-1,i) !smallest nc (column counter), for vy_m(j1-1,i)
685  k = k+1
686  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
687  lgs_a_value(k) = -4.0_dp*inv_deta*vis_int_g(j1,i)
688  lgs_a_index(k) = nc
689 
690  nc = 2*nn(j1,i-1)-1 ! next nc (column counter), for vx_m(j1,i-1)
691  k = k+1
692  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
693  lgs_a_value(k) = -2.0_dp*inv_dxi*vis_int_g(j1,i)
694  lgs_a_index(k) = nc
695 
696  nc = 2*nn(j1,i)-1 ! next nc (column counter), for vx_m(j1,i)
697  k = k+1
698  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
699  lgs_a_value(k) = 2.0_dp*inv_dxi*vis_int_g(j1,i)
700  lgs_a_index(k) = nc
701 
702  nc = 2*nn(j1,i) ! largest nc (column counter), for vy_m(j1,i)
703  k = k+1
704  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
705  lgs_a_value(k) = 4.0_dp*inv_deta*vis_int_g(j1,i)
706  lgs_a_index(k) = nc
707 
708  lgs_b_value(nr) = factor_rhs_2 &
709  *(h_c(j1,i)+h_t(j1,i))*(h_c(j1,i)+h_t(j1,i))
710 
711  lgs_x_value(nr) = vy_m(j,i)
712 
713  else ! (maske(j1-1,i)==2_i2b).and.(maske(j1+1,i)==2_i2b);
714  ! velocity assumed to be zero
715 
716  k = k+1
717  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
718  lgs_a_value(k) = 1.0_dp ! diagonal element only
719  lgs_a_index(k) = nr
720 
721  lgs_b_value(nr) = 0.0_dp
722 
723  lgs_x_value(nr) = 0.0_dp
724 
725  end if
726 
727  else if ( (maske(j,i)==0_i2b).or.(maske(j+1,i)==0_i2b) ) then
728  ! neither neighbour is floating ice, but at least one neighbour is
729  ! grounded ice; velocity taken from the solution for grounded ice
730 
731  k = k+1
732  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
733  lgs_a_value(k) = 1.0_dp ! diagonal element only
734  lgs_a_index(k) = nr
735 
736  lgs_b_value(nr) = vy_m(j,i)
737 
738  lgs_x_value(nr) = vy_m(j,i)
739 
740  else ! neither neighbour is floating or grounded ice,
741  ! velocity assumed to be zero
742 
743  k = k+1
744  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
745  lgs_a_value(k) = 1.0_dp ! diagonal element only
746  lgs_a_index(k) = nr
747 
748  lgs_b_value(nr) = 0.0_dp
749  lgs_x_value(nr) = 0.0_dp
750 
751  end if
752 
753  else ! boundary condition, velocity assumed to be zero
754 
755  k = k+1
756  if (k > n_sprs) stop ' calc_vxy_ssa_matrix: n_sprs too small!'
757  lgs_a_value(k) = 1.0_dp ! diagonal element only
758  lgs_a_index(k) = nr
759 
760  lgs_b_value(nr) = 0.0_dp
761  lgs_x_value(nr) = 0.0_dp
762 
763  end if
764 
765  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
766 
767 end do
768 
769 !-------- Settings for Lis --------
770 
771 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
772 call lis_vector_create(lis_comm_world, lgs_b, ierr)
773 call lis_vector_create(lis_comm_world, lgs_x, ierr)
774 
775 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
776 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
777 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
778 
779 do nr=1, nmax
780 
781  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
782  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
783  lgs_a_value(nc), lgs_a, ierr)
784  end do
785 
786  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
787  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
788 
789 end do
790 
791 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
792 call lis_matrix_assemble(lgs_a, ierr)
793 
794 !-------- Solution of the system of linear equations with Lis --------
795 
796 call lis_solver_create(solver, ierr)
797 
798 ch_solver_set_option = '-i bicgstabl -p ilu '// &
799  '-maxiter 1000 -tol 1.0e-11 -initx_zeros false'
800 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
801 
802 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
803 
804 ! call lis_solver_get_iters(solver, iter, ierr)
805 ! write(6,'(11x,a,i5)') 'calc_vxy_ssa_matrix: iter = ', iter
806 
807 lgs_x_value = 0.0_dp
808 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
809 call lis_matrix_destroy(lgs_a, ierr)
810 call lis_vector_destroy(lgs_b, ierr)
811 call lis_vector_destroy(lgs_x, ierr)
812 call lis_solver_destroy(solver, ierr)
813 
814 do n=1, nmax-1, 2
815 
816  i = ii((n+1)/2)
817  j = jj((n+1)/2)
818 
819  nr = n
820  vx_m(j,i) = lgs_x_value(nr)
821 
822  nr = n+1
823  vy_m(j,i) = lgs_x_value(nr)
824 
825 end do
826 
827 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
828 deallocate(lgs_b_value, lgs_x_value)
829 
830 end subroutine calc_vxy_ssa_matrix
831 !
Generated on Thu Sep 12 2013 09:45:13 for SICOPOLIS V3.1 by   doxygen 1.8.2