SICOPOLIS V3.3
calc_vz_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ v z _ m
4 !
5 !> @file
6 !!
7 !! Computation of the vertical velocity vz.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 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 vertical velocity vz.
34 !<------------------------------------------------------------------------------
35 module calc_vz_m
36 
37  use sico_types_m
39  use sico_vars_m
40 
41  implicit none
42 
43  private
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Computation of the vertical velocity vz for grounded ice.
50 !<------------------------------------------------------------------------------
51 subroutine calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
52 
53 implicit none
54 
55 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t
56 
57 integer(i4b) :: i, j, kc, kt
58 real(dp) :: avz3(0:kcmax)
59 real(dp) :: cvz00, cvz0(0:ktmax), cvz1(0:ktmax), cvz2(0:ktmax), &
60  cvz3(0:kcmax), cvz4(0:kcmax), cvz5(0:kcmax)
61 real(dp) :: dxi_inv, deta_inv
62 
63 !-------- Term abbreviations --------
64 
65 dxi_inv = 1.0_dp/dxi
66 deta_inv = 1.0_dp/deta
67 
68 do kc=0, kcmax
69  if (flag_aa_nonzero) then
70  avz3(kc) = aa*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
71  else
72  avz3(kc) = dzeta_c
73  end if
74 end do
75 
76 do i=1, imax-1
77 do j=1, jmax-1
78 
79  if (maske(j,i)==0) then ! grounded ice
80 
81 !-------- Abbreviations --------
82 
83  cvz00 = 0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1))*dzb_dxi_g(j,i) &
84  +0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i))*dzb_deta_g(j,i) &
85  +dzb_dtau(j,i)-q_b_tot(j,i)
86 
87  kt=0
88  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
89  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
90  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
91  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
92  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
93  *dzeta_t
94  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
95  *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
96  -vx_t(kt,j,i) -vx_t(kt,j,i-1))*0.5_dp
97  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
98  *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
99  -vy_t(kt,j,i) -vy_t(kt,j-1,i))*0.5_dp
100 
101  do kt=1, ktmax-1
102  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
103  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
104  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
105  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
106  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
107  *dzeta_t
108  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
109  *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
110  -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.25_dp
111  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
112  *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
113  -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.25_dp
114  end do
115 
116  kt=ktmax
117  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
118  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
119  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
120  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
121  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
122  *dzeta_t
123  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
124  *(vx_t(kt,j,i) +vx_t(kt,j,i-1) &
125  -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.5_dp
126  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
127  *(vy_t(kt,j,i) +vy_t(kt,j-1,i) &
128  -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.5_dp
129 
130  kc=0
131  cvz3(kc) = avz3(kc)*h_c(j,i) &
132  *insq_g11_g(j,i)*insq_g22_g(j,i) &
133  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
134  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
135  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
136  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
137  cvz4(kc) = (dzm_dxi_g(j,i) &
138  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
139  *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
140  -vx_c(kc,j,i) -vx_c(kc,j,i-1))*0.5_dp
141  cvz5(kc) = (dzm_deta_g(j,i) &
142  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
143  *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
144  -vy_c(kc,j,i) -vy_c(kc,j-1,i))*0.5_dp
145 
146  do kc=1, kcmax-1
147  cvz3(kc) = avz3(kc)*h_c(j,i) &
148  *insq_g11_g(j,i)*insq_g22_g(j,i) &
149  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
150  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
151  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
152  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
153  cvz4(kc) = (dzm_dxi_g(j,i) &
154  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
155  *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
156  -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.25_dp
157  cvz5(kc) = (dzm_deta_g(j,i) &
158  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
159  *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
160  -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.25_dp
161  end do
162 
163  kc=kcmax
164  cvz3(kc) = avz3(kc)*h_c(j,i) &
165  *insq_g11_g(j,i)*insq_g22_g(j,i) &
166  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
167  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
168  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
169  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
170  cvz4(kc) = (dzm_dxi_g(j,i) &
171  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
172  *(vx_c(kc,j,i) +vx_c(kc,j,i-1) &
173  -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.5_dp
174  cvz5(kc) = (dzm_deta_g(j,i) &
175  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
176  *(vy_c(kc,j,i) +vy_c(kc,j-1,i) &
177  -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.5_dp
178 
179 !-------- Computation of vz_b --------
180 
181  vz_b(j,i) = cvz00
182 
183 !-------- Computation of vz --------
184 
185  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
186  ! cold ice base, temperate ice base
187 
188  do kt=0, ktmax-1
189  vz_t(kt,j,i) = vz_b(j,i)
190  end do
191 
192  vz_m(j,i) = vz_b(j,i)
193 
194  vz_c(0,j,i) = vz_m(j,i) &
195  +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))
196 
197  do kc=1, kcmax-1
198  vz_c(kc,j,i) = vz_c(kc-1,j,i) &
199  -cvz3(kc)+cvz4(kc)+cvz5(kc)
200  end do
201 
202  else ! n_cts(j,i) == 1, temperate ice layer
203 
204  vz_t(0,j,i) = vz_b(j,i) &
205  +0.5_dp*(-cvz0(0)+cvz1(0)+cvz2(0))
206 
207  do kt=1, ktmax-1
208  vz_t(kt,j,i) = vz_t(kt-1,j,i) &
209  -cvz0(kt)+cvz1(kt)+cvz2(kt)
210  end do
211 
212  vz_m(j,i) = vz_t(ktmax-1,j,i) &
213  +0.5_dp*(-cvz0(ktmax)+cvz1(ktmax)+cvz2(ktmax))
214 
215  vz_c(0,j,i) = vz_m(j,i) &
216  +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))
217 
218  do kc=1, kcmax-1
219  vz_c(kc,j,i) = vz_c(kc-1,j,i) &
220  -cvz3(kc)+cvz4(kc)+cvz5(kc)
221  end do
222 
223  end if
224 
225 !-------- Computation of vz_s --------
226 
227  kc=kcmax
228 
229  vz_s(j,i) = vz_c(kc-1,j,i) &
230  +0.5_dp*(-cvz3(kc)+cvz4(kc)+cvz5(kc))
231 
232  else ! maske(j,i) /= 0 (not grounded ice)
233 
234  vz_b(j,i) = 0.0_dp
235 
236  do kt=0, ktmax-1
237  vz_t(kt,j,i) = 0.0_dp
238  end do
239 
240  vz_m(j,i) = 0.0_dp
241 
242  do kc=0, kcmax-1
243  vz_c(kc,j,i) = 0.0_dp
244  end do
245 
246  vz_s(j,i) = 0.0_dp
247 
248  end if
249 
250 end do
251 end do
252 
253 end subroutine calc_vz_grounded
254 
255 !-------------------------------------------------------------------------------
256 !> Computation of the vertical velocity vz for floating ice.
257 !<------------------------------------------------------------------------------
258 subroutine calc_vz_floating(z_sl, dxi, deta, dzeta_c)
260 implicit none
261 
262 real(dp), intent(in) :: z_sl
263 real(dp), intent(in) :: dxi, deta, dzeta_c
264 
265 integer(i4b) :: i, j, kt, kc
266 real(dp), dimension(0:JMAX,0:IMAX) :: vz_sl
267 real(dp), dimension(0:KCMAX) :: zeta_c_sgz, eaz_c_sgz, eaz_c_quotient_sgz
268 real(dp) :: dxi_inv, deta_inv
269 real(dp) :: dvx_dxi, dvy_deta
270 
271 dxi_inv = 1.0_dp/dxi
272 deta_inv = 1.0_dp/deta
273 
274 !-------- Abbreviations --------
275 
276 do kc=0, kcmax-1
277 
278  zeta_c_sgz(kc) = (kc+0.5_dp)*dzeta_c
279 
280  if (flag_aa_nonzero) then
281  eaz_c_sgz(kc) = exp(aa*zeta_c_sgz(kc))
282  eaz_c_quotient_sgz(kc) = (eaz_c_sgz(kc)-1.0_dp)/(ea-1.0_dp)
283  else
284  eaz_c_sgz(kc) = 1.0_dp
285  eaz_c_quotient_sgz(kc) = zeta_c_sgz(kc)
286  end if
287 
288 end do
289 
290 !-------- Computation of vz --------
291 
292 do i=1, imax-1
293 do j=1, jmax-1
294 
295  if (maske(j,i)==3_i2b) then ! floating ice
296 
297 ! ------ Derivatives of the horizontal velocity
298 
299  dvx_dxi = insq_g11_g(j,i)*insq_g22_g(j,i) &
300  *(vx_m(j,i)*sq_g22_sgx(j,i)-vx_m(j,i-1)*sq_g22_sgx(j,i-1)) &
301  *dxi_inv
302 
303  dvy_deta = insq_g11_g(j,i)*insq_g22_g(j,i) &
304  *(vy_m(j,i)*sq_g11_sgy(j,i)-vy_m(j-1,i)*sq_g11_sgy(j-1,i)) &
305  *deta_inv
306 
307 ! ------ Basal velocity vz_b
308 
309  vz_b(j,i) = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1))*dzb_dxi_g(j,i) &
310  +0.5_dp*(vy_m(j,i)+vy_m(j-1,i))*dzb_deta_g(j,i) &
311  +dzb_dtau(j,i)-q_b_tot(j,i)
312  ! kinematic boundary condition at the ice base
313 
314 ! ------ Velocity at sea level vz_sl
315 
316  vz_sl(j,i) = vz_b(j,i) - (z_sl-zb(j,i))*(dvx_dxi+dvy_deta)
317 
318 ! ------ Surface velocity vz_s
319 
320  vz_s(j,i) = vz_sl(j,i) - (zs(j,i)-z_sl)*(dvx_dxi+dvy_deta)
321 
322 ! ------ Velocity vz_m at the interface between
323 ! the upper (kc) and the lower (kt) domain
324 
325  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
326  ! cold ice base, temperate ice base
327 
328  vz_m(j,i) = vz_b(j,i)
329 
330  else ! n_cts(j,i) == 1, temperate ice layer
331 
332  vz_m(j,i) = vz_b(j,i) - h_t(j,i)*(dvx_dxi+dvy_deta)
333 
334  end if
335 
336 ! ------ 3-D velocity vz_c and vz_t
337 
338  do kc=0, kcmax-1
339  vz_c(kc,j,i) = vz_sl(j,i) &
340  -(zm(j,i)+eaz_c_quotient_sgz(kc)*h_c(j,i)-z_sl) &
341  *(dvx_dxi+dvy_deta)
342  end do
343 
344  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
345  ! cold ice base, temperate ice base
346 
347  do kt=0, ktmax-1
348  vz_t(kt,j,i) = vz_b(j,i)
349  end do
350 
351  else ! n_cts(j,i) == 1, temperate ice layer
352 
353  do kt=0, ktmax-1
354  vz_t(kt,j,i) = vz_sl(j,i) &
355  -(zb(j,i) &
356  +0.5_dp*(zeta_t(kt)+zeta_t(kt+1))*h_t(j,i) &
357  -z_sl) &
358  *(dvx_dxi+dvy_deta)
359  end do
360 
361  end if
362 
363  end if
364 
365 end do
366 end do
367 
368 end subroutine calc_vz_floating
369 
370 !-------------------------------------------------------------------------------
371 !> Computation of the vertical velocity vz for static ice.
372 !<------------------------------------------------------------------------------
373 subroutine calc_vz_static()
375 implicit none
376 
377 integer(i4b) :: i, j, kc, kt
378 
379 do i=0, imax
380 do j=0, jmax
381 
382  if ((maske(j,i)==0).or.(maske(j,i)==3)) then ! grounded or floating ice
383 
384  vz_b(j,i) = dzb_dtau(j,i)-q_b_tot(j,i)
385  ! kinematic boundary condition at the ice base
386 
387  do kt=0, ktmax-1
388  vz_t(kt,j,i) = vz_b(j,i)
389  end do
390 
391  vz_m(j,i) = vz_b(j,i)
392 
393  do kc=0, kcmax-1
394  vz_c(kc,j,i) = vz_b(j,i)
395  end do
396 
397  vz_s(j,i) = vz_b(j,i)
398 
399  else ! maske(j,i) == (1 or 2)
400 
401  vz_b(j,i) = 0.0_dp
402 
403  do kt=0, ktmax-1
404  vz_t(kt,j,i) = 0.0_dp
405  end do
406 
407  vz_m(j,i) = 0.0_dp
408 
409  do kc=0, kcmax-1
410  vz_c(kc,j,i) = 0.0_dp
411  end do
412 
413  vz_s(j,i) = 0.0_dp
414 
415  end if
416 
417 end do
418 end do
419 
420 end subroutine calc_vz_static
421 
422 !-------------------------------------------------------------------------------
423 
424 end module calc_vz_m
425 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:jmax, 0:imax) vz_m
vz_m(j,i): Velocity in z-direction at the position z=zm (interface between the upper (kc) and the low...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) vx_m
vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
subroutine, public calc_vz_floating(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz for floating ice.
Definition: calc_vz_m.F90:259
real(dp), dimension(0:jmax, 0:imax) vz_b
vz_b(j,i): Velocity in z-direction at the ice base, at (i,j)
real(dp), dimension(0:jmax, 0:imax) vz_s
vz_s(j,i): Velocity in z-direction at the ice surface, at (i,j)
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) sq_g22_sgx
sq_g22_sgx(j,i): Square root of g22, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
real(dp), dimension(0:jmax, 0:imax) vy_m
vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp) ea
ea: Abbreviation for exp(aa)
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp), dimension(0:jmax, 0:imax) dh_c_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
real(dp), dimension(0:jmax, 0:imax) sq_g11_sgy
sq_g11_sgy(j,i): Square root of g11, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) dzb_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
subroutine, public calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_m.F90:374
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(0:jmax, 0:imax) dh_t_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
real(dp), dimension(0:jmax, 0:imax) dzb_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)