SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_temp2.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p 2
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age for an ice column with a temperate base
8 !! overlain by cold ice.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-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 temperature and age for an ice column with a temperate base
35 !! overlain by cold ice.
36 !<------------------------------------------------------------------------------
37 subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
38  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
39  ai1, ai2, &
40  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
41 
42 use sico_types
43 use sico_variables
44 use sico_vars
45 use ice_material_quantities, only : ratefac_c, kappa_val, c_val
46 use sico_sle_solvers, only : tri_sle
47 
48 implicit none
49 
50 integer(i4b), intent(in) :: i, j
51 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
52  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
53  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
54  ai1(0:kcmax), ai2(0:kcmax), &
55  atr1, alb1
56 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
57 
58 integer(i4b) :: kc, kt, kr
59 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
60  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
61 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
62  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
63 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
64 real(dp) :: temp_c_help(0:kcmax)
65 real(dp) :: vx_c_help, vy_c_help
66 real(dp) :: adv_vert_help
67 real(dp) :: dtt_dxi, dtt_deta
68 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
69  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
70  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
71  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
72  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
73 real(dp), parameter :: zero=0.0_dp
74 
75 !-------- Check for boundary points --------
76 
77 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
78  stop ' calc_temp2: Boundary points not allowed.'
79 
80 !-------- Abbreviations --------
81 
82 ctr1 = atr1
83 clb1 = alb1*q_geo(j,i)
84 
85 #if (ADV_VERT==1)
86 
87 do kc=1, kcmax-1
88  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
89 end do
90 
91 kc=0
92 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
93  ! only needed for kc=0 ...
94 kc=kcmax-1
95 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
96  ! ... and kc=KCMAX-1
97 
98 #elif (ADV_VERT==2 || ADV_VERT==3)
99 
100 do kc=0, kcmax-1
101  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
102 end do
103 
104 #endif
105 
106 do kc=0, kcmax
107 
108  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
109  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
110  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
111  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
112  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
113  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
114  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
115  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
116  ct5(kc) = at5(kc) &
117  /c_val(temp_c(kc,j,i)) &
118  /h_c(j,i)
119 
120 #if (DYNAMICS==2)
121  if (.not.flag_shelfy_stream(j,i)) then
122 #endif
123  ct7(kc) = at7 &
124  /c_val(temp_c(kc,j,i)) &
125  *enh_c(kc,j,i) &
126  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
127  *creep(sigma_c(kc,j,i)) &
128  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
129 #if (DYNAMICS==2)
130  else
131  ct7(kc) = 2.0_dp*at7 &
132  /c_val(temp_c(kc,j,i)) &
133  *viscosity(de_c(kc,j,i), &
134  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
135  enh_c(kc,j,i), 0_i2b) &
136  *de_c(kc,j,i)**2
137  end if
138 #endif
139 
140  ci1(kc) = ai1(kc)/h_c(j,i)
141 
142 end do
143 
144 #if (ADV_VERT==1)
145 
146 kc=0
147 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
148 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
149 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
150 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
151 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
152 kc=kcmax-1
153 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
154 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
155 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
156 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
157 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
158 
159 #elif (ADV_VERT==2 || ADV_VERT==3)
160 
161 do kc=0, kcmax-1
162  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
163  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
164  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
165  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
166  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
167 end do
168 
169 #endif
170 
171 do kc=0, kcmax-1
172  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
173  ct6(kc) = at6(kc) &
174  *kappa_val(temp_c_help(kc)) &
175  /h_c(j,i)
176  ci2(kc) = ai2(kc)/h_c(j,i)
177 end do
178 
179 #if (ADV_HOR==3)
180 dtt_dxi = 2.0_dp*dtt_2dxi
181 dtt_deta = 2.0_dp*dtt_2deta
182 #endif
183 
184 !-------- Set up the equations for the bedrock temperature --------
185 
186 kr=0
187 lgs_a1(kr) = 1.0_dp
188 lgs_a2(kr) = -1.0_dp
189 lgs_b(kr) = clb1
190 
191 #if (Q_LITHO==1)
192 ! (coupled heat-conducting bedrock)
193 
194 do kr=1, krmax-1
195  lgs_a0(kr) = - ctr1
196  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
197  lgs_a2(kr) = - ctr1
198  lgs_b(kr) = temp_r(kr,j,i)
199 end do
200 
201 #elif (Q_LITHO==0)
202 ! (no coupled heat-conducting bedrock)
203 
204 do kr=1, krmax-1
205  lgs_a0(kr) = 1.0_dp
206  lgs_a1(kr) = 0.0_dp
207  lgs_a2(kr) = -1.0_dp
208  lgs_b(kr) = 2.0_dp*clb1
209 end do
210 
211 #endif
212 
213 kr=krmax
214 lgs_a0(kr) = 0.0_dp
215 lgs_a1(kr) = 1.0_dp
216 lgs_b(kr) = temp_t_m(0,j,i)
217 
218 !-------- Solve system of linear equations --------
219 
220 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
221 
222 !-------- Assign the result --------
223 
224 do kr=0, krmax
225  temp_r_neu(kr,j,i) = lgs_x(kr)
226 end do
227 
228 !-------- Set up the equations for the ice temperature --------
229 
230 kc=0
231 lgs_a1(kc) = 1.0_dp
232 lgs_a2(kc) = 0.0_dp
233 lgs_b(kc) = temp_c_m(0,j,i)
234 
235 do kc=1, kcmax-1
236 
237 #if (ADV_VERT==1)
238 
239  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
240  -ct5(kc)*ct6(kc-1)
241  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
242  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
243  -ct5(kc)*ct6(kc)
244 
245 #elif (ADV_VERT==2)
246 
247  lgs_a0(kc) &
248  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
249  -ct5(kc)*ct6(kc-1)
250  lgs_a1(kc) &
251  = 1.0_dp &
252  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
253  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
254  +ct5(kc)*(ct6(kc)+ct6(kc-1))
255  lgs_a2(kc) &
256  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
257  -ct5(kc)*ct6(kc)
258 
259 #elif (ADV_VERT==3)
260 
261  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
262 
263  lgs_a0(kc) &
264  = -max(adv_vert_help, 0.0_dp) &
265  -ct5(kc)*ct6(kc-1)
266  lgs_a1(kc) &
267  = 1.0_dp &
268  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
269  +ct5(kc)*(ct6(kc)+ct6(kc-1))
270  lgs_a2(kc) &
271  = min(adv_vert_help, 0.0_dp) &
272  -ct5(kc)*ct6(kc)
273 
274 #endif
275 
276 #if (ADV_HOR==2)
277 
278  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
279  -dtt_2dxi* &
280  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
281  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
282  *insq_g11_sgx(j,i) &
283  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
284  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
285  *insq_g11_sgx(j,i-1) ) &
286  -dtt_2deta* &
287  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
288  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
289  *insq_g22_sgy(j,i) &
290  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
291  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
292  *insq_g22_sgy(j-1,i) )
293 
294 #elif (ADV_HOR==3)
295 
296  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
297  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
298 
299  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
300  -dtt_dxi* &
301  ( min(vx_c_help, 0.0_dp) &
302  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
303  *insq_g11_sgx(j,i) &
304  +max(vx_c_help, 0.0_dp) &
305  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
306  *insq_g11_sgx(j,i-1) ) &
307  -dtt_deta* &
308  ( min(vy_c_help, 0.0_dp) &
309  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
310  *insq_g22_sgy(j,i) &
311  +max(vy_c_help, 0.0_dp) &
312  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
313  *insq_g22_sgy(j-1,i) )
314 
315 #endif
316 
317 end do
318 
319 kc=kcmax
320 lgs_a0(kc) = 0.0_dp
321 lgs_a1(kc) = 1.0_dp
322 lgs_b(kc) = temp_s(j,i)
323 
324 !-------- Solve system of linear equations --------
325 
326 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
327 
328 !-------- Assign the result --------
329 
330 do kc=0, kcmax
331  temp_c_neu(kc,j,i) = lgs_x(kc)
332 end do
333 
334 !-------- Set water content in the non-existing temperate layer
335 ! to zero --------
336 
337 do kt=0, ktmax
338  omega_t_neu(kt,j,i) = 0.0_dp
339 end do
340 
341 !-------- Water drainage from the non-existing temperate layer --------
342 
343 q_tld(j,i) = 0.0_dp
344 
345 !-------- Set up the equations for the age of cold ice --------
346 
347 kc=0 ! adv_vert_sg(0) <= 0
348 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
349 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
350 
351 #if (ADV_HOR==2)
352 
353 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
354  -dtt_2dxi* &
355  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
356  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
357  *insq_g11_sgx(j,i) &
358  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
359  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
360  *insq_g11_sgx(j,i-1) ) &
361  -dtt_2deta* &
362  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
363  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
364  *insq_g22_sgy(j,i) &
365  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
366  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
367  *insq_g22_sgy(j-1,i) )
368 
369 #elif (ADV_HOR==3)
370 
371 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
372 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
373 
374 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
375  -dtt_dxi* &
376  ( min(vx_c_help, 0.0_dp) &
377  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
378  *insq_g11_sgx(j,i) &
379  +max(vx_c_help, 0.0_dp) &
380  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
381  *insq_g11_sgx(j,i-1) ) &
382  -dtt_deta* &
383  ( min(vy_c_help, 0.0_dp) &
384  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
385  *insq_g22_sgy(j,i) &
386  +max(vy_c_help, 0.0_dp) &
387  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
388  *insq_g22_sgy(j-1,i) )
389 
390 #endif
391 
392 do kc=1, kcmax-1
393 
394 #if (ADV_VERT==1)
395 
396  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
397  -ci1(kc)*ci2(kc-1)
398  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
399  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
400  -ci1(kc)*ci2(kc)
401 
402 #elif (ADV_VERT==2)
403 
404  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
405  lgs_a1(kc) = 1.0_dp &
406  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
407  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
408  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
409 
410 #elif (ADV_VERT==3)
411 
412  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
413 
414  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
415  lgs_a1(kc) = 1.0_dp &
416  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
417  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
418 
419 #endif
420 
421 #if (ADV_HOR==2)
422 
423  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
424  -dtt_2dxi* &
425  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
426  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
427  *insq_g11_sgx(j,i) &
428  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
429  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
430  *insq_g11_sgx(j,i-1) ) &
431  -dtt_2deta* &
432  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
433  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
434  *insq_g22_sgy(j,i) &
435  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
436  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
437  *insq_g22_sgy(j-1,i) )
438 
439 #elif (ADV_HOR==3)
440 
441  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
442  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
443 
444  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
445  -dtt_dxi* &
446  ( min(vx_c_help, 0.0_dp) &
447  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
448  *insq_g11_sgx(j,i) &
449  +max(vx_c_help, 0.0_dp) &
450  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
451  *insq_g11_sgx(j,i-1) ) &
452  -dtt_deta* &
453  ( min(vy_c_help, 0.0_dp) &
454  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
455  *insq_g22_sgy(j,i) &
456  +max(vy_c_help, 0.0_dp) &
457  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
458  *insq_g22_sgy(j-1,i) )
459 
460 #endif
461 
462 end do
463 
464 kc=kcmax
465 if (as_perp(j,i) >= zero) then
466  lgs_a0(kc) = 0.0_dp
467  lgs_a1(kc) = 1.0_dp
468  lgs_b(kc) = 0.0_dp
469 else
470  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
471  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
472  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
473  ! assumed/enforced
474 #if (ADV_HOR==2)
475 
476  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
477  -dtt_2dxi* &
478  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
479  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
480  *insq_g11_sgx(j,i) &
481  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
482  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
483  *insq_g11_sgx(j,i-1) ) &
484  -dtt_2deta* &
485  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
486  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
487  *insq_g22_sgy(j,i) &
488  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
489  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
490  *insq_g22_sgy(j-1,i) )
491 
492 #elif (ADV_HOR==3)
493 
494  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
495  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
496 
497  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
498  -dtt_dxi* &
499  ( min(vx_c_help, 0.0_dp) &
500  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
501  *insq_g11_sgx(j,i) &
502  +max(vx_c_help, 0.0_dp) &
503  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
504  *insq_g11_sgx(j,i-1) ) &
505  -dtt_deta* &
506  ( min(vy_c_help, 0.0_dp) &
507  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
508  *insq_g22_sgy(j,i) &
509  +max(vy_c_help, 0.0_dp) &
510  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
511  *insq_g22_sgy(j-1,i) )
512 
513 #endif
514 
515 end if
516 
517 !-------- Solve system of linear equations --------
518 
519 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
520 
521 !-------- Assign the result,
522 ! restriction to interval [0, AGE_MAX yr] --------
523 
524 do kc=0, kcmax
525 
526  age_c_neu(kc,j,i) = lgs_x(kc)
527 
528  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
529  age_c_neu(kc,j,i) = 0.0_dp
530  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
531  age_c_neu(kc,j,i) = age_max*year_sec
532 
533 end do
534 
535 !-------- Age of the ice in the non-existing temperate layer --------
536 
537 do kt=0, ktmax
538  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
539 end do
540 
541 end subroutine calc_temp2
542 !
sico_sle_solvers::tri_sle
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
Definition: sico_sle_solvers.F90:122
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
viscosity
real(dp) function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, i_flag_cold_temp)
Ice viscosity as a function of the effective strain rate and the temperature (in cold ice) or the wat...
Definition: viscosity.F90:39
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
ice_material_quantities::ratefac_c
real(dp) function, public ratefac_c(temp_val, temp_m_val)
Rate factor for cold ice: Linear interpolation of tabulated values in RF(.).
Definition: ice_material_quantities.F90:151
sico_sle_solvers
Solvers for systems of linear equations used by SICOPOLIS.
Definition: sico_sle_solvers.F90:35
ice_material_quantities::kappa_val
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
Definition: ice_material_quantities.F90:217
ice_material_quantities::c_val
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
Definition: ice_material_quantities.F90:254
ice_material_quantities
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Definition: ice_material_quantities.F90:37
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
creep
real(dp) function creep(sigma_val)
Creep response function for ice.
Definition: creep.F90:35
calc_temp2
subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, alb1, ai1, ai2, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature and age for an ice column with a temperate base overlain by cold ice...
Definition: calc_temp2.F90:37
Generated on Wed Mar 2 2016 11:26:54 for SICOPOLIS V3.2 by   doxygen 1.8.5