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