SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_top.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t o p
4 !
5 !> @file
6 !!
7 !! Computation of the topography (including gradients).
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve, Tatsuru Sato
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of the topography (including gradients).
34 !<------------------------------------------------------------------------------
35 subroutine calc_top(time, time_init, &
36  z_sl, dzsl_dtau, z_mar, &
37  dtime, dxi, deta, dzeta_t, &
38  mean_accum, &
39  ii, jj, nn, itercount, iter_wss)
40 
41 use sico_types
42 use sico_variables
43 use sico_sle_solvers, only : sor_sprs
44 
45 implicit none
46 
47 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
48 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
49 integer(i4b), intent(in) :: itercount, iter_wss
50 real(dp), intent(in) :: time, time_init
51 real(dp), intent(in) :: z_sl, dzsl_dtau, z_mar
52 real(dp), intent(in) :: dtime, dxi, deta, dzeta_t
53 real(dp), intent(in) :: mean_accum
54 
55 integer(i4b) :: i, j, kc, kt
56 integer(i4b) :: k, nnz
57 real(dp) :: eps_sor
58 real(dp) :: czs2(0:jmax,0:imax), czs3(0:jmax,0:imax)
59 real(dp) :: year_sec_inv
60 real(dp) :: tldt_inv
61 real(dp) :: dtime_inv, dxi_inv, deta_inv
62 real(dp) :: dt_dxi, dt_deta
63 real(dp) :: azs2, azs3
64 real(dp), dimension(0:JMAX,0:IMAX) :: h_ice, h_sea
65 real(dp) :: rho_rhoa_ratio, rhosw_rhoa_ratio
66 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio, rho_sw_rho_sw_minus_rho_ratio
67 real(dp) :: upvx_p, upvy_p, upvx_d, upvy_d, h_balance
68 real(dp) :: calv_uw_coeff, r1_calv_uw, r2_calv_uw
69 real(dp) :: target_topo_time_lag, target_topo_weight_lin
70 logical :: flag_calving_event
71 
72 real(dp), parameter :: h_isol_max = 1.0e+03_dp
73  ! Thickness limit of isolated glaciated points
74 
75 #if CALCZS==3
76 
77 integer(i4b) :: ierr
78 integer(i4b) :: iter
79 integer(i4b) :: nc, nr
80 integer(i4b), parameter :: nmax = (imax+1)*(jmax+1)
81 integer(i4b), parameter :: n_sprs = 10*(imax+1)*(jmax+1)
82 integer(i4b), allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
83 integer(i4b), allocatable, dimension(:) :: lgs_a_diag_index
84 integer(i4b), allocatable, dimension(:) :: lgs_a_index_trim
85 real(dp), allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
86 real(dp), allocatable, dimension(:) :: lgs_a_value_trim
87 
88 #elif CALCZS==4
89 
90 lis_integer :: ierr
91 lis_integer :: iter
92 lis_integer :: nc, nr
93 lis_integer, parameter :: nmax = (imax+1)*(jmax+1)
94 lis_integer, parameter :: n_sprs = 10*(imax+1)*(jmax+1)
95 lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
96 lis_integer, allocatable, dimension(:) :: lgs_a_diag_index
97 lis_matrix :: lgs_a
98 lis_vector :: lgs_b, lgs_x
99 lis_scalar, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
100 lis_solver :: solver
101 character(len=256) :: ch_solver_set_option
102 
103 #endif
104 
105 !-------- Term abbreviations --------
106 
107 year_sec_inv = 1.0_dp/year_sec
108 
109 tldt_inv = 1.0_dp/(time_lag+dtime)
110 
111 rho_rhoa_ratio = rho/rho_a
112 rhosw_rhoa_ratio = rho_sw/rho_a
113 
114 rho_rhosw_ratio = rho/rho_sw
115 rhosw_rho_ratio = rho_sw/rho
116 rho_sw_rho_sw_minus_rho_ratio = rho_sw/(rho_sw-rho)
117 
118 dtime_inv = 1.0_dp/dtime
119 dxi_inv = 1.0_dp/dxi
120 deta_inv = 1.0_dp/deta
121 
122 dt_dxi = dtime/dxi
123 dt_deta = dtime/deta
124 
125 azs2 = dtime/(dxi*dxi)
126 azs3 = dtime/(deta*deta)
127 
128 !-------- Computation of zl_neu and its time derivative --------
129 
130 #if REBOUND==0
131 
132 do i=0, imax
133 do j=0, jmax
134  zl_neu(j,i) = zl(j,i)
135  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
136 end do
137 end do
138 
139 #elif REBOUND==1
140 
141 do i=0, imax
142 do j=0, jmax
143 
144  if (maske(j,i) <= 1_i2b) then
145 
146  zl_neu(j,i) = tldt_inv*( time_lag*zl(j,i) &
147  + dtime*(zl0(j,i) &
148  -frac_llra*rho_rhoa_ratio*(h_c(j,i)+h_t(j,i))) )
149 
150  else ! (maske(j,i) >= 2_i2b)
151 
152  zl_neu(j,i) = tldt_inv*( time_lag*zl(j,i) &
153  + dtime*(zl0(j,i) &
154  -frac_llra*rhosw_rhoa_ratio*z_sl) )
155 
156  ! Water load relative to the present sea-level stand (0 m)
157  ! -> can be positive or negative
158 
159  end if
160 
161  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
162 
163 end do
164 end do
165 
166 #elif REBOUND==2
167 
168 if ((mod(itercount, iter_wss)==0).or.(itercount==1)) then
169  write(6,*) ' -------- Computation of wss'
170  call calc_elra(z_sl, dxi, deta) ! Update of the steady-state displacement
171  ! of the lithosphere (wss)
172 end if
173 
174 do i=0, imax
175 do j=0, jmax
176  zl_neu(j,i) = tldt_inv * ( time_lag*zl(j,i) + dtime*(zl0(j,i)-wss(j,i)) )
177  dzl_dtau(j,i) = (zl_neu(j,i)-zl(j,i))*dtime_inv
178 end do
179 end do
180 
181 #endif
182 
183 ! ------ Adjustment due to prescribed target topography
184 
185 #if ZS_EVOL==2
186 
187 if (time >= time_target_topo_final) then
188 
189  zl_neu = zl_target
190  dzl_dtau = 0.0_dp ! previously this was (zl_neu-zl)*dtime_inv
191 
192 else if (time >= time_target_topo_init) then
193 
194  target_topo_time_lag = (time_target_topo_final-time)/3.0_dp
195 
196  zl_neu = ( target_topo_time_lag*zl_neu + dtime*zl_target ) &
197  / ( target_topo_time_lag + dtime )
198 
199  dzl_dtau = (zl_neu-zl)*dtime_inv
200 
201 end if
202 
203 #endif
204 
205 !-------- Computation of zb_neu and its time derivative --------
206 
207 do i=0, imax
208 do j=0, jmax
209 
210  if (maske(j,i) <= 1_i2b) then
211 
212  zb_neu(j,i) = zl_neu(j,i)
213  dzb_dtau(j,i) = dzl_dtau(j,i)
214 
215  else ! (maske(j,i) >= 2_i2b)
216 
217 #if ( MARGIN==1 || MARGIN==2 )
218  zb_neu(j,i) = zl_neu(j,i)
219  dzb_dtau(j,i) = dzl_dtau(j,i)
220 #elif ( MARGIN==3 )
221  zb_neu(j,i) = zb(j,i) ! this is only preliminary
222  dzb_dtau(j,i) = 0.0_dp ! and will be corrected later
223 #endif
224 
225  end if
226 
227 end do
228 end do
229 
230 !-------- Computation of zs_neu --------
231 
232 ! ------ Explicit scheme (without or with ice shelves)
233 
234 #if CALCZS==0
235 
236 ! ---- Abbreviations
237 
238 do i=0, imax-1
239 do j=0, jmax
240  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
241  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
242 end do
243 end do
244 
245 do i=0, imax
246 do j=0, jmax-1
247  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
248  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
249 end do
250 end do
251 
252 #if ( MARGIN==1 || MARGIN==2 )
253 
254 do i=1, imax-1
255 do j=1, jmax-1
256 
257 #if MARGIN==1
258 
259  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
260 
261 #endif
262 
263  zs_neu(j,i) = zs(j,i) &
264  + dtime*(as_perp(j,i)+dzb_dtau(j,i)-q_b_tot(j,i)) &
265  + ( czs2(j,i) *(zs(j,i+1)-zs(j,i) ) &
266  -czs2(j,i-1)*(zs(j,i) -zs(j,i-1)) &
267  +czs3(j,i) *(zs(j+1,i)-zs(j,i) ) &
268  -czs3(j-1,i)*(zs(j,i) -zs(j-1,i)) ) &
269  *insq_g11_g(j,i)*insq_g22_g(j,i)
270 
271 #if MARGIN==1
272 
273  else ! maske(j,i) == 2_i2b, sea
274  zs_neu(j,i) = zb_neu(j,i)
275  end if
276 
277 #endif
278 
279 end do
280 end do
281 
282 #elif MARGIN==3
283 
284 do i=1, imax-1
285 do j=1, jmax-1
286 
287  if (.not.flag_grounding_line(j,i).and.maske(j,i)<2_i2b) then
288  ! inland ice/no_ice conditions
289 
290  h_ice(j,i) = (h_c(j,i)+h_t(j,i)) &
291  + dtime*(as_perp(j,i)-q_b_tot(j,i)) &
292  + ( czs2(j,i) *(zs(j,i+1)-zs(j,i) ) &
293  -czs2(j,i-1)*(zs(j,i) -zs(j,i-1)) &
294  +czs3(j,i) *(zs(j+1,i)-zs(j,i) ) &
295  -czs3(j-1,i)*(zs(j,i) -zs(j-1,i)) ) &
296  *insq_g11_g(j,i)*insq_g22_g(j,i)
297 
298  else ! ice shelves, ocean and on grounding lines
299 
300  upvx_p = 0.25_dp*( (vx_m(j,i)+vx_m(j,i-1)) &
301  -dabs(vx_m(j,i)+vx_m(j,i-1)) )*dt_dxi
302  upvx_d = 0.25_dp*( (vx_m(j,i)+vx_m(j,i-1)) &
303  +dabs(vx_m(j,i)+vx_m(j,i-1)) )*dt_dxi
304  upvy_p = 0.25_dp*( (vy_m(j,i)+vy_m(j-1,i)) &
305  -dabs(vy_m(j,i)+vy_m(j-1,i)) )*dt_deta
306  upvy_d = 0.25_dp*( (vy_m(j,i)+vy_m(j-1,i)) &
307  +dabs(vy_m(j,i)+vy_m(j-1,i)) )*dt_deta
308 
309  h_ice(j,i) = h_c(j,i)+h_t(j,i) &
310  + dtime*(as_perp(j,i)-q_b_tot(j,i) ) &
311  - ( upvx_p &
312  *((h_c(j,i+1)+h_t(j,i+1))-(h_c(j,i) +h_t(j,i) )) &
313  +upvx_d &
314  *((h_c(j,i) +h_t(j,i)) -(h_c(j,i-1)+h_c(j,i-1))) &
315  +upvy_p &
316  *((h_c(j+1,i)+h_t(j+1,i))-(h_c(j,i) +h_t(j,i) )) &
317  +upvy_d &
318  *((h_c(j,i) +h_t(j,i)) -(h_c(j-1,i)-h_t(j-1,i))) &
319  +dt_dxi &
320  *(vx_m(j,i)-vx_m(j,i-1))*(h_c(j,i)+h_t(j,i)) &
321  +dt_deta &
322  *(vy_m(j,i)-vy_m(j-1,i))*(h_c(j,i)+h_t(j,i)) )
323 
324  end if
325 
326 end do
327 end do
328 
329 #endif
330 
331 do i=0, imax
332 #if ( MARGIN==1 || MARGIN==2 )
333  zs_neu(0,i) = zb_neu(0,i)
334  zs_neu(jmax,i) = zb_neu(jmax,i)
335 #elif MARGIN==3
336  zs_neu(0,i) = z_sl
337  zs_neu(jmax,i) = z_sl
338 #endif
339 end do
340 
341 do j=0, jmax
342 #if ( MARGIN==1 || MARGIN==2 )
343  zs_neu(j,0) = zb_neu(j,0)
344  zs_neu(j,imax) = zb_neu(j,imax)
345 #elif MARGIN==3
346  zs_neu(j,0) = z_sl
347  zs_neu(j,imax) = z_sl
348 #endif
349 end do
350 
351 ! ------ ADI scheme / over-implicit ADI scheme (DELETED)
352 
353 #elif ( CALCZS==1 || CALCZS==2 )
354 
355 stop ' calc_top: ADI and over-implicit ADI schemes not available any more!'
356 
357 ! ------ Over-implicit scheme (without ice shelves)
358 
359 #elif ( ( CALCZS==3 || CALCZS==4 ) && ( MARGIN==1 || MARGIN==2 ) )
360 
361 ! ---- Abbreviations
362 
363 do i=0, imax-1
364 do j=0, jmax
365  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
366  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
367 end do
368 end do
369 
370 do i=0, imax
371 do j=0, jmax-1
372  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
373  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
374 end do
375 end do
376 
377 ! ---- Assembly of the system of linear equations
378 ! (matrix storage: compressed sparse row CSR)
379 
380 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
381 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
382 
383 lgs_a_value = 0.0_dp
384 lgs_a_index = 0
385 lgs_a_ptr = 0
386 
387 lgs_b_value = 0.0_dp
388 lgs_x_value = 0.0_dp
389 
390 lgs_a_ptr(1) = 1
391 
392 k = 0
393 
394 do nr=1, nmax ! loop over rows
395 
396  i = ii(nr)
397  j = jj(nr)
398 
399  if ( &
400  (i /= 0).and.(i /= imax).and.(j /= 0).and.(j /= jmax) & ! inner point
401 #if MARGIN==1
402  .and.(maske(j,i) <= 1_i2b) & ! no sea
403 #endif
404  ) then
405 
406  nc = nn(j-1,i) ! smallest nc (column counter)
407  k = k+1
408  if (k > n_sprs) stop ' calc_top: n_sprs too small!'
409  lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
410  *insq_g11_g(j,i)*insq_g22_g(j,i)
411  lgs_a_index(k) = nc
412 
413  nc = nn(j,i-1) ! next nc (column counter)
414  k = k+1
415  if (k > n_sprs) stop ' calc_top: n_sprs too small!'
416  lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
417  *insq_g11_g(j,i)*insq_g22_g(j,i)
418  lgs_a_index(k) = nc
419 
420  nc = nn(j,i) ! next nc (column counter)
421  if (nc /= nr) & ! (diagonal element)
422  stop ' calc_top: Check for diagonal element failed!'
423  k = k+1
424  if (k > n_sprs) stop ' calc_top: n_sprs too small!'
425  lgs_a_value(k) = 1.0_dp &
426  + (czs2(j,i)+czs2(j,i-1)+czs3(j,i)+czs3(j-1,i)) &
427  *ovi_weight &
428  *insq_g11_g(j,i)*insq_g22_g(j,i)
429  lgs_a_diag_index(nr) = k
430  lgs_a_index(k) = nc
431 
432  nc = nn(j,i+1) ! next nc (column counter)
433  k = k+1
434  if (k > n_sprs) stop ' calc_top: n_sprs too small!'
435  lgs_a_value(k) = -czs2(j,i)*ovi_weight &
436  *insq_g11_g(j,i)*insq_g22_g(j,i)
437  lgs_a_index(k) = nc
438 
439  nc = nn(j+1,i) ! largest nc (column counter)
440  k = k+1
441  if (k > n_sprs) stop ' calc_top: n_sprs too small!'
442  lgs_a_value(k) = -czs3(j,i)*ovi_weight &
443  *insq_g11_g(j,i)*insq_g22_g(j,i)
444  lgs_a_index(k) = nc
445 
446  lgs_b_value(nr) = zs(j,i) &
447  + dtime*(as_perp(j,i)+dzb_dtau(j,i)-q_b_tot(j,i)) &
448  + ( czs2(j,i)*(zs(j,i+1)-zs(j,i)) &
449  -czs2(j,i-1)*(zs(j,i)-zs(j,i-1)) &
450  +czs3(j,i)*(zs(j+1,i)-zs(j,i)) &
451  -czs3(j-1,i)*(zs(j,i)-zs(j-1,i)) ) &
452  *(1.0_dp-ovi_weight) &
453  *insq_g11_g(j,i)*insq_g22_g(j,i)
454 
455  else ! boundary condition
456 
457  k = k+1
458  if (k > n_sprs) stop ' calc_top: n_sprs too small!'
459  lgs_a_value(k) = 1.0_dp ! diagonal element only
460  lgs_a_diag_index(nr) = k
461  lgs_a_index(k) = nr
462 
463  lgs_b_value(nr) = zb_neu(j,i)
464 
465  end if
466 
467  lgs_x_value(nr) = zs(j,i) ! old surface topography,
468  ! initial guess for solution vector
469 
470  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
471 
472 end do
473 
474 nnz = k ! number of non-zero elements of the matrix
475 
476 #if CALCZS==3
477 
478 ! ---- Solution of the system of linear equations
479 
480 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
481 
482 do k=1, nnz ! relocate matrix to trimmed arrays
483  lgs_a_value_trim(k) = lgs_a_value(k)
484  lgs_a_index_trim(k) = lgs_a_index(k)
485 end do
486 
487 deallocate(lgs_a_value, lgs_a_index)
488 
489 eps_sor = 1.0e-05_dp*mean_accum*dtime ! convergence parameter
490 
491 call sor_sprs(lgs_a_value_trim, &
492  lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
493  lgs_b_value, &
494  nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
495 
496 do nr=1, nmax
497  i = ii(nr)
498  j = jj(nr)
499  zs_neu(j,i) = lgs_x_value(nr)
500 end do
501 
502 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
503 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
504 
505 #elif CALCZS==4
506 
507 ! ---- Settings for Lis
508 
509 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
510 call lis_vector_create(lis_comm_world, lgs_b, ierr)
511 call lis_vector_create(lis_comm_world, lgs_x, ierr)
512 
513 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
514 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
515 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
516 
517 do nr=1, nmax
518 
519  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
520  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
521  lgs_a_value(nc), lgs_a, ierr)
522  end do
523 
524  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
525  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
526 
527 end do
528 
529 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
530 call lis_matrix_assemble(lgs_a, ierr)
531 
532 ! ---- Solution of the system of linear equations with Lis
533 
534 call lis_solver_create(solver, ierr)
535 
536 ch_solver_set_option = '-i bicg -p ilu '// &
537  '-maxiter 1000 -tol 1.0e-12 -initx_zeros false'
538 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
539 
540 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
541 
542 call lis_solver_get_iters(solver, iter, ierr)
543 write(6,'(11x,a,i5)') 'calc_top: iter = ', iter
544 
545 lgs_x_value = 0.0_dp
546 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
547 call lis_matrix_destroy(lgs_a, ierr)
548 call lis_vector_destroy(lgs_b, ierr)
549 call lis_vector_destroy(lgs_x, ierr)
550 call lis_solver_destroy(solver, ierr)
551 
552 do nr=1, nmax
553  i = ii(nr)
554  j = jj(nr)
555  zs_neu(j,i) = lgs_x_value(nr)
556 end do
557 
558 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
559 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
560 
561 #endif
562 
563 ! ------ Over-implicit ice thickness scheme (with ice shelves)
564 
565 #elif ( CALCZS==3 && MARGIN==3 )
566 
567 do i=0, imax
568 do j=0, jmax
569  h_ice(j,i) =(h_c(j,i)+h_t(j,i))
570 end do
571 end do
572 
573 do i=0, imax-1
574 do j=0, jmax
575  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
576  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
577 end do
578 end do
579 
580 do i=0, imax
581 do j=0, jmax-1
582  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
583  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
584 end do
585 end do
586 
587 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
588 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
589 
590 ! init of a*x=b, a is sparse matrix of A
591 lgs_a_value = 1.0_dp ; lgs_b_value=0.0_dp ; lgs_x_value=0.0_dp
592 ! total term number of a and ija
593 k=0
594 
595 ! main loop for making matrix
596 do nr=1, nmax
597 
598 ! set numbering
599  i = ii(nr) ; j = jj(nr)
600 
601 ! non_edge_points
602  if ( (i /= 0).and.(i /= imax).and.(j /= 0).and.(j /= jmax) ) then
603 
604 !#if SHS==0
605  if (.not.flag_grounding_line(j,i).and.maske(j,i)<2_i2b) then
606 !#else
607  ! if ((.not.flag_grounding_line(j,i).or..not.flag_sf(j,i)).and.maske(j,i)<2_i2b) then
608 !#endif
609 
610  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=h_ice(j,i)
611 ! right hand side
612  lgs_b_value(nr) = h_ice(j,i) &
613  + dtime*(as_perp(j,i)-q_b_tot(j,i)) &
614  + ( czs2(j,i)*(h_ice(j,i+1)-h_ice(j,i)) &
615  -czs2(j,i-1)*(h_ice(j,i)-h_ice(j,i-1)) &
616  +czs3(j,i)*(h_ice(j+1,i)-h_ice(j,i)) &
617  -czs3(j-1,i)*(h_ice(j,i)-h_ice(j-1,i)) ) &
618  *(1.0_dp-ovi_weight) &
619  *insq_g11_g(j,i)*insq_g22_g(j,i) &
620  + ( czs2(j,i)*(zb_neu(j,i+1)-zb_neu(j,i)) &
621  -czs2(j,i-1)*(zb_neu(j,i)-zb_neu(j,i-1)) &
622  +czs3(j,i)*(zb_neu(j+1,i)-zb_neu(j,i)) &
623  -czs3(j-1,i)*(zb_neu(j,i)-zb_neu(j-1,i)) ) &
624  *insq_g11_g(j,i)*insq_g22_g(j,i)
625 
626 ! H(j-1,i)
627  k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc
628  lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
629  *insq_g11_g(j,i)*insq_g22_g(j,i)
630 ! zs(j,i-1)
631  k=k+1 ; nc = nn(j,i-1) ; lgs_a_index(k)=nc
632  lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
633  *insq_g11_g(j,i)*insq_g22_g(j,i)
634 ! zs(j,i) diagonal
635  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
636  lgs_a_value(k) = 1.0_dp + (czs2(j,i)+czs2(j,i-1) &
637  +czs3(j,i)+czs3(j-1,i)) &
638  *ovi_weight &
639  *insq_g11_g(j,i)*insq_g22_g(j,i)
640 ! zs(j,i+1)
641  k=k+1 ; nc = nn(j,i+1) ; lgs_a_index(k)=nc
642  lgs_a_value(k) = -czs2(j,i)*ovi_weight &
643  *insq_g11_g(j,i)*insq_g22_g(j,i)
644 ! zs(j+1,i)
645  k=k+1 ; nc = nn(j+1,i) ; lgs_a_index(k)=nc
646  lgs_a_value(k) = -czs3(j,i)*ovi_weight &
647  *insq_g11_g(j,i)*insq_g22_g(j,i)
648 
649 #if ( !defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1 )
650 
651  else
652 
653 #elif ZS_EVOL_SSA==0
654 
655  else if (flag_grounding_line(j,i)) then
656 
657 #else
658  stop ' calc_top: ZS_EVOL_SSA must be either 1 or 0!'
659 #endif
660 
661 ! ab
662  upvx_p = 0.25_dp*( (vx_m(j,i)+vx_m(j,i-1)) &
663  -dabs(vx_m(j,i)+vx_m(j,i-1)) )*dt_dxi
664  upvx_d = 0.25_dp*( (vx_m(j,i)+vx_m(j,i-1)) &
665  +dabs(vx_m(j,i)+vx_m(j,i-1)) )*dt_dxi
666  upvy_p = 0.25_dp*( (vy_m(j,i)+vy_m(j-1,i)) &
667  -dabs(vy_m(j,i)+vy_m(j-1,i)) )*dt_deta
668  upvy_d = 0.25_dp*( (vy_m(j,i)+vy_m(j-1,i)) &
669  +dabs(vy_m(j,i)+vy_m(j-1,i)) )*dt_deta
670 
671  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=h_ice(j,i)
672  lgs_b_value(nr)= h_ice(j,i) &
673  +dtime*(as_perp(j,i)-q_b_tot(j,i)) &
674  -(1.0_dp-ovi_weight) &
675  *( upvx_p*(h_ice(j,i+1)-h_ice(j,i)) &
676  +upvx_d*(h_ice(j,i)-h_ice(j,i-1)) &
677  +upvy_p*(h_ice(j+1,i)-h_ice(j,i)) &
678  +upvy_d*(h_ice(j,i)-h_ice(j-1,i)) &
679  +dt_dxi*(vx_m(j,i)-vx_m(j,i-1))*h_ice(j,i) &
680  +dt_deta*(vy_m(j,i)-vy_m(j-1,i))*h_ice(j,i) )
681 
682 ! zs(j-1,i)
683  k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc
684  lgs_a_value(k) = -upvy_d *ovi_weight !&
685  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
686 ! zs(j,i-1)
687  k=k+1 ; nc = nn(j,i-1) ; lgs_a_index(k)=nc
688  lgs_a_value(k) = -upvx_d *ovi_weight !&
689  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
690 ! zs(j,i) diagonal
691  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
692  lgs_a_value(k) = 1.0_dp &
693  +ovi_weight &
694  *( -upvx_p+upvx_d-upvy_p+upvy_d &
695  +dt_dxi*(vx_m(j,i)-vx_m(j,i-1)) &
696  +dt_deta*(vy_m(j,i)-vy_m(j-1,i)) ) !&
697  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
698 ! zs(j,i+1)
699  k=k+1 ; nc = nn(j,i+1) ; lgs_a_index(k)=nc
700  lgs_a_value(k) = upvx_p * ovi_weight !&
701  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
702 ! zs(j+1,i)
703  k=k+1 ; nc = nn(j+1,i) ; lgs_a_index(k)=nc
704  lgs_a_value(k) = upvy_p * ovi_weight !&
705  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
706 
707 #if ( !defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1 )
708 
709  continue ! do nothing
710 
711 #elif ZS_EVOL_SSA==0
712 
713  else if (maske(j,i)==3_i2b) then
714 
715  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=h_ice(j,i)
716  lgs_b_value(nr) = h_ice(j,i)
717 
718  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
719  lgs_a_value(k) = 1.0_dp
720 
721  else if (maske(j,i)==2_i2b) then
722 
723  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=0.0_dp
724  lgs_b_value(nr) = 0.0_dp ! zb_neu(j,i)
725 
726  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
727  lgs_a_value(k) = 1.0_dp
728 
729 #else
730  stop ' calc_top: ZS_EVOL_SSA must be either 1 or 0!'
731 #endif
732 
733  end if
734 
735  else ! boundary condition, I=0 or IMAX, J=0 or JMAX
736 
737  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=0.0_dp
738  lgs_b_value(nr) = 0.0_dp !zb_neu(j,i)
739 
740  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
741  lgs_a_value(k) = 1.0_dp
742 
743  end if
744 
745 end do
746 
747 ! set csr's last index
748 lgs_a_ptr(nmax+1)=k+1
749 
750 nnz = k ! number of non-zero elements of the matrix
751 
752 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
753 
754 do k=1, nnz ! relocate matrix to trimmed arrays
755  lgs_a_value_trim(k) = lgs_a_value(k)
756  lgs_a_index_trim(k) = lgs_a_index(k)
757 end do
758 
759 deallocate(lgs_a_value, lgs_a_index)
760 
761 eps_sor = 1.0e-05_dp*mean_accum*dtime ! convergence parameter
762 
763 call sor_sprs(lgs_a_value_trim, &
764  lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
765  lgs_b_value, &
766  nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
767 
768 do nr=1, nmax
769  i = ii(nr)
770  j = jj(nr)
771  h_ice(j,i) = lgs_x_value(nr)
772 end do
773 
774 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
775 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
776 
777 ! #elif CALCZS == 5
778 !
779 ! #if MARGIN==3
780 !
781 ! do i=1,IMAX-1
782 ! do j=1,JMAX-1
783 !
784 ! if (maske(j,i)/=2) then
785 ! H_ice(j,i)= dtime*(as_perp(j,i)-Q_b_tot(j,i))+(H_c(j,i)+ H_t(j,i))
786 ! else
787 ! H_ice(j,i) = 0.0_dp
788 ! end if
789 !
790 ! end do
791 ! end do
792 !
793 ! H_ice(:,0)=0.0_dp ; H_ice(:,IMAX)=0.0_dp
794 ! H_ice(0,:)=0.0_dp ; H_ice(JMAX,:)=0.0_dp
795 !
796 ! #endif
797 
798 #endif
799 
800 ! ------ Adjustment due to prescribed target topography
801 
802 #if ZS_EVOL==2
803 
804 #if ( MARGIN==1 || MARGIN==2 )
805 
806 if (time >= time_target_topo_final) then
807  zs_neu = zs_target
808 else if (time >= time_target_topo_init) then
809  zs_neu = ( target_topo_time_lag*zs_neu + dtime*zs_target ) &
810  / ( target_topo_time_lag + dtime )
811 end if
812 
813 #elif ( MARGIN==3 )
814 
815 if (time >= time_target_topo_final) then
816  h_ice = h_target
817 else if (time >= time_target_topo_init) then
818  h_ice = ( target_topo_time_lag*h_ice + dtime*h_target ) &
819  / ( target_topo_time_lag + dtime )
820 end if
821 
822 #endif
823 
824 #endif
825 
826 !-------- Update of the mask --------
827 
828 #if ( ( MARGIN==2 ) \
829  && ( marine_ice_formation==2 ) \
830  && ( marine_ice_calving==9 ) )
831 
832 calv_uw_coeff = calv_uw_coeff * year_sec_inv
833 r1_calv_uw = r1_calv_uw
834 r2_calv_uw = r2_calv_uw
835 
836 #if ZS_EVOL==2
837 if (time >= time_target_topo_final) then
838  calv_uw_coeff = 0.0_dp
839  ! adjustment due to prescribed target topography
840 else if (time >= time_target_topo_init) then
841  target_topo_weight_lin = (time-time_target_topo_init) &
842  /(time_target_topo_final-time_target_topo_init)
843  calv_uw_coeff = (1.0_dp-target_topo_weight_lin)*calv_uw_coeff
844  ! adjustment due to prescribed target topography
845 end if
846 #endif
847 
848 #endif
849 
850 #if ( ZS_EVOL>=1 && ( MARGIN==1 || MARGIN==2 ) )
851 
852 do i=0, imax
853 do j=0, jmax
854 
855  if (maske(j,i) <= 1_i2b) then ! land
856 
857  h_ice(j,i) = zs_neu(j,i)-zb_neu(j,i) ! ice thickness
858  h_sea(j,i) = z_sl -zl_neu(j,i) ! sea depth (only makes sense for
859  ! marine ice and "underwater ice",
860  ! otherwise meaningless)
861 
862  if (h_ice(j,i) >= h_min) then
863  maske(j,i) = 0_i2b ! grounded ice
864 #if ( ( MARGIN==2 ) \
865  && ( marine_ice_formation==2 ) \
866  && ( marine_ice_calving==9 ) )
867  if ( h_ice(j,i) < (rhosw_rho_ratio*h_sea(j,i)) ) then
868  ! "underwater ice"
869  h_ice(j,i) = h_ice(j,i) - ( calv_uw_coeff &
870  *h_ice(j,i)**r1_calv_uw &
871  *h_sea(j,i)**r2_calv_uw ) * dtime
872  ! calving of "underwater ice"
873  zs_neu(j,i) = zb_neu(j,i) + h_ice(j,i)
874  if (h_ice(j,i) < h_min) maske(j,i) = 2_i2b ! sea
875  end if
876 #endif
877  else
878  maske(j,i) = 1_i2b ! ice-free land
879  end if
880 
881 #if MARGIN==2
882 
883  else ! (maske(j,i) == 2_i2b); sea
884 
885  h_ice(j,i) = zs_neu(j,i)-zb_neu(j,i) ! ice thickness
886  h_sea(j,i) = z_sl -zl_neu(j,i) ! sea depth
887 
888  if (h_ice(j,i) >= h_min) then
889 
890  if ( h_ice(j,i) < (rhosw_rho_ratio*h_sea(j,i)) ) then
891 
892 #if MARINE_ICE_FORMATION==1
893  maske(j,i) = 2_i2b ! floating ice cut off -> sea
894 #elif MARINE_ICE_FORMATION==2
895  maske(j,i) = 0_i2b ! "underwater ice"
896 #if MARINE_ICE_CALVING==9
897  h_ice(j,i) = h_ice(j,i) - ( calv_uw_coeff &
898  *h_ice(j,i)**r1_calv_uw &
899  *h_sea(j,i)**r2_calv_uw ) * dtime
900  ! calving of "underwater ice"
901  zs_neu(j,i) = zb_neu(j,i) + h_ice(j,i)
902  if (h_ice(j,i) < h_min) maske(j,i) = 2_i2b ! sea
903 #endif
904 #endif
905  else
906  maske(j,i) = 0_i2b ! grounded ice
907  end if
908 
909 #if ( MARINE_ICE_CALVING==2 \
910  || marine_ice_calving==4 \
911  || marine_ice_calving==6 )
912  if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
913 #elif ( MARINE_ICE_CALVING==3 \
914  || marine_ice_calving==5 \
915  || marine_ice_calving==7 )
916  if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
917 #endif
918 
919  else
920  maske(j,i) = 2_i2b ! sea
921  end if
922 
923 #endif
924 
925  end if
926 
927 end do
928 end do
929 
930 #elif ( ZS_EVOL>=1 && MARGIN==3 )
931 
932 do i=1, imax-1
933 do j=1, jmax-1
934 
935  h_balance = (z_sl-zl_neu(j,i))*rhosw_rho_ratio
936 
937 ! grounding_line migration check
938 #if ( !defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1 )
939 
940  if ( ( maske(j,i)<2_i2b &
941  .and. (maske(j,i+1)>1_i2b.or.maske(j,i-1)>1_i2b &
942  .or.maske(j+1,i)>1_i2b.or.maske(j-1,i)>1_i2b) ) &
943  .or. &
944  ( maske(j,i)>=2_i2b &
945  .and. (maske(j,i+1)==0_i2b.or.maske(j,i-1)==0_i2b &
946  .or.maske(j+1,i)==0_i2b.or.maske(j-1,i)==0_i2b) ) ) then
947 
948  if (h_ice(j,i)>=h_min) then
949 
950  if (h_ice(j,i)<h_balance.and.zl_neu(j,i)<z_sl) then
951  maske(j,i) = 3_i2b
952  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_ice(j,i)
953  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
954  zs_neu(j,i) = zb_neu(j,i)+h_ice(j,i)
955  else
956  maske(j,i) = 0_i2b
957  zb_neu(j,i) = zl_neu(j,i)
958  dzb_dtau(j,i) = dzl_dtau(j,i)
959  zs_neu(j,i) = zb_neu(j,i)+h_ice(j,i)
960  end if
961 
962  else ! if (H_ice(j,i)<=H_MIN) then
963 
964  if (zl_neu(j,i)>z_sl) then
965  maske(j,i) = 1_i2b
966  zb_neu(j,i) = zl_neu(j,i)
967  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
968  zs_neu(j,i) = zb_neu(j,i)
969  else
970  maske(j,i) = 2_i2b
971  zb_neu(j,i) = z_sl
972  dzb_dtau(j,i) = 0.0_dp
973  zs_neu(j,i) = z_sl
974  end if
975 
976  end if
977 
978 #elif ZS_EVOL_SSA==0
979 
980  if (maske(j,i)==2_i2b) then
981  maske(j,i) = 2_i2b
982  zb_neu(j,i) = z_sl
983  dzb_dtau(j,i) = 0.0_dp
984  zs_neu(j,i) = z_sl
985  else if (maske(j,i)==3_i2b) then
986  maske(j,i) = 3_i2b
987  zb_neu(j,i) = zb(j,i)
988  dzb_dtau(j,i) = 0.0_dp
989  zs_neu(j,i) = zs(j,i)
990 
991 #else
992  stop ' calc_top: ZS_EVOL_SSA must be either 1 or 0!'
993 #endif
994 
995  else if (maske(j,i)<2_i2b) then
996 
997 ! can change maske 0 or 1
998  if (h_ice(j,i)>=h_min) then
999  maske(j,i) = 0_i2b
1000  zb_neu(j,i) = zl_neu(j,i)
1001  dzb_dtau(j,i) = dzl_dtau(j,i)
1002  zs_neu(j,i) = zb_neu(j,i)+h_ice(j,i)
1003  else
1004  maske(j,i) = 1_i2b
1005  zb_neu(j,i) = zl_neu(j,i)
1006  dzb_dtau(j,i) = dzl_dtau(j,i)
1007  zs_neu(j,i) = zl_neu(j,i)
1008  end if
1009 
1010 #if ( !defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1 )
1011 
1012  else ! if (maske(j,i)==2.or.maske(j,i)==3) then
1013 
1014  if (h_ice(j,i)>h_min) then
1015 
1016  if (h_ice(j,i)<h_balance.and.zl_neu(j,i)<z_sl) then
1017  maske(j,i) = 3_i2b
1018  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_ice(j,i)
1019  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
1020  zs_neu(j,i) = zb_neu(j,i)+h_ice(j,i)
1021  else
1022  maske(j,i) = 0_i2b
1023  zb_neu(j,i) = zl_neu(j,i)
1024  dzb_dtau(j,i) = dzl_dtau(j,i)
1025  zs_neu(j,i) = zb_neu(j,i)+h_ice(j,i)
1026  end if
1027 
1028  else ! if (H_ice(j,i)<=H_MIN) then
1029 
1030  if (zl_neu(j,i)>z_sl) then
1031  maske(j,i) = 1_i2b
1032  zb_neu(j,i) = zl_neu(j,i)
1033  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
1034  zs_neu(j,i) = zb_neu(j,i)
1035  else
1036  maske(j,i) = 2_i2b
1037  zb_neu(j,i) = z_sl
1038  dzb_dtau(j,i) = 0.0_dp
1039  zs_neu(j,i) = z_sl
1040  end if
1041 
1042  end if
1043 
1044 #endif
1045 
1046  end if
1047 
1048 end do
1049 end do
1050 
1051 #if ICE_SHELF_CALVING==2
1052 
1053 do
1054 
1055  flag_calving_front = .false.
1056  flag_calving_event = .false.
1057 
1058  do i=1, imax-1
1059  do j=1, jmax-1
1060 
1061  if ( (maske(j,i)==3_i2b) & ! floating ice
1062  .and. &
1063  ( (maske(j,i+1)==2_i2b) & ! with
1064  .or.(maske(j,i-1)==2_i2b) & ! one
1065  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1066  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1067  ) &
1068  flag_calving_front(j,i) = .true. ! preliminary detection
1069  ! of the calving front
1070 
1071  if ( flag_calving_front(j,i).and.(h_ice(j,i) < h_calv) ) then
1072  flag_calving_event = .true. ! calving event,
1073  maske(j,i) = 2_i2b ! floating ice point changes to sea point
1074  end if
1075 
1076  end do
1077  end do
1078 
1079  if (.not.flag_calving_event) exit
1080 
1081 end do
1082 
1083 #endif
1084 
1085 #endif
1086 
1087 ! ------ Adjustment due to prescribed target topography
1088 
1089 #if ZS_EVOL==2
1090 if (time >= time_target_topo_final) maske = maske_target
1091 #endif
1092 
1093 ! ------ Adjustment due to prescribed maximum ice extent
1094 
1095 #if ZS_EVOL==3
1096 
1097 do i=0, imax
1098 do j=0, jmax
1099 
1100  if (maske_maxextent(j,i)==0_i2b) then ! not allowed to glaciate
1101 
1102  if (maske(j,i)==0_i2b) then
1103  maske(j,i) = 1_i2b ! grounded ice -> ice-free land
1104  else if (maske(j,i)==3_i2b) then
1105  maske(j,i) = 2_i2b ! floating ice -> sea
1106  end if
1107 
1108  end if
1109 
1110 end do
1111 end do
1112 
1113 #endif
1114 
1115 !-------- Detection of the grounding line and the calving front --------
1116 
1117 flag_grounding_line = .false.
1118 flag_calving_front = .false.
1119 
1120 #if MARGIN==3
1121 
1122 do i=1, imax-1
1123 do j=1, jmax-1
1124 
1125  if ( (maske(j,i)==0_i2b) & ! grounded ice
1126  .and. &
1127  ( (maske(j,i+1)==3_i2b) & ! with
1128  .or.(maske(j,i-1)==3_i2b) & ! one
1129  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1130  .or.(maske(j-1,i)==3_i2b) ) & ! grounded ice point
1131  ) &
1132  flag_grounding_line(j,i) = .true.
1133 
1134  if ( (maske(j,i)==3_i2b) & ! floating ice
1135  .and. &
1136  ( (maske(j,i+1)==2_i2b) & ! with
1137  .or.(maske(j,i-1)==2_i2b) & ! one
1138  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1139  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1140  ) &
1141  flag_calving_front(j,i) = .true.
1142 
1143 end do
1144 end do
1145 
1146 #endif
1147 
1148 !-------- Correction of zs_neu and zb_neu
1149 ! for ice-free land, sea and floating ice --------
1150 
1151 do i=0, imax
1152 do j=0, jmax
1153 
1154  if (maske(j,i) == 1_i2b) then ! ice-free land
1155 
1156  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1157  ! from being below zb_neu(j,i)
1158 
1159  else if (maske(j,i) == 2_i2b) then ! sea
1160 
1161 #if ( MARGIN==1 || MARGIN==2 )
1162  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
1163  ! from being below zb_neu(j,i)
1164 #elif ( MARGIN==3 )
1165  zs_neu(j,i) = z_sl ! both zs_neu(j,i) and zb_neu(j,i)
1166  zb_neu(j,i) = z_sl ! set to the current sea level stand
1167 #endif
1168 
1169  else if (maske(j,i) == 3_i2b) then ! floating ice
1170 
1171 #if ( MARGIN==1 || MARGIN==2 )
1172 
1173  stop ' calc_top: maske(j,i)==3 only allowed for MARGIN==3!'
1174 
1175 #elif ( MARGIN==3 )
1176 
1177  h_ice(j,i) = zs_neu(j,i)-zb_neu(j,i) ! ice thickness
1178 
1179  zb_neu(j,i) = z_sl - rho_rhosw_ratio*h_ice(j,i) ! floating condition
1180  zs_neu(j,i) = zb_neu(j,i) + h_ice(j,i)
1181 
1182 #endif
1183 
1184  end if
1185 
1186 end do
1187 end do
1188 
1189 !-------- Computation of further quantities --------
1190 
1191 #if ZS_EVOL==0
1192 
1193 do i=0, imax
1194 do j=0, jmax
1195 
1196  if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
1197  ! grounded or floating ice
1198 
1199  as_perp(j,i) = as_perp(j,i) - (zs_neu(j,i)-zs(j,i))*dtime_inv
1200  ! correction of the accumulation-ablation function
1201  ! such that it is consistent with a steady surface
1202 
1203  else ! maske(j,i)==0_i2b or 3_i2b, ice-free land or sea
1204 
1205  as_perp(j,i) = 0.0_dp
1206 
1207  end if
1208 
1209  zs_neu(j,i) = zs(j,i) ! newly computed surface topography is discarded
1210 
1211 end do
1212 end do
1213 
1214 #endif
1215 
1216 do i=0, imax
1217 do j=0, jmax
1218 
1219  if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
1220  ! grounded or floating ice
1221 
1222 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
1223 
1224  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
1225  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
1226  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
1227  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
1228  ) then ! all nearest neighbours ice-free
1229  zs_neu(j,i) = min(zs_neu(j,i),(zb_neu(j,i)+h_isol_max))
1230  end if
1231 
1232 ! ------ Further stuff
1233 
1234  if (n_cts(j,i) == 1) then
1235  h_t_neu(j,i) = h_t(j,i)*(zs_neu(j,i)-zb_neu(j,i)) &
1236  /(zs(j,i)-zb(j,i))
1237  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
1238  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1239  else
1240  h_t_neu(j,i) = 0.0_dp
1241  zm_neu(j,i) = zb_neu(j,i)
1242  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
1243  end if
1244 
1245  else ! maske(j,i) == 1_i2b, 2_i2b
1246 
1247  zs_neu(j,i) = zb_neu(j,i)
1248  zm_neu(j,i) = zb_neu(j,i)
1249  h_c_neu(j,i) = 0.0_dp
1250  h_t_neu(j,i) = 0.0_dp
1251 
1252  end if
1253 
1254  dzs_dtau(j,i) = (zs_neu(j,i)-zs(j,i))*dtime_inv
1255  dzb_dtau(j,i) = (zb_neu(j,i)-zb(j,i))*dtime_inv
1256  dzm_dtau(j,i) = dh_t_dtau(j,i)+dzb_dtau(j,i)
1257  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
1258 
1259 #if ZS_EVOL==2
1260  if ( abs((time-time_target_topo_final)*year_sec_inv) < eps ) then
1261  dzs_dtau = 0.0_dp
1262  dzb_dtau = 0.0_dp
1263  dzm_dtau = 0.0_dp
1264  dh_c_dtau(j,i) = 0.0_dp
1265  dh_t_dtau(j,i) = 0.0_dp
1266  end if
1267 #endif
1268 
1269 end do
1270 end do
1271 
1272 !-------- New gradients --------
1273 
1274 #if TOPOGRAD==0
1275 call topograd_1(dxi, deta, 2)
1276 #elif TOPOGRAD==1
1277 call topograd_2(dxi, deta, 2)
1278 #endif
1279 
1280 !-------- Compute the volume flux across the CTS, am_perp --------
1281 
1282 #if CALCMOD==1
1283 
1284 do i=1, imax-1
1285 do j=1, jmax-1
1286 
1287  if ( ((maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b)) &
1288  .and.(n_cts(j,i) == 1_i2b)) then
1289 
1290  am_perp_st(j,i) = &
1291  0.5_dp*(vx_c(0,j,i)+vx_c(0,j,i-1))*dzm_dxi_g(j,i) &
1292  + 0.5_dp*(vy_c(0,j,i)+vy_c(0,j-1,i))*dzm_deta_g(j,i) &
1293  - 0.5_dp*(vz_c(0,j,i)+vz_t(ktmax-1,j,i))
1294  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
1295 
1296  else
1297  am_perp_st(j,i) = 0.0_dp
1298  am_perp(j,i) = 0.0_dp
1299  end if
1300 
1301 end do
1302 end do
1303 
1304 #endif
1305 
1306 end subroutine calc_top
1307 !
Generated on Thu Sep 12 2013 09:45:12 for SICOPOLIS V3.1 by   doxygen 1.8.2