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