SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_dxyz.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ d x y z
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-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 all components of the strain-rate tensor, the full
35 !! effective strain rate and the shear fraction.
36 !<------------------------------------------------------------------------------
37 subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
38 
39 use sico_types
40 use sico_variables
41 use sico_vars
42 
43 implicit none
44 
45 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t
46 
47 integer(i4b) :: i, j, kc, kt
48 real(dp) :: dxi_inv, deta_inv
49 real(dp), dimension(0:JMAX,0:IMAX) :: fact_x, fact_y
50 real(dp), dimension(0:KCMAX) :: fact_z_c
51 real(dp) :: fact_z_t
52 real(dp) :: h_c_inv, h_t_inv
53 real(dp), dimension(0:KCMAX) :: lxy_c, lyx_c, lxz_c, lzx_c, lyz_c, lzy_c
54 real(dp), dimension(0:KTMAX) :: lxy_t, lyx_t, lxz_t, lzx_t, lyz_t, lzy_t
55 real(dp), dimension(0:KCMAX) :: shear_c_squared
56 real(dp), dimension(0:KTMAX) :: shear_t_squared
57 real(dp) :: abs_v_ssa_inv, nx, ny
58 real(dp) :: shear_x_help, shear_y_help
59 real(dp) :: lambda_shear_help
60 
61 !-------- Term abbreviations --------
62 
63 dxi_inv = 1.0_dp/dxi
64 deta_inv = 1.0_dp/deta
65 
66 fact_x = dxi_inv *insq_g11_g
67 fact_y = deta_inv*insq_g22_g
68 
69 do kc=0, kcmax
70  if (flag_aa_nonzero) then
71  fact_z_c(kc) = (ea-1.0_dp)/(aa*eaz_c(kc)*dzeta_c)
72  else
73  fact_z_c(kc) = 1.0_dp/dzeta_c
74  end if
75 end do
76 
77 fact_z_t = 1.0_dp/dzeta_t
78 
79 !-------- Initialisation --------
80 
81 dxx_c = 0.0_dp
82 dyy_c = 0.0_dp
83 dxy_c = 0.0_dp
84 dxz_c = 0.0_dp
85 dyz_c = 0.0_dp
86 de_c = 0.0_dp
87 lambda_shear_c = 0.0_dp
88 
89 dxx_t = 0.0_dp
90 dyy_t = 0.0_dp
91 dxy_t = 0.0_dp
92 dxz_t = 0.0_dp
93 dyz_t = 0.0_dp
94 de_t = 0.0_dp
95 lambda_shear_t = 0.0_dp
96 
97 !-------- Computation --------
98 
99 do i=1, imax-1
100 do j=1, jmax-1
101 
102  if ((maske(j,i) == 0).or.(maske(j,i) == 3)) then ! grounded or floating ice
103 
104  h_c_inv = 1.0_dp/(abs(h_c(j,i))+epsi)
105 
106  kc=0
107 
108  dxx_c(kc,j,i) = (vx_c(kc,j,i)-vx_c(kc,j,i-1))*fact_x(j,i)
109  dyy_c(kc,j,i) = (vy_c(kc,j,i)-vy_c(kc,j-1,i))*fact_y(j,i)
110 
111  lxy_c(kc) = ( (vx_c(kc,j+1,i)+vx_c(kc,j+1,i-1)) &
112  - (vx_c(kc,j-1,i)+vx_c(kc,j-1,i-1)) ) &
113  *0.25_dp*fact_y(j,i)
114 
115  lyx_c(kc) = ( (vy_c(kc,j,i+1)+vy_c(kc,j-1,i+1)) &
116  - (vy_c(kc,j,i-1)+vy_c(kc,j-1,i-1)) ) &
117  *0.25_dp*fact_x(j,i)
118 
119  lzx_c(kc) = (vz_m(j,i+1)-vz_m(j,i-1))*0.5_dp*fact_x(j,i)
120  lzy_c(kc) = (vz_m(j+1,i)-vz_m(j-1,i))*0.5_dp*fact_y(j,i)
121 
122  lxz_c(kc) = ( (vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1)) &
123  - (vx_c(kc ,j,i)+vx_c(kc ,j,i-1)) ) &
124  *0.5_dp*fact_z_c(kc)*h_c_inv
125 
126  lyz_c(kc) = ( (vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i)) &
127  - (vy_c(kc ,j,i)+vy_c(kc ,j-1,i)) ) &
128  *0.5_dp*fact_z_c(kc)*h_c_inv
129 
130  continue ! end kc=0
131 
132  do kc=1, kcmax-1
133 
134  dxx_c(kc,j,i) = (vx_c(kc,j,i)-vx_c(kc,j,i-1))*fact_x(j,i)
135  dyy_c(kc,j,i) = (vy_c(kc,j,i)-vy_c(kc,j-1,i))*fact_y(j,i)
136 
137  lxy_c(kc) = ( (vx_c(kc,j+1,i)+vx_c(kc,j+1,i-1)) &
138  - (vx_c(kc,j-1,i)+vx_c(kc,j-1,i-1)) ) &
139  *0.25_dp*fact_y(j,i)
140 
141  lyx_c(kc) = ( (vy_c(kc,j,i+1)+vy_c(kc,j-1,i+1)) &
142  - (vy_c(kc,j,i-1)+vy_c(kc,j-1,i-1)) ) &
143  *0.25_dp*fact_x(j,i)
144 
145  lzx_c(kc) = ( (vz_c(kc,j,i+1)+vz_c(kc-1,j,i+1)) &
146  - (vz_c(kc,j,i-1)+vz_c(kc-1,j,i-1)) ) &
147  *0.25_dp*fact_x(j,i)
148 
149  lzy_c(kc) = ( (vz_c(kc,j+1,i)+vz_c(kc-1,j+1,i)) &
150  - (vz_c(kc,j-1,i)+vz_c(kc-1,j-1,i)) ) &
151  *0.25_dp*fact_y(j,i)
152 
153  lxz_c(kc) = ( (vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1)) &
154  - (vx_c(kc-1,j,i)+vx_c(kc-1,j,i-1)) ) &
155  *0.25_dp*fact_z_c(kc)*h_c_inv
156 
157  lyz_c(kc) = ( (vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i)) &
158  - (vy_c(kc-1,j,i)+vy_c(kc-1,j-1,i)) ) &
159  *0.25_dp*fact_z_c(kc)*h_c_inv
160 
161  end do
162 
163  kc=kcmax
164 
165  dxx_c(kc,j,i) = (vx_c(kc,j,i)-vx_c(kc,j,i-1))*fact_x(j,i)
166  dyy_c(kc,j,i) = (vy_c(kc,j,i)-vy_c(kc,j-1,i))*fact_y(j,i)
167 
168  lxy_c(kc) = ( (vx_c(kc,j+1,i)+vx_c(kc,j+1,i-1)) &
169  - (vx_c(kc,j-1,i)+vx_c(kc,j-1,i-1)) ) &
170  *0.25_dp*fact_y(j,i)
171 
172  lyx_c(kc) = ( (vy_c(kc,j,i+1)+vy_c(kc,j-1,i+1)) &
173  - (vy_c(kc,j,i-1)+vy_c(kc,j-1,i-1)) ) &
174  *0.25_dp*fact_x(j,i)
175 
176  lzx_c(kc) = (vz_s(j,i+1)-vz_s(j,i-1))*0.5_dp*fact_x(j,i)
177  lzy_c(kc) = (vz_s(j+1,i)-vz_s(j-1,i))*0.5_dp*fact_y(j,i)
178 
179  lxz_c(kc) = ( (vx_c(kc ,j,i)+vx_c(kc ,j,i-1)) &
180  - (vx_c(kc-1,j,i)+vx_c(kc-1,j,i-1)) ) &
181  *0.5_dp*fact_z_c(kc)*h_c_inv
182 
183  lyz_c(kc) = ( (vy_c(kc ,j,i)+vy_c(kc ,j-1,i)) &
184  - (vy_c(kc-1,j,i)+vy_c(kc-1,j-1,i)) ) &
185  *0.5_dp*fact_z_c(kc)*h_c_inv
186 
187  continue ! end kc=KCMAX
188 
189  dxy_c(:,j,i) = 0.5_dp*(lxy_c+lyx_c)
190  dxz_c(:,j,i) = 0.5_dp*(lxz_c+lzx_c)
191  dyz_c(:,j,i) = 0.5_dp*(lyz_c+lzy_c)
192 
193  do kc=0, kcmax
194 
195  shear_c_squared(kc) = dxz_c(kc,j,i)*dxz_c(kc,j,i) &
196  + dyz_c(kc,j,i)*dyz_c(kc,j,i)
197 
198  de_c(kc,j,i) = sqrt( dxx_c(kc,j,i)*dxx_c(kc,j,i) &
199  + dyy_c(kc,j,i)*dyy_c(kc,j,i) &
200  + dxx_c(kc,j,i)*dyy_c(kc,j,i) &
201  + dxy_c(kc,j,i)*dxy_c(kc,j,i) &
202  + shear_c_squared(kc) )
203 
204  lambda_shear_c(kc,j,i) = sqrt(shear_c_squared(kc))/de_c(kc,j,i)
205 
206  end do
207 
208 #if (CALCMOD==1)
209 
210  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
211  ! cold ice base, temperate ice base
212 
213  dxx_t(:,j,i) = dxx_c(0,j,i)
214  dyy_t(:,j,i) = dyy_c(0,j,i)
215  dxy_t(:,j,i) = dxy_c(0,j,i)
216  dxz_t(:,j,i) = dxz_c(0,j,i)
217  dyz_t(:,j,i) = dyz_c(0,j,i)
218  de_t(:,j,i) = de_c(0,j,i)
219  lambda_shear_t(:,j,i) = lambda_shear_c(0,j,i)
220 
221  else ! n_cts(j,i) == 1, temperate ice layer
222 
223  h_t_inv = 1.0_dp/(abs(h_t(j,i))+epsi)
224 
225  kt=0
226 
227  dxx_t(kt,j,i) = (vx_t(kt,j,i)-vx_t(kt,j,i-1))*fact_x(j,i)
228  dyy_t(kt,j,i) = (vy_t(kt,j,i)-vy_t(kt,j-1,i))*fact_y(j,i)
229 
230  lxy_t(kt) = ( (vx_t(kt,j+1,i)+vx_t(kt,j+1,i-1)) &
231  - (vx_t(kt,j-1,i)+vx_t(kt,j-1,i-1)) ) &
232  *0.25_dp*fact_y(j,i)
233 
234  lyx_t(kt) = ( (vy_t(kt,j,i+1)+vy_t(kt,j-1,i+1)) &
235  - (vy_t(kt,j,i-1)+vy_t(kt,j-1,i-1)) ) &
236  *0.25_dp*fact_x(j,i)
237 
238  lzx_t(kt) = (vz_b(j,i+1)-vz_b(j,i-1))*0.5_dp*fact_x(j,i)
239  lzy_t(kt) = (vz_b(j+1,i)-vz_b(j-1,i))*0.5_dp*fact_y(j,i)
240 
241  lxz_t(kt) = ( (vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1)) &
242  - (vx_t(kt ,j,i)+vx_t(kt ,j,i-1)) ) &
243  *0.5_dp*fact_z_t*h_t_inv
244 
245  lyz_t(kt) = ( (vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i)) &
246  - (vy_t(kt ,j,i)+vy_t(kt ,j-1,i)) ) &
247  *0.5_dp*fact_z_t*h_t_inv
248 
249  continue ! end kt=0
250 
251  do kt=1, ktmax-1
252 
253  dxx_t(kt,j,i) = (vx_t(kt,j,i)-vx_t(kt,j,i-1))*fact_x(j,i)
254  dyy_t(kt,j,i) = (vy_t(kt,j,i)-vy_t(kt,j-1,i))*fact_y(j,i)
255 
256  lxy_t(kt) = ( (vx_t(kt,j+1,i)+vx_t(kt,j+1,i-1)) &
257  - (vx_t(kt,j-1,i)+vx_t(kt,j-1,i-1)) ) &
258  *0.25_dp*fact_y(j,i)
259 
260  lyx_t(kt) = ( (vy_t(kt,j,i+1)+vy_t(kt,j-1,i+1)) &
261  - (vy_t(kt,j,i-1)+vy_t(kt,j-1,i-1)) ) &
262  *0.25_dp*fact_x(j,i)
263 
264  lzx_t(kt) = ( (vz_t(kt,j,i+1)+vz_t(kt-1,j,i+1)) &
265  - (vz_t(kt,j,i-1)+vz_t(kt-1,j,i-1)) ) &
266  *0.25_dp*fact_x(j,i)
267 
268  lzy_t(kt) = ( (vz_t(kt,j+1,i)+vz_t(kt-1,j+1,i)) &
269  - (vz_t(kt,j-1,i)+vz_t(kt-1,j-1,i)) ) &
270  *0.25_dp*fact_y(j,i)
271 
272  lxz_t(kt) = ( (vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1)) &
273  - (vx_t(kt-1,j,i)+vx_t(kt-1,j,i-1)) ) &
274  *0.25_dp*fact_z_t*h_t_inv
275 
276  lyz_t(kt) = ( (vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i)) &
277  - (vy_t(kt-1,j,i)+vy_t(kt-1,j-1,i)) ) &
278  *0.25_dp*fact_z_t*h_t_inv
279 
280  end do
281 
282  kt=ktmax
283 
284  dxx_t(kt,j,i) = (vx_t(kt,j,i)-vx_t(kt,j,i-1))*fact_x(j,i)
285  dyy_t(kt,j,i) = (vy_t(kt,j,i)-vy_t(kt,j-1,i))*fact_y(j,i)
286 
287  lxy_t(kt) = ( (vx_t(kt,j+1,i)+vx_t(kt,j+1,i-1)) &
288  - (vx_t(kt,j-1,i)+vx_t(kt,j-1,i-1)) ) &
289  *0.25_dp*fact_y(j,i)
290 
291  lyx_t(kt) = ( (vy_t(kt,j,i+1)+vy_t(kt,j-1,i+1)) &
292  - (vy_t(kt,j,i-1)+vy_t(kt,j-1,i-1)) ) &
293  *0.25_dp*fact_x(j,i)
294 
295  lzx_t(kt) = (vz_m(j,i+1)-vz_m(j,i-1))*0.5_dp*fact_x(j,i)
296  lzy_t(kt) = (vz_m(j+1,i)-vz_m(j-1,i))*0.5_dp*fact_y(j,i)
297 
298  lxz_t(kt) = ( (vx_t(kt ,j,i)+vx_t(kt ,j,i-1)) &
299  - (vx_t(kt-1,j,i)+vx_t(kt-1,j,i-1)) ) &
300  *0.5_dp*fact_z_t*h_t_inv
301 
302  lyz_t(kt) = ( (vy_t(kt ,j,i)+vy_t(kt ,j-1,i)) &
303  - (vy_t(kt-1,j,i)+vy_t(kt-1,j-1,i)) ) &
304  *0.5_dp*fact_z_t*h_t_inv
305 
306  continue ! end kt=KTMAX
307 
308  dxy_t(:,j,i) = 0.5_dp*(lxy_t+lyx_t)
309  dxz_t(:,j,i) = 0.5_dp*(lxz_t+lzx_t)
310  dyz_t(:,j,i) = 0.5_dp*(lyz_t+lzy_t)
311 
312  do kt=0, ktmax
313 
314  shear_t_squared(kt) = dxz_t(kt,j,i)*dxz_t(kt,j,i) &
315  + dyz_t(kt,j,i)*dyz_t(kt,j,i)
316 
317  de_t(kt,j,i) = sqrt( dxx_t(kt,j,i)*dxx_t(kt,j,i) &
318  + dyy_t(kt,j,i)*dyy_t(kt,j,i) &
319  + dxx_t(kt,j,i)*dyy_t(kt,j,i) &
320  + dxy_t(kt,j,i)*dxy_t(kt,j,i) &
321  + shear_t_squared(kt) )
322 
323  lambda_shear_t(kt,j,i) = sqrt(shear_t_squared(kt))/de_t(kt,j,i)
324 
325  end do
326 
327  end if
328 
329 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
330 
331  dxx_t(:,j,i) = dxx_c(0,j,i)
332  dyy_t(:,j,i) = dyy_c(0,j,i)
333  dxy_t(:,j,i) = dxy_c(0,j,i)
334  dxz_t(:,j,i) = dxz_c(0,j,i)
335  dyz_t(:,j,i) = dyz_c(0,j,i)
336  de_t(:,j,i) = de_c(0,j,i)
337  lambda_shear_t(:,j,i) = lambda_shear_c(0,j,i)
338 
339 #else
340  stop ' calc_dxyz: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
341 #endif
342 
343 ! ------ Modification of the shear fraction for floating ice (ice shelves)
344 
345  if (maske(j,i) == 3) then ! floating ice
346 
347  abs_v_ssa_inv = 1.0_dp / &
348  sqrt( 0.25_dp*(vx_m(j,i)+vx_m(j,i-1))**2 &
349  +0.25_dp*(vy_m(j,i)+vy_m(j-1,i))**2 )
350 
351  nx = -0.5_dp*(vy_m(j,i)+vy_m(j-1,i)) * abs_v_ssa_inv
352  ny = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1)) * abs_v_ssa_inv
353 
354  shear_x_help = dxx_c(kcmax,j,i)*nx &
355  + dxy_c(kcmax,j,i)*ny &
356  - dxx_c(kcmax,j,i)*nx**3 &
357  - 2.0_dp*dxy_c(kcmax,j,i)*nx**2*ny &
358  - dyy_c(kcmax,j,i)*nx*ny**2
359  ! strain rate for ice shelves independent of depth,
360  ! thus surface values used here
361 
362  shear_y_help = dxy_c(kcmax,j,i)*nx &
363  + dyy_c(kcmax,j,i)*ny &
364  - dxx_c(kcmax,j,i)*nx**2*ny &
365  - 2.0_dp*dxy_c(kcmax,j,i)*nx*ny**2 &
366  - dyy_c(kcmax,j,i)*ny**3
367  ! strain rate for ice shelves independent of depth,
368  ! thus surface values used here
369 
370  lambda_shear_help = sqrt(shear_x_help**2 + shear_y_help**2) &
371  / de_ssa(j,i)
372 
373  lambda_shear_c(:,j,i) = lambda_shear_help
374  lambda_shear_t(:,j,i) = lambda_shear_help
375 
376  end if
377 
378 ! ------ Constrain the shear fraction to reasonable [0,1] interval
379 
380  lambda_shear_c(:,j,i) = min(max(lambda_shear_c(:,j,i), 0.0_dp), 1.0_dp)
381  lambda_shear_t(:,j,i) = min(max(lambda_shear_t(:,j,i), 0.0_dp), 1.0_dp)
382 
383  else ! maske(j,i) == 1 or 2; ice-free land or ocean
384 
385  dxx_c(:,j,i) = 0.0_dp
386  dyy_c(:,j,i) = 0.0_dp
387  dxy_c(:,j,i) = 0.0_dp
388  dxz_c(:,j,i) = 0.0_dp
389  dyz_c(:,j,i) = 0.0_dp
390  de_c(:,j,i) = 0.0_dp
391  lambda_shear_c(:,j,i) = 0.0_dp
392 
393  dxx_t(:,j,i) = 0.0_dp
394  dyy_t(:,j,i) = 0.0_dp
395  dxy_t(:,j,i) = 0.0_dp
396  dxz_t(:,j,i) = 0.0_dp
397  dyz_t(:,j,i) = 0.0_dp
398  de_t(:,j,i) = 0.0_dp
399  lambda_shear_t(:,j,i) = 0.0_dp
400 
401  end if
402 
403 end do
404 end do
405 
406 end subroutine calc_dxyz
407 !
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
calc_dxyz
subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz.F90:37
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
Generated on Wed Mar 2 2016 11:26:53 for SICOPOLIS V3.2 by   doxygen 1.8.5