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