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