SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_top_3.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t o p _ 3
4 !
5 !> @file
6 !!
7 !! Computation of the ice topography (including gradients) for
8 !! coupled SIA/SSA or SIA/SStA/SSA dynamics of ice sheets with 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 !! coupled SIA/SSA or SIA/SStA/SSA dynamics of ice sheets with ice shelves.
36 !<------------------------------------------------------------------------------
37 subroutine calc_top_3(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) :: dt_dxi, dt_deta
61 real(dp) :: azs2, azs3
62 real(dp), dimension(0:JMAX,0:IMAX) :: h, h_neu
63 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
64 real(dp) :: vx_m_help, vy_m_help
65 real(dp) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
66 real(dp) :: uph_x_1, uph_x_2, uph_y_1, uph_y_2
67 real(dp) :: h_balance
68 real(dp) :: target_topo_time_lag, target_topo_weight_lin
69 logical :: flag_calving_event
70 
71 real(dp), parameter :: h_isol_max = 1.0e+03_dp
72  ! Thickness limit of isolated glaciated points
73 
74 #if (CALCZS==3)
75 
76 integer(i4b) :: ierr
77 integer(i4b) :: iter
78 integer(i4b) :: nc, nr
79 integer(i4b), parameter :: nmax = (imax+1)*(jmax+1)
80 integer(i4b), parameter :: n_sprs = 10*(imax+1)*(jmax+1)
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 
87 #elif (CALCZS==4)
88 
89 ! LIS_INTEGER :: ierr
90 ! LIS_INTEGER :: iter
91 ! LIS_INTEGER :: nc, nr
92 ! LIS_INTEGER, parameter :: nmax = (IMAX+1)*(JMAX+1)
93 ! LIS_INTEGER, parameter :: n_sprs = 10*(IMAX+1)*(JMAX+1)
94 ! LIS_INTEGER, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
95 ! LIS_INTEGER, allocatable, dimension(:) :: lgs_a_diag_index
96 ! LIS_MATRIX :: lgs_a
97 ! LIS_VECTOR :: lgs_b, lgs_x
98 ! LIS_SCALAR, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
99 ! LIS_SOLVER :: solver
100 ! character(len=256) :: ch_solver_set_option
101 
102 stop ' calc_top_3: Option CALCZS==4 not supported by this routine!'
103 
104 #endif
105 
106 !-------- Term abbreviations --------
107 
108 year_sec_inv = 1.0_dp/year_sec
109 
110 rho_rhosw_ratio = rho/rho_sw
111 rhosw_rho_ratio = rho_sw/rho
112 
113 dtime_inv = 1.0_dp/dtime
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 zb_neu (updated ice base topography)
122 ! and its time derivative --------
123 
124 do i=0, imax
125 do j=0, jmax
126 
127  if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
128 
129  zb_neu(j,i) = zl_neu(j,i)
130  dzb_dtau(j,i) = dzl_dtau(j,i)
131 
132  else ! (maske(j,i) >= 2_i2b)
133 
134  zb_neu(j,i) = zb(j,i) ! this is only preliminary
135  dzb_dtau(j,i) = 0.0_dp ! and will be corrected later
136 
137  end if
138 
139 end do
140 end do
141 
142 !-------- Computation of H_neu (updated ice thickness) --------
143 
144 h = h_c + h_t
145 
146 zs_neu = zs ! initialisation, will be overwritten later
147 h_neu = h ! initialisation, will be overwritten later
148 
149 ! ------ Explicit scheme
150 
151 #if (CALCZS==0)
152 
153 ! ---- Abbreviations
154 
155 do i=0, imax-1
156 do j=0, jmax
157  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
158  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
159 end do
160 end do
161 
162 do i=0, imax
163 do j=0, jmax-1
164  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
165  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
166 end do
167 end do
168 
169 do i=1, imax-1
170 do j=1, jmax-1
171 
172  if ( (maske(j,i) <= 1_i2b) &
173  .and. &
174  (.not.flag_grounding_line_1(j,i)) &
175  .and. &
176  (.not.flag_shelfy_stream(j,i)) &
177  ) then
178  ! grounded ice or ice-free land,
179  ! but neither grounding line nor shelfy stream
180 
181  h_neu(j,i) = h(j,i) &
182  + dtime*(as_perp(j,i)-q_b_tot(j,i)) &
183  + ( czs2(j,i) *(zs(j,i+1)-zs(j,i) ) &
184  -czs2(j,i-1)*(zs(j,i) -zs(j,i-1)) &
185  +czs3(j,i) *(zs(j+1,i)-zs(j,i) ) &
186  -czs3(j-1,i)*(zs(j,i) -zs(j-1,i)) ) &
187  *insq_g11_g(j,i)*insq_g22_g(j,i)
188 
189  else ! floating ice, grounding line, shelfy stream or sea
190 
191  vx_m_help = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1))
192  vy_m_help = 0.5_dp*(vy_m(j,i)+vy_m(j-1,i))
193 
194  vx_m_1 = min(vx_m_help, 0.0_dp)
195  vx_m_2 = max(vx_m_help, 0.0_dp)
196  vy_m_1 = min(vy_m_help, 0.0_dp)
197  vy_m_2 = max(vy_m_help, 0.0_dp)
198 
199  uph_x_1 = h(j,i+1)-h(j,i )
200  uph_x_2 = h(j,i )-h(j,i-1)
201  uph_y_1 = h(j+1,i)-h(j ,i)
202  uph_y_2 = h(j ,i)-h(j-1,i)
203 
204  h_neu(j,i) = h(j,i) &
205  + dtime*(as_perp(j,i)-q_b_tot(j,i)) &
206  - dt_dxi &
207  * ( vx_m_1*uph_x_1 + vx_m_2*uph_x_2 &
208  + (vx_m(j,i)-vx_m(j,i-1))*h(j,i) ) &
209  - dt_deta &
210  * ( vy_m_1*uph_y_1 + vy_m_2*uph_y_2 &
211  + (vy_m(j,i)-vy_m(j-1,i))*h(j,i) )
212  end if
213 
214 end do
215 end do
216 
217 do i=0, imax
218  h_neu(0,i) = 0.0_dp
219  h_neu(jmax,i) = 0.0_dp
220 end do
221 
222 do j=0, jmax
223  h_neu(j,0) = 0.0_dp
224  h_neu(j,imax) = 0.0_dp
225 end do
226 
227 ! ------ ADI scheme / over-implicit ADI scheme (DELETED)
228 
229 #elif (CALCZS==1 || CALCZS==2)
230 
231 stop ' calc_top_3: Options CALCZS==1 or 2 not supported by this routine!'
232 
233 ! ------ Over-implicit ice thickness scheme
234 
235 #elif (CALCZS==3)
236 
237 do i=0, imax-1
238 do j=0, jmax
239  czs2(j,i) = azs2*0.5_dp*(h_diff(j,i)+h_diff(j,i+1)) &
240  *sq_g22_sgx(j,i)*insq_g11_sgx(j,i)
241 end do
242 end do
243 
244 do i=0, imax
245 do j=0, jmax-1
246  czs3(j,i) = azs3*0.5_dp*(h_diff(j,i)+h_diff(j+1,i)) &
247  *sq_g11_sgy(j,i)*insq_g22_sgy(j,i)
248 end do
249 end do
250 
251 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
252 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
253 
254 lgs_a_value = 1.0_dp
255 lgs_b_value = 0.0_dp
256 lgs_x_value = 0.0_dp
257 
258 k = 0
259 
260 ! ---- Main loop for matrix assembly
261 
262 do nr=1, nmax
263 
264  i = ii(nr) ; j = jj(nr) ! set numbering
265 
266  if ( (i /= 0).and.(i /= imax).and.(j /= 0).and.(j /= jmax) ) then
267  ! non-edge points
268  if ( (maske(j,i) <= 1_i2b) &
269  .and. &
270  (.not.flag_grounding_line_1(j,i)) &
271  .and. &
272  (.not.flag_shelfy_stream(j,i)) &
273  ) then
274  ! grounded ice or ice-free land,
275  ! but neither grounding line nor shelfy stream
276 
277  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=h(j,i)
278 
279  k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc ! for H(j-1,i)
280  lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
281  *insq_g11_g(j,i)*insq_g22_g(j,i)
282 
283  k=k+1 ; nc=nn(j,i-1) ; lgs_a_index(k)=nc ! for H(j,i-1)
284  lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
285  *insq_g11_g(j,i)*insq_g22_g(j,i)
286 
287  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k ! for H(j,i)
288  lgs_a_value(k) = 1.0_dp + (czs2(j,i)+czs2(j,i-1) & ! (diagonal
289  +czs3(j,i)+czs3(j-1,i)) & ! element)
290  *ovi_weight &
291  *insq_g11_g(j,i)*insq_g22_g(j,i)
292 
293  k=k+1 ; nc=nn(j,i+1) ; lgs_a_index(k)=nc ! for H(j,i+1)
294  lgs_a_value(k) = -czs2(j,i)*ovi_weight &
295  *insq_g11_g(j,i)*insq_g22_g(j,i)
296 
297  k=k+1 ; nc=nn(j+1,i) ; lgs_a_index(k)=nc ! for H(j+1,i)
298  lgs_a_value(k) = -czs3(j,i)*ovi_weight &
299  *insq_g11_g(j,i)*insq_g22_g(j,i)
300 
301  lgs_b_value(nr) = h(j,i) &
302  + dtime*(as_perp(j,i)-q_b_tot(j,i)) &
303  + ( czs2(j,i) *(h(j,i+1)-h(j,i )) &
304  -czs2(j,i-1)*(h(j,i )-h(j,i-1)) &
305  +czs3(j,i) *(h(j+1,i)-h(j,i )) &
306  -czs3(j-1,i)*(h(j,i )-h(j-1,i)) ) &
307  *(1.0_dp-ovi_weight) &
308  *insq_g11_g(j,i)*insq_g22_g(j,i) &
309  + ( czs2(j,i) *(zb_neu(j,i+1)-zb_neu(j,i )) &
310  -czs2(j,i-1)*(zb_neu(j,i )-zb_neu(j,i-1)) &
311  +czs3(j,i) *(zb_neu(j+1,i)-zb_neu(j,i )) &
312  -czs3(j-1,i)*(zb_neu(j,i )-zb_neu(j-1,i)) ) &
313  *insq_g11_g(j,i)*insq_g22_g(j,i)
314  ! right-hand side
315 
316 #if (!defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1)
317 
318  else ! floating ice, grounding line, shelfy stream or sea
319 
320 #elif (ZS_EVOL_SSA==0)
321 
322  else if ( flag_grounding_line_1(j,i) &
323  .or. &
324  flag_shelfy_stream(j,i) &
325  ) then
326  ! grounding line or shelfy stream
327 #else
328  else
329  stop ' calc_top_3: ZS_EVOL_SSA must be either 1 or 0!'
330 #endif
331 
332  vx_m_help = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1))
333  vy_m_help = 0.5_dp*(vy_m(j,i)+vy_m(j-1,i))
334 
335  vx_m_1 = min(vx_m_help, 0.0_dp)
336  vx_m_2 = max(vx_m_help, 0.0_dp)
337  vy_m_1 = min(vy_m_help, 0.0_dp)
338  vy_m_2 = max(vy_m_help, 0.0_dp)
339 
340  uph_x_1 = h(j,i+1)-h(j,i )
341  uph_x_2 = h(j,i )-h(j,i-1)
342  uph_y_1 = h(j+1,i)-h(j ,i)
343  uph_y_2 = h(j ,i)-h(j-1,i)
344 
345  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=h(j,i)
346 
347  k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc ! for H(j-1,i)
348  lgs_a_value(k) = -vy_m_2*dt_deta*ovi_weight !&
349  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
350 
351  k=k+1 ; nc=nn(j,i-1) ; lgs_a_index(k)=nc ! for H(j,i-1)
352  lgs_a_value(k) = -vx_m_2*dt_dxi *ovi_weight !&
353  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
354 
355  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k ! for H(j,i)
356  lgs_a_value(k) = 1.0_dp & ! (diagonal
357  +ovi_weight & ! element)
358  *( dt_dxi &
359  * ( -vx_m_1+vx_m_2 &
360  +(vx_m(j,i)-vx_m(j,i-1)) ) &
361  +dt_deta &
362  * ( -vy_m_1+vy_m_2 &
363  +(vy_m(j,i)-vy_m(j-1,i)) ) ) !&
364  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
365 
366  k=k+1 ; nc=nn(j,i+1) ; lgs_a_index(k)=nc ! for H(j,i+1)
367  lgs_a_value(k) = vx_m_1*dt_dxi *ovi_weight !&
368  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
369 
370  k=k+1 ; nc=nn(j+1,i) ; lgs_a_index(k)=nc ! for H(j+1,i)
371  lgs_a_value(k) = vy_m_1*dt_deta*ovi_weight !&
372  ! *insq_g11_g(j,i)*insq_g22_g(j,i)
373 
374  lgs_b_value(nr)= h(j,i) &
375  +dtime*(as_perp(j,i)-q_b_tot(j,i)) &
376  -(1.0_dp-ovi_weight) &
377  *( dt_dxi &
378  * ( vx_m_1*uph_x_1 + vx_m_2*uph_x_2 &
379  + (vx_m(j,i)-vx_m(j,i-1))*h(j,i) ) &
380  +dt_deta &
381  * ( vy_m_1*uph_y_1 + vy_m_2*uph_y_2 &
382  + (vy_m(j,i)-vy_m(j-1,i))*h(j,i) ) )
383  ! right-hand side
384 
385 #if (!defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1)
386 
387  continue ! do nothing
388 
389 #elif (ZS_EVOL_SSA==0)
390 
391  else if (maske(j,i)==3_i2b) then ! floating ice
392 
393  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=h(j,i)
394  lgs_b_value(nr) = h(j,i)
395 
396  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
397  lgs_a_value(k) = 1.0_dp
398 
399  else if (maske(j,i)==2_i2b) then ! sea
400 
401  lgs_a_ptr(nr)=k+1 ; lgs_x_value(nr)=0.0_dp
402  lgs_b_value(nr) = 0.0_dp
403 
404  k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
405  lgs_a_value(k) = 1.0_dp
406 
407 #else
408  stop ' calc_top_3: ZS_EVOL_SSA must be either 1 or 0!'
409 #endif
410 
411  end if
412 
413  else ! zero-ice-thickness boundary condition
414 
415  lgs_a_ptr(nr) = k+1
416 
417  k = k+1
418  lgs_a_value(k) = 1.0_dp ! diagonal element only
419  lgs_a_diag_index(nr) = k
420  lgs_a_index(k) = nr
421  lgs_x_value(nr) = 0.0_dp
422  lgs_b_value(nr) = 0.0_dp
423 
424  end if
425 
426 end do
427 
428 lgs_a_ptr(nmax+1)=k+1 ! set csr's last index
429 
430 nnz = k ! number of non-zero elements of the matrix
431 
432 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
433 
434 do k=1, nnz ! relocate matrix to trimmed arrays
435  lgs_a_value_trim(k) = lgs_a_value(k)
436  lgs_a_index_trim(k) = lgs_a_index(k)
437 end do
438 
439 deallocate(lgs_a_value, lgs_a_index)
440 
441 eps_sor = 1.0e-05_dp*mean_accum*dtime ! convergence parameter
442 
443 call sor_sprs(lgs_a_value_trim, &
444  lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
445  lgs_b_value, &
446  nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
447 
448 do nr=1, nmax
449  i = ii(nr)
450  j = jj(nr)
451  h_neu(j,i) = lgs_x_value(nr)
452 end do
453 
454 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
455 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
456 
457 #endif
458 
459 ! ------ Adjustment due to prescribed target topography
460 
461 #if (ZS_EVOL==2)
462 
463 if (time >= time_target_topo_final) then
464  h_neu = h_target
465 else if (time >= time_target_topo_init) then
466  h_neu = ( target_topo_time_lag*h_neu + dtime*h_target ) &
467  / ( target_topo_time_lag + dtime )
468 end if
469 
470 #endif
471 
472 !-------- Update of the mask --------
473 
474 ! zs_neu = zb_neu + H_neu ! ice surface topography
475 !
476 ! For ZS_EVOL==0 (no evolution of the free surface,
477 ! like in run v32_ant40_sr_spinup01_fixtopo1), this causes troubles.
478 ! Therefore commented out for the time being.
479 
480 #if (ZS_EVOL>=1)
481 
482 do i=1, imax-1
483 do j=1, jmax-1
484 
485  h_balance = (z_sl-zl_neu(j,i))*rhosw_rho_ratio
486 
487 #if (!defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1)
488 
489  ! grounding_line migration check
490 
491  if ( ( maske(j,i) <= 1_i2b &
492  .and. (maske(j,i+1)>1_i2b.or.maske(j,i-1)>1_i2b &
493  .or.maske(j+1,i)>1_i2b.or.maske(j-1,i)>1_i2b) ) &
494  .or. &
495  ( maske(j,i)>=2_i2b &
496  .and. (maske(j,i+1)==0_i2b.or.maske(j,i-1)==0_i2b &
497  .or.maske(j+1,i)==0_i2b.or.maske(j-1,i)==0_i2b) ) ) then
498 
499  if (h_neu(j,i)>=h_min) then
500 
501  if (h_neu(j,i)<h_balance.and.zl_neu(j,i)<z_sl) then
502  maske(j,i) = 3_i2b
503  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_neu(j,i)
504  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
505  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
506  else
507  maske(j,i) = 0_i2b
508  zb_neu(j,i) = zl_neu(j,i)
509  dzb_dtau(j,i) = dzl_dtau(j,i)
510  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
511  end if
512 
513  else ! if (H_neu(j,i)<=H_MIN) then
514 
515  if (zl_neu(j,i)>z_sl) then
516  maske(j,i) = 1_i2b
517  zb_neu(j,i) = zl_neu(j,i)
518  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
519  zs_neu(j,i) = zb_neu(j,i)
520  else
521  maske(j,i) = 2_i2b
522  zb_neu(j,i) = z_sl
523  dzb_dtau(j,i) = 0.0_dp
524  zs_neu(j,i) = z_sl
525  end if
526 
527  end if
528 
529 #elif (ZS_EVOL_SSA==0)
530 
531  if (maske(j,i)==2_i2b) then
532  maske(j,i) = 2_i2b
533  zb_neu(j,i) = z_sl
534  dzb_dtau(j,i) = 0.0_dp
535  zs_neu(j,i) = z_sl
536  else if (maske(j,i)==3_i2b) then
537  maske(j,i) = 3_i2b
538  zb_neu(j,i) = zb(j,i)
539  dzb_dtau(j,i) = 0.0_dp
540  zs_neu(j,i) = zs(j,i)
541 
542 #else
543  stop ' calc_top_3: ZS_EVOL_SSA must be either 1 or 0!'
544 #endif
545 
546  else if (maske(j,i) <= 1_i2b) then ! grounded ice or ice-free land
547 
548  if (h_neu(j,i)>=h_min) then ! can change maske 0 or 1
549  maske(j,i) = 0_i2b
550  zb_neu(j,i) = zl_neu(j,i)
551  dzb_dtau(j,i) = dzl_dtau(j,i)
552  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
553  else
554  maske(j,i) = 1_i2b
555  zb_neu(j,i) = zl_neu(j,i)
556  dzb_dtau(j,i) = dzl_dtau(j,i)
557  zs_neu(j,i) = zl_neu(j,i)
558  end if
559 
560 #if (!defined(ZS_EVOL_SSA) || ZS_EVOL_SSA==1)
561 
562  else ! if (maske(j,i)==2.or.maske(j,i)==3) then
563 
564  if (h_neu(j,i)>h_min) then
565 
566  if (h_neu(j,i)<h_balance.and.zl_neu(j,i)<z_sl) then
567  maske(j,i) = 3_i2b
568  zb_neu(j,i) = z_sl-rho_rhosw_ratio*h_neu(j,i)
569  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
570  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
571  else
572  maske(j,i) = 0_i2b
573  zb_neu(j,i) = zl_neu(j,i)
574  dzb_dtau(j,i) = dzl_dtau(j,i)
575  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
576  end if
577 
578  else ! if (H_neu(j,i)<=H_MIN) then
579 
580  if (zl_neu(j,i)>z_sl) then
581  maske(j,i) = 1_i2b
582  zb_neu(j,i) = zl_neu(j,i)
583  dzb_dtau(j,i) = dtime_inv*(zb_neu(j,i)-zb(j,i))
584  zs_neu(j,i) = zb_neu(j,i)
585  else
586  maske(j,i) = 2_i2b
587  zb_neu(j,i) = z_sl
588  dzb_dtau(j,i) = 0.0_dp
589  zs_neu(j,i) = z_sl
590  end if
591 
592  end if
593 
594 #endif
595 
596  end if
597 
598 end do
599 end do
600 
601 #if (ICE_SHELF_CALVING==2)
602 
603 do
604 
605  flag_calving_front_1 = .false.
606  flag_calving_event = .false.
607 
608  do i=1, imax-1
609  do j=1, jmax-1
610 
611  if ( (maske(j,i)==3_i2b) & ! floating ice
612  .and. &
613  ( (maske(j,i+1)==2_i2b) & ! with
614  .or.(maske(j,i-1)==2_i2b) & ! one
615  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
616  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
617  ) &
618  flag_calving_front_1(j,i) = .true. ! preliminary detection
619  ! of the calving front
620 
621  if ( flag_calving_front_1(j,i).and.(h_neu(j,i) < h_calv) ) then
622  flag_calving_event = .true. ! calving event,
623  maske(j,i) = 2_i2b ! floating ice point changes to sea point
624  end if
625 
626  end do
627  end do
628 
629  if (.not.flag_calving_event) exit
630 
631 end do
632 
633 #endif
634 
635 #endif
636 
637 ! ------ Adjustment due to prescribed target topography
638 
639 #if (ZS_EVOL==2)
640 if (time >= time_target_topo_final) maske = maske_target
641 #endif
642 
643 ! ------ Adjustment due to prescribed maximum ice extent
644 
645 #if (ZS_EVOL==3)
646 
647 do i=0, imax
648 do j=0, jmax
649 
650  if (maske_maxextent(j,i)==0_i2b) then ! not allowed to glaciate
651 
652  if (maske(j,i)==0_i2b) then
653  maske(j,i) = 1_i2b ! grounded ice -> ice-free land
654  else if (maske(j,i)==3_i2b) then
655  maske(j,i) = 2_i2b ! floating ice -> sea
656  end if
657 
658  end if
659 
660 end do
661 end do
662 
663 #endif
664 
665 !-------- Detection of the grounding line and the calving front --------
666 
667 flag_grounding_line_1 = .false.
668 flag_grounding_line_2 = .false.
669 
670 flag_calving_front_1 = .false.
671 flag_calving_front_2 = .false.
672 
673 do i=1, imax-1
674 do j=1, jmax-1
675 
676  if ( (maske(j,i)==0_i2b) & ! grounded ice
677  .and. &
678  ( (maske(j,i+1)==3_i2b) & ! with
679  .or.(maske(j,i-1)==3_i2b) & ! one
680  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
681  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
682  ) &
683  flag_grounding_line_1(j,i) = .true.
684 
685  if ( (maske(j,i)==3_i2b) & ! floating ice
686  .and. &
687  ( (maske(j,i+1)==0_i2b) & ! with
688  .or.(maske(j,i-1)==0_i2b) & ! one
689  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
690  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
691  ) &
692  flag_grounding_line_2(j,i) = .true.
693 
694  if ( (maske(j,i)==3_i2b) & ! floating ice
695  .and. &
696  ( (maske(j,i+1)==2_i2b) & ! with
697  .or.(maske(j,i-1)==2_i2b) & ! one
698  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
699  .or.(maske(j-1,i)==2_i2b) ) & ! sea point
700  ) &
701  flag_calving_front_1(j,i) = .true.
702 
703  if ( (maske(j,i)==2_i2b) & ! sea
704  .and. &
705  ( (maske(j,i+1)==3_i2b) & ! with
706  .or.(maske(j,i-1)==3_i2b) & ! one
707  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
708  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
709  ) &
710  flag_calving_front_2(j,i) = .true.
711 
712 end do
713 end do
714 
715 do i=1, imax-1
716 
717  j=0
718  if ( (maske(j,i)==2_i2b) & ! sea
719  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
720  ) & ! floating ice point
721  flag_calving_front_2(j,i) = .true.
722 
723  j=jmax
724  if ( (maske(j,i)==2_i2b) & ! sea
725  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
726  ) & ! floating ice point
727  flag_calving_front_2(j,i) = .true.
728 
729 end do
730 
731 do j=1, jmax-1
732 
733  i=0
734  if ( (maske(j,i)==2_i2b) & ! sea
735  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
736  ) & ! floating ice point
737  flag_calving_front_2(j,i) = .true.
738 
739  i=imax
740  if ( (maske(j,i)==2_i2b) & ! sea
741  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
742  ) & ! floating ice point
743  flag_calving_front_2(j,i) = .true.
744 
745 end do
746 
747 !-------- Correction of zs_neu and zb_neu
748 ! for ice-free land, sea and floating ice --------
749 
750 do i=0, imax
751 do j=0, jmax
752 
753  if (maske(j,i) == 1_i2b) then ! ice-free land
754 
755  zs_neu(j,i) = zb_neu(j,i) ! this prevents zs_neu(j,i)
756  h_neu(j,i) = 0.0_dp ! from being below zb_neu(j,i)
757 
758  else if (maske(j,i) == 2_i2b) then ! sea
759 
760  zs_neu(j,i) = z_sl ! both zs_neu(j,i) and zb_neu(j,i)
761  zb_neu(j,i) = z_sl ! set to the current sea level stand
762  h_neu(j,i) = 0.0_dp
763 
764  else if (maske(j,i) == 3_i2b) then ! floating ice
765 
766  h_neu(j,i) = zs_neu(j,i)-zb_neu(j,i) ! ice thickness
767 
768  zb_neu(j,i) = z_sl - rho_rhosw_ratio*h_neu(j,i) ! floating condition
769  zs_neu(j,i) = zb_neu(j,i) + h_neu(j,i)
770 
771  end if
772 
773 end do
774 end do
775 
776 !-------- Computation of further quantities --------
777 
778 #if (ZS_EVOL==0)
779 
780 do i=0, imax
781 do j=0, jmax
782 
783  if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
784  ! grounded or floating ice
785 
786  as_perp(j,i) = as_perp(j,i) - (zs_neu(j,i)-zs(j,i))*dtime_inv
787  ! correction of the accumulation-ablation function
788  ! such that it is consistent with a steady surface
789 
790  else ! maske(j,i)==1_i2b or 2_i2b, ice-free land or sea
791 
792  as_perp(j,i) = 0.0_dp
793 
794  end if
795 
796  zs_neu(j,i) = zs(j,i) ! newly computed surface topography
797  h_neu(j,i) = zs_neu(j,i)-zb_neu(j,i) ! is discarded
798 
799 end do
800 end do
801 
802 #endif
803 
804 do i=0, imax
805 do j=0, jmax
806 
807  if ( (maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b) ) then
808  ! grounded or floating ice
809 
810 ! ------ Limit the thickness of isolated glaciated points to H_isol_max
811 
812  if ( ((maske(j,i+1) == 1_i2b).or.(maske(j,i+1) == 2_i2b)) &
813  .and. ((maske(j,i-1) == 1_i2b).or.(maske(j,i-1) == 2_i2b)) &
814  .and. ((maske(j+1,i) == 1_i2b).or.(maske(j+1,i) == 2_i2b)) &
815  .and. ((maske(j-1,i) == 1_i2b).or.(maske(j-1,i) == 2_i2b)) &
816  ) then ! all nearest neighbours ice-free
817  h_neu(j,i) = min(h_neu(j,i), h_isol_max)
818  zs_neu(j,i) = zb_neu(j,i)+h_neu(j,i)
819  end if
820 
821 ! ------ Further stuff
822 
823  if (n_cts(j,i) == 1) then
824  h_t_neu(j,i) = h_t(j,i) * h_neu(j,i)/h(j,i)
825  zm_neu(j,i) = zb_neu(j,i)+h_t_neu(j,i)
826  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
827  else
828  h_t_neu(j,i) = 0.0_dp
829  zm_neu(j,i) = zb_neu(j,i)
830  h_c_neu(j,i) = zs_neu(j,i)-zm_neu(j,i)
831  end if
832 
833  else ! maske(j,i) == 1_i2b or 2_i2b, ice-free land or sea
834 
835  zs_neu(j,i) = zb_neu(j,i)
836  zm_neu(j,i) = zb_neu(j,i)
837  h_c_neu(j,i) = 0.0_dp
838  h_t_neu(j,i) = 0.0_dp
839  h_neu(j,i) = 0.0_dp
840 
841  end if
842 
843  dzs_dtau(j,i) = (zs_neu(j,i)-zs(j,i))*dtime_inv
844  dzb_dtau(j,i) = (zb_neu(j,i)-zb(j,i))*dtime_inv
845  dzm_dtau(j,i) = dh_t_dtau(j,i)+dzb_dtau(j,i)
846  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
847 
848 #if (ZS_EVOL==2)
849  if ( abs((time-time_target_topo_final)*year_sec_inv) < eps ) then
850  dzs_dtau = 0.0_dp
851  dzb_dtau = 0.0_dp
852  dzm_dtau = 0.0_dp
853  dh_c_dtau(j,i) = 0.0_dp
854  dh_t_dtau(j,i) = 0.0_dp
855  end if
856 #endif
857 
858 end do
859 end do
860 
861 !-------- New gradients --------
862 
863 #if (TOPOGRAD==0)
864 call topograd_1(dxi, deta, 2)
865 #elif (TOPOGRAD==1)
866 call topograd_2(dxi, deta, 2)
867 #endif
868 
869 !-------- Compute the volume flux across the CTS, am_perp --------
870 
871 #if (CALCMOD==1)
872 
873 do i=1, imax-1
874 do j=1, jmax-1
875 
876  if ( ((maske(j,i) == 0_i2b).or.(maske(j,i) == 3_i2b)) &
877  .and.(n_cts(j,i) == 1_i2b)) then
878 
879  am_perp_st(j,i) = &
880  0.5_dp*(vx_c(0,j,i)+vx_c(0,j,i-1))*dzm_dxi_g(j,i) &
881  + 0.5_dp*(vy_c(0,j,i)+vy_c(0,j-1,i))*dzm_deta_g(j,i) &
882  - 0.5_dp*(vz_c(0,j,i)+vz_t(ktmax-1,j,i))
883  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
884 
885  else
886  am_perp_st(j,i) = 0.0_dp
887  am_perp(j,i) = 0.0_dp
888  end if
889 
890 end do
891 end do
892 
893 #endif
894 
895 end subroutine calc_top_3
896 !
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
calc_top_3
subroutine calc_top_3(time, z_sl, z_mar, dtime, dxi, deta, mean_accum, ii, jj, nn)
Computation of the ice topography (including gradients) for coupled SIA/SSA or SIA/SStA/SSA dynamics ...
Definition: calc_top_3.F90:37
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
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