SICOPOLIS V3.0
 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 
44 implicit none
45 integer(i4b) :: i, j, kc, kt, kr
46 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
47  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
48  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
49  ai1(0:kcmax), ai2(0:kcmax), ai4, &
50  atr1, acb1, acb2, acb3, acb4, alb1
51 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
52  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, &
53  ccb1, ccb2, ccb3, ccb4, 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_temp1: Boundary points not allowed.'
77 
78 !-------- Abbreviations --------
79 
80 ctr1 = atr1
81 
82 ccb1 = acb1 &
83  *kappa_val(temp_c(0,j,i)) &
84  /h_c(j,i)
85 ccb2 = acb2
86 ccb3 = acb3*0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
87  *h_c(j,i)*dzs_dxi_g(j,i)
88 ccb4 = acb4*0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
89  *h_c(j,i)*dzs_deta_g(j,i)
90 
91 clb1 = alb1*q_geo(j,i)
92 
93 #if ADV_VERT==1
94 
95 do kc=1, kcmax
96  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
97 end do
98 
99 #elif (ADV_VERT==2 || ADV_VERT==3)
100 
101 do kc=0, kcmax-1
102  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
103 end do
104 
105 #endif
106 
107 do kc=0, kcmax
108 
109  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
110  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
111  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
112  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
113  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i) ! new factor !!!
114  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
115  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
116  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i) ! new factor !!!
117  ct5(kc) = at5(kc) &
118  /c_val(temp_c(kc,j,i)) &
119  /h_c(j,i)
120  ct7(kc) = at7 &
121  /c_val(temp_c(kc,j,i)) &
122  *enh_c(kc,j,i) &
123  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
124  *creep(sigma_c(kc,j,i)) &
125  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
126  ci1(kc) = ai1(kc)/h_c(j,i)
127 
128 ! ------ Fluxes for horizontal advection of temperature and age
129 
130 #if ADV_HOR==4
131  if (vx_c(kc,j,i-1).ge.zero) then
132  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
133  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
134  else
135  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
136  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
137  end if
138 
139  if (vx_c(kc,j,i).ge.zero) then
140  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
141  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
142  else
143  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
144  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
145  end if
146 
147  if (vy_c(kc,j-1,i).ge.zero) then
148  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
149  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
150  else
151  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
152  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
153  end if
154 
155  if (vy_c(kc,j,i).ge.zero) then
156  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
157  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
158  else
159  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
160  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
161  end if
162 #endif
163 
164 end do
165 
166 #if (ADV_VERT==2 || ADV_VERT==3)
167 
168 do kc=0, kcmax-1
169  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
170  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
171  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
172  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
173  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
174 end do
175 
176 #endif
177 
178 do kc=0, kcmax-1
179  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
180  ct6(kc) = at6(kc) &
181  *kappa_val(temp_c_help(kc)) &
182  /h_c(j,i)
183  ci2(kc) = ai2(kc)/h_c(j,i)
184 end do
185 
186 #if ADV_HOR==4
187 dtt_dxi = 2.0_dp*dtt_2dxi
188 dtt_deta = 2.0_dp*dtt_2deta
189 #endif
190 
191 !-------- Set up the temperature equations (ice and bedrock
192 ! simultaneously) --------
193 
194 kr=0
195 lgs_a1(kr) = 1.0_dp
196 lgs_a2(kr) = -1.0_dp
197 lgs_b(kr) = clb1
198 
199 #if Q_LITHO==1
200 ! (coupled heat-conducting bedrock)
201 
202 do kr=1, krmax-1
203  lgs_a0(kr) = - ctr1
204  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
205  lgs_a2(kr) = - ctr1
206  lgs_b(kr) = temp_r(kr,j,i)
207 end do
208 
209 #elif Q_LITHO==0
210 ! (no coupled heat-conducting bedrock)
211 
212 do kr=1, krmax-1
213  lgs_a0(kr) = 1.0_dp
214  lgs_a1(kr) = 0.0_dp
215  lgs_a2(kr) = -1.0_dp
216  lgs_b(kr) = 2.0_dp*clb1
217 end do
218 
219 #endif
220 
221 kr=krmax
222 kc=0
223 lgs_a0(kr) = ccb2
224 lgs_a1(kr) = -(ccb1+ccb2)
225 lgs_a2(kr) = ccb1
226 lgs_b(kr) = ccb3+ccb4
227 
228 do kc=1, kcmax-1
229 
230 #if ADV_VERT==1
231 
232  lgs_a0(krmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
233  -ct5(kc)*ct6(kc-1)
234  lgs_a1(krmax+kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
235  lgs_a2(krmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
236  -ct5(kc)*ct6(kc)
237 
238 #elif (ADV_VERT==2 || ADV_VERT==3)
239 
240  lgs_a0(krmax+kc) &
241  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
242  -ct5(kc)*ct6(kc-1)
243  lgs_a1(krmax+kc) &
244  = 1.0_dp &
245  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
246  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
247  +ct5(kc)*(ct6(kc)+ct6(kc-1))
248  lgs_a2(krmax+kc) &
249  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
250  -ct5(kc)*ct6(kc)
251 
252 #endif
253 
254 #if (ADV_HOR==2 || ADV_HOR==3)
255 
256  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
257  -dtt_2dxi* &
258  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
259  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
260  *insq_g11_sgx(j,i) &
261  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
262  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
263  *insq_g11_sgx(j,i-1) ) &
264  -dtt_2deta* &
265  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
266  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
267  *insq_g22_sgy(j,i) &
268  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
269  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
270  *insq_g22_sgy(j-1,i) )
271 
272 #elif ADV_HOR==4
273 
274  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
275  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
276  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
277 
278 #endif
279 
280 end do
281 
282 kc=kcmax
283 lgs_a0(krmax+kc) = 0.0_dp
284 lgs_a1(krmax+kc) = 1.0_dp
285 lgs_b(krmax+kc) = temp_s(j,i)
286 
287 !-------- Solve system of linear equations --------
288 
289 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+krmax)
290 
291 !-------- Assign the result --------
292 
293 do kr=0, krmax
294  temp_r_neu(kr,j,i) = lgs_x(kr)
295 end do
296 
297 do kc=0, kcmax
298  temp_c_neu(kc,j,i) = lgs_x(krmax+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 the 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 
441 !-------- Solve system of linear equations --------
442 
443 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
444 
445 !-------- Assign the result,
446 ! restriction to interval [0, AGE_MAX yr] --------
447 
448 do kc=0, kcmax
449 
450  age_c_neu(kc,j,i) = lgs_x(kc)
451 
452  if (age_c_neu(kc,j,i).lt.(age_min*year_sec)) &
453  age_c_neu(kc,j,i) = 0.0_dp
454  if (age_c_neu(kc,j,i).gt.(age_max*year_sec)) &
455  age_c_neu(kc,j,i) = age_max*year_sec
456 
457 end do
458 
459 !-------- Age of the ice in the non-existing temperate layer --------
460 
461 do kt=0, ktmax
462  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
463 end do
464 
465 #endif
466 
467 end subroutine calc_temp1
468 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1