SICOPOLIS V3.0
 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 
45 implicit none
46 integer(i4b) :: i, j, kc, kt, kr
47 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
48  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
49  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
50  ai1(0:kcmax), ai2(0:kcmax), ai4, &
51  atr1, alb1
52 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
53  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
54 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
55  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
56 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
57 real(dp) :: ftx_c_l(0:kcmax), ftx_c_r(0:kcmax), &
58  fty_c_l(0:kcmax), fty_c_r(0:kcmax), &
59  fax_c_l(0:kcmax), fax_c_r(0:kcmax), &
60  fay_c_l(0:kcmax), fay_c_r(0:kcmax)
61 real(dp) :: temp_c_help(0:kcmax)
62 real(dp) :: mean_accum_inv
63 real(dp) :: dtt_2dxi, dtt_2deta, dtime_temp
64 real(dp) :: dtt_dxi, dtt_deta
65 real(dp) :: ratefac, creep, kappa_val, c_val
66 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
67  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
68  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
69  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
70  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
71 real(dp), parameter :: zero=0.0_dp
72 
73 !-------- Check for boundary points --------
74 
75 if ((i.eq.0).or.(i.eq.imax).or.(j.eq.0).or.(j.eq.jmax)) &
76  stop ' calc_temp2: Boundary points not allowed.'
77 
78 !-------- Abbreviations --------
79 
80 ctr1 = atr1
81 clb1 = alb1*q_geo(j,i)
82 
83 #if ADV_VERT==1
84 
85 do kc=1, kcmax
86  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
87 end do
88 
89 #elif (ADV_VERT==2 || ADV_VERT==3)
90 
91 do kc=0, kcmax-1
92  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
93 end do
94 
95 #endif
96 
97 do kc=0, kcmax
98 
99  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
100  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
101  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
102  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
103  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i) ! new factor !!!
104  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
105  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
106  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i) ! new factor !!!
107  ct5(kc) = at5(kc) &
108  /c_val(temp_c(kc,j,i)) &
109  /h_c(j,i)
110  ct7(kc) = at7 &
111  /c_val(temp_c(kc,j,i)) &
112  *enh_c(kc,j,i) &
113  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
114  *creep(sigma_c(kc,j,i)) &
115  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
116  ci1(kc) = ai1(kc)/h_c(j,i)
117 
118 ! ------ Fluxes for horizontal advection of temperature and age
119 
120 #if ADV_HOR==4
121  if (vx_c(kc,j,i-1).ge.zero) then
122  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
123  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
124  else
125  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
126  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
127  end if
128 
129  if (vx_c(kc,j,i).ge.zero) then
130  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
131  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
132  else
133  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
134  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
135  end if
136 
137  if (vy_c(kc,j-1,i).ge.zero) then
138  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
139  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
140  else
141  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
142  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
143  end if
144 
145  if (vy_c(kc,j,i).ge.zero) then
146  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
147  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
148  else
149  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
150  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
151  end if
152 #endif
153 
154 end do
155 
156 #if (ADV_VERT==2 || ADV_VERT==3)
157 
158 do kc=0, kcmax-1
159  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
160  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
161  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
162  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
163  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
164 end do
165 
166 #endif
167 
168 do kc=0, kcmax-1
169  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
170  ct6(kc) = at6(kc) &
171  *kappa_val(temp_c_help(kc)) &
172  /h_c(j,i)
173  ci2(kc) = ai2(kc)/h_c(j,i)
174 end do
175 
176 #if ADV_HOR==4
177 dtt_dxi = 2.0_dp*dtt_2dxi
178 dtt_deta = 2.0_dp*dtt_2deta
179 #endif
180 
181 !-------- Set up the equations for the bedrock temperature --------
182 
183 kr=0
184 lgs_a1(kr) = 1.0_dp
185 lgs_a2(kr) = -1.0_dp
186 lgs_b(kr) = clb1
187 
188 #if Q_LITHO==1
189 ! (coupled heat-conducting bedrock)
190 
191 do kr=1, krmax-1
192  lgs_a0(kr) = - ctr1
193  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
194  lgs_a2(kr) = - ctr1
195  lgs_b(kr) = temp_r(kr,j,i)
196 end do
197 
198 #elif Q_LITHO==0
199 ! (no coupled heat-conducting bedrock)
200 
201 do kr=1, krmax-1
202  lgs_a0(kr) = 1.0_dp
203  lgs_a1(kr) = 0.0_dp
204  lgs_a2(kr) = -1.0_dp
205  lgs_b(kr) = 2.0_dp*clb1
206 end do
207 
208 #endif
209 
210 kr=krmax
211 lgs_a0(kr) = 0.0_dp
212 lgs_a1(kr) = 1.0_dp
213 lgs_b(kr) = temp_t_m(0,j,i)
214 
215 !-------- Solve system of linear equations --------
216 
217 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
218 
219 !-------- Assign the result --------
220 
221 do kr=0, krmax
222  temp_r_neu(kr,j,i) = lgs_x(kr)
223 end do
224 
225 !-------- Set up the equations for the ice temperature --------
226 
227 kc=0
228 lgs_a1(kc) = 1.0_dp
229 lgs_a2(kc) = 0.0_dp
230 lgs_b(kc) = temp_c_m(0,j,i)
231 
232 do kc=1, kcmax-1
233 
234 #if ADV_VERT==1
235 
236  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
237  -ct5(kc)*ct6(kc-1)
238  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
239  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
240  -ct5(kc)*ct6(kc)
241 
242 #elif (ADV_VERT==2 || ADV_VERT==3)
243 
244  lgs_a0(kc) &
245  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
246  -ct5(kc)*ct6(kc-1)
247  lgs_a1(kc) &
248  = 1.0_dp &
249  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
250  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
251  +ct5(kc)*(ct6(kc)+ct6(kc-1))
252  lgs_a2(kc) &
253  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
254  -ct5(kc)*ct6(kc)
255 
256 #endif
257 
258 #if (ADV_HOR==2 || ADV_HOR==3)
259 
260  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
261  -dtt_2dxi* &
262  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
263  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
264  *insq_g11_sgx(j,i) &
265  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
266  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
267  *insq_g11_sgx(j,i-1) ) &
268  -dtt_2deta* &
269  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
270  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
271  *insq_g22_sgy(j,i) &
272  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
273  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
274  *insq_g22_sgy(j-1,i) )
275 
276 #elif ADV_HOR==4
277 
278  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
279  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
280  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
281 
282 #endif
283 
284 end do
285 
286 kc=kcmax
287 lgs_a0(kc) = 0.0_dp
288 lgs_a1(kc) = 1.0_dp
289 lgs_b(kc) = temp_s(j,i)
290 
291 !-------- Solve system of linear equations --------
292 
293 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
294 
295 !-------- Assign the result --------
296 
297 do kc=0, kcmax
298  temp_c_neu(kc,j,i) = lgs_x(kc)
299 end do
300 
301 !-------- Set water content in the non-existing temperate layer
302 ! to zero --------
303 
304 do kt=0, ktmax
305  omega_t_neu(kt,j,i) = 0.0_dp
306 end do
307 
308 !-------- Water drainage from the non-existing temperate layer --------
309 
310 q_tld(j,i) = 0.0_dp
311 
312 !-------- Set up the equations for the age of cold ice --------
313 
314 #if (ADV_HOR==3 && ADV_VERT==3)
315 
316 call calc_age_c(dtime_temp, i, j) ! Algorithm by B. Muegge
317 
318 #else
319 
320 #if ADV_VERT==1
321 
322 kc=0
323 lgs_a1(kc) = 1.0_dp
324 lgs_a2(kc) = -1.0_dp
325 lgs_b(kc) = m_age*h_c(j,i)*ai4*mean_accum_inv
326 
327 #elif ADV_VERT==2
328 
329 kc=0
330 lgs_a1(kc) = 1.0_dp - adv_vert_sg(kc)
331 lgs_a2(kc) = adv_vert_sg(kc)
332 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
333  -dtt_2dxi* &
334  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
335  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
336  *insq_g11_sgx(j,i) &
337  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
338  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
339  *insq_g11_sgx(j,i-1) ) &
340  -dtt_2deta* &
341  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
342  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
343  *insq_g22_sgy(j,i) &
344  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
345  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
346  *insq_g22_sgy(j-1,i) )
347 
348 #endif
349 
350 do kc=1, kcmax-1
351 
352 #if ADV_VERT==1
353 
354  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
355  -ci1(kc)*ci2(kc-1)
356  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
357  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
358  -ci1(kc)*ci2(kc)
359 
360 #elif ADV_VERT==2
361 
362  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
363  lgs_a1(kc) = 1.0_dp &
364  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
365  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
366  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
367 
368 #endif
369 
370 #if ADV_HOR==2
371 
372  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
373  -dtt_2dxi* &
374  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
375  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
376  *insq_g11_sgx(j,i) &
377  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
378  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
379  *insq_g11_sgx(j,i-1) ) &
380  -dtt_2deta* &
381  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
382  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
383  *insq_g22_sgy(j,i) &
384  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
385  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
386  *insq_g22_sgy(j-1,i) )
387 
388 #elif ADV_HOR==4
389 
390  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
391  -dtt_dxi *(fax_c_r(kc)-fax_c_l(kc)) &
392  -dtt_deta*(fay_c_r(kc)-fay_c_l(kc))
393 
394 #endif
395 
396 end do
397 
398 #if ADV_VERT==1
399 
400 kc=kcmax
401 if (as_perp(j,i).ge.zero) then
402  lgs_a0(kc) = 0.0_dp
403  lgs_a1(kc) = 1.0_dp
404  lgs_b(kc) = 0.0_dp
405 else
406  lgs_a0(kc) = -1.0_dp
407  lgs_a1(kc) = 1.0_dp
408  lgs_b(kc) = 0.0_dp
409 end if
410 
411 #elif ADV_VERT==2
412 
413 kc=kcmax
414 if (as_perp(j,i).ge.zero) then
415  lgs_a0(kc) = 0.0_dp
416  lgs_a1(kc) = 1.0_dp
417  lgs_b(kc) = 0.0_dp
418 else
419  lgs_a0(kc) = -adv_vert_sg(kc-1)
420  lgs_a1(kc) = 1.0_dp + adv_vert_sg(kc-1)
421  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
422  -dtt_2dxi* &
423  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
424  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
425  *insq_g11_sgx(j,i) &
426  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
427  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
428  *insq_g11_sgx(j,i-1) ) &
429  -dtt_2deta* &
430  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
431  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
432  *insq_g22_sgy(j,i) &
433  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
434  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
435  *insq_g22_sgy(j-1,i) )
436 end if
437 
438 #endif
439 
440 !-------- Solve system of linear equations --------
441 
442 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
443 
444 !-------- Assign the result,
445 ! restriction to interval [0, AGE_MAX yr] --------
446 
447 do kc=0, kcmax
448 
449  age_c_neu(kc,j,i) = lgs_x(kc)
450 
451  if (age_c_neu(kc,j,i).lt.(age_min*year_sec)) &
452  age_c_neu(kc,j,i) = 0.0_dp
453  if (age_c_neu(kc,j,i).gt.(age_max*year_sec)) &
454  age_c_neu(kc,j,i) = age_max*year_sec
455 
456 end do
457 
458 !-------- Age of the ice in the non-existing temperate layer --------
459 
460 do kt=0, ktmax
461  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
462 end do
463 
464 #endif
465 
466 end subroutine calc_temp2
467 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1