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