SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_temp3.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p 3
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age for an ice column with a
8 !! temperate base overlain by a temperate-ice layer.
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, water content and age for an ice column with a
35 !! temperate base overlain by a temperate-ice layer.
36 !<------------------------------------------------------------------------------
37 subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
38  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
39  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
40  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
41  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
42 
43 use sico_types
44 use sico_variables
45 
46 implicit none
47 integer(i4b) :: i, j, kc, kt, kr
48 real(dp) :: 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), ai3, &
52  atr1, am1, am2, alb1
53 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
54  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, cm1, cm2, &
55  clb1
56 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
57  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
58 real(dp) :: ci1(0:kcmax), ci2(0:kcmax), ci3
59 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
60 real(dp) :: cw1(0:ktmax), cw2(0:ktmax), cw3(0:ktmax), cw4(0:ktmax), &
61  cw5, cw7(0:ktmax), cw8, cw9(0:ktmax)
62 real(dp) :: cw1_sg(0:ktmax), cw2_sg(0:ktmax), cw3_sg(0:ktmax), &
63  cw4_sg(0:ktmax), adv_vert_w_sg(0:ktmax), abs_adv_vert_w_sg(0:ktmax)
64 real(dp) :: ftx_c_l(0:kcmax), ftx_c_r(0:kcmax), &
65  fty_c_l(0:kcmax), fty_c_r(0:kcmax), &
66  fax_c_l(0:kcmax), fax_c_r(0:kcmax), &
67  fay_c_l(0:kcmax), fay_c_r(0:kcmax)
68 real(dp) :: fwx_t_l(0:ktmax), fwx_t_r(0:ktmax), &
69  fwy_t_l(0:ktmax), fwy_t_r(0:ktmax), &
70  fax_t_l(0:ktmax), fax_t_r(0:ktmax), &
71  fay_t_l(0:ktmax), fay_t_r(0:ktmax)
72 real(dp) :: sigma_c_help(0:kcmax), sigma_t_help(0:ktmax), &
73  temp_c_help(0:kcmax)
74 real(dp) :: dtt_2dxi, dtt_2deta, dtime_temp, dzeta_t
75 real(dp) :: dtt_dxi, dtt_deta
76 real(dp) :: mean_accum_inv
77 real(dp) :: ratefac, ratefac_t, creep, kappa_val, c_val
78 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
79  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
80  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
81  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
82  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
83 real(dp), parameter :: zero=0.0_dp
84 
85 !-------- Check for boundary points --------
86 
87 if ((i.eq.0).or.(i.eq.imax).or.(j.eq.0).or.(j.eq.jmax)) &
88  stop ' calc_temp3: Boundary points not allowed.'
89 
90 !-------- Abbreviations --------
91 
92 ctr1 = atr1
93 cm1 = am1*h_c_neu(j,i)
94 clb1 = alb1*q_geo(j,i)
95 
96 #if ADV_VERT==1
97 
98 do kc=1, kcmax
99  ct1(kc) = at1(kc)/h_c_neu(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
100 end do
101 
102 #elif (ADV_VERT==2 || ADV_VERT==3)
103 
104 do kc=0, kcmax-1
105  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
106 end do
107 
108 #endif
109 
110 do kc=0, kcmax
111 
112  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
113  +at2_2(kc)*dh_c_dtau(j,i) )/h_c_neu(j,i)
114  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
115  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c_neu(j,i) &
116  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i) ! new factor !!!
117  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
118  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c_neu(j,i) &
119  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i) ! new factor !!!
120  ct5(kc) = at5(kc) &
121  /c_val(temp_c(kc,j,i)) &
122  /h_c_neu(j,i)
123 
124  sigma_c_help(kc) &
125  = rho*g*h_c_neu(j,i)*(1.0_dp-eaz_c_quotient(kc)) &
126  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
127  ct7(kc) = at7 &
128  /c_val(temp_c(kc,j,i)) &
129  *enh_c(kc,j,i) &
130  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
131  *creep(sigma_c_help(kc)) &
132  *sigma_c_help(kc)*sigma_c_help(kc)
133 
134  ci1(kc) = ai1(kc)/h_c_neu(j,i)
135 
136 ! ------ Fluxes for horizontal advection of temperature and age
137 
138 #if ADV_HOR==4
139  if (vx_c(kc,j,i-1).ge.zero) then
140  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
141  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
142  else
143  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
144  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
145  end if
146 
147  if (vx_c(kc,j,i).ge.zero) then
148  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
149  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
150  else
151  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
152  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
153  end if
154 
155  if (vy_c(kc,j-1,i).ge.zero) then
156  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
157  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
158  else
159  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
160  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
161  end if
162 
163  if (vy_c(kc,j,i).ge.zero) then
164  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
165  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
166  else
167  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
168  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
169  end if
170 #endif
171 
172 end do
173 
174 #if (ADV_VERT==2 || ADV_VERT==3)
175 
176 do kc=0, 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))
182 end do
183 
184 #endif
185 
186 do kc=0, kcmax-1
187  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
188  ct6(kc) = at6(kc) &
189  *kappa_val(temp_c_help(kc)) &
190  /h_c_neu(j,i)
191  ci2(kc) = ai2(kc)/h_c_neu(j,i)
192 end do
193 
194 cw5 = aw5/(h_t_neu(j,i)**2)
195 cw8 = aw8
196 ci3 = ai3/(h_t_neu(j,i)**2)
197 
198 #if ADV_VERT==1
199 
200 do kt=1, ktmax
201  cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i))
202 end do
203 
204 #elif (ADV_VERT==2 || ADV_VERT==3)
205 
206 do kt=0, ktmax-1
207  cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
208 end do
209 
210 #endif
211 
212 do kt=1, ktmax
213  cw9(kt) = aw9 &
214  *c_val(temp_t_m(kt,j,i)) &
215  *( dzs_dtau(j,i) &
216  +0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))*dzs_dxi_g(j,i) &
217  +0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))*dzs_deta_g(j,i) &
218  -0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i)) )
219 end do
220 
221 do kt=0, ktmax
222 
223  cw2(kt) = aw2*(dzb_dtau(j,i)+zeta_t(kt)*dh_t_dtau(j,i)) &
224  /h_t_neu(j,i)
225  cw3(kt) = aw3*(dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
226  /h_t_neu(j,i) &
227  *0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1)) *insq_g11_g(j,i) ! new factor !!!
228  cw4(kt) = aw4*(dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
229  /h_t_neu(j,i) &
230  *0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i)) *insq_g22_g(j,i) ! new factor !!!
231 
232  sigma_t_help(kt) &
233  = sigma_c_help(0) &
234  + rho*g*h_t_neu(j,i)*(1.0_dp-zeta_t(kt)) &
235  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
236  cw7(kt) = aw7 &
237  *enh_t(kt,j,i) &
238  *ratefac_t(omega_t(kt,j,i)) &
239  *creep(sigma_t_help(kt)) &
240  *sigma_t_help(kt)*sigma_t_help(kt)
241 
242 ! ------ Fluxes for horizontal advection of water content and age
243 
244 #if ADV_HOR==4
245  if (vx_t(kt,j,i-1).ge.zero) then
246  fwx_t_l(kt) = omega_t(kt,j,i-1)*vx_t(kt,j,i-1)
247  fax_t_l(kt) = age_t(kt,j,i-1)*vx_t(kt,j,i-1)
248  else
249  fwx_t_l(kt) = omega_t(kt,j,i)*vx_t(kt,j,i-1)
250  fax_t_l(kt) = age_t(kt,j,i)*vx_t(kt,j,i-1)
251  end if
252 
253  if (vx_t(kt,j,i).ge.zero) then
254  fwx_t_r(kt) = omega_t(kt,j,i)*vx_t(kt,j,i)
255  fax_t_r(kt) = age_t(kt,j,i)*vx_t(kt,j,i)
256  else
257  fwx_t_r(kt) = omega_t(kt,j,i+1)*vx_t(kt,j,i)
258  fax_t_r(kt) = age_t(kt,j,i+1)*vx_t(kt,j,i)
259  end if
260 
261  if (vy_t(kt,j-1,i).ge.zero) then
262  fwy_t_l(kt) = omega_t(kt,j-1,i)*vy_t(kt,j-1,i)
263  fay_t_l(kt) = age_t(kt,j-1,i)*vy_t(kt,j-1,i)
264  else
265  fwy_t_l(kt) = omega_t(kt,j,i)*vy_t(kt,j-1,i)
266  fay_t_l(kt) = age_t(kt,j,i)*vy_t(kt,j-1,i)
267  end if
268 
269  if (vy_t(kt,j,i).ge.zero) then
270  fwy_t_r(kt) = omega_t(kt,j,i)*vy_t(kt,j,i)
271  fay_t_r(kt) = age_t(kt,j,i)*vy_t(kt,j,i)
272  else
273  fwy_t_r(kt) = omega_t(kt,j+1,i)*vy_t(kt,j,i)
274  fay_t_r(kt) = age_t(kt,j+1,i)*vy_t(kt,j,i)
275  end if
276 #endif
277 
278 end do
279 
280 #if (ADV_VERT==2 || ADV_VERT==3)
281 
282 do kt=0, ktmax-1
283  cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
284  cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
285  cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
286  adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
287  abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
288 end do
289 
290 #endif
291 
292 #if ADV_HOR==4
293 dtt_dxi = 2.0_dp*dtt_2dxi
294 dtt_deta = 2.0_dp*dtt_2deta
295 #endif
296 
297 !-------- Set up the equations for the bedrock temperature --------
298 
299 kr=0
300 lgs_a1(kr) = 1.0_dp
301 lgs_a2(kr) = -1.0_dp
302 lgs_b(kr) = clb1
303 
304 #if Q_LITHO==1
305 ! (coupled heat-conducting bedrock)
306 
307 do kr=1, krmax-1
308  lgs_a0(kr) = - ctr1
309  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
310  lgs_a2(kr) = - ctr1
311  lgs_b(kr) = temp_r(kr,j,i)
312 end do
313 
314 #elif Q_LITHO==0
315 ! (no coupled heat-conducting bedrock)
316 
317 do kr=1, krmax-1
318  lgs_a0(kr) = 1.0_dp
319  lgs_a1(kr) = 0.0_dp
320  lgs_a2(kr) = -1.0_dp
321  lgs_b(kr) = 2.0_dp*clb1
322 end do
323 
324 #endif
325 
326 kr=krmax
327 lgs_a0(kr) = 0.0_dp
328 lgs_a1(kr) = 1.0_dp
329 lgs_b(kr) = temp_t_m(0,j,i)
330 
331 !-------- Solve system of linear equations --------
332 
333 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
334 
335 !-------- Assign the result --------
336 
337 do kr=0, krmax
338  temp_r_neu(kr,j,i) = lgs_x(kr)
339 end do
340 
341 !-------- Set up the equations for the water content in
342 ! temperate ice --------
343 
344 kt=0
345 lgs_a1(kt) = 1.0_dp
346 lgs_a2(kt) = -1.0_dp
347 lgs_b(kt) = 0.0_dp
348 
349 do kt=1, ktmax-1
350 
351 #if ADV_VERT==1
352 
353  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
354  lgs_a1(kt) = 1.0_dp + 2.0_dp*cw5
355  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
356 
357 #elif (ADV_VERT==2 || ADV_VERT==3)
358 
359  lgs_a0(kt) &
360  = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
361  -cw5
362  lgs_a1(kt) &
363  = 1.0_dp &
364  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
365  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
366  +2.0_dp*cw5
367  lgs_a2(kt) &
368  = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
369  -cw5
370 
371 #endif
372 
373 #if (ADV_HOR==2 || ADV_HOR==3)
374 
375  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
376  -dtt_2dxi* &
377  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
378  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
379  *insq_g11_sgx(j,i) &
380  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
381  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
382  *insq_g11_sgx(j,i-1) ) &
383  -dtt_2deta* &
384  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
385  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
386  *insq_g22_sgy(j,i) &
387  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
388  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
389  *insq_g22_sgy(j-1,i) )
390 
391 #elif ADV_HOR==4
392 
393  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
394  -dtt_dxi *(fwx_t_r(kt)-fwx_t_l(kt)) &
395  -dtt_deta*(fwy_t_r(kt)-fwy_t_l(kt))
396 
397 #endif
398 
399 end do
400 
401 kt=ktmax
402 
403 if (am_perp(j,i).ge.zero) then ! melting condition
404  lgs_a0(kt) = 0.0_dp
405  lgs_a1(kt) = 1.0_dp
406  lgs_b(kt) = 0.0_dp
407 else ! am_perp(j,i).lt.0.0, freezing condition
408  lgs_a0(kt) = -1.0_dp
409  lgs_a1(kt) = 1.0_dp
410  lgs_b(kt) = 0.0_dp
411 end if
412 
413 !-------- Solve system of linear equations --------
414 
415 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, ktmax)
416 
417 !-------- Assign the result, compute the water drainage --------
418 
419 q_tld(j,i) = 0.0_dp
420 
421 do kt=0, ktmax
422 
423  if (lgs_x(kt).lt.zero) then
424  omega_t_neu(kt,j,i) = 0.0_dp ! (as a precaution)
425  else if (lgs_x(kt).lt.omega_max) then
426  omega_t_neu(kt,j,i) = lgs_x(kt)
427  else
428  omega_t_neu(kt,j,i) = omega_max
429  q_tld(j,i) = q_tld(j,i) &
430  +aqtld*h_t_neu(j,i)*(lgs_x(kt)-omega_max)
431  end if
432 
433 end do
434 
435 !-------- Set up the equations for the ice temperature --------
436 
437 ! ------ Abbreviation for the jump of the temperature gradient with
438 ! the new omega
439 
440 if (am_perp(j,i).ge.zero) then ! melting condition
441  cm2 = 0.0_dp
442 else ! am_perp(j,i).lt.0.0, freezing condition
443  cm2 = am2*h_c_neu(j,i)*omega_t_neu(ktmax,j,i)*am_perp(j,i) &
444  /kappa_val(temp_c(0,j,i))
445 end if
446 
447 kc=0
448 lgs_a1(kc) = 1.0_dp
449 lgs_a2(kc) = -1.0_dp
450 lgs_b(kc) = -cm1-cm2
451 
452 do kc=1, kcmax-1
453 
454 #if ADV_VERT==1
455 
456  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
457  -ct5(kc)*ct6(kc-1)
458  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
459  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
460  -ct5(kc)*ct6(kc)
461 
462 #elif (ADV_VERT==2 || ADV_VERT==3)
463 
464  lgs_a0(kc) &
465  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
466  -ct5(kc)*ct6(kc-1)
467  lgs_a1(kc) &
468  = 1.0_dp &
469  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
470  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
471  +ct5(kc)*(ct6(kc)+ct6(kc-1))
472  lgs_a2(kc) &
473  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
474  -ct5(kc)*ct6(kc)
475 
476 #endif
477 
478 #if (ADV_HOR==2 || ADV_HOR==3)
479 
480  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
481  -dtt_2dxi* &
482  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
483  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
484  *insq_g11_sgx(j,i) &
485  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
486  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
487  *insq_g11_sgx(j,i-1) ) &
488  -dtt_2deta* &
489  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
490  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
491  *insq_g22_sgy(j,i) &
492  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
493  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
494  *insq_g22_sgy(j-1,i) )
495 
496 #elif ADV_HOR==4
497 
498  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
499  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
500  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
501 
502 #endif
503 
504 end do
505 
506 kc=kcmax
507 lgs_a0(kc) = 0.0_dp
508 lgs_a1(kc) = 1.0_dp
509 lgs_b(kc) = temp_s(j,i)
510 
511 !-------- Solve system of linear equations --------
512 
513 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
514 
515 !-------- Assign the result --------
516 
517 do kc=0, kcmax
518  temp_c_neu(kc,j,i) = lgs_x(kc)
519 end do
520 
521 !-------- Set up the equations for the age (cold and temperate ice
522 ! simultaneously) --------
523 
524 #if (ADV_HOR==3 && ADV_VERT==3)
525 
526 call calc_age_t(dtime_temp, i, j) ! Algorithm by B. Muegge
527 
528 #else
529 
530 #if ADV_VERT==1
531 
532 kt=0
533 lgs_a1(kt) = 1.0_dp
534 lgs_a2(kt) = -1.0_dp
535 lgs_b(kt) = m_age*h_t(j,i)*dzeta_t*mean_accum_inv
536 
537 #elif ADV_VERT==2
538 
539 kt=0
540 lgs_a1(kt) = 1.0_dp - adv_vert_w_sg(kt)
541 lgs_a2(kt) = adv_vert_w_sg(kt)
542 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
543  -dtt_2dxi* &
544  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
545  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
546  *insq_g11_sgx(j,i) &
547  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
548  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
549  *insq_g11_sgx(j,i-1) ) &
550  -dtt_2deta* &
551  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
552  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
553  *insq_g22_sgy(j,i) &
554  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
555  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
556  *insq_g22_sgy(j-1,i) )
557 
558 #endif
559 
560 do kt=1, ktmax-1
561 
562 #if ADV_VERT==1
563 
564  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
565  lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
566  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
567 
568 #elif ADV_VERT==2
569 
570  lgs_a0(kt) = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1))
571  lgs_a1(kt) = 1.0_dp &
572  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
573  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
574  lgs_a2(kt) = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
575 
576 #endif
577 
578 #if ADV_HOR==2
579 
580  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
581  -dtt_2dxi* &
582  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
583  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
584  *insq_g11_sgx(j,i) &
585  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
586  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
587  *insq_g11_sgx(j,i-1) ) &
588  -dtt_2deta* &
589  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
590  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
591  *insq_g22_sgy(j,i) &
592  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
593  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
594  *insq_g22_sgy(j-1,i) )
595 
596 #elif ADV_HOR==4
597 
598  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
599  -dtt_dxi *(fax_t_r(kt)-fax_t_l(kt)) &
600  -dtt_deta*(fay_t_r(kt)-fay_t_l(kt))
601 
602 #endif
603 
604 end do
605 
606 #if ADV_VERT==1
607 
608 kt=ktmax
609 kc=0
610 
611 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
612 lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
613 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
614 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
615  -dtt_2dxi* &
616  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
617  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
618  *insq_g11_sgx(j,i) &
619  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
620  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
621  *insq_g11_sgx(j,i-1) ) &
622  -dtt_2deta* &
623  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
624  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
625  *insq_g22_sgy(j,i) &
626  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
627  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
628  *insq_g22_sgy(j-1,i) )
629 
630 #elif ADV_VERT==2
631 
632 kt=ktmax
633 kc=0
634 
635 if (adv_vert_sg(kc).le.zero) then
636 
637  lgs_a0(ktmax+kc) = 0.0_dp
638  lgs_a1(ktmax+kc) = 1.0_dp - adv_vert_sg(kc)
639  lgs_a2(ktmax+kc) = adv_vert_sg(kc)
640  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
641  -dtt_2dxi* &
642  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
643  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
644  *insq_g11_sgx(j,i) &
645  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
646  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
647  *insq_g11_sgx(j,i-1) ) &
648  -dtt_2deta* &
649  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
650  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
651  *insq_g22_sgy(j,i) &
652  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
653  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
654  *insq_g22_sgy(j-1,i) )
655 
656 else if (adv_vert_w_sg(kt-1).ge.zero) then
657 
658  lgs_a0(kt) = -adv_vert_w_sg(kt-1)
659  lgs_a1(kt) = 1.0_dp + adv_vert_w_sg(kt-1)
660  lgs_a2(kt) = 0.0_dp
661  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
662  -dtt_2dxi* &
663  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
664  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
665  *insq_g11_sgx(j,i) &
666  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
667  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
668  *insq_g11_sgx(j,i-1) ) &
669  -dtt_2deta* &
670  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
671  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
672  *insq_g22_sgy(j,i) &
673  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
674  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
675  *insq_g22_sgy(j-1,i) )
676 
677 else
678 
679  lgs_a0(kt) = -0.5_dp
680  lgs_a1(kt) = 1.0_dp
681  lgs_a2(kt) = -0.5_dp
682  lgs_b(kt) = 0.0_dp
683  ! Makeshift: Average of age_c(kc=1) and age_t(kt=KTMAX-1)
684 
685 end if
686 
687 #endif
688 
689 do kc=1, kcmax-1
690 
691 #if ADV_VERT==1
692 
693  lgs_a0(ktmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
694  -ci1(kc)*ci2(kc-1)
695  lgs_a1(ktmax+kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
696  lgs_a2(ktmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
697  -ci1(kc)*ci2(kc)
698 
699 #elif ADV_VERT==2
700 
701  lgs_a0(ktmax+kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
702  lgs_a1(ktmax+kc) = 1.0_dp &
703  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
704  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
705  lgs_a2(ktmax+kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
706 
707 #endif
708 
709 #if ADV_HOR==2
710 
711  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
712  -dtt_2dxi* &
713  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
714  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
715  *insq_g11_sgx(j,i) &
716  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
717  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
718  *insq_g11_sgx(j,i-1) ) &
719  -dtt_2deta* &
720  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
721  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
722  *insq_g22_sgy(j,i) &
723  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
724  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
725  *insq_g22_sgy(j-1,i) )
726 
727 #elif ADV_HOR==4
728 
729  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
730  -dtt_dxi *(fax_c_r(kc)-fax_c_l(kc)) &
731  -dtt_deta*(fay_c_r(kc)-fay_c_l(kc))
732 
733 #endif
734 
735 end do
736 
737 #if ADV_VERT==1
738 
739 kc=kcmax
740 if (as_perp(j,i).ge.zero) then
741  lgs_a0(ktmax+kc) = 0.0_dp
742  lgs_a1(ktmax+kc) = 1.0_dp
743  lgs_b(ktmax+kc) = 0.0_dp
744 else
745  lgs_a0(ktmax+kc) = -1.0_dp
746  lgs_a1(ktmax+kc) = 1.0_dp
747  lgs_b(ktmax+kc) = 0.0_dp
748 end if
749 
750 #elif ADV_VERT==2
751 
752 kc=kcmax
753 if (as_perp(j,i).ge.zero) then
754  lgs_a0(ktmax+kc) = 0.0_dp
755  lgs_a1(ktmax+kc) = 1.0_dp
756  lgs_b(ktmax+kc) = 0.0_dp
757 else
758  lgs_a0(ktmax+kc) = -adv_vert_sg(kc-1)
759  lgs_a1(ktmax+kc) = 1.0_dp + adv_vert_sg(kc-1)
760  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
761  -dtt_2dxi* &
762  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
763  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
764  *insq_g11_sgx(j,i) &
765  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
766  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
767  *insq_g11_sgx(j,i-1) ) &
768  -dtt_2deta* &
769  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
770  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
771  *insq_g22_sgy(j,i) &
772  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
773  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
774  *insq_g22_sgy(j-1,i) )
775 end if
776 
777 #endif
778 
779 !-------- Solve system of linear equations --------
780 
781 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+ktmax)
782 
783 !-------- Assign the result,
784 ! restriction to interval [0, AGE_MAX yr] --------
785 
786 do kt=0, ktmax
787 
788  age_t_neu(kt,j,i) = lgs_x(kt)
789 
790  if (age_t_neu(kt,j,i).lt.(age_min*year_sec)) &
791  age_t_neu(kt,j,i) = 0.0_dp
792  if (age_t_neu(kt,j,i).gt.(age_max*year_sec)) &
793  age_t_neu(kt,j,i) = age_max*year_sec
794 
795 end do
796 
797 do kc=0, kcmax
798 
799  age_c_neu(kc,j,i) = lgs_x(ktmax+kc)
800 
801  if (age_c_neu(kc,j,i).lt.(age_min*year_sec)) &
802  age_c_neu(kc,j,i) = 0.0_dp
803  if (age_c_neu(kc,j,i).gt.(age_max*year_sec)) &
804  age_c_neu(kc,j,i) = age_max*year_sec
805 
806 end do
807 
808 #endif
809 
810 end subroutine calc_temp3
811 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1