SICOPOLIS V5-dev  Revision 1288
calc_dxyz_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ d x y z _ m
4 !
5 !> @file
6 !!
7 !! Computation of all components of the strain-rate tensor, the full
8 !! effective strain rate and the shear fraction.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2014-2018 Ralf Greve, Fuyuki Saito
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 all components of the strain-rate tensor, the full
35 !! effective strain rate and the shear fraction.
36 !<------------------------------------------------------------------------------
38 
39  use sico_types_m
41  use sico_vars_m
42  use error_m
43 
44  implicit none
45 
46  private
47  public :: calc_dxyz
48 
49 contains
50 
51 !-------------------------------------------------------------------------------
52 !> Main subroutine of calc_dxyz_m:
53 !! Computation of all components of the strain-rate tensor, the full
54 !! effective strain rate and the shear fraction.
55 !<------------------------------------------------------------------------------
56  subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
57 
58  implicit none
59 
60  real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t
61 
62  integer(i4b) :: i, j, kc, kt
63  real(dp) :: dxi_inv, deta_inv
64  real(dp), dimension(0:JMAX,0:IMAX) :: fact_x, fact_y
65  real(dp), dimension(0:KCMAX) :: fact_z_c
66  real(dp) :: fact_z_t
67  real(dp) :: H_c_inv, H_t_inv
68  real(dp), dimension(0:KCMAX) :: lxy_c, lyx_c, lxz_c, lzx_c, lyz_c, lzy_c
69  real(dp), dimension(0:KTMAX) :: lxy_t, lyx_t, lxz_t, lzx_t, lyz_t, lzy_t
70  real(dp), dimension(0:KCMAX) :: shear_c_squared
71  real(dp), dimension(0:KTMAX) :: shear_t_squared
72  real(dp) :: abs_v_ssa_inv, nx, ny
73  real(dp) :: shear_x_help, shear_y_help
74  real(dp) :: lambda_shear_help
75 
76 !-------- Term abbreviations --------
77 
78  dxi_inv = 1.0_dp/dxi
79  deta_inv = 1.0_dp/deta
80 
81  fact_x = dxi_inv *insq_g11_g
82  fact_y = deta_inv*insq_g22_g
83 
84  do kc=0, kcmax
85  if (flag_aa_nonzero) then
86  fact_z_c(kc) = (ea-1.0_dp)/(aa*eaz_c(kc)*dzeta_c)
87  else
88  fact_z_c(kc) = 1.0_dp/dzeta_c
89  end if
90  end do
91 
92  fact_z_t = 1.0_dp/dzeta_t
93 
94 !-------- Initialisation --------
95 
96  dxx_c = 0.0_dp
97  dyy_c = 0.0_dp
98  dxy_c = 0.0_dp
99  dxz_c = 0.0_dp
100  dyz_c = 0.0_dp
101  de_c = 0.0_dp
102  lambda_shear_c = 0.0_dp
103 
104  dxx_t = 0.0_dp
105  dyy_t = 0.0_dp
106  dxy_t = 0.0_dp
107  dxz_t = 0.0_dp
108  dyz_t = 0.0_dp
109  de_t = 0.0_dp
110  lambda_shear_t = 0.0_dp
111 
112 !-------- Computation --------
113 
114  do i=1, imax-1
115  do j=1, jmax-1
116 
117  if ((maske(j,i) == 0).or.(maske(j,i) == 3)) then
118  ! grounded or floating ice
119 
120  h_c_inv = 1.0_dp/(abs(h_c(j,i))+eps_dp)
121 
122  kc=0
123 
124  dxx_c(kc,j,i) = (vx_c(kc,j,i)-vx_c(kc,j,i-1))*fact_x(j,i)
125  dyy_c(kc,j,i) = (vy_c(kc,j,i)-vy_c(kc,j-1,i))*fact_y(j,i)
126 
127  lxy_c(kc) = ( (vx_c(kc,j+1,i)+vx_c(kc,j+1,i-1)) &
128  - (vx_c(kc,j-1,i)+vx_c(kc,j-1,i-1)) ) &
129  *0.25_dp*fact_y(j,i)
130 
131  lyx_c(kc) = ( (vy_c(kc,j,i+1)+vy_c(kc,j-1,i+1)) &
132  - (vy_c(kc,j,i-1)+vy_c(kc,j-1,i-1)) ) &
133  *0.25_dp*fact_x(j,i)
134 
135  lzx_c(kc) = (vz_m(j,i+1)-vz_m(j,i-1))*0.5_dp*fact_x(j,i)
136  lzy_c(kc) = (vz_m(j+1,i)-vz_m(j-1,i))*0.5_dp*fact_y(j,i)
137 
138  lxz_c(kc) = ( (vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1)) &
139  - (vx_c(kc ,j,i)+vx_c(kc ,j,i-1)) ) &
140  *0.5_dp*fact_z_c(kc)*h_c_inv
141 
142  lyz_c(kc) = ( (vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i)) &
143  - (vy_c(kc ,j,i)+vy_c(kc ,j-1,i)) ) &
144  *0.5_dp*fact_z_c(kc)*h_c_inv
145 
146  ! end kc=0
147 
148  do kc=1, kcmax-1
149 
150  dxx_c(kc,j,i) = (vx_c(kc,j,i)-vx_c(kc,j,i-1))*fact_x(j,i)
151  dyy_c(kc,j,i) = (vy_c(kc,j,i)-vy_c(kc,j-1,i))*fact_y(j,i)
152 
153  lxy_c(kc) = ( (vx_c(kc,j+1,i)+vx_c(kc,j+1,i-1)) &
154  - (vx_c(kc,j-1,i)+vx_c(kc,j-1,i-1)) ) &
155  *0.25_dp*fact_y(j,i)
156 
157  lyx_c(kc) = ( (vy_c(kc,j,i+1)+vy_c(kc,j-1,i+1)) &
158  - (vy_c(kc,j,i-1)+vy_c(kc,j-1,i-1)) ) &
159  *0.25_dp*fact_x(j,i)
160 
161  lzx_c(kc) = ( (vz_c(kc,j,i+1)+vz_c(kc-1,j,i+1)) &
162  - (vz_c(kc,j,i-1)+vz_c(kc-1,j,i-1)) ) &
163  *0.25_dp*fact_x(j,i)
164 
165  lzy_c(kc) = ( (vz_c(kc,j+1,i)+vz_c(kc-1,j+1,i)) &
166  - (vz_c(kc,j-1,i)+vz_c(kc-1,j-1,i)) ) &
167  *0.25_dp*fact_y(j,i)
168 
169  lxz_c(kc) = ( (vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1)) &
170  - (vx_c(kc-1,j,i)+vx_c(kc-1,j,i-1)) ) &
171  *0.25_dp*fact_z_c(kc)*h_c_inv
172 
173  lyz_c(kc) = ( (vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i)) &
174  - (vy_c(kc-1,j,i)+vy_c(kc-1,j-1,i)) ) &
175  *0.25_dp*fact_z_c(kc)*h_c_inv
176 
177  end do
178 
179  kc=kcmax
180 
181  dxx_c(kc,j,i) = (vx_c(kc,j,i)-vx_c(kc,j,i-1))*fact_x(j,i)
182  dyy_c(kc,j,i) = (vy_c(kc,j,i)-vy_c(kc,j-1,i))*fact_y(j,i)
183 
184  lxy_c(kc) = ( (vx_c(kc,j+1,i)+vx_c(kc,j+1,i-1)) &
185  - (vx_c(kc,j-1,i)+vx_c(kc,j-1,i-1)) ) &
186  *0.25_dp*fact_y(j,i)
187 
188  lyx_c(kc) = ( (vy_c(kc,j,i+1)+vy_c(kc,j-1,i+1)) &
189  - (vy_c(kc,j,i-1)+vy_c(kc,j-1,i-1)) ) &
190  *0.25_dp*fact_x(j,i)
191 
192  lzx_c(kc) = (vz_s(j,i+1)-vz_s(j,i-1))*0.5_dp*fact_x(j,i)
193  lzy_c(kc) = (vz_s(j+1,i)-vz_s(j-1,i))*0.5_dp*fact_y(j,i)
194 
195  lxz_c(kc) = ( (vx_c(kc ,j,i)+vx_c(kc ,j,i-1)) &
196  - (vx_c(kc-1,j,i)+vx_c(kc-1,j,i-1)) ) &
197  *0.5_dp*fact_z_c(kc)*h_c_inv
198 
199  lyz_c(kc) = ( (vy_c(kc ,j,i)+vy_c(kc ,j-1,i)) &
200  - (vy_c(kc-1,j,i)+vy_c(kc-1,j-1,i)) ) &
201  *0.5_dp*fact_z_c(kc)*h_c_inv
202 
203  ! end kc=KCMAX
204 
205 #ifndef ALLOW_OPENAD /* Normal */
206  dxy_c(:,j,i) = 0.5_dp*(lxy_c+lyx_c)
207  dxz_c(:,j,i) = 0.5_dp*(lxz_c+lzx_c)
208  dyz_c(:,j,i) = 0.5_dp*(lyz_c+lzy_c)
209 #else /* OpenAD: cannot parse var(:) syntax */
210  do kc=0, kcmax
211  dxy_c(kc,j,i) = 0.5_dp*(lxy_c(kc)+lyx_c(kc))
212  dxz_c(kc,j,i) = 0.5_dp*(lxz_c(kc)+lzx_c(kc))
213  dyz_c(kc,j,i) = 0.5_dp*(lyz_c(kc)+lzy_c(kc))
214  end do
215 #endif /* Normal vs. OpenAD */
216 
217  do kc=0, kcmax
218 
219  shear_c_squared(kc) = dxz_c(kc,j,i)*dxz_c(kc,j,i) &
220  + dyz_c(kc,j,i)*dyz_c(kc,j,i)
221 
222  de_c(kc,j,i) = sqrt( ( dxx_c(kc,j,i)*dxx_c(kc,j,i) &
223  + dyy_c(kc,j,i)*dyy_c(kc,j,i) &
224  + dxx_c(kc,j,i)*dyy_c(kc,j,i) &
225  + dxy_c(kc,j,i)*dxy_c(kc,j,i) &
226  + shear_c_squared(kc) ) + eps_dp*eps_dp )
227 
228  lambda_shear_c(kc,j,i) = sqrt(shear_c_squared(kc)) &
229  /(de_c(kc,j,i)+eps_dp)
230 
231  end do
232 
233 #if (CALCMOD==1)
234 
235  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
236  ! cold ice base, temperate ice base
237 
238 #ifndef ALLOW_OPENAD /* Normal */
239  dxx_t(:,j,i) = dxx_c(0,j,i)
240  dyy_t(:,j,i) = dyy_c(0,j,i)
241  dxy_t(:,j,i) = dxy_c(0,j,i)
242  dxz_t(:,j,i) = dxz_c(0,j,i)
243  dyz_t(:,j,i) = dyz_c(0,j,i)
244  de_t(:,j,i) = de_c(0,j,i)
245  lambda_shear_t(:,j,i) = lambda_shear_c(0,j,i)
246 #else /* OpenAD: cannot parse var(:) syntax */
247  do kc=0, kcmax
248  dxx_t(kc,j,i) = dxx_c(0,j,i)
249  dyy_t(kc,j,i) = dyy_c(0,j,i)
250  dxy_t(kc,j,i) = dxy_c(0,j,i)
251  dxz_t(kc,j,i) = dxz_c(0,j,i)
252  dyz_t(kc,j,i) = dyz_c(0,j,i)
253  de_t(kc,j,i) = de_c(0,j,i)
254  lambda_shear_t(kc,j,i) = lambda_shear_c(0,j,i)
255  end do
256 #endif /* Normal vs. OpenAD */
257 
258  else ! n_cts(j,i) == 1, temperate ice layer
259 
260  h_t_inv = 1.0_dp/(abs(h_t(j,i))+eps_dp)
261 
262  kt=0
263 
264  dxx_t(kt,j,i) = (vx_t(kt,j,i)-vx_t(kt,j,i-1))*fact_x(j,i)
265  dyy_t(kt,j,i) = (vy_t(kt,j,i)-vy_t(kt,j-1,i))*fact_y(j,i)
266 
267  lxy_t(kt) = ( (vx_t(kt,j+1,i)+vx_t(kt,j+1,i-1)) &
268  - (vx_t(kt,j-1,i)+vx_t(kt,j-1,i-1)) ) &
269  *0.25_dp*fact_y(j,i)
270 
271  lyx_t(kt) = ( (vy_t(kt,j,i+1)+vy_t(kt,j-1,i+1)) &
272  - (vy_t(kt,j,i-1)+vy_t(kt,j-1,i-1)) ) &
273  *0.25_dp*fact_x(j,i)
274 
275  lzx_t(kt) = (vz_b(j,i+1)-vz_b(j,i-1))*0.5_dp*fact_x(j,i)
276  lzy_t(kt) = (vz_b(j+1,i)-vz_b(j-1,i))*0.5_dp*fact_y(j,i)
277 
278  lxz_t(kt) = ( (vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1)) &
279  - (vx_t(kt ,j,i)+vx_t(kt ,j,i-1)) ) &
280  *0.5_dp*fact_z_t*h_t_inv
281 
282  lyz_t(kt) = ( (vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i)) &
283  - (vy_t(kt ,j,i)+vy_t(kt ,j-1,i)) ) &
284  *0.5_dp*fact_z_t*h_t_inv
285 
286  ! end kt=0
287 
288  do kt=1, ktmax-1
289 
290  dxx_t(kt,j,i) = (vx_t(kt,j,i)-vx_t(kt,j,i-1))*fact_x(j,i)
291  dyy_t(kt,j,i) = (vy_t(kt,j,i)-vy_t(kt,j-1,i))*fact_y(j,i)
292 
293  lxy_t(kt) = ( (vx_t(kt,j+1,i)+vx_t(kt,j+1,i-1)) &
294  - (vx_t(kt,j-1,i)+vx_t(kt,j-1,i-1)) ) &
295  *0.25_dp*fact_y(j,i)
296 
297  lyx_t(kt) = ( (vy_t(kt,j,i+1)+vy_t(kt,j-1,i+1)) &
298  - (vy_t(kt,j,i-1)+vy_t(kt,j-1,i-1)) ) &
299  *0.25_dp*fact_x(j,i)
300 
301  lzx_t(kt) = ( (vz_t(kt,j,i+1)+vz_t(kt-1,j,i+1)) &
302  - (vz_t(kt,j,i-1)+vz_t(kt-1,j,i-1)) ) &
303  *0.25_dp*fact_x(j,i)
304 
305  lzy_t(kt) = ( (vz_t(kt,j+1,i)+vz_t(kt-1,j+1,i)) &
306  - (vz_t(kt,j-1,i)+vz_t(kt-1,j-1,i)) ) &
307  *0.25_dp*fact_y(j,i)
308 
309  lxz_t(kt) = ( (vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1)) &
310  - (vx_t(kt-1,j,i)+vx_t(kt-1,j,i-1)) ) &
311  *0.25_dp*fact_z_t*h_t_inv
312 
313  lyz_t(kt) = ( (vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i)) &
314  - (vy_t(kt-1,j,i)+vy_t(kt-1,j-1,i)) ) &
315  *0.25_dp*fact_z_t*h_t_inv
316 
317  end do
318 
319  kt=ktmax
320 
321  dxx_t(kt,j,i) = (vx_t(kt,j,i)-vx_t(kt,j,i-1))*fact_x(j,i)
322  dyy_t(kt,j,i) = (vy_t(kt,j,i)-vy_t(kt,j-1,i))*fact_y(j,i)
323 
324  lxy_t(kt) = ( (vx_t(kt,j+1,i)+vx_t(kt,j+1,i-1)) &
325  - (vx_t(kt,j-1,i)+vx_t(kt,j-1,i-1)) ) &
326  *0.25_dp*fact_y(j,i)
327 
328  lyx_t(kt) = ( (vy_t(kt,j,i+1)+vy_t(kt,j-1,i+1)) &
329  - (vy_t(kt,j,i-1)+vy_t(kt,j-1,i-1)) ) &
330  *0.25_dp*fact_x(j,i)
331 
332  lzx_t(kt) = (vz_m(j,i+1)-vz_m(j,i-1))*0.5_dp*fact_x(j,i)
333  lzy_t(kt) = (vz_m(j+1,i)-vz_m(j-1,i))*0.5_dp*fact_y(j,i)
334 
335  lxz_t(kt) = ( (vx_t(kt ,j,i)+vx_t(kt ,j,i-1)) &
336  - (vx_t(kt-1,j,i)+vx_t(kt-1,j,i-1)) ) &
337  *0.5_dp*fact_z_t*h_t_inv
338 
339  lyz_t(kt) = ( (vy_t(kt ,j,i)+vy_t(kt ,j-1,i)) &
340  - (vy_t(kt-1,j,i)+vy_t(kt-1,j-1,i)) ) &
341  *0.5_dp*fact_z_t*h_t_inv
342 
343  ! end kt=KTMAX
344 
345 #ifndef ALLOW_OPENAD /* Normal */
346  dxy_t(:,j,i) = 0.5_dp*(lxy_t+lyx_t)
347  dxz_t(:,j,i) = 0.5_dp*(lxz_t+lzx_t)
348  dyz_t(:,j,i) = 0.5_dp*(lyz_t+lzy_t)
349 #else /* OpenAD: cannot parse var(:) syntax */
350  do kt=0, ktmax
351  dxy_t(kt,j,i) = 0.5_dp*(lxy_t(kt)+lyx_t(kt))
352  dxz_t(kt,j,i) = 0.5_dp*(lxz_t(kt)+lzx_t(kt))
353  dyz_t(kt,j,i) = 0.5_dp*(lyz_t(kt)+lzy_t(kt))
354  end do
355 #endif /* Normal vs. OpenAD */
356 
357  do kt=0, ktmax
358 
359  shear_t_squared(kt) = dxz_t(kt,j,i)*dxz_t(kt,j,i) &
360  + dyz_t(kt,j,i)*dyz_t(kt,j,i)
361 
362  de_t(kt,j,i) = sqrt( ( dxx_t(kt,j,i)*dxx_t(kt,j,i) &
363  + dyy_t(kt,j,i)*dyy_t(kt,j,i) &
364  + dxx_t(kt,j,i)*dyy_t(kt,j,i) &
365  + dxy_t(kt,j,i)*dxy_t(kt,j,i) &
366  + shear_t_squared(kt) ) + eps_dp*eps_dp )
367 
368  lambda_shear_t(kt,j,i) = sqrt(shear_t_squared(kt)) &
369  /(de_t(kt,j,i)+eps_dp)
370 
371  end do
372 
373  end if
374 
375 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
376 
377 #ifndef ALLOW_OPENAD /* Normal */
378  dxx_t(:,j,i) = dxx_c(0,j,i)
379  dyy_t(:,j,i) = dyy_c(0,j,i)
380  dxy_t(:,j,i) = dxy_c(0,j,i)
381  dxz_t(:,j,i) = dxz_c(0,j,i)
382  dyz_t(:,j,i) = dyz_c(0,j,i)
383  de_t(:,j,i) = de_c(0,j,i)
384  lambda_shear_t(:,j,i) = lambda_shear_c(0,j,i)
385 #else /* OpenAD: cannot parse var(:) syntax */
386  do kt=0, ktmax
387  dxx_t(kt,j,i) = dxx_c(0,j,i)
388  dyy_t(kt,j,i) = dyy_c(0,j,i)
389  dxy_t(kt,j,i) = dxy_c(0,j,i)
390  dxz_t(kt,j,i) = dxz_c(0,j,i)
391  dyz_t(kt,j,i) = dyz_c(0,j,i)
392  de_t(kt,j,i) = de_c(0,j,i)
393  lambda_shear_t(kt,j,i) = lambda_shear_c(0,j,i)
394  end do
395 #endif /* Normal vs. OpenAD */
396 
397 #else
398  errormsg = ' >>> calc_dxyz: Parameter CALCMOD must be -1, 0, 1, 2 or 3!'
399  call error(errormsg)
400 #endif
401 
402 ! ------ Modification of the shear fraction for floating ice (ice shelves)
403 
404  if (maske(j,i) == 3) then ! floating ice
405 
406  abs_v_ssa_inv = 1.0_dp / &
407  sqrt( ( 0.25_dp*(vx_m(j,i)+vx_m(j,i-1))**2 &
408  +0.25_dp*(vy_m(j,i)+vy_m(j-1,i))**2 ) &
409  + eps_dp**2 )
410 
411  nx = -0.5_dp*(vy_m(j,i)+vy_m(j-1,i)) * abs_v_ssa_inv
412  ny = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1)) * abs_v_ssa_inv
413 
414  shear_x_help = dxx_c(kcmax,j,i)*nx &
415  + dxy_c(kcmax,j,i)*ny &
416  - dxx_c(kcmax,j,i)*nx**3 &
417  - 2.0_dp*dxy_c(kcmax,j,i)*nx**2*ny &
418  - dyy_c(kcmax,j,i)*nx*ny**2
419  ! strain rate for ice shelves independent of depth,
420  ! thus surface values used here
421 
422  shear_y_help = dxy_c(kcmax,j,i)*nx &
423  + dyy_c(kcmax,j,i)*ny &
424  - dxx_c(kcmax,j,i)*nx**2*ny &
425  - 2.0_dp*dxy_c(kcmax,j,i)*nx*ny**2 &
426  - dyy_c(kcmax,j,i)*ny**3
427  ! strain rate for ice shelves independent of depth,
428  ! thus surface values used here
429 
430  lambda_shear_help = sqrt( ( shear_x_help**2 + shear_y_help**2 ) &
431  + eps_dp**2 ) &
432  / (de_ssa(j,i)+eps_dp)
433 
434 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
435 
436  lambda_shear_c(:,j,i) = lambda_shear_help
437  lambda_shear_t(:,j,i) = lambda_shear_help
438 
439 #else /* OpenAD: cannot parse var(:) syntax */
440 
441  do kc=0, kcmax
442  lambda_shear_c(kc,j,i) = lambda_shear_help
443  end do
444 
445  do kt=0, ktmax
446  lambda_shear_t(kt,j,i) = lambda_shear_help
447  end do
448 
449 #endif /* Normal vs. OpenAD */
450 
451  end if
452 
453 ! ------ Constrain the shear fraction to reasonable [0,1] interval
454 
455 #ifndef ALLOW_OPENAD /* Normal */
456 
457  lambda_shear_c(:,j,i) = min(max(lambda_shear_c(:,j,i), 0.0_dp), 1.0_dp)
458  lambda_shear_t(:,j,i) = min(max(lambda_shear_t(:,j,i), 0.0_dp), 1.0_dp)
459 
460 #else /* OpenAD: cannot parse var(:) syntax */
461 
462  do kc=0, kcmax
463  lambda_shear_c(kc,j,i) = min(max(lambda_shear_c(kc,j,i),0.0_dp), &
464  1.0_dp)
465  end do
466 
467  do kt=0, ktmax
468  lambda_shear_t(kt,j,i) = min(max(lambda_shear_t(kt,j,i),0.0_dp), &
469  1.0_dp)
470  end do
471 
472 #endif /* Normal vs. OpenAD */
473 
474  else ! maske(j,i) == 1 or 2; ice-free land or ocean
475 
476 #ifndef ALLOW_OPENAD /* Normal */
477 
478  dxx_c(:,j,i) = 0.0_dp
479  dyy_c(:,j,i) = 0.0_dp
480  dxy_c(:,j,i) = 0.0_dp
481  dxz_c(:,j,i) = 0.0_dp
482  dyz_c(:,j,i) = 0.0_dp
483  de_c(:,j,i) = 0.0_dp
484  lambda_shear_c(:,j,i) = 0.0_dp
485 
486  dxx_t(:,j,i) = 0.0_dp
487  dyy_t(:,j,i) = 0.0_dp
488  dxy_t(:,j,i) = 0.0_dp
489  dxz_t(:,j,i) = 0.0_dp
490  dyz_t(:,j,i) = 0.0_dp
491  de_t(:,j,i) = 0.0_dp
492  lambda_shear_t(:,j,i) = 0.0_dp
493 
494 #else /* OpenAD: cannot parse var(:) syntax */
495 
496  do kc=0, kcmax
497  dxx_c(kc,j,i) = 0.0_dp
498  dyy_c(kc,j,i) = 0.0_dp
499  dxy_c(kc,j,i) = 0.0_dp
500  dxz_c(kc,j,i) = 0.0_dp
501  dyz_c(kc,j,i) = 0.0_dp
502  de_c(kc,j,i) = 0.0_dp
503  lambda_shear_c(kc,j,i) = 0.0_dp
504  end do
505 
506  do kt=0, ktmax
507  dxx_t(kt,j,i) = 0.0_dp
508  dyy_t(kt,j,i) = 0.0_dp
509  dxy_t(kt,j,i) = 0.0_dp
510  dxz_t(kt,j,i) = 0.0_dp
511  dyz_t(kt,j,i) = 0.0_dp
512  de_t(kt,j,i) = 0.0_dp
513  lambda_shear_t(kt,j,i) = 0.0_dp
514  end do
515 
516 #endif /* Normal vs. OpenAD */
517 
518  end if
519 
520  end do
521  end do
522 
523  end subroutine calc_dxyz
524 
525 !-------------------------------------------------------------------------------
526 
527 end module calc_dxyz_m
528 !
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)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) de_c
de_c(kc,j,i): Full effective strain rate in the upper (kc) ice domain
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:ktmax, 0:jmax, 0:imax) lambda_shear_t
lambda_shear_t(kt,j,i): Shear fraction in the lower (kt) ice domain
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dxy_t
dxy_t(kt,j,i): Strain rate dxy in the lower (kt) ice domain
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) de_t
de_t(kt,j,i): Full effective strain rate in the lower (kt) ice domain
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dyy_t
dyy_t(kt,j,i): Strain rate dyy in the lower (kt) ice domain
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 error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
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)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dyz_c
dyz_c(kc,j,i): Strain rate dyz in the upper (kc) ice domain
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 ...
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz_m.F90:37
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (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)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) lambda_shear_c
lambda_shear_c(kc,j,i): Shear fraction in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dyz_t
dyz_t(kt,j,i): Strain rate dyz in the lower (kt) ice domain
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dxz_t
dxz_t(kt,j,i): Strain rate dxz in the lower (kt) ice domain
Declarations of kind types for SICOPOLIS.
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...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dxx_c
dxx_c(kc,j,i): Strain rate dxx in the upper (kc) ice domain
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dxx_t
dxx_t(kt,j,i): Strain rate dxx in the lower (kt) ice domain
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:kcmax, 0:jmax, 0:imax) dyy_c
dyy_c(kc,j,i): Strain rate dyy in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dxz_c
dxz_c(kc,j,i): Strain rate dxz in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dxy_c
dxy_c(kc,j,i): Strain rate dxy in the upper (kc) ice domain
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) 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)) ...
subroutine, public calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Main subroutine of calc_dxyz_m: Computation of all components of the strain-rate tensor, the full effective strain rate and the shear fraction.
Definition: calc_dxyz_m.F90:57
character(len=256) errormsg
errormsg: Error-message string
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)) ...