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