SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_temp1.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p 1
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age for a cold ice column.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve
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 temperature and age for a cold ice column.
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
36  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
37  acb3, acb4, alb1, ai1, ai2, ai4, &
38  mean_accum_inv, &
39  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
40 
41 use sico_types
42 use sico_variables
43 use sico_sle_solvers, only : tri_sle
44 
45 implicit none
46 
47 integer(i4b), intent(in) :: i, j
48 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
49  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
50  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
51  ai1(0:kcmax), ai2(0:kcmax), ai4, &
52  atr1, acb1, acb2, acb3, acb4, alb1
53 real(dp), intent(in) :: mean_accum_inv
54 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
55 
56 integer(i4b) :: kc, kt, kr
57 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
58  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, &
59  ccb1, ccb2, ccb3, ccb4, 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_temp1: Boundary points not allowed.'
82 
83 !-------- Abbreviations --------
84 
85 ctr1 = atr1
86 
87 ccb1 = acb1 &
88  *kappa_val(temp_c(0,j,i)) &
89  /h_c(j,i)
90 ccb2 = acb2
91 ccb3 = acb3*0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
92  *h_c(j,i)*dzs_dxi_g(j,i)
93 ccb4 = acb4*0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
94  *h_c(j,i)*dzs_deta_g(j,i)
95 
96 clb1 = alb1*q_geo(j,i)
97 
98 #if ADV_VERT==1
99 
100 do kc=1, kcmax-1
101  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
102 end do
103 
104 kc=0
105 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
106  ! only needed for kc=0 ...
107 kc=kcmax-1
108 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
109  ! ... and kc=KCMAX-1
110 
111 #elif ( ADV_VERT==2 || ADV_VERT==3 )
112 
113 do kc=0, kcmax-1
114  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
115 end do
116 
117 #endif
118 
119 do kc=0, kcmax
120 
121  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
122  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
123  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
124  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
125  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
126  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
127  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
128  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
129  ct5(kc) = at5(kc) &
130  /c_val(temp_c(kc,j,i)) &
131  /h_c(j,i)
132  ct7(kc) = at7 &
133  /c_val(temp_c(kc,j,i)) &
134  *enh_c(kc,j,i) &
135  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
136  *creep(sigma_c(kc,j,i)) &
137  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
138  ci1(kc) = ai1(kc)/h_c(j,i)
139 
140 ! ------ Fluxes for horizontal advection of temperature and age
141 
142 #if ADV_HOR==4
143  if (vx_c(kc,j,i-1) >= zero) then
144  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
145  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
146  else
147  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
148  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
149  end if
150 
151  if (vx_c(kc,j,i) >= zero) then
152  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
153  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
154  else
155  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
156  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
157  end if
158 
159  if (vy_c(kc,j-1,i) >= zero) then
160  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
161  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
162  else
163  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
164  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
165  end if
166 
167  if (vy_c(kc,j,i) >= zero) then
168  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
169  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
170  else
171  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
172  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
173  end if
174 #endif
175 
176 end do
177 
178 #if ADV_VERT==1
179 
180 kc=0
181 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
182 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
183 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
184 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
185 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
186 kc=kcmax-1
187 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
188 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
189 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
190 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
191 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
192 
193 #elif ( ADV_VERT==2 || ADV_VERT==3 )
194 
195 do kc=0, kcmax-1
196  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
197  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
198  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
199  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
200  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
201 end do
202 
203 #endif
204 
205 do kc=0, kcmax-1
206  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
207  ct6(kc) = at6(kc) &
208  *kappa_val(temp_c_help(kc)) &
209  /h_c(j,i)
210  ci2(kc) = ai2(kc)/h_c(j,i)
211 end do
212 
213 #if ( ADV_HOR==3 || ADV_HOR==4 )
214 dtt_dxi = 2.0_dp*dtt_2dxi
215 dtt_deta = 2.0_dp*dtt_2deta
216 #endif
217 
218 !-------- Set up the temperature equations (ice and bedrock
219 ! simultaneously) --------
220 
221 kr=0
222 lgs_a1(kr) = 1.0_dp
223 lgs_a2(kr) = -1.0_dp
224 lgs_b(kr) = clb1
225 
226 #if Q_LITHO==1
227 ! (coupled heat-conducting bedrock)
228 
229 do kr=1, krmax-1
230  lgs_a0(kr) = - ctr1
231  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
232  lgs_a2(kr) = - ctr1
233  lgs_b(kr) = temp_r(kr,j,i)
234 end do
235 
236 #elif Q_LITHO==0
237 ! (no coupled heat-conducting bedrock)
238 
239 do kr=1, krmax-1
240  lgs_a0(kr) = 1.0_dp
241  lgs_a1(kr) = 0.0_dp
242  lgs_a2(kr) = -1.0_dp
243  lgs_b(kr) = 2.0_dp*clb1
244 end do
245 
246 #endif
247 
248 kr=krmax
249 kc=0
250 lgs_a0(kr) = ccb2
251 lgs_a1(kr) = -(ccb1+ccb2)
252 lgs_a2(kr) = ccb1
253 lgs_b(kr) = ccb3+ccb4
254 
255 do kc=1, kcmax-1
256 
257 #if ADV_VERT==1
258 
259  lgs_a0(krmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
260  -ct5(kc)*ct6(kc-1)
261  lgs_a1(krmax+kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
262  lgs_a2(krmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
263  -ct5(kc)*ct6(kc)
264 
265 #elif ADV_VERT==2
266 
267  lgs_a0(krmax+kc) &
268  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
269  -ct5(kc)*ct6(kc-1)
270  lgs_a1(krmax+kc) &
271  = 1.0_dp &
272  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
273  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
274  +ct5(kc)*(ct6(kc)+ct6(kc-1))
275  lgs_a2(krmax+kc) &
276  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
277  -ct5(kc)*ct6(kc)
278 
279 #elif ADV_VERT==3
280 
281  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
282 
283  lgs_a0(krmax+kc) &
284  = -max(adv_vert_help, 0.0_dp) &
285  -ct5(kc)*ct6(kc-1)
286  lgs_a1(krmax+kc) &
287  = 1.0_dp &
288  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
289  +ct5(kc)*(ct6(kc)+ct6(kc-1))
290  lgs_a2(krmax+kc) &
291  = min(adv_vert_help, 0.0_dp) &
292  -ct5(kc)*ct6(kc)
293 
294 #endif
295 
296 #if ADV_HOR==2
297 
298  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
299  -dtt_2dxi* &
300  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
301  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
302  *insq_g11_sgx(j,i) &
303  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
304  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
305  *insq_g11_sgx(j,i-1) ) &
306  -dtt_2deta* &
307  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
308  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
309  *insq_g22_sgy(j,i) &
310  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
311  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
312  *insq_g22_sgy(j-1,i) )
313 
314 #elif ADV_HOR==3
315 
316  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
317  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
318 
319  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
320  -dtt_dxi* &
321  ( min(vx_c_help, 0.0_dp) &
322  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
323  *insq_g11_sgx(j,i) &
324  +max(vx_c_help, 0.0_dp) &
325  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
326  *insq_g11_sgx(j,i-1) ) &
327  -dtt_deta* &
328  ( min(vy_c_help, 0.0_dp) &
329  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
330  *insq_g22_sgy(j,i) &
331  +max(vy_c_help, 0.0_dp) &
332  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
333  *insq_g22_sgy(j-1,i) )
334 
335 #elif ADV_HOR==4
336 
337  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
338  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
339  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
340 
341 #endif
342 
343 end do
344 
345 kc=kcmax
346 lgs_a0(krmax+kc) = 0.0_dp
347 lgs_a1(krmax+kc) = 1.0_dp
348 lgs_b(krmax+kc) = temp_s(j,i)
349 
350 !-------- Solve system of linear equations --------
351 
352 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+krmax)
353 
354 !-------- Assign the result --------
355 
356 do kr=0, krmax
357  temp_r_neu(kr,j,i) = lgs_x(kr)
358 end do
359 
360 do kc=0, kcmax
361  temp_c_neu(kc,j,i) = lgs_x(krmax+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_temp1
586 !
Generated on Thu Sep 12 2013 09:45:09 for SICOPOLIS V3.1 by   doxygen 1.8.2