SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_temp_enth_2.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ e n t h _ 2
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age for an ice column with a temperate base
8 !! with the enthalpy method.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2013-2016 Ralf Greve, Heinz Blatter
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 !! with the enthalpy method.
36 !<------------------------------------------------------------------------------
37 subroutine calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
38  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
39  ai1, ai2, aqtlde, am3, &
40  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
41 
42 use sico_types
44 use sico_vars
45 use sico_sle_solvers, only : tri_sle
47 
48 implicit none
49 
50 integer(i4b), intent(in) :: i, j
51 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
52  at3_1(0:kcmax), at3_2(0:kcmax), &
53  at4_1(0:kcmax), at4_2(0:kcmax), &
54  at5(0:kcmax), at6(0:kcmax), at7, &
55  ai1(0:kcmax), ai2(0:kcmax), &
56  atr1, alb1, aqtlde(0:kcmax), am3(0:kcmax)
57 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
58 
59 integer(i4b) :: kc, kt, kr
60 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
61  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
62 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
63  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
64 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
65 real(dp) :: cqtlde(0:kcmax), cm3(0:kcmax)
66 real(dp) :: dtt_dxi, dtt_deta
67 real(dp) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
68 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
69  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
70  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
71  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
72  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
73 
74 real(dp), parameter :: eps_omega=1.0e-12_dp
75 
76 !-------- Check for boundary points --------
77 
78 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
79  stop ' calc_temp_enth_2: Boundary points not allowed.'
80 
81 !-------- Abbreviations --------
82 
83 call calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
84  at4_1, at4_2, at5, atr1, alb1, &
85  ai1, aqtlde, &
86  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
87  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
88  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
89  adv_vert_sg, abs_adv_vert_sg, &
90  ci1, cqtlde, dtt_dxi, dtt_deta)
91 
92 do kc=0, kcmax
93  temp_c_val(kc) = temp_c(kc,j,i)
94  omega_c_val(kc) = omega_c(kc,j,i)
95 end do
96 
97 call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
98  i, j, ce6, ce7, ci2, cm3)
99 
100 !-------- Computation of the bedrock temperature --------
101 
102 ! ------ Set-up of the the equations
103 
104 call calc_temp_enth_2_b(ctr1, clb1, i, j, &
105  lgs_a0, lgs_a1, lgs_a2, lgs_b)
106 
107 ! ------ Solution of the system of linear equations
108 
109 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
110 
111 ! ------ Assignment of the result
112 
113 do kr=0, krmax
114  temp_r_neu(kr,j,i) = lgs_x(kr)
115 end do
116 
117 !-------- Computation of the ice enthalpy (predictor step) --------
118 
119 ! ------ Set-up of the equations
120 
121 call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
122  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
123  adv_vert_sg, abs_adv_vert_sg, &
124  dtime_temp, dtt_dxi, dtt_deta, &
125  dtt_2dxi, dtt_2deta, &
126  i, j, 0_i2b, &
127  lgs_a0, lgs_a1, lgs_a2, lgs_b)
128 
129 ! ------ Solution of the system of linear equations
130 
131 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
132 
133 ! ------ Assignment of the result
134 
135 do kc=0, kcmax
136  enth_c_neu(kc,j,i) = lgs_x(kc)
137  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
138  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
139 end do
140 
141 ! ------ Determination of the CTS
142 
143 kc_cts_neu(j,i) = 0
144 
145 do kc=1, kcmax-1
146  if (omega_c_neu(kc,j,i) > eps_omega) then
147  kc_cts_neu(j,i) = kc
148  else
149  exit
150  end if
151 end do
152 
153 !-------- Computation of the ice enthalpy
154 ! (corrector step for the cold-ice domain only
155 ! in order to fulfull the transition condition at the CTS) --------
156 
157 #if (CALCMOD==3) /* ENTM scheme */
158 
159 if (kc_cts_neu(j,i) > 0) then
160 
161 ! ------ Update of the abbreviations where needed
162 
163  do kc=0, kcmax
164  temp_c_val(kc) = temp_c_neu(kc,j,i)
165  omega_c_val(kc) = omega_c_neu(kc,j,i)
166  end do
167 
168  call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
169  i, j, ce6, ce7, ci2, cm3)
170 
171 ! ------ Set-up of the equations
172 
173  call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
174  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
175  adv_vert_sg, abs_adv_vert_sg, &
176  dtime_temp, dtt_dxi, dtt_deta, &
177  dtt_2dxi, dtt_2deta, &
178  i, j, kc_cts_neu(j,i), &
179  lgs_a0, lgs_a1, lgs_a2, lgs_b)
180 
181 ! ------ Solution of the system of linear equations
182 
183  call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
184 
185 ! ------ Assignment of the result
186 
187  kc=kc_cts_neu(j,i)
188  enth_c_neu(kc,j,i) = lgs_x(kc)
189  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
190  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
191 
192  do kc=kc_cts_neu(j,i)+1, kcmax
193  enth_c_neu(kc,j,i) = lgs_x(kc)
194  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
195  omega_c_neu(kc,j,i) = 0.0_dp ! cold-ice domain
196  end do
197 
198 end if
199 
200 #elif (CALCMOD==2) /* ENTC scheme */
201 
202 continue ! no corrector step
203 
204 #else
205 stop ' calc_temp_enth_2: CALCMOD must be either 2 or 3!'
206 #endif
207 
208 !-------- Water drainage from temperate ice (if existing) --------
209 
210 q_tld(j,i) = 0.0_dp
211 
212 do kc=0, kc_cts_neu(j,i)
213 
214  if (omega_c_neu(kc,j,i) > omega_max) then
215 
216  q_tld(j,i) = q_tld(j,i) + cqtlde(kc)*(omega_c_neu(kc,j,i)-omega_max)
217 
218  omega_c_neu(kc,j,i) = omega_max
219  enth_c_neu(kc,j,i) = enth_fct_temp_omega(temp_c_neu(kc,j,i), omega_max)
220 
221  end if
222 
223 end do
224 
225 !-------- Set enthalpy and water content in the redundant,
226 ! lower (kt) ice layer to the value at the ice base --------
227 
228 do kt=0, ktmax
229  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
230  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
231 end do
232 
233 !-------- Computation of the age of ice --------
234 
235 ! ------ Set-up of the the equations
236 
237 call calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
238  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
239  adv_vert_sg, abs_adv_vert_sg, &
240  dtime_temp, dtt_dxi, dtt_deta, &
241  dtt_2dxi, dtt_2deta, &
242  i, j, &
243  lgs_a0, lgs_a1, lgs_a2, lgs_b)
244 
245 ! ------ Solution of the system of linear equations
246 
247 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
248 
249 ! ------ Assignment of the result
250 ! restriction to interval [0, AGE_MAX yr]
251 
252 do kc=0, kcmax
253 
254  age_c_neu(kc,j,i) = lgs_x(kc)
255 
256  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
257  age_c_neu(kc,j,i) = 0.0_dp
258  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
259  age_c_neu(kc,j,i) = age_max*year_sec
260 
261 end do
262 
263 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
264 
265 do kt=0, ktmax
266  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
267 end do
268 
269 end subroutine calc_temp_enth_2
270 
271 !-------------------------------------------------------------------------------
272 !> Computation of temperature and age for an ice column with a temperate base
273 !! with the enthalpy method: Abbreviations I.
274 !<------------------------------------------------------------------------------
275 subroutine calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
276  at4_1, at4_2, at5, atr1, alb1, &
277  ai1, aqtlde, &
278  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
279  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
280  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
281  adv_vert_sg, abs_adv_vert_sg, &
282  ci1, cqtlde, dtt_dxi, dtt_deta)
283 
284 use sico_types
285 use sico_variables
286 use sico_vars
287 
288 implicit none
289 
290 integer(i4b), intent(in) :: i, j
291 real(dp), intent(in) :: at1(0:kcmax), &
292  at2_1(0:kcmax), at2_2(0:kcmax), &
293  at3_1(0:kcmax), at3_2(0:kcmax), &
294  at4_1(0:kcmax), at4_2(0:kcmax), &
295  at5(0:kcmax), ai1(0:kcmax), &
296  atr1, alb1, aqtlde(0:kcmax)
297 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
298 
299 real(dp), intent(out) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
300  ct4(0:kcmax), ce5(0:kcmax), &
301  ctr1, clb1
302 real(dp), intent(out) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
303  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
304  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
305 real(dp), intent(out) :: ci1(0:kcmax), cqtlde(0:kcmax)
306 real(dp), intent(out) :: dtt_dxi, dtt_deta
307 
308 integer(i4b) :: kc
309 
310 !-------- Initialisation --------
311 
312 ct1 = 0.0_dp
313 ct2 = 0.0_dp
314 ct3 = 0.0_dp
315 ct4 = 0.0_dp
316 ce5 = 0.0_dp
317 ctr1 = 0.0_dp
318 clb1 = 0.0_dp
319 ct1_sg = 0.0_dp
320 ct2_sg = 0.0_dp
321 ct3_sg = 0.0_dp
322 ct4_sg = 0.0_dp
323 adv_vert_sg = 0.0_dp
324 abs_adv_vert_sg = 0.0_dp
325 ci1 = 0.0_dp
326 cqtlde = 0.0_dp
327 dtt_dxi = 0.0_dp
328 dtt_deta = 0.0_dp
329 
330 !-------- Actual computation --------
331 
332 ctr1 = atr1
333 clb1 = alb1*q_geo(j,i)
334 
335 #if (ADV_VERT==1)
336 
337 do kc=1, kcmax-1
338  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
339 end do
340 
341 kc=0
342 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
343  ! only needed for kc=0 ...
344 kc=kcmax-1
345 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
346  ! ... and kc=KCMAX-1
347 
348 #elif (ADV_VERT==2 || ADV_VERT==3)
349 
350 do kc=0, kcmax-1
351  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
352 end do
353 
354 #endif
355 
356 do kc=0, kcmax
357 
358  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
359  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
360  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
361  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
362  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
363  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
364  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
365  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
366  ce5(kc) = at5(kc)/h_c(j,i)
367  ci1(kc) = ai1(kc)/h_c(j,i)
368  cqtlde(kc) = aqtlde(kc)*h_c(j,i)
369 
370 end do
371 
372 #if (ADV_VERT==1)
373 
374 kc=0
375 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
376 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
377 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
378 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
379 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
380 kc=kcmax-1
381 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
382 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
383 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
384 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
385 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
386 
387 #elif (ADV_VERT==2 || ADV_VERT==3)
388 
389 do kc=0, kcmax-1
390  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
391  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
392  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
393  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
394  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
395 end do
396 
397 #endif
398 
399 #if (ADV_HOR==3)
400 dtt_dxi = 2.0_dp*dtt_2dxi
401 dtt_deta = 2.0_dp*dtt_2deta
402 #endif
403 
404 end subroutine calc_temp_enth_2_a1
405 
406 !-------------------------------------------------------------------------------
407 !> Computation of temperature and age for an ice column with a temperate base
408 !! with the enthalpy method: Abbreviations II.
409 !<------------------------------------------------------------------------------
410 subroutine calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
411  i, j, ce6, ce7, ci2, cm3)
412 
413 use sico_types
414 use sico_variables
415 use sico_vars
417 
418 implicit none
419 
420 integer(i4b), intent(in) :: i, j
421 real(dp), intent(in) :: at6(0:kcmax), at7, ai2(0:kcmax), am3(0:kcmax)
422 real(dp), intent(in) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
423 
424 real(dp), intent(out) :: ce6(0:kcmax), ce7(0:kcmax), ci2(0:kcmax), &
425  cm3(0:kcmax)
426 
427 integer(i4b) :: kc
428 real(dp) :: temp_c_help(0:kcmax)
429 
430 !-------- Initialisation --------
431 
432 ce6 = 0.0_dp
433 ce7 = 0.0_dp
434 ci2 = 0.0_dp
435 cm3 = 0.0_dp
436 
437 !-------- Actual computation --------
438 
439 do kc=0, kcmax
440 
441 #if (DYNAMICS==2)
442  if (.not.flag_shelfy_stream(j,i)) then
443 #endif
444  ce7(kc) = at7 &
445  *enh_c(kc,j,i) &
446  *ratefac_c_t(temp_c_val(kc), omega_c_val(kc), temp_c_m(kc,j,i)) &
447  *creep(sigma_c(kc,j,i)) &
448  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
449 #if (DYNAMICS==2)
450  else
451  ce7(kc) = 2.0_dp*at7 &
452  *viscosity(de_c(kc,j,i), &
453  temp_c_val(kc), temp_c_m(kc,j,i), omega_c_val(kc), &
454  enh_c(kc,j,i), 2_i2b) &
455  *de_c(kc,j,i)**2
456  end if
457 #endif
458 
459  cm3(kc) = am3(kc)*h_c(j,i)*c_val(temp_c_val(kc))
460 
461 end do
462 
463 do kc=0, kc_cts_neu(j,i)-1 ! temperate layer
464  ce6(kc) = at6(kc) &
465  *nue/h_c(j,i)
466  ci2(kc) = ai2(kc)/h_c(j,i)
467 end do
468 
469 do kc=kc_cts_neu(j,i), kcmax-1 ! cold layer
470  temp_c_help(kc) = 0.5_dp*(temp_c_val(kc)+temp_c_val(kc+1))
471  ce6(kc) = at6(kc) &
472  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
473  ci2(kc) = ai2(kc)/h_c(j,i)
474 end do
475 
476 end subroutine calc_temp_enth_2_a2
477 
478 !-------------------------------------------------------------------------------
479 !> Computation of temperature and age for an ice column with a temperate base
480 !! with the enthalpy method:
481 !! Set-up of the equations for the bedrock temperature.
482 !<------------------------------------------------------------------------------
483 subroutine calc_temp_enth_2_b(ctr1, clb1, i, j, &
484  lgs_a0, lgs_a1, lgs_a2, lgs_b)
485 
486 use sico_types
487 use sico_variables
488 use sico_vars
489 
490 implicit none
491 
492 integer(i4b), intent(in) :: i, j
493 real(dp), intent(in) :: ctr1, clb1
494 
495 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
496  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
497  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
498  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
499 
500 integer(i4b) :: kr
501 
502 !-------- Initialisation --------
503 
504 lgs_a0 = 0.0_dp
505 lgs_a1 = 0.0_dp
506 lgs_a2 = 0.0_dp
507 lgs_b = 0.0_dp
508 
509 !-------- Actual computation --------
510 
511 kr=0
512 lgs_a1(kr) = 1.0_dp
513 lgs_a2(kr) = -1.0_dp
514 lgs_b(kr) = clb1
515 
516 #if (Q_LITHO==1)
517 ! (coupled heat-conducting bedrock)
518 
519 do kr=1, krmax-1
520  lgs_a0(kr) = - ctr1
521  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
522  lgs_a2(kr) = - ctr1
523  lgs_b(kr) = temp_r(kr,j,i)
524 end do
525 
526 #elif (Q_LITHO==0)
527 ! (no coupled heat-conducting bedrock)
528 
529 do kr=1, krmax-1
530  lgs_a0(kr) = 1.0_dp
531  lgs_a1(kr) = 0.0_dp
532  lgs_a2(kr) = -1.0_dp
533  lgs_b(kr) = 2.0_dp*clb1
534 end do
535 
536 #endif
537 
538 kr=krmax
539 lgs_a0(kr) = 0.0_dp
540 lgs_a1(kr) = 1.0_dp
541 lgs_b(kr) = temp_t_m(0,j,i)
542 
543 end subroutine calc_temp_enth_2_b
544 
545 !-------------------------------------------------------------------------------
546 !> Computation of temperature and age for an ice column with a temperate base
547 !! with the enthalpy method:
548 !! Set-up of the equations for the ice enthalpy.
549 !<------------------------------------------------------------------------------
550 subroutine calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
551  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
552  adv_vert_sg, abs_adv_vert_sg, &
553  dtime_temp, dtt_dxi, dtt_deta, &
554  dtt_2dxi, dtt_2deta, &
555  i, j, kcmin, &
556  lgs_a0, lgs_a1, lgs_a2, lgs_b)
557 
558 use sico_types
559 use sico_variables
560 use sico_vars
562 
563 implicit none
564 
565 integer(i4b), intent(in) :: i, j
566 integer(i2b), intent(in) :: kcmin
567 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
568  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
569  ce7(0:kcmax)
570 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
571  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
572  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
573 real(dp), intent(in) :: cm3(0:kcmax)
574 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
575 
576 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
577  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
578  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
579  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
580 
581 integer(i4b) :: kc
582 real(dp) :: vx_c_help, vy_c_help
583 real(dp) :: adv_vert_help
584 
585 !-------- Initialisation --------
586 
587 lgs_a0 = 0.0_dp
588 lgs_a1 = 0.0_dp
589 lgs_a2 = 0.0_dp
590 lgs_b = 0.0_dp
591 
592 !-------- Actual computation --------
593 
594 if (kcmin == 0) then ! predictor step
595 
596  kc=0
597 
598  if (kc_cts_neu(j,i) == 0) then ! temperate base without temperate layer
599 
600  lgs_a1(kc) = 1.0_dp
601  lgs_a2(kc) = 0.0_dp
602  lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i), 0.0_dp)
603 
604  else ! kc_cts_neu(j,i) > 0, temperate base with temperate layer
605 
606  lgs_a1(kc) = 1.0_dp
607  lgs_a2(kc) = -1.0_dp
608  lgs_b(kc) = 0.0_dp
609 
610  end if
611 
612 else ! kcmin > 0, corrector step
613 
614  kc=0
615 
616  lgs_a1(kc) = 1.0_dp ! dummy
617  lgs_a2(kc) = 0.0_dp ! setting,
618  lgs_b(kc) = 0.0_dp ! not needed
619 
620  do kc=1, kcmin-1
621 
622  lgs_a0(kc) = 0.0_dp ! dummy
623  lgs_a1(kc) = 1.0_dp ! setting,
624  lgs_a2(kc) = 0.0_dp ! not
625  lgs_b(kc) = 0.0_dp ! needed
626 
627  end do
628 
629  kc=kcmin
630 
631  lgs_a0(kc) = 0.0_dp
632  lgs_a1(kc) = 1.0_dp
633  lgs_a2(kc) = -1.0_dp
634  lgs_b(kc) = -cm3(kc)
635 
636 end if
637 
638 do kc=kcmin+1, kcmax-1
639 
640 #if (ADV_VERT==1)
641 
642  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
643  -ce5(kc)*ce6(kc-1)
644  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
645  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
646  -ce5(kc)*ce6(kc)
647 
648 #elif (ADV_VERT==2)
649 
650  lgs_a0(kc) &
651  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
652  -ce5(kc)*ce6(kc-1)
653  lgs_a1(kc) &
654  = 1.0_dp &
655  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
656  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
657  +ce5(kc)*(ce6(kc)+ce6(kc-1))
658  lgs_a2(kc) &
659  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
660  -ce5(kc)*ce6(kc)
661 
662 #elif (ADV_VERT==3)
663 
664  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
665 
666  lgs_a0(kc) &
667  = -max(adv_vert_help, 0.0_dp) &
668  -ce5(kc)*ce6(kc-1)
669  lgs_a1(kc) &
670  = 1.0_dp &
671  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
672  +ce5(kc)*(ce6(kc)+ce6(kc-1))
673  lgs_a2(kc) &
674  = min(adv_vert_help, 0.0_dp) &
675  -ce5(kc)*ce6(kc)
676 
677 #endif
678 
679 #if (ADV_HOR==2)
680 
681  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
682  -dtt_2dxi* &
683  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
684  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
685  *insq_g11_sgx(j,i) &
686  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
687  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
688  *insq_g11_sgx(j,i-1) ) &
689  -dtt_2deta* &
690  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
691  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
692  *insq_g22_sgy(j,i) &
693  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
694  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
695  *insq_g22_sgy(j-1,i) )
696 
697 #elif (ADV_HOR==3)
698 
699  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
700  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
701 
702  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
703  -dtt_dxi* &
704  ( min(vx_c_help, 0.0_dp) &
705  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
706  *insq_g11_sgx(j,i) &
707  +max(vx_c_help, 0.0_dp) &
708  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
709  *insq_g11_sgx(j,i-1) ) &
710  -dtt_deta* &
711  ( min(vy_c_help, 0.0_dp) &
712  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
713  *insq_g22_sgy(j,i) &
714  +max(vy_c_help, 0.0_dp) &
715  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
716  *insq_g22_sgy(j-1,i) )
717 
718 #endif
719 
720 end do
721 
722 kc=kcmax
723 lgs_a0(kc) = 0.0_dp
724 lgs_a1(kc) = 1.0_dp
725 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
726  ! zero water content at the ice surface
727 
728 end subroutine calc_temp_enth_2_c
729 
730 !-------------------------------------------------------------------------------
731 !> Computation of temperature and age for an ice column with a temperate base
732 !! with the enthalpy method:
733 !! Set-up of the equations for the age of ice.
734 !<------------------------------------------------------------------------------
735 subroutine calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
736  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
737  adv_vert_sg, abs_adv_vert_sg, &
738  dtime_temp, dtt_dxi, dtt_deta, &
739  dtt_2dxi, dtt_2deta, &
740  i, j, &
741  lgs_a0, lgs_a1, lgs_a2, lgs_b)
742 
743 use sico_types
744 use sico_variables
745 use sico_vars
746 
747 implicit none
748 
749 integer(i4b), intent(in) :: i, j
750 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
751  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
752 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
753  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
754  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
755 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
756 
757 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
758  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
759  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
760  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
761 
762 integer(i4b) :: kc
763 real(dp) :: vx_c_help, vy_c_help
764 real(dp) :: adv_vert_help
765 
766 !-------- Initialisation --------
767 
768 lgs_a0 = 0.0_dp
769 lgs_a1 = 0.0_dp
770 lgs_a2 = 0.0_dp
771 lgs_b = 0.0_dp
772 
773 !-------- Actual computation --------
774 
775 kc=0 ! adv_vert_sg(0) <= 0
776 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
777 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
778 
779 #if (ADV_HOR==2)
780 
781 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
782  -dtt_2dxi* &
783  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
784  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
785  *insq_g11_sgx(j,i) &
786  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
787  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
788  *insq_g11_sgx(j,i-1) ) &
789  -dtt_2deta* &
790  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
791  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
792  *insq_g22_sgy(j,i) &
793  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
794  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
795  *insq_g22_sgy(j-1,i) )
796 
797 #elif (ADV_HOR==3)
798 
799 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
800 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
801 
802 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
803  -dtt_dxi* &
804  ( min(vx_c_help, 0.0_dp) &
805  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
806  *insq_g11_sgx(j,i) &
807  +max(vx_c_help, 0.0_dp) &
808  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
809  *insq_g11_sgx(j,i-1) ) &
810  -dtt_deta* &
811  ( min(vy_c_help, 0.0_dp) &
812  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
813  *insq_g22_sgy(j,i) &
814  +max(vy_c_help, 0.0_dp) &
815  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
816  *insq_g22_sgy(j-1,i) )
817 
818 #endif
819 
820 do kc=1, kcmax-1
821 
822 #if (ADV_VERT==1)
823 
824  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
825  -ci1(kc)*ci2(kc-1)
826  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
827  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
828  -ci1(kc)*ci2(kc)
829 
830 #elif (ADV_VERT==2)
831 
832  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
833  lgs_a1(kc) = 1.0_dp &
834  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
835  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
836  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
837 
838 #elif (ADV_VERT==3)
839 
840  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
841 
842  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
843  lgs_a1(kc) = 1.0_dp &
844  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
845  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
846 
847 #endif
848 
849 #if (ADV_HOR==2)
850 
851  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
852  -dtt_2dxi* &
853  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
854  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
855  *insq_g11_sgx(j,i) &
856  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
857  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
858  *insq_g11_sgx(j,i-1) ) &
859  -dtt_2deta* &
860  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
861  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
862  *insq_g22_sgy(j,i) &
863  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
864  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
865  *insq_g22_sgy(j-1,i) )
866 
867 #elif (ADV_HOR==3)
868 
869  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
870  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
871 
872  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
873  -dtt_dxi* &
874  ( min(vx_c_help, 0.0_dp) &
875  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
876  *insq_g11_sgx(j,i) &
877  +max(vx_c_help, 0.0_dp) &
878  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
879  *insq_g11_sgx(j,i-1) ) &
880  -dtt_deta* &
881  ( min(vy_c_help, 0.0_dp) &
882  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
883  *insq_g22_sgy(j,i) &
884  +max(vy_c_help, 0.0_dp) &
885  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
886  *insq_g22_sgy(j-1,i) )
887 
888 #endif
889 
890 end do
891 
892 kc=kcmax
893 if (as_perp(j,i) >= 0.0_dp) then
894  lgs_a0(kc) = 0.0_dp
895  lgs_a1(kc) = 1.0_dp
896  lgs_b(kc) = 0.0_dp
897 else
898  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
899  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
900  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
901  ! assumed/enforced
902 #if (ADV_HOR==2)
903 
904  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
905  -dtt_2dxi* &
906  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
907  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
908  *insq_g11_sgx(j,i) &
909  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
910  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
911  *insq_g11_sgx(j,i-1) ) &
912  -dtt_2deta* &
913  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
914  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
915  *insq_g22_sgy(j,i) &
916  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
917  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
918  *insq_g22_sgy(j-1,i) )
919 
920 #elif (ADV_HOR==3)
921 
922  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
923  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
924 
925  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
926  -dtt_dxi* &
927  ( min(vx_c_help, 0.0_dp) &
928  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
929  *insq_g11_sgx(j,i) &
930  +max(vx_c_help, 0.0_dp) &
931  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
932  *insq_g11_sgx(j,i-1) ) &
933  -dtt_deta* &
934  ( min(vy_c_help, 0.0_dp) &
935  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
936  *insq_g22_sgy(j,i) &
937  +max(vy_c_help, 0.0_dp) &
938  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
939  *insq_g22_sgy(j-1,i) )
940 
941 #endif
942 
943 end if
944 
945 end subroutine calc_temp_enth_2_d
946 !
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
subroutine calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, atr1, alb1, ai1, aqtlde, dtime_temp, dtt_2dxi, dtt_2deta, i, j, ct1, ct2, ct3, ct4, ce5, ctr1, clb1, ct1_sg, ct2_sg, ct3_sg, ct4_sg, adv_vert_sg, abs_adv_vert_sg, ci1, cqtlde, dtt_dxi, dtt_deta)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
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
real(dp) function, public ratefac_c_t(temp_val, omega_val, temp_m_val)
Rate factor for cold and temperate ice: Combination of ratefac_c and ratefac_t (only for the enthalpy...
real(dp) function, public temp_fct_enth(enth_val, temp_m_val)
Temperature as a function of enthalpy.
subroutine calc_temp_enth_2_b(ctr1, clb1, i, j, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
subroutine calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, alb1, ai1, ai2, aqtlde, am3, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method...
real(dp) function, public omega_fct_enth(enth_val, temp_m_val)
Water content as a function of enthalpy.
subroutine calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, adv_vert_sg, abs_adv_vert_sg, dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta, i, j, kcmin, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
Solvers for systems of linear equations used by SICOPOLIS.
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Declarations of global variables for SICOPOLIS.
real(dp) function creep(sigma_val)
Creep response function for ice.
Definition: creep.F90:35
subroutine calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, ct1_sg, ct2_sg, ct3_sg, ct4_sg, adv_vert_sg, abs_adv_vert_sg, dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta, i, j, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, i, j, ce6, ce7, ci2, cm3)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.