SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_top_1.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t o p _ 1
4 !
5 !> @file
6 !!
7 !! Computation of the ice topography (including gradients) for
8 !! SIA-only dynamics of ice sheets without ice shelves.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2016 Ralf Greve
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of the ice topography (including gradients) for
35 !! SIA-only dynamics of ice sheets without ice shelves.
36 !<------------------------------------------------------------------------------
37 subroutine calc_top_1(time, z_sl, z_mar, dtime, dxi, deta, mean_accum, &
38  ii, jj, nn)
39 
40 use sico_types
41 use sico_variables
42 use sico_vars
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 real(dp), intent(in) :: time
50 real(dp), intent(in) :: z_sl, z_mar
51 real(dp), intent(in) :: dtime, dxi, deta
52 real(dp), intent(in) :: mean_accum
53 
54 integer(i4b) :: i, j, kc, kt
55 integer(i4b) :: k, nnz
56 real(dp) :: eps_sor
57 real(dp) :: czs2(0:jmax,0:imax), czs3(0:jmax,0:imax)
58 real(dp) :: year_sec_inv
59 real(dp) :: dtime_inv
60 real(dp) :: azs2, azs3
61 real(dp), dimension(0:JMAX,0:IMAX) :: h, h_neu, h_sea_neu
62 real(dp) :: rhosw_rho_ratio
63 real(dp) :: calv_uw_coeff, r1_calv_uw, r2_calv_uw
64 real(dp) :: target_topo_time_lag, target_topo_weight_lin
65 logical :: flag_calving_event
66 
67 real(dp), parameter :: h_isol_max = 1.0e+03_dp
68  ! Thickness limit of isolated glaciated points
69 
70 #if (CALCZS==3)
71 
72 integer(i4b) :: ierr
73 integer(i4b) :: iter
74 integer(i4b) :: nc, nr
75 integer(i4b), parameter :: nmax = (imax+1)*(jmax+1)
76 integer(i4b), parameter :: n_sprs = 10*(imax+1)*(jmax+1)
77 integer(i4b), allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
78 integer(i4b), allocatable, dimension(:) :: lgs_a_diag_index
79 integer(i4b), allocatable, dimension(:) :: lgs_a_index_trim
80 real(dp), allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
81 real(dp), allocatable, dimension(:) :: lgs_a_value_trim
82 
83 #elif (CALCZS==4)
84 
85 lis_integer :: ierr
86 lis_integer :: iter
87 lis_integer :: nc, nr
88 lis_integer, parameter :: nmax = (imax+1)*(jmax+1)
89 lis_integer, parameter :: n_sprs = 10*(imax+1)*(jmax+1)
90 lis_integer, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
91 lis_integer, allocatable, dimension(:) :: lgs_a_diag_index
92 lis_matrix :: lgs_a
93 lis_vector :: lgs_b, lgs_x
94 lis_scalar, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
95 lis_solver :: solver
96 character(len=256) :: ch_solver_set_option
97 
98 #endif
99 
100 !-------- Term abbreviations --------
101 
102 year_sec_inv = 1.0_dp/year_sec
103 
104 rhosw_rho_ratio = rho_sw/rho
105 
106 dtime_inv = 1.0_dp/dtime
107 
108 azs2 = dtime/(dxi*dxi)
109 azs3 = dtime/(deta*deta)
110 
111 !-------- Computation of zb_neu (updated ice base topography)
112 ! and its time derivative --------
113 
114 zb_neu = zl_neu
115 dzb_dtau = dzl_dtau
116 
117 !-------- Computation of zs_neu (updated ice surface topography) --------
118 
119 h = h_c + h_t
120 
121 zs_neu = zs ! initialisation, will be overwritten later
122 h_neu = h ! initialisation, will be overwritten later
123 
124 ! ------ Explicit scheme
125 
126 #if (CALCZS==0)
127 
128 ! ---- Abbreviations
129 
130 do i=0, imax-1
131 do j=0, jmax
132  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
133  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
134 end do
135 end do
136 
137 do i=0, imax
138 do j=0, jmax-1
139  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
140  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
141 end do
142 end do
143 
144 do i=1, imax-1
145 do j=1, jmax-1
146 
147 #if (MARGIN==1)
148 
149  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
150 
151 #endif
152 
153  zs_neu(j,i) = zs(j,i) &
154  + dtime*(as_perp(j,i)+dzb_dtau(j,i)-q_b_tot(j,i)) &
155  + ( czs2(j,i) *(zs(j,i+1)-zs(j,i) ) &
156  -czs2(j,i-1)*(zs(j,i) -zs(j,i-1)) &
157  +czs3(j,i) *(zs(j+1,i)-zs(j,i) ) &
158  -czs3(j-1,i)*(zs(j,i) -zs(j-1,i)) ) &
159  *insq_g11_g(j,i)*insq_g22_g(j,i)
160 
161 #if (MARGIN==1)
162 
163  else ! maske(j,i) == 2_i2b, sea
164  zs_neu(j,i) = zb_neu(j,i)
165  end if
166 
167 #endif
168 
169 end do
170 end do
171 
172 do i=0, imax
173  zs_neu(0,i) = zb_neu(0,i)
174  zs_neu(jmax,i) = zb_neu(jmax,i)
175 end do
176 
177 do j=0, jmax
178  zs_neu(j,0) = zb_neu(j,0)
179  zs_neu(j,imax) = zb_neu(j,imax)
180 end do
181 
182 ! ------ ADI scheme / over-implicit ADI scheme (DELETED)
183 
184 #elif (CALCZS==1 || CALCZS==2)
185 
186 stop ' calc_top_1: ADI and over-implicit ADI schemes not available any more!'
187 
188 ! ------ Over-implicit scheme
189 
190 #elif (CALCZS==3 || CALCZS==4)
191 
192 ! ---- Abbreviations
193 
194 do i=0, imax-1
195 do j=0, jmax
196  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
197  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
198 end do
199 end do
200 
201 do i=0, imax
202 do j=0, jmax-1
203  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
204  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
205 end do
206 end do
207 
208 ! ---- Assembly of the system of linear equations
209 ! (matrix storage: compressed sparse row CSR)
210 
211 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
212 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
213 
214 lgs_a_value = 0.0_dp
215 lgs_a_index = 0
216 lgs_a_ptr = 0
217 
218 lgs_b_value = 0.0_dp
219 lgs_x_value = 0.0_dp
220 
221 lgs_a_ptr(1) = 1
222 
223 k = 0
224 
225 do nr=1, nmax ! loop over rows
226 
227  i = ii(nr)
228  j = jj(nr)
229 
230  if ( &
231  (i /= 0).and.(i /= imax).and.(j /= 0).and.(j /= jmax) & ! inner point
232 #if (MARGIN==1)
233  .and.(maske(j,i) <= 1_i2b) & ! grounded ice or ice-free land
234 #endif
235  ) then
236 
237  nc = nn(j-1,i) ! smallest nc (column counter), for zs(j-1,i)
238  k = k+1
239  ! if (k > n_sprs) stop ' calc_top_1: n_sprs too small!'
240  lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
241  *insq_g11_g(j,i)*insq_g22_g(j,i)
242  lgs_a_index(k) = nc
243 
244  nc = nn(j,i-1) ! next nc (column counter), for zs(j,i-1)
245  k = k+1
246  ! if (k > n_sprs) stop ' calc_top_1: n_sprs too small!'
247  lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
248  *insq_g11_g(j,i)*insq_g22_g(j,i)
249  lgs_a_index(k) = nc
250 
251  nc = nn(j,i) ! next nc (column counter), for zs(j,i)
252  ! if (nc /= nr) & (diagonal element)
253  ! stop ' calc_top_1: Check for diagonal element failed!'
254  k = k+1
255  ! if (k > n_sprs) stop ' calc_top_1: n_sprs too small!'
256  lgs_a_value(k) = 1.0_dp &
257  + (czs2(j,i)+czs2(j,i-1)+czs3(j,i)+czs3(j-1,i)) &
258  *ovi_weight &
259  *insq_g11_g(j,i)*insq_g22_g(j,i)
260  lgs_a_diag_index(nr) = k
261  lgs_a_index(k) = nc
262 
263  nc = nn(j,i+1) ! next nc (column counter), for zs(j,i+1)
264  k = k+1
265  ! if (k > n_sprs) stop ' calc_top_1: n_sprs too small!'
266  lgs_a_value(k) = -czs2(j,i)*ovi_weight &
267  *insq_g11_g(j,i)*insq_g22_g(j,i)
268  lgs_a_index(k) = nc
269 
270  nc = nn(j+1,i) ! largest nc (column counter), for zs(j+1,i)
271  k = k+1
272  ! if (k > n_sprs) stop ' calc_top_1: n_sprs too small!'
273  lgs_a_value(k) = -czs3(j,i)*ovi_weight &
274  *insq_g11_g(j,i)*insq_g22_g(j,i)
275  lgs_a_index(k) = nc
276 
277  lgs_b_value(nr) = zs(j,i) &
278  + dtime*(as_perp(j,i)+dzb_dtau(j,i)-q_b_tot(j,i)) &
279  + ( czs2(j,i)*(zs(j,i+1)-zs(j,i)) &
280  -czs2(j,i-1)*(zs(j,i)-zs(j,i-1)) &
281  +czs3(j,i)*(zs(j+1,i)-zs(j,i)) &
282  -czs3(j-1,i)*(zs(j,i)-zs(j-1,i)) ) &
283  *(1.0_dp-ovi_weight) &
284  *insq_g11_g(j,i)*insq_g22_g(j,i)
285  ! right-hand side
286 
287  else ! zero-ice-thickness boundary condition
288 
289  k = k+1
290  ! if (k > n_sprs) stop ' calc_top_1: n_sprs too small!'
291  lgs_a_value(k) = 1.0_dp ! diagonal element only
292  lgs_a_diag_index(nr) = k
293  lgs_a_index(k) = nr
294  lgs_b_value(nr) = zb_neu(j,i)
295 
296  end if
297 
298  lgs_x_value(nr) = zs(j,i) ! old surface topography,
299  ! initial guess for solution vector
300 
301  lgs_a_ptr(nr+1) = k+1 ! row is completed, store index to next row
302 
303 end do
304 
305 nnz = k ! number of non-zero elements of the matrix
306 
307 #if (CALCZS==3)
308 
309 ! ---- Solution of the system of linear equations
310 
311 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
312 
313 do k=1, nnz ! relocate matrix to trimmed arrays
314  lgs_a_value_trim(k) = lgs_a_value(k)
315  lgs_a_index_trim(k) = lgs_a_index(k)
316 end do
317 
318 deallocate(lgs_a_value, lgs_a_index)
319 
320 eps_sor = 1.0e-05_dp*mean_accum*dtime ! convergence parameter
321 
322 call sor_sprs(lgs_a_value_trim, &
323  lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
324  lgs_b_value, &
325  nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
326 
327 do nr=1, nmax
328  i = ii(nr)
329  j = jj(nr)
330  zs_neu(j,i) = lgs_x_value(nr)
331 end do
332 
333 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
334 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
335 
336 #elif (CALCZS==4)
337 
338 ! ---- Settings for Lis
339 
340 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
341 call lis_vector_create(lis_comm_world, lgs_b, ierr)
342 call lis_vector_create(lis_comm_world, lgs_x, ierr)
343 
344 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
345 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
346 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
347 
348 do nr=1, nmax
349 
350  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
351  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
352  lgs_a_value(nc), lgs_a, ierr)
353  end do
354 
355  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
356  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
357 
358 end do
359 
360 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
361 call lis_matrix_assemble(lgs_a, ierr)
362 
363 ! ---- Solution of the system of linear equations with Lis
364 
365 call lis_solver_create(solver, ierr)
366 
367 ch_solver_set_option = '-i bicg -p ilu '// &
368  '-maxiter 1000 -tol 1.0e-12 -initx_zeros false'
369 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
370 
371 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
372 
373 call lis_solver_get_iters(solver, iter, ierr)
374 write(6,'(11x,a,i5)') 'calc_top_1: iter = ', iter
375 
376 lgs_x_value = 0.0_dp
377 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
378 call lis_matrix_destroy(lgs_a, ierr)
379 call lis_vector_destroy(lgs_b, ierr)
380 call lis_vector_destroy(lgs_x, ierr)
381 call lis_solver_destroy(solver, ierr)
382 
383 do nr=1, nmax
384  i = ii(nr)
385  j = jj(nr)
386  zs_neu(j,i) = lgs_x_value(nr)
387 end do
388 
389 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
390 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
391 
392 #endif
393 
394 #endif
395 
396 ! ------ Adjustment due to prescribed target topography
397 
398 #if (ZS_EVOL==2)
399 
400 if (time >= time_target_topo_final) then
401  zs_neu = zs_target
402 else if (time >= time_target_topo_init) then
403  zs_neu = ( target_topo_time_lag*zs_neu + dtime*zs_target ) &
404  / ( target_topo_time_lag + dtime )
405 end if
406 
407 #endif
408 
409 !-------- Update of the mask --------
410 
411 h_neu = zs_neu - zb_neu ! ice thickness
412 
413 #if ((MARGIN==2) \
414  && (marine_ice_formation==2) \
415  && (marine_ice_calving==9))
416 
417 calv_uw_coeff = calv_uw_coeff * year_sec_inv
418 r1_calv_uw = r1_calv_uw
419 r2_calv_uw = r2_calv_uw
420 
421 #if (ZS_EVOL==2)
422 if (time >= time_target_topo_final) then
423  calv_uw_coeff = 0.0_dp
424  ! adjustment due to prescribed target topography
425 else if (time >= time_target_topo_init) then
426  target_topo_weight_lin = (time-time_target_topo_init) &
427  /(time_target_topo_final-time_target_topo_init)
428  calv_uw_coeff = (1.0_dp-target_topo_weight_lin)*calv_uw_coeff
429  ! adjustment due to prescribed target topography
430 end if
431 #endif
432 
433 #endif
434 
435 #if (ZS_EVOL>=1)
436 
437 do i=0, imax
438 do j=0, jmax
439 
440  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
441 
442  h_sea_neu(j,i) = z_sl-zl_neu(j,i) ! sea depth (only makes sense for
443  ! marine ice and "underwater ice",
444  ! otherwise meaningless)
445 
446  if (h_neu(j,i) >= h_min) then
447  maske(j,i) = 0_i2b ! grounded ice
448 #if ((MARGIN==2) \
449  && (marine_ice_formation==2) \
450  && (marine_ice_calving==9))
451  if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then
452  ! "underwater ice"
453  h_neu(j,i) = h_neu(j,i) - ( calv_uw_coeff &
454  *h_neu(j,i)**r1_calv_uw &
455  *h_sea_neu(j,i)**r2_calv_uw ) * dtime
456  ! calving of "underwater ice"
457  zs_neu(j,i) = zb_neu(j,i) + h_neu(j,i)
458  if (h_neu(j,i) < h_min) maske(j,i) = 2_i2b ! sea
459  end if
460 #endif
461  else
462  maske(j,i) = 1_i2b ! ice-free land
463  end if
464 
465 #if (MARGIN==2)
466 
467  else ! (maske(j,i) == 2_i2b); sea
468 
469  h_sea_neu(j,i) = z_sl-zl_neu(j,i) ! sea depth
470 
471  if (h_neu(j,i) >= h_min) then
472 
473  if ( h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) ) then
474 
475 #if (MARINE_ICE_FORMATION==1)
476  maske(j,i) = 2_i2b ! floating ice cut off -> sea
477 #elif (MARINE_ICE_FORMATION==2)
478  maske(j,i) = 0_i2b ! "underwater ice"
479 #if (MARINE_ICE_CALVING==9)
480  h_neu(j,i) = h_neu(j,i) - ( calv_uw_coeff &
481  *h_neu(j,i)**r1_calv_uw &
482  *h_sea_neu(j,i)**r2_calv_uw ) * dtime
483  ! calving of "underwater ice"
484  zs_neu(j,i) = zb_neu(j,i) + h_neu(j,i)
485  if (h_neu(j,i) < h_min) maske(j,i) = 2_i2b ! sea
486 #endif
487 #endif
488  else
489  maske(j,i) = 0_i2b ! grounded ice
490  end if
491 
492 #if (MARINE_ICE_CALVING==2 \
493  || marine_ice_calving==4 \
494  || marine_ice_calving==6)
495  if (zl0(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
496 #elif (MARINE_ICE_CALVING==3 \
497  || marine_ice_calving==5 \
498  || marine_ice_calving==7)
499  if (zl_neu(j,i) < z_mar) maske(j,i) = 2_i2b ! sea
500 #endif
501 
502  else
503  maske(j,i) = 2_i2b ! sea
504  end if
505 
506 #endif
507 
508  end if
509 
510 end do
511 end do
512 
513 #endif
514 
515 ! ------ Adjustment due to prescribed target topography
516 
517 #if (ZS_EVOL==2)
518 if (time >= time_target_topo_final) maske = maske_target
519 #endif
520 
521 ! ------ Adjustment due to prescribed maximum ice extent
522 
523 #if (ZS_EVOL==3)
524 
525 do i=0, imax
526 do j=0, jmax
527 
528  if (maske_maxextent(j,i)==0_i2b) then ! not allowed to glaciate
529 
530  if (maske(j,i)==0_i2b) then
531  maske(j,i) = 1_i2b ! grounded ice -> ice-free land
532  else if (maske(j,i)==3_i2b) then
533  maske(j,i) = 2_i2b ! floating ice -> sea
534  end if
535 
536  end if
537 
538 end do
539 end do
540 
541 #endif
542 
543 !-------- Detection of the grounding line and the calving front
544 ! (not relevant for SIA-only dynamics) --------
545 
546 flag_grounding_line_1 = .false.
547 flag_grounding_line_2 = .false.
548 
549 flag_calving_front_1 = .false.
550 flag_calving_front_2 = .false.
551 
552 !-------- Correction of zs_neu and zb_neu for ice-free land and sea --------
553 
554 do i=0, imax
555 do j=0, jmax
556 
557  if (maske(j,i) == 1_i2b) then ! ice-free land
558 
559  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
560  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
561 
562  else if (maske(j,i) == 2_i2b) then ! sea
563 
564  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
565  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
566 
567  else if (maske(j,i) == 3_i2b) then ! floating ice
568 
569  stop ' calc_top_1: maske(j,i)==3 not allowed for SIA-only dynamics!'
570 
571  end if
572 
573 end do
574 end do
575 
576 !-------- Computation of further quantities --------
577 
578 #if (ZS_EVOL==0)
579 
580 do i=0, imax
581 do j=0, jmax
582 
583  if (maske(j,i) == 0_i2b) then ! grounded ice
584 
585  as_perp(j,i) = as_perp(j,i) - (zs_neu(j,i)-zs(j,i))*dtime_inv
586  ! correction of the accumulation-ablation function
587  ! such that it is consistent with a steady surface
588 
589  else ! maske(j,i)==1_i2b or 2_i2b, ice-free land or sea
590 
591  as_perp(j,i) = 0.0_dp
592 
593  end if
594 
595  zs_neu(j,i) = zs(j,i) ! newly computed surface topography
596  h_neu(j,i) = zs_neu(j,i)-zb_neu(j,i) ! is discarded
597 
598 end do
599 end do
600 
601 #endif
602 
603 do i=0, imax
604 do j=0, jmax
605 
606  if (maske(j,i) == 0_i2b) then ! grounded ice
607 
608 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
609 
610  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
611  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
612  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
613  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
614  ) then ! all nearest neighbours ice-free
615  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
616  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
617  end if
618 
619 ! ------ Further stuff
620 
621  if (n_cts(j,i) == 1) then
622  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
623  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
624  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
625  else
626  h_t_neu(j,i) = 0.0_dp
627  zm_neu(j,i) = zb_neu(j,i)
628  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
629  end if
630 
631  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
632 
633  zs_neu(j,i) = zb_neu(j,i)
634  zm_neu(j,i) = zb_neu(j,i)
635  h_c_neu(j,i) = 0.0_dp
636  h_t_neu(j,i) = 0.0_dp
637  h_neu(j,i) = 0.0_dp
638 
639  end if
640 
641  dzs_dtau(j,i) = (zs_neu(j,i)-zs(j,i))*dtime_inv
642  dzb_dtau(j,i) = (zb_neu(j,i)-zb(j,i))*dtime_inv
643  dzm_dtau(j,i) = dh_t_dtau(j,i)+dzb_dtau(j,i)
644  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
645 
646 #if (ZS_EVOL==2)
647  if ( abs((time-time_target_topo_final)*year_sec_inv) < eps ) then
648  dzs_dtau = 0.0_dp
649  dzb_dtau = 0.0_dp
650  dzm_dtau = 0.0_dp
651  dh_c_dtau(j,i) = 0.0_dp
652  dh_t_dtau(j,i) = 0.0_dp
653  end if
654 #endif
655 
656 end do
657 end do
658 
659 !-------- New gradients --------
660 
661 #if (TOPOGRAD==0)
662 call topograd_1(dxi, deta, 2)
663 #elif (TOPOGRAD==1)
664 call topograd_2(dxi, deta, 2)
665 #endif
666 
667 !-------- Compute the volume flux across the CTS, am_perp --------
668 
669 #if (CALCMOD==1)
670 
671 do i=1, imax-1
672 do j=1, jmax-1
673 
674  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
675 
676  am_perp_st(j,i) = &
677  0.5_dp*(vx_c(0,j,i)+vx_c(0,j,i-1))*dzm_dxi_g(j,i) &
678  + 0.5_dp*(vy_c(0,j,i)+vy_c(0,j-1,i))*dzm_deta_g(j,i) &
679  - 0.5_dp*(vz_c(0,j,i)+vz_t(ktmax-1,j,i))
680  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
681 
682  else
683  am_perp_st(j,i) = 0.0_dp
684  am_perp(j,i) = 0.0_dp
685  end if
686 
687 end do
688 end do
689 
690 #endif
691 
692 end subroutine calc_top_1
693 !
topograd_2
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
Definition: topograd_2.F90:41
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
topograd_1
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_1.F90:41
sico_sle_solvers::sor_sprs
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, lgs_b_value, nnz, nmax, omega, eps_sor, lgs_x_value, ierr)
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
Definition: sico_sle_solvers.F90:47
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
sico_sle_solvers
Solvers for systems of linear equations used by SICOPOLIS.
Definition: sico_sle_solvers.F90:35
calc_top_1
subroutine calc_top_1(time, z_sl, z_mar, dtime, dxi, deta, mean_accum, ii, jj, nn)
Computation of the ice topography (including gradients) for SIA-only dynamics of ice sheets without i...
Definition: calc_top_1.F90:37
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
Generated on Wed Mar 2 2016 11:26:55 for SICOPOLIS V3.2 by   doxygen 1.8.5