SICOPOLIS V3.3
calc_temp_enth_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ t e m p _ e n t h _ m
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age with the enthalpy method.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2013-2017 Ralf Greve, Heinz Blatter
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, water content and age with the enthalpy method.
34 !<------------------------------------------------------------------------------
35 module calc_temp_enth_m
36 
37  use sico_types_m
38  use sico_variables_m
39  use sico_vars_m
40 
41  implicit none
42 
43  private
44  public :: calc_temp_enth
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Main subroutine of calc_temp_enth_m:
50 !! Computation of temperature, water content and age with the enthalpy method.
51 !<------------------------------------------------------------------------------
52 subroutine calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
53  dtime_temp)
54 
55 implicit none
56 
57 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
58 real(dp), intent(in) :: dtime_temp
59 
60 integer(i4b) :: i, j, kc, kt, kr, ii, jj
61 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
62  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
63  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
64  ai1(0:kcmax), ai2(0:kcmax), &
65  atr1, acb1, acb2, acb3, acb4, alb1, aqtlde(0:kcmax), &
66  am1, am3(0:kcmax)
67 real(dp) :: dtt_2dxi, dtt_2deta
68 
69 !-------- Term abbreviations
70 
71 at7 = 2.0_dp/rho*dtime_temp
72 
73 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
74 
75 if (flag_aa_nonzero) then
76  am1 = aa*beta*dzeta_c/(ea-1.0_dp)
77 else
78  am1 = beta*dzeta_c
79 end if
80 
81 if (flag_aa_nonzero) then
82  acb1 = (ea-1.0_dp)/aa/dzeta_c
83 else
84  acb1 = 1.0_dp/dzeta_c
85 end if
86 
87 acb2 = kappa_r/h_r/dzeta_r
88 acb3 = rho*g
89 acb4 = rho*g
90 
91 alb1 = h_r/kappa_r*dzeta_r
92 
93 dtt_2dxi = 0.5_dp*dtime_temp/dxi
94 dtt_2deta = 0.5_dp*dtime_temp/deta
95 
96 do kc=0, kcmax
97 
98  if (flag_aa_nonzero) then
99 
100  at1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
101  at2_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
102  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
103  *dtime_temp/dzeta_c
104  at3_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
105  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
106  *dtime_temp/dzeta_c
107  at4_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
108  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
109  *dtime_temp/dzeta_c
110  at5(kc) = (ea-1.0_dp)/(rho*aa*eaz_c(kc)) &
111  *dtime_temp/dzeta_c
112  if (kc /= kcmax) then
113  at6(kc) = (ea-1.0_dp) &
114  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
115  /dzeta_c
116  else
117  at6(kc) = 0.0_dp
118  end if
119  ai1(kc) = agediff*(ea-1.0_dp)/(aa*eaz_c(kc)) &
120  *dtime_temp/dzeta_c
121  if (kc /= kcmax) then
122  ai2(kc) = (ea-1.0_dp) &
123  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
124  /dzeta_c
125  else
126  ai2(kc) = 0.0_dp
127  end if
128  aqtlde(kc) = (aa*eaz_c(kc))/(ea-1.0_dp)*dzeta_c/dtime_temp
129  am3(kc) = (aa*eaz_c(kc))/(ea-1.0_dp)*dzeta_c*beta
130 
131  else
132 
133  at1(kc) = dtime_temp/dzeta_c
134  at2_1(kc) = dtime_temp/dzeta_c
135  at2_2(kc) = zeta_c(kc) &
136  *dtime_temp/dzeta_c
137  at3_1(kc) = dtime_temp/dzeta_c
138  at3_2(kc) = zeta_c(kc) &
139  *dtime_temp/dzeta_c
140  at4_1(kc) = dtime_temp/dzeta_c
141  at4_2(kc) = zeta_c(kc) &
142  *dtime_temp/dzeta_c
143  at5(kc) = 1.0_dp/rho &
144  *dtime_temp/dzeta_c
145  if (kc /= kcmax) then
146  at6(kc) = 1.0_dp &
147  /dzeta_c
148  else
149  at6(kc) = 0.0_dp
150  end if
151  ai1(kc) = agediff &
152  *dtime_temp/dzeta_c
153  if (kc /= kcmax) then
154  ai2(kc) = 1.0_dp &
155  /dzeta_c
156  else
157  ai2(kc) = 0.0_dp
158  end if
159  aqtlde(kc) = dzeta_c/dtime_temp
160  am3(kc) = dzeta_c*beta
161 
162  end if
163 
164 end do
165 
166 !-------- Computation loop --------
167 
168 do i=1, imax-1 ! skipping domain margins
169 do j=1, jmax-1 ! skipping domain margins
170 
171  if (maske(j,i)==0) then ! glaciated land
172 
173 ! ------ Old vertical column cold
174 
175  if (n_cts(j,i) == -1) then
176 
177  n_cts_neu(j,i) = -1
178  kc_cts_neu(j,i) = 0
179  zm_neu(j,i) = zb(j,i)
180  h_c_neu(j,i) = h_c(j,i)
181  h_t_neu(j,i) = 0.0_dp
182 
183  call calc_temp_enth_1(at1, at2_1, at2_2, at3_1, at3_2, &
184  at4_1, at4_2, at5, at6, at7, &
185  atr1, acb1, acb2, acb3, acb4, alb1, &
186  ai1, ai2, &
187  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
188 
189 ! ---- Check whether base has become temperate
190 
191  if (temp_c_neu(0,j,i) > temp_c_m(0,j,i)-eps) then
192 
193  n_cts_neu(j,i) = 0
194  kc_cts_neu(j,i) = 0
195 
196  call calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
197  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
198  ai1, ai2, aqtlde, am3, &
199  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
200 
201  end if
202 
203 ! ------ Old vertical column with temperate base
204 
205  else if (n_cts(j,i) == 0) then
206 
207  n_cts_neu(j,i) = 0
208  kc_cts_neu(j,i) = kc_cts(j,i)
209  zm_neu(j,i) = zb(j,i)
210  h_c_neu(j,i) = h_c(j,i)
211  h_t_neu(j,i) = h_t(j,i)
212 
213  call calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
214  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
215  ai1, ai2, aqtlde, am3, &
216  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
217 
218 ! ---- Check whether temperate base becomes cold
219 
220  if ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) < (am1*h_c(j,i)) ) then
221 
222  n_cts_neu(j,i) = -1
223  kc_cts_neu(j,i) = 0
224 
225  call calc_temp_enth_1(at1, at2_1, at2_2, at3_1, at3_2, &
226  at4_1, at4_2, at5, at6, at7, &
227  atr1, acb1, acb2, acb3, acb4, alb1, &
228  ai1, ai2, &
229  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
230 
231  if (temp_c_neu(0,j,i) > temp_c_m(0,j,i)-eps) then
232 
233  n_cts_neu(j,i) = 0
234  kc_cts_neu(j,i) = 0
235 
236  call calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
237  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
238  ai1, ai2, aqtlde, am3, &
239  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
240 
241  end if
242 
243  end if
244 
245  end if
246 
247 #if (MARGIN==3)
248 
249  else if (maske(j,i)==3) then ! floating ice
250 
251  n_cts_neu(j,i) = -1
252  kc_cts_neu(j,i) = 0
253  zm_neu(j,i) = zb(j,i)
254  h_c_neu(j,i) = h_c(j,i)
255  h_t_neu(j,i) = 0.0_dp
256 
257  call calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
258  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
259  ai1, ai2, &
260  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
261 
262 ! ------ Reset temperatures above melting to the melting point
263 ! and water contents above zero to zero
264 ! (should not occur, but just in case)
265 
266  do kc=0, kcmax
267  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
268  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
269  if (omega_c_neu(kc,j,i) > 0.0_dp) &
270  omega_c_neu(kc,j,i) = 0.0_dp
271  end do
272 
273 #endif
274 
275  else ! maske(j,i) == 1,2 (ice-free land or sea point)
276 
277  n_cts_neu(j,i) = -1
278  kc_cts_neu(j,i) = 0
279  zm_neu(j,i) = zb(j,i)
280  h_c_neu(j,i) = h_c(j,i)
281  h_t_neu(j,i) = 0.0_dp
282 
283  call calc_temp_enth_r(atr1, alb1, i, j)
284 
285  end if
286 
287 end do
288 end do
289 
290 !-------- Extrapolate values on margins --------
291 
292 ! ------ Lower left corner
293 
294 i=0
295 j=0
296 
297 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
298  ! glaciated land or floating ice
299  ii=i+1
300  jj=j+1
301 
302  do kc=0, kcmax
303  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
304  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
305  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
306  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
307  end do
308 
309  do kt=0, ktmax
310  ! redundant, lower (kt) ice layer
311  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
312  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
313  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
314  end do
315 
316  do kr=0, krmax
317  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
318  end do
319 
320  n_cts_neu(j,i) = n_cts_neu(jj,ii)
321  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
322  zm_neu(j,i) = zb(j,i)
323  h_c_neu(j,i) = h_c(j,i)
324  h_t_neu(j,i) = h_t(j,i)
325 
326 else ! maske(j,i) == 1,2 (ice-free land or sea point)
327 
328  n_cts_neu(j,i) = -1
329  kc_cts_neu(j,i) = 0
330  zm_neu(j,i) = zb(j,i)
331  h_c_neu(j,i) = h_c(j,i)
332  h_t_neu(j,i) = h_t(j,i)
333 
334  call calc_temp_enth_r(atr1, alb1, i, j)
335 
336 end if
337 
338 ! ------ Lower right corner
339 
340 i=imax
341 j=0
342 
343 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
344  ! glaciated land or floating ice
345  ii=i-1
346  jj=j+1
347 
348  do kc=0, kcmax
349  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
350  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
351  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
352  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
353  end do
354 
355  do kt=0, ktmax
356  ! redundant, lower (kt) ice layer
357  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
358  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
359  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
360  end do
361 
362  do kr=0, krmax
363  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
364  end do
365 
366  n_cts_neu(j,i) = n_cts_neu(jj,ii)
367  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
368  zm_neu(j,i) = zb(j,i)
369  h_c_neu(j,i) = h_c(j,i)
370  h_t_neu(j,i) = h_t(j,i)
371 
372 else ! maske(j,i) == 1,2 (ice-free land or sea point)
373 
374  n_cts_neu(j,i) = -1
375  kc_cts_neu(j,i) = 0
376  zm_neu(j,i) = zb(j,i)
377  h_c_neu(j,i) = h_c(j,i)
378  h_t_neu(j,i) = h_t(j,i)
379 
380  call calc_temp_enth_r(atr1, alb1, i, j)
381 
382 end if
383 
384 ! ------ Upper left corner
385 
386 i=0
387 j=jmax
388 
389 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
390  ! glaciated land or floating ice
391  ii=i+1
392  jj=j-1
393 
394  do kc=0, kcmax
395  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
396  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
397  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
398  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
399  end do
400 
401  do kt=0, ktmax
402  ! redundant, lower (kt) ice layer
403  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
404  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
405  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
406  end do
407 
408  do kr=0, krmax
409  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
410  end do
411 
412  n_cts_neu(j,i) = n_cts_neu(jj,ii)
413  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
414  zm_neu(j,i) = zb(j,i)
415  h_c_neu(j,i) = h_c(j,i)
416  h_t_neu(j,i) = h_t(j,i)
417 
418 else ! maske(j,i) == 1,2 (ice-free land or sea point)
419 
420  n_cts_neu(j,i) = -1
421  kc_cts_neu(j,i) = 0
422  zm_neu(j,i) = zb(j,i)
423  h_c_neu(j,i) = h_c(j,i)
424  h_t_neu(j,i) = h_t(j,i)
425 
426  call calc_temp_enth_r(atr1, alb1, i, j)
427 
428 end if
429 
430 ! ------ Upper right corner
431 
432 i=imax
433 j=jmax
434 
435 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
436  ! glaciated land or floating ice
437  ii=i-1
438  jj=j-1
439 
440  do kc=0, kcmax
441  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
442  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
443  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
444  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
445  end do
446 
447  do kt=0, ktmax
448  ! redundant, lower (kt) ice layer
449  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
450  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
451  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
452  end do
453 
454  do kr=0, krmax
455  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
456  end do
457 
458  n_cts_neu(j,i) = n_cts_neu(jj,ii)
459  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
460  zm_neu(j,i) = zb(j,i)
461  h_c_neu(j,i) = h_c(j,i)
462  h_t_neu(j,i) = h_t(j,i)
463 
464 else ! maske(j,i) == 1,2 (ice-free land or sea point)
465 
466  n_cts_neu(j,i) = -1
467  kc_cts_neu(j,i) = 0
468  zm_neu(j,i) = zb(j,i)
469  h_c_neu(j,i) = h_c(j,i)
470  h_t_neu(j,i) = h_t(j,i)
471 
472  call calc_temp_enth_r(atr1, alb1, i, j)
473 
474 end if
475 
476 ! ------ Lower and upper margins
477 
478 do i=1, imax-1
479 
480 ! ---- Lower margin
481 
482  j=0
483 
484  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
485  ! glaciated land or floating ice
486  ii=i
487  jj=j+1
488 
489  do kc=0, kcmax
490  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
491  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
492  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
493  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
494  end do
495 
496  do kt=0, ktmax
497  ! redundant, lower (kt) ice layer
498  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
499  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
500  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
501  end do
502 
503  do kr=0, krmax
504  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
505  end do
506 
507  n_cts_neu(j,i) = n_cts_neu(jj,ii)
508  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
509  zm_neu(j,i) = zb(j,i)
510  h_c_neu(j,i) = h_c(j,i)
511  h_t_neu(j,i) = h_t(j,i)
512 
513  else ! maske(j,i) == 1,2 (ice-free land or sea point)
514 
515  n_cts_neu(j,i) = -1
516  kc_cts_neu(j,i) = 0
517  zm_neu(j,i) = zb(j,i)
518  h_c_neu(j,i) = h_c(j,i)
519  h_t_neu(j,i) = h_t(j,i)
520 
521  call calc_temp_enth_r(atr1, alb1, i, j)
522 
523  end if
524 
525 ! ---- Upper margin
526 
527  j=jmax
528 
529  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
530  ! glaciated land or floating ice
531  ii=i
532  jj=j-1
533 
534  do kc=0, kcmax
535  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
536  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
537  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
538  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
539  end do
540 
541  do kt=0, ktmax
542  ! redundant, lower (kt) ice layer
543  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
544  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
545  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
546  end do
547 
548  do kr=0, krmax
549  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
550  end do
551 
552  n_cts_neu(j,i) = n_cts_neu(jj,ii)
553  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
554  zm_neu(j,i) = zb(j,i)
555  h_c_neu(j,i) = h_c(j,i)
556  h_t_neu(j,i) = h_t(j,i)
557 
558  else ! maske(j,i) == 1,2 (ice-free land or sea point)
559 
560  n_cts_neu(j,i) = -1
561  kc_cts_neu(j,i) = 0
562  zm_neu(j,i) = zb(j,i)
563  h_c_neu(j,i) = h_c(j,i)
564  h_t_neu(j,i) = h_t(j,i)
565 
566  call calc_temp_enth_r(atr1, alb1, i, j)
567 
568  end if
569 
570 end do
571 
572 ! ------ Left and right margins
573 
574 do j=1, jmax-1
575 
576 ! ---- Left margin
577 
578  i=0
579 
580  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
581  ! glaciated land or floating ice
582  ii=i+1
583  jj=j
584 
585  do kc=0, kcmax
586  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
587  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
588  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
589  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
590  end do
591 
592  do kt=0, ktmax
593  ! redundant, lower (kt) ice layer
594  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
595  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
596  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
597  end do
598 
599  do kr=0, krmax
600  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
601  end do
602 
603  n_cts_neu(j,i) = n_cts_neu(jj,ii)
604  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
605  zm_neu(j,i) = zb(j,i)
606  h_c_neu(j,i) = h_c(j,i)
607  h_t_neu(j,i) = h_t(j,i)
608 
609  else ! maske(j,i) == 1,2 (ice-free land or sea point)
610 
611  n_cts_neu(j,i) = -1
612  kc_cts_neu(j,i) = 0
613  zm_neu(j,i) = zb(j,i)
614  h_c_neu(j,i) = h_c(j,i)
615  h_t_neu(j,i) = h_t(j,i)
616 
617  call calc_temp_enth_r(atr1, alb1, i, j)
618 
619  end if
620 
621 ! ---- Right margin
622 
623  i=imax
624 
625  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
626  ! glaciated land or floating ice
627  ii=i-1
628  jj=j
629 
630  do kc=0, kcmax
631  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
632  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
633  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
634  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
635  end do
636 
637  do kt=0, ktmax
638  ! redundant, lower (kt) ice layer
639  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
640  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
641  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
642  end do
643 
644  do kr=0, krmax
645  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
646  end do
647 
648  n_cts_neu(j,i) = n_cts_neu(jj,ii)
649  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
650  zm_neu(j,i) = zb(j,i)
651  h_c_neu(j,i) = h_c(j,i)
652  h_t_neu(j,i) = h_t(j,i)
653 
654  else ! maske(j,i) == 1,2 (ice-free land or sea point)
655 
656  n_cts_neu(j,i) = -1
657  kc_cts_neu(j,i) = 0
658  zm_neu(j,i) = zb(j,i)
659  h_c_neu(j,i) = h_c(j,i)
660  h_t_neu(j,i) = h_t(j,i)
661 
662  call calc_temp_enth_r(atr1, alb1, i, j)
663 
664  end if
665 
666 end do
667 
668 end subroutine calc_temp_enth
669 
670 !-------------------------------------------------------------------------------
671 !> Computation of temperature and age for a cold ice column with the
672 !! enthalpy method.
673 !<------------------------------------------------------------------------------
674 subroutine calc_temp_enth_1(at1, at2_1, at2_2, at3_1, at3_2, &
675  at4_1, at4_2, at5, at6, at7, &
676  atr1, acb1, acb2, acb3, acb4, alb1, &
677  ai1, ai2, &
678  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
679 
680 use sico_maths_m, only : tri_sle
681 use enth_temp_omega_m, only : temp_fct_enth
682 
683 implicit none
684 
685 integer(i4b), intent(in) :: i, j
686 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
687  at3_1(0:kcmax), at3_2(0:kcmax), &
688  at4_1(0:kcmax), at4_2(0:kcmax), &
689  at5(0:kcmax), at6(0:kcmax), at7, &
690  ai1(0:kcmax), ai2(0:kcmax), &
691  atr1, acb1, acb2, acb3, acb4, alb1
692 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
693 
694 integer(i4b) :: kc, kt, kr
695 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
696  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), &
697  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1
698 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
699  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
700 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
701 real(dp) :: dtt_dxi, dtt_deta
702 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
703  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
704  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
705  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
706  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
707 
708 !-------- Check for boundary points --------
709 
710 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
711  stop ' >>> calc_temp_enth_1: Boundary points not allowed.'
712 
713 !-------- Abbreviations --------
714 
715 call calc_temp_enth_1_a(at1, at2_1, at2_2, at3_1, at3_2, &
716  at4_1, at4_2, at5, at6, at7, &
717  atr1, acb1, acb2, acb3, acb4, alb1, &
718  ai1, ai2, &
719  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
720  ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
721  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1, &
722  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
723  adv_vert_sg, abs_adv_vert_sg, &
724  ci1, ci2, dtt_dxi, dtt_deta)
725 
726 !-------- Computation of the bedrock temperature
727 ! (upper boundary condition: old temperature at the ice base) --------
728 
729 ! ------ Set-up of the the equations
730 
731 call calc_temp_enth_1_b(ctr1, clb1, i, j, &
732  lgs_a0, lgs_a1, lgs_a2, lgs_b)
733 
734 ! ------ Solution of the system of linear equations
735 
736 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
737 
738 ! ------ Assignment of the result (predictor values)
739 
740 do kr=0, krmax
741  temp_r_neu(kr,j,i) = lgs_x(kr)
742 end do
743 
744 !-------- Computation of the ice enthalpy --------
745 
746 ! ------ Set-up of the the equations
747 
748 call calc_temp_enth_1_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
749  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
750  ccbe1, ccb2, ccb3, ccb4, &
751  adv_vert_sg, abs_adv_vert_sg, &
752  dtime_temp, dtt_dxi, dtt_deta, &
753  dtt_2dxi, dtt_2deta, &
754  i, j, &
755  lgs_a0, lgs_a1, lgs_a2, lgs_b)
756 
757 ! ------ Solution of the system of linear equations
758 
759 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
760 
761 ! ------ Assignment of the result
762 
763 do kc=0, kcmax
764  enth_c_neu(kc,j,i) = lgs_x(kc)
765  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
766  omega_c_neu(kc,j,i) = 0.0_dp ! solution is supposed to be for cold ice
767 end do
768 
769 !-------- Water drainage from the non-existing temperate ice --------
770 
771 q_tld(j,i) = 0.0_dp
772 
773 !-------- Set enthalpy and water content in the redundant,
774 ! lower (kt) ice layer to the value at the ice base --------
775 
776 do kt=0, ktmax
777  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
778  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
779 end do
780 
781 !-------- Computation of the age of ice --------
782 
783 ! ------ Set-up of the the equations
784 
785 call calc_temp_enth_1_d(ct1, ct2, ct3, ct4, ci1, ci2, &
786  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
787  adv_vert_sg, abs_adv_vert_sg, &
788  dtime_temp, dtt_dxi, dtt_deta, &
789  dtt_2dxi, dtt_2deta, &
790  i, j, &
791  lgs_a0, lgs_a1, lgs_a2, lgs_b)
792 
793 ! ------ Solution of the system of linear equations
794 
795 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
796 
797 ! ------ Assignment of the result,
798 ! restriction to interval [0, AGE_MAX yr] --------
799 
800 do kc=0, kcmax
801 
802  age_c_neu(kc,j,i) = lgs_x(kc)
803 
804  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
805  age_c_neu(kc,j,i) = 0.0_dp
806  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
807  age_c_neu(kc,j,i) = age_max*year_sec
808 
809 end do
810 
811 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
812 
813 do kt=0, ktmax
814  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
815 end do
816 
817 end subroutine calc_temp_enth_1
818 
819 !-------------------------------------------------------------------------------
820 !> Computation of temperature and age for a cold ice column with the
821 !! enthalpy method: Abbreviations.
822 !<------------------------------------------------------------------------------
823 subroutine calc_temp_enth_1_a(at1, at2_1, at2_2, at3_1, at3_2, &
824  at4_1, at4_2, at5, at6, at7, &
825  atr1, acb1, acb2, acb3, acb4, alb1, &
826  ai1, ai2, &
827  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
828  ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
829  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1, &
830  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
831  adv_vert_sg, abs_adv_vert_sg, &
832  ci1, ci2, dtt_dxi, dtt_deta)
833 
834 use ice_material_properties_m, only : ratefac_c_t, kappa_val, c_val, &
835  creep, viscosity
836 
837 implicit none
838 
839 integer(i4b), intent(in) :: i, j
840 real(dp), intent(in) :: at1(0:kcmax), &
841  at2_1(0:kcmax), at2_2(0:kcmax), &
842  at3_1(0:kcmax), at3_2(0:kcmax), &
843  at4_1(0:kcmax), at4_2(0:kcmax), &
844  at5(0:kcmax), at6(0:kcmax), at7, &
845  ai1(0:kcmax), ai2(0:kcmax), &
846  atr1, acb1, acb2, acb3, acb4, alb1
847 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
848 
849 real(dp), intent(out) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
850  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
851  ce7(0:kcmax), &
852  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1
853 real(dp), intent(out) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
854  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
855  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
856 real(dp), intent(out) :: ci1(0:kcmax), ci2(0:kcmax)
857 real(dp), intent(out) :: dtt_dxi, dtt_deta
858 
859 integer(i4b) :: kc
860 real(dp) :: temp_c_help(0:kcmax)
861 
862 !-------- Initialisation --------
863 
864 ct1 = 0.0_dp
865 ct2 = 0.0_dp
866 ct3 = 0.0_dp
867 ct4 = 0.0_dp
868 ce5 = 0.0_dp
869 ce6 = 0.0_dp
870 ce7 = 0.0_dp
871 ctr1 = 0.0_dp
872 clb1 = 0.0_dp
873 ct1_sg = 0.0_dp
874 ct2_sg = 0.0_dp
875 ct3_sg = 0.0_dp
876 ct4_sg = 0.0_dp
877 adv_vert_sg = 0.0_dp
878 abs_adv_vert_sg = 0.0_dp
879 ci1 = 0.0_dp
880 ci2 = 0.0_dp
881 dtt_dxi = 0.0_dp
882 dtt_deta = 0.0_dp
883 
884 !-------- Actual computation --------
885 
886 ctr1 = atr1
887 
888 ccbe1 = acb1 &
889  *kappa_val(temp_c(0,j,i)) &
890  /(c_val(temp_c(0,j,i))*h_c(j,i))
891 ccb2 = acb2
892 
893 #if (DYNAMICS==2)
894 if (.not.flag_shelfy_stream(j,i)) then
895 #endif
896 
897  ccb3 = acb3*0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
898  *h_c(j,i)*dzs_dxi_g(j,i)
899  ccb4 = acb4*0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
900  *h_c(j,i)*dzs_deta_g(j,i)
901 
902 #if (DYNAMICS==2)
903 else ! flag_shelfy_stream(j,i) == .true.
904 
905  ccb3 = -c_drag(j,i) &
906  * sqrt(vx_b_g(j,i)**2 &
907  +vy_b_g(j,i)**2) &
908  **(1.0_dp+p_weert_inv(j,i))
909  ccb4 = 0.0_dp
910 
911 end if
912 #endif
913 
914 clb1 = alb1*q_geo(j,i)
915 
916 #if (ADV_VERT==1)
917 
918 do kc=1, kcmax-1
919  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
920 end do
921 
922 kc=0
923 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
924  ! only needed for kc=0 ...
925 kc=kcmax-1
926 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
927  ! ... and kc=KCMAX-1
928 
929 #elif (ADV_VERT==2 || ADV_VERT==3)
930 
931 do kc=0, kcmax-1
932  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
933 end do
934 
935 #endif
936 
937 do kc=0, kcmax
938 
939  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
940  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
941  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
942  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
943  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
944  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
945  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
946  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
947  ce5(kc) = at5(kc)/h_c(j,i)
948 
949 #if (DYNAMICS==2)
950  if (.not.flag_shelfy_stream(j,i)) then
951 #endif
952  ce7(kc) = at7 &
953  *enh_c(kc,j,i) &
954  *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), temp_c_m(kc,j,i)) &
955  *creep(sigma_c(kc,j,i)) &
956  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
957 #if (DYNAMICS==2)
958  else
959  ce7(kc) = 2.0_dp*at7 &
960  *viscosity(de_c(kc,j,i), &
961  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
962  enh_c(kc,j,i), 0_i2b) &
963  *de_c(kc,j,i)**2
964  end if
965 #endif
966 
967  ci1(kc) = ai1(kc)/h_c(j,i)
968 
969 end do
970 
971 #if (ADV_VERT==1)
972 
973 kc=0
974 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
975 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
976 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
977 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
978 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
979 kc=kcmax-1
980 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
981 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
982 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
983 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
984 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
985 
986 #elif (ADV_VERT==2 || ADV_VERT==3)
987 
988 do kc=0, kcmax-1
989  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
990  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
991  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
992  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
993  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
994 end do
995 
996 #endif
997 
998 do kc=0, kcmax-1
999  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
1000  ce6(kc) = at6(kc) &
1001  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
1002  ci2(kc) = ai2(kc)/h_c(j,i)
1003 end do
1004 
1005 #if (ADV_HOR==3)
1006 dtt_dxi = 2.0_dp*dtt_2dxi
1007 dtt_deta = 2.0_dp*dtt_2deta
1008 #endif
1009 
1010 end subroutine calc_temp_enth_1_a
1011 
1012 !-------------------------------------------------------------------------------
1013 !> Computation of temperature and age for a cold ice column with the
1014 !! enthalpy method:
1015 !! Set-up of the equations for the bedrock temperature.
1016 !<------------------------------------------------------------------------------
1017 subroutine calc_temp_enth_1_b(ctr1, clb1, i, j, &
1018  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1019 
1020 implicit none
1021 
1022 integer(i4b), intent(in) :: i, j
1023 real(dp), intent(in) :: ctr1, clb1
1024 
1025 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1026  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1027  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1028  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1029 
1030 integer(i4b) :: kr
1031 
1032 !-------- Initialisation --------
1033 
1034 lgs_a0 = 0.0_dp
1035 lgs_a1 = 0.0_dp
1036 lgs_a2 = 0.0_dp
1037 lgs_b = 0.0_dp
1038 
1039 !-------- Actual computation --------
1040 
1041 kr=0
1042 lgs_a1(kr) = 1.0_dp
1043 lgs_a2(kr) = -1.0_dp
1044 lgs_b(kr) = clb1
1045 
1046 #if (Q_LITHO==1)
1047 ! (coupled heat-conducting bedrock)
1048 
1049 do kr=1, krmax-1
1050  lgs_a0(kr) = -ctr1
1051  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1052  lgs_a2(kr) = -ctr1
1053  lgs_b(kr) = temp_r(kr,j,i)
1054 end do
1055 
1056 #elif (Q_LITHO==0)
1057 ! (no coupled heat-conducting bedrock)
1058 
1059 do kr=1, krmax-1
1060  lgs_a0(kr) = 1.0_dp
1061  lgs_a1(kr) = 0.0_dp
1062  lgs_a2(kr) = -1.0_dp
1063  lgs_b(kr) = 2.0_dp*clb1
1064 end do
1065 
1066 #endif
1067 
1068 kr=krmax
1069 lgs_a0(kr) = 0.0_dp
1070 lgs_a1(kr) = 1.0_dp
1071 lgs_b(kr) = temp_c(0,j,i)
1072 
1073 end subroutine calc_temp_enth_1_b
1074 
1075 !-------------------------------------------------------------------------------
1076 !> Computation of temperature and age for a cold ice column with the
1077 !! enthalpy method:
1078 !! Set-up of the equations for the ice enthalpy.
1079 !<------------------------------------------------------------------------------
1080 subroutine calc_temp_enth_1_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1081  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1082  ccbe1, ccb2, ccb3, ccb4, &
1083  adv_vert_sg, abs_adv_vert_sg, &
1084  dtime_temp, dtt_dxi, dtt_deta, &
1085  dtt_2dxi, dtt_2deta, &
1086  i, j, &
1087  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1088 
1089 use enth_temp_omega_m, only : enth_fct_temp_omega
1090 
1091 implicit none
1092 
1093 integer(i4b), intent(in) :: i, j
1094 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1095  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
1096  ce7(0:kcmax)
1097 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1098  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1099  ccbe1, ccb2, ccb3, ccb4, &
1100  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1101 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1102 
1103 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1104  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1105  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1106  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1107 
1108 integer(i4b) :: kc, kr
1109 real(dp) :: vx_c_help, vy_c_help
1110 real(dp) :: adv_vert_help
1111 
1112 !-------- Initialisation --------
1113 
1114 lgs_a0 = 0.0_dp
1115 lgs_a1 = 0.0_dp
1116 lgs_a2 = 0.0_dp
1117 lgs_b = 0.0_dp
1118 
1119 !-------- Actual computation --------
1120 
1121 kr=krmax
1122 kc=0
1123 lgs_a1(kc) = -ccbe1
1124 lgs_a2(kc) = ccbe1
1125 lgs_b(kc) = ccb2*(temp_r_neu(kr,j,i)-temp_r_neu(kr-1,j,i)) + ccb3 + ccb4
1126 
1127 do kc=1, kcmax-1
1128 
1129 #if (ADV_VERT==1)
1130 
1131  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1132  -ce5(kc)*ce6(kc-1)
1133  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
1134  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1135  -ce5(kc)*ce6(kc)
1136 
1137 #elif (ADV_VERT==2)
1138 
1139  lgs_a0(kc) &
1140  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1141  -ce5(kc)*ce6(kc-1)
1142  lgs_a1(kc) &
1143  = 1.0_dp &
1144  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1145  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1146  +ce5(kc)*(ce6(kc)+ce6(kc-1))
1147  lgs_a2(kc) &
1148  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1149  -ce5(kc)*ce6(kc)
1150 
1151 #elif (ADV_VERT==3)
1152 
1153  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1154 
1155  lgs_a0(kc) &
1156  = -max(adv_vert_help, 0.0_dp) &
1157  -ce5(kc)*ce6(kc-1)
1158  lgs_a1(kc) &
1159  = 1.0_dp &
1160  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
1161  +ce5(kc)*(ce6(kc)+ce6(kc-1))
1162  lgs_a2(kc) &
1163  = min(adv_vert_help, 0.0_dp) &
1164  -ce5(kc)*ce6(kc)
1165 
1166 #endif
1167 
1168 #if (ADV_HOR==2)
1169 
1170  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
1171  -dtt_2dxi* &
1172  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1173  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
1174  *insq_g11_sgx(j,i) &
1175  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1176  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
1177  *insq_g11_sgx(j,i-1) ) &
1178  -dtt_2deta* &
1179  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1180  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
1181  *insq_g22_sgy(j,i) &
1182  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1183  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
1184  *insq_g22_sgy(j-1,i) )
1185 
1186 #elif (ADV_HOR==3)
1187 
1188  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1189  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1190 
1191  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
1192  -dtt_dxi* &
1193  ( min(vx_c_help, 0.0_dp) &
1194  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
1195  *insq_g11_sgx(j,i) &
1196  +max(vx_c_help, 0.0_dp) &
1197  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
1198  *insq_g11_sgx(j,i-1) ) &
1199  -dtt_deta* &
1200  ( min(vy_c_help, 0.0_dp) &
1201  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
1202  *insq_g22_sgy(j,i) &
1203  +max(vy_c_help, 0.0_dp) &
1204  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
1205  *insq_g22_sgy(j-1,i) )
1206 
1207 #endif
1208 
1209 end do
1210 
1211 kc=kcmax
1212 lgs_a0(kc) = 0.0_dp
1213 lgs_a1(kc) = 1.0_dp
1214 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
1215  ! zero water content at the ice surface
1216 
1217 end subroutine calc_temp_enth_1_c
1218 
1219 !-------------------------------------------------------------------------------
1220 !> Computation of temperature and age for a cold ice column with the
1221 !! enthalpy method:
1222 !! Set-up of the equations for the age of ice.
1223 !<------------------------------------------------------------------------------
1224 subroutine calc_temp_enth_1_d(ct1, ct2, ct3, ct4, ci1, ci2, &
1225  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1226  adv_vert_sg, abs_adv_vert_sg, &
1227  dtime_temp, dtt_dxi, dtt_deta, &
1228  dtt_2dxi, dtt_2deta, &
1229  i, j, &
1230  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1231 
1232 implicit none
1233 
1234 integer(i4b), intent(in) :: i, j
1235 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1236  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
1237 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1238  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1239  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1240 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1241 
1242 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1243  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1244  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1245  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1246 
1247 integer(i4b) :: kc
1248 real(dp) :: vx_c_help, vy_c_help
1249 real(dp) :: adv_vert_help
1250 
1251 !-------- Initialisation --------
1252 
1253 lgs_a0 = 0.0_dp
1254 lgs_a1 = 0.0_dp
1255 lgs_a2 = 0.0_dp
1256 lgs_b = 0.0_dp
1257 
1258 !-------- Actual computation --------
1259 
1260 kc=0 ! adv_vert_sg(0) <= 0
1261 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
1262 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
1263 
1264 #if (ADV_HOR==2)
1265 
1266 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1267  -dtt_2dxi* &
1268  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1269  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1270  *insq_g11_sgx(j,i) &
1271  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1272  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1273  *insq_g11_sgx(j,i-1) ) &
1274  -dtt_2deta* &
1275  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1276  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1277  *insq_g22_sgy(j,i) &
1278  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1279  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1280  *insq_g22_sgy(j-1,i) )
1281 
1282 #elif (ADV_HOR==3)
1283 
1284 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1285 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1286 
1287 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1288  -dtt_dxi* &
1289  ( min(vx_c_help, 0.0_dp) &
1290  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1291  *insq_g11_sgx(j,i) &
1292  +max(vx_c_help, 0.0_dp) &
1293  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1294  *insq_g11_sgx(j,i-1) ) &
1295  -dtt_deta* &
1296  ( min(vy_c_help, 0.0_dp) &
1297  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1298  *insq_g22_sgy(j,i) &
1299  +max(vy_c_help, 0.0_dp) &
1300  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1301  *insq_g22_sgy(j-1,i) )
1302 
1303 #endif
1304 
1305 do kc=1, kcmax-1
1306 
1307 #if (ADV_VERT==1)
1308 
1309  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1310  -ci1(kc)*ci2(kc-1)
1311  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
1312  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1313  -ci1(kc)*ci2(kc)
1314 
1315 #elif (ADV_VERT==2)
1316 
1317  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
1318  lgs_a1(kc) = 1.0_dp &
1319  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1320  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1321  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1322 
1323 #elif (ADV_VERT==3)
1324 
1325  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1326 
1327  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
1328  lgs_a1(kc) = 1.0_dp &
1329  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
1330  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
1331 
1332 #endif
1333 
1334 #if (ADV_HOR==2)
1335 
1336  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1337  -dtt_2dxi* &
1338  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1339  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1340  *insq_g11_sgx(j,i) &
1341  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1342  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1343  *insq_g11_sgx(j,i-1) ) &
1344  -dtt_2deta* &
1345  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1346  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1347  *insq_g22_sgy(j,i) &
1348  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1349  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1350  *insq_g22_sgy(j-1,i) )
1351 
1352 #elif (ADV_HOR==3)
1353 
1354  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1355  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1356 
1357  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1358  -dtt_dxi* &
1359  ( min(vx_c_help, 0.0_dp) &
1360  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1361  *insq_g11_sgx(j,i) &
1362  +max(vx_c_help, 0.0_dp) &
1363  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1364  *insq_g11_sgx(j,i-1) ) &
1365  -dtt_deta* &
1366  ( min(vy_c_help, 0.0_dp) &
1367  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1368  *insq_g22_sgy(j,i) &
1369  +max(vy_c_help, 0.0_dp) &
1370  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1371  *insq_g22_sgy(j-1,i) )
1372 
1373 #endif
1374 
1375 end do
1376 
1377 kc=kcmax
1378 if (as_perp(j,i) >= 0.0_dp) then
1379  lgs_a0(kc) = 0.0_dp
1380  lgs_a1(kc) = 1.0_dp
1381  lgs_b(kc) = 0.0_dp
1382 else
1383  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1384  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1385  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
1386  ! assumed/enforced
1387 #if (ADV_HOR==2)
1388 
1389  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1390  -dtt_2dxi* &
1391  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1392  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1393  *insq_g11_sgx(j,i) &
1394  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1395  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1396  *insq_g11_sgx(j,i-1) ) &
1397  -dtt_2deta* &
1398  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1399  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1400  *insq_g22_sgy(j,i) &
1401  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1402  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1403  *insq_g22_sgy(j-1,i) )
1404 
1405 #elif (ADV_HOR==3)
1406 
1407  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1408  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1409 
1410  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1411  -dtt_dxi* &
1412  ( min(vx_c_help, 0.0_dp) &
1413  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1414  *insq_g11_sgx(j,i) &
1415  +max(vx_c_help, 0.0_dp) &
1416  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1417  *insq_g11_sgx(j,i-1) ) &
1418  -dtt_deta* &
1419  ( min(vy_c_help, 0.0_dp) &
1420  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1421  *insq_g22_sgy(j,i) &
1422  +max(vy_c_help, 0.0_dp) &
1423  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1424  *insq_g22_sgy(j-1,i) )
1425 
1426 #endif
1427 
1428 end if
1429 
1430 end subroutine calc_temp_enth_1_d
1431 
1432 !-------------------------------------------------------------------------------
1433 !> Computation of temperature and age for an ice column with a temperate base
1434 !! with the enthalpy method.
1435 !<------------------------------------------------------------------------------
1436 subroutine calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
1437  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1438  ai1, ai2, aqtlde, am3, &
1439  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1440 
1441 use sico_maths_m, only : tri_sle
1442 use enth_temp_omega_m, only : enth_fct_temp_omega, &
1443  temp_fct_enth, omega_fct_enth
1444 
1445 implicit none
1446 
1447 integer(i4b), intent(in) :: i, j
1448 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1449  at3_1(0:kcmax), at3_2(0:kcmax), &
1450  at4_1(0:kcmax), at4_2(0:kcmax), &
1451  at5(0:kcmax), at6(0:kcmax), at7, &
1452  ai1(0:kcmax), ai2(0:kcmax), &
1453  atr1, alb1, aqtlde(0:kcmax), am3(0:kcmax)
1454 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1455 
1456 integer(i4b) :: kc, kt, kr
1457 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1458  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
1459 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1460  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1461 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1462 real(dp) :: cqtlde(0:kcmax), cm3(0:kcmax)
1463 real(dp) :: dtt_dxi, dtt_deta
1464 real(dp) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
1465 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1466  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1467  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1468  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1469  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1470 
1471 real(dp), parameter :: eps_omega=1.0e-12_dp
1472 
1473 !-------- Check for boundary points --------
1474 
1475 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
1476  stop ' >>> calc_temp_enth_2: Boundary points not allowed.'
1477 
1478 !-------- Abbreviations --------
1479 
1480 call calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
1481  at4_1, at4_2, at5, atr1, alb1, &
1482  ai1, aqtlde, &
1483  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
1484  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
1485  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1486  adv_vert_sg, abs_adv_vert_sg, &
1487  ci1, cqtlde, dtt_dxi, dtt_deta)
1488 
1489 do kc=0, kcmax
1490  temp_c_val(kc) = temp_c(kc,j,i)
1491  omega_c_val(kc) = omega_c(kc,j,i)
1492 end do
1493 
1494 call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1495  i, j, ce6, ce7, ci2, cm3)
1496 
1497 !-------- Computation of the bedrock temperature --------
1498 
1499 ! ------ Set-up of the the equations
1500 
1501 call calc_temp_enth_2_b(ctr1, clb1, i, j, &
1502  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1503 
1504 ! ------ Solution of the system of linear equations
1505 
1506 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
1507 
1508 ! ------ Assignment of the result
1509 
1510 do kr=0, krmax
1511  temp_r_neu(kr,j,i) = lgs_x(kr)
1512 end do
1513 
1514 !-------- Computation of the ice enthalpy (predictor step) --------
1515 
1516 ! ------ Set-up of the equations
1517 
1518 call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1519  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1520  adv_vert_sg, abs_adv_vert_sg, &
1521  dtime_temp, dtt_dxi, dtt_deta, &
1522  dtt_2dxi, dtt_2deta, &
1523  i, j, 0_i2b, &
1524  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1525 
1526 ! ------ Solution of the system of linear equations
1527 
1528 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1529 
1530 ! ------ Assignment of the result
1531 
1532 do kc=0, kcmax
1533  enth_c_neu(kc,j,i) = lgs_x(kc)
1534  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1535  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1536 end do
1537 
1538 ! ------ Determination of the CTS
1539 
1540 kc_cts_neu(j,i) = 0
1541 
1542 do kc=1, kcmax-1
1543  if (omega_c_neu(kc,j,i) > eps_omega) then
1544  kc_cts_neu(j,i) = kc
1545  else
1546  exit
1547  end if
1548 end do
1549 
1550 !-------- Computation of the ice enthalpy
1551 ! (corrector step for the cold-ice domain only
1552 ! in order to fulfull the transition condition at the CTS) --------
1553 
1554 #if (CALCMOD==3) /* ENTM scheme */
1555 
1556 if (kc_cts_neu(j,i) > 0) then
1557 
1558 ! ------ Update of the abbreviations where needed
1559 
1560  do kc=0, kcmax
1561  temp_c_val(kc) = temp_c_neu(kc,j,i)
1562  omega_c_val(kc) = omega_c_neu(kc,j,i)
1563  end do
1564 
1565  call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1566  i, j, ce6, ce7, ci2, cm3)
1567 
1568 ! ------ Set-up of the equations
1569 
1570  call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1571  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1572  adv_vert_sg, abs_adv_vert_sg, &
1573  dtime_temp, dtt_dxi, dtt_deta, &
1574  dtt_2dxi, dtt_2deta, &
1575  i, j, kc_cts_neu(j,i), &
1576  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1577 
1578 ! ------ Solution of the system of linear equations
1579 
1580  call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1581 
1582 ! ------ Assignment of the result
1583 
1584  kc=kc_cts_neu(j,i)
1585  enth_c_neu(kc,j,i) = lgs_x(kc)
1586  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1587  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1588 
1589  do kc=kc_cts_neu(j,i)+1, kcmax
1590  enth_c_neu(kc,j,i) = lgs_x(kc)
1591  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1592  omega_c_neu(kc,j,i) = 0.0_dp ! cold-ice domain
1593  end do
1594 
1595 end if
1596 
1597 #elif (CALCMOD==2) /* ENTC scheme */
1598 
1599 !!! continue ! no corrector step
1600 
1601 #else
1602 stop ' >>> calc_temp_enth_2: CALCMOD must be either 2 or 3!'
1603 #endif
1604 
1605 !-------- Water drainage from temperate ice (if existing) --------
1606 
1607 q_tld(j,i) = 0.0_dp
1608 
1609 do kc=0, kc_cts_neu(j,i)
1610 
1611  if (omega_c_neu(kc,j,i) > omega_max) then
1612 
1613  q_tld(j,i) = q_tld(j,i) + cqtlde(kc)*(omega_c_neu(kc,j,i)-omega_max)
1614 
1615  omega_c_neu(kc,j,i) = omega_max
1616  enth_c_neu(kc,j,i) = enth_fct_temp_omega(temp_c_neu(kc,j,i), omega_max)
1617 
1618  end if
1619 
1620 end do
1621 
1622 !-------- Set enthalpy and water content in the redundant,
1623 ! lower (kt) ice layer to the value at the ice base --------
1624 
1625 do kt=0, ktmax
1626  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
1627  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
1628 end do
1629 
1630 !-------- Computation of the age of ice --------
1631 
1632 ! ------ Set-up of the the equations
1633 
1634 call calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
1635  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1636  adv_vert_sg, abs_adv_vert_sg, &
1637  dtime_temp, dtt_dxi, dtt_deta, &
1638  dtt_2dxi, dtt_2deta, &
1639  i, j, &
1640  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1641 
1642 ! ------ Solution of the system of linear equations
1643 
1644 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1645 
1646 ! ------ Assignment of the result
1647 ! restriction to interval [0, AGE_MAX yr]
1648 
1649 do kc=0, kcmax
1650 
1651  age_c_neu(kc,j,i) = lgs_x(kc)
1652 
1653  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
1654  age_c_neu(kc,j,i) = 0.0_dp
1655  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
1656  age_c_neu(kc,j,i) = age_max*year_sec
1657 
1658 end do
1659 
1660 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
1661 
1662 do kt=0, ktmax
1663  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
1664 end do
1665 
1666 end subroutine calc_temp_enth_2
1667 
1668 !-------------------------------------------------------------------------------
1669 !> Computation of temperature and age for an ice column with a temperate base
1670 !! with the enthalpy method: Abbreviations I.
1671 !<------------------------------------------------------------------------------
1672 subroutine calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
1673  at4_1, at4_2, at5, atr1, alb1, &
1674  ai1, aqtlde, &
1675  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
1676  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
1677  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1678  adv_vert_sg, abs_adv_vert_sg, &
1679  ci1, cqtlde, dtt_dxi, dtt_deta)
1680 
1681 implicit none
1682 
1683 integer(i4b), intent(in) :: i, j
1684 real(dp), intent(in) :: at1(0:kcmax), &
1685  at2_1(0:kcmax), at2_2(0:kcmax), &
1686  at3_1(0:kcmax), at3_2(0:kcmax), &
1687  at4_1(0:kcmax), at4_2(0:kcmax), &
1688  at5(0:kcmax), ai1(0:kcmax), &
1689  atr1, alb1, aqtlde(0:kcmax)
1690 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1691 
1692 real(dp), intent(out) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1693  ct4(0:kcmax), ce5(0:kcmax), &
1694  ctr1, clb1
1695 real(dp), intent(out) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1696  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1697  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1698 real(dp), intent(out) :: ci1(0:kcmax), cqtlde(0:kcmax)
1699 real(dp), intent(out) :: dtt_dxi, dtt_deta
1700 
1701 integer(i4b) :: kc
1702 
1703 !-------- Initialisation --------
1704 
1705 ct1 = 0.0_dp
1706 ct2 = 0.0_dp
1707 ct3 = 0.0_dp
1708 ct4 = 0.0_dp
1709 ce5 = 0.0_dp
1710 ctr1 = 0.0_dp
1711 clb1 = 0.0_dp
1712 ct1_sg = 0.0_dp
1713 ct2_sg = 0.0_dp
1714 ct3_sg = 0.0_dp
1715 ct4_sg = 0.0_dp
1716 adv_vert_sg = 0.0_dp
1717 abs_adv_vert_sg = 0.0_dp
1718 ci1 = 0.0_dp
1719 cqtlde = 0.0_dp
1720 dtt_dxi = 0.0_dp
1721 dtt_deta = 0.0_dp
1722 
1723 !-------- Actual computation --------
1724 
1725 ctr1 = atr1
1726 clb1 = alb1*q_geo(j,i)
1727 
1728 #if (ADV_VERT==1)
1729 
1730 do kc=1, kcmax-1
1731  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
1732 end do
1733 
1734 kc=0
1735 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1736  ! only needed for kc=0 ...
1737 kc=kcmax-1
1738 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1739  ! ... and kc=KCMAX-1
1740 
1741 #elif (ADV_VERT==2 || ADV_VERT==3)
1742 
1743 do kc=0, kcmax-1
1744  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1745 end do
1746 
1747 #endif
1748 
1749 do kc=0, kcmax
1750 
1751  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
1752  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
1753  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
1754  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
1755  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
1756  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
1757  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
1758  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
1759  ce5(kc) = at5(kc)/h_c(j,i)
1760  ci1(kc) = ai1(kc)/h_c(j,i)
1761  cqtlde(kc) = aqtlde(kc)*h_c(j,i)
1762 
1763 end do
1764 
1765 #if (ADV_VERT==1)
1766 
1767 kc=0
1768 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1769 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1770 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1771 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1772 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
1773 kc=kcmax-1
1774 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1775 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1776 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1777 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1778 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1779 
1780 #elif (ADV_VERT==2 || ADV_VERT==3)
1781 
1782 do kc=0, kcmax-1
1783  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1784  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1785  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1786  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1787  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1788 end do
1789 
1790 #endif
1791 
1792 #if (ADV_HOR==3)
1793 dtt_dxi = 2.0_dp*dtt_2dxi
1794 dtt_deta = 2.0_dp*dtt_2deta
1795 #endif
1796 
1797 end subroutine calc_temp_enth_2_a1
1798 
1799 !-------------------------------------------------------------------------------
1800 !> Computation of temperature and age for an ice column with a temperate base
1801 !! with the enthalpy method: Abbreviations II.
1802 !<------------------------------------------------------------------------------
1803 subroutine calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1804  i, j, ce6, ce7, ci2, cm3)
1805 
1806 use ice_material_properties_m, only : ratefac_c_t, kappa_val, c_val, &
1807  creep, viscosity
1808 
1809 implicit none
1810 
1811 integer(i4b), intent(in) :: i, j
1812 real(dp), intent(in) :: at6(0:kcmax), at7, ai2(0:kcmax), am3(0:kcmax)
1813 real(dp), intent(in) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
1814 
1815 real(dp), intent(out) :: ce6(0:kcmax), ce7(0:kcmax), ci2(0:kcmax), &
1816  cm3(0:kcmax)
1817 
1818 integer(i4b) :: kc
1819 real(dp) :: temp_c_help(0:kcmax)
1820 
1821 !-------- Initialisation --------
1822 
1823 ce6 = 0.0_dp
1824 ce7 = 0.0_dp
1825 ci2 = 0.0_dp
1826 cm3 = 0.0_dp
1827 
1828 !-------- Actual computation --------
1829 
1830 do kc=0, kcmax
1831 
1832 #if (DYNAMICS==2)
1833  if (.not.flag_shelfy_stream(j,i)) then
1834 #endif
1835  ce7(kc) = at7 &
1836  *enh_c(kc,j,i) &
1837  *ratefac_c_t(temp_c_val(kc), omega_c_val(kc), temp_c_m(kc,j,i)) &
1838  *creep(sigma_c(kc,j,i)) &
1839  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
1840 #if (DYNAMICS==2)
1841  else
1842  ce7(kc) = 2.0_dp*at7 &
1843  *viscosity(de_c(kc,j,i), &
1844  temp_c_val(kc), temp_c_m(kc,j,i), omega_c_val(kc), &
1845  enh_c(kc,j,i), 2_i2b) &
1846  *de_c(kc,j,i)**2
1847  end if
1848 #endif
1849 
1850  cm3(kc) = am3(kc)*h_c(j,i)*c_val(temp_c_val(kc))
1851 
1852 end do
1853 
1854 do kc=0, kc_cts_neu(j,i)-1 ! temperate layer
1855  ce6(kc) = at6(kc) &
1856  *nue/h_c(j,i)
1857  ci2(kc) = ai2(kc)/h_c(j,i)
1858 end do
1859 
1860 do kc=kc_cts_neu(j,i), kcmax-1 ! cold layer
1861  temp_c_help(kc) = 0.5_dp*(temp_c_val(kc)+temp_c_val(kc+1))
1862  ce6(kc) = at6(kc) &
1863  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
1864  ci2(kc) = ai2(kc)/h_c(j,i)
1865 end do
1866 
1867 end subroutine calc_temp_enth_2_a2
1868 
1869 !-------------------------------------------------------------------------------
1870 !> Computation of temperature and age for an ice column with a temperate base
1871 !! with the enthalpy method:
1872 !! Set-up of the equations for the bedrock temperature.
1873 !<------------------------------------------------------------------------------
1874 subroutine calc_temp_enth_2_b(ctr1, clb1, i, j, &
1875  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1876 
1877 implicit none
1878 
1879 integer(i4b), intent(in) :: i, j
1880 real(dp), intent(in) :: ctr1, clb1
1881 
1882 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1883  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1884  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1885  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1886 
1887 integer(i4b) :: kr
1888 
1889 !-------- Initialisation --------
1890 
1891 lgs_a0 = 0.0_dp
1892 lgs_a1 = 0.0_dp
1893 lgs_a2 = 0.0_dp
1894 lgs_b = 0.0_dp
1895 
1896 !-------- Actual computation --------
1897 
1898 kr=0
1899 lgs_a1(kr) = 1.0_dp
1900 lgs_a2(kr) = -1.0_dp
1901 lgs_b(kr) = clb1
1902 
1903 #if (Q_LITHO==1)
1904 ! (coupled heat-conducting bedrock)
1905 
1906 do kr=1, krmax-1
1907  lgs_a0(kr) = - ctr1
1908  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1909  lgs_a2(kr) = - ctr1
1910  lgs_b(kr) = temp_r(kr,j,i)
1911 end do
1912 
1913 #elif (Q_LITHO==0)
1914 ! (no coupled heat-conducting bedrock)
1915 
1916 do kr=1, krmax-1
1917  lgs_a0(kr) = 1.0_dp
1918  lgs_a1(kr) = 0.0_dp
1919  lgs_a2(kr) = -1.0_dp
1920  lgs_b(kr) = 2.0_dp*clb1
1921 end do
1922 
1923 #endif
1924 
1925 kr=krmax
1926 lgs_a0(kr) = 0.0_dp
1927 lgs_a1(kr) = 1.0_dp
1928 lgs_b(kr) = temp_t_m(0,j,i)
1929 
1930 end subroutine calc_temp_enth_2_b
1931 
1932 !-------------------------------------------------------------------------------
1933 !> Computation of temperature and age for an ice column with a temperate base
1934 !! with the enthalpy method:
1935 !! Set-up of the equations for the ice enthalpy.
1936 !<------------------------------------------------------------------------------
1937 subroutine calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1938  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1939  adv_vert_sg, abs_adv_vert_sg, &
1940  dtime_temp, dtt_dxi, dtt_deta, &
1941  dtt_2dxi, dtt_2deta, &
1942  i, j, kcmin, &
1943  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1944 
1945 use enth_temp_omega_m, only : enth_fct_temp_omega
1946 
1947 implicit none
1948 
1949 integer(i4b), intent(in) :: i, j
1950 integer(i2b), intent(in) :: kcmin
1951 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1952  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
1953  ce7(0:kcmax)
1954 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1955  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1956  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1957 real(dp), intent(in) :: cm3(0:kcmax)
1958 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1959 
1960 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1961  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1962  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1963  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1964 
1965 integer(i4b) :: kc
1966 real(dp) :: vx_c_help, vy_c_help
1967 real(dp) :: adv_vert_help
1968 
1969 !-------- Initialisation --------
1970 
1971 lgs_a0 = 0.0_dp
1972 lgs_a1 = 0.0_dp
1973 lgs_a2 = 0.0_dp
1974 lgs_b = 0.0_dp
1975 
1976 !-------- Actual computation --------
1977 
1978 if (kcmin == 0) then ! predictor step
1979 
1980  kc=0
1981 
1982  if (kc_cts_neu(j,i) == 0) then ! temperate base without temperate layer
1983 
1984  lgs_a1(kc) = 1.0_dp
1985  lgs_a2(kc) = 0.0_dp
1986  lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i), 0.0_dp)
1987 
1988  else ! kc_cts_neu(j,i) > 0, temperate base with temperate layer
1989 
1990  lgs_a1(kc) = 1.0_dp
1991  lgs_a2(kc) = -1.0_dp
1992  lgs_b(kc) = 0.0_dp
1993 
1994  end if
1995 
1996 else ! kcmin > 0, corrector step
1997 
1998  kc=0
1999 
2000  lgs_a1(kc) = 1.0_dp ! dummy
2001  lgs_a2(kc) = 0.0_dp ! setting,
2002  lgs_b(kc) = 0.0_dp ! not needed
2003 
2004  do kc=1, kcmin-1
2005 
2006  lgs_a0(kc) = 0.0_dp ! dummy
2007  lgs_a1(kc) = 1.0_dp ! setting,
2008  lgs_a2(kc) = 0.0_dp ! not
2009  lgs_b(kc) = 0.0_dp ! needed
2010 
2011  end do
2012 
2013  kc=kcmin
2014 
2015  lgs_a0(kc) = 0.0_dp
2016  lgs_a1(kc) = 1.0_dp
2017  lgs_a2(kc) = -1.0_dp
2018  lgs_b(kc) = -cm3(kc)
2019 
2020 end if
2021 
2022 do kc=kcmin+1, kcmax-1
2023 
2024 #if (ADV_VERT==1)
2025 
2026  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2027  -ce5(kc)*ce6(kc-1)
2028  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
2029  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2030  -ce5(kc)*ce6(kc)
2031 
2032 #elif (ADV_VERT==2)
2033 
2034  lgs_a0(kc) &
2035  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2036  -ce5(kc)*ce6(kc-1)
2037  lgs_a1(kc) &
2038  = 1.0_dp &
2039  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2040  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2041  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2042  lgs_a2(kc) &
2043  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2044  -ce5(kc)*ce6(kc)
2045 
2046 #elif (ADV_VERT==3)
2047 
2048  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2049 
2050  lgs_a0(kc) &
2051  = -max(adv_vert_help, 0.0_dp) &
2052  -ce5(kc)*ce6(kc-1)
2053  lgs_a1(kc) &
2054  = 1.0_dp &
2055  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2056  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2057  lgs_a2(kc) &
2058  = min(adv_vert_help, 0.0_dp) &
2059  -ce5(kc)*ce6(kc)
2060 
2061 #endif
2062 
2063 #if (ADV_HOR==2)
2064 
2065  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2066  -dtt_2dxi* &
2067  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2068  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2069  *insq_g11_sgx(j,i) &
2070  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2071  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2072  *insq_g11_sgx(j,i-1) ) &
2073  -dtt_2deta* &
2074  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2075  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2076  *insq_g22_sgy(j,i) &
2077  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2078  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2079  *insq_g22_sgy(j-1,i) )
2080 
2081 #elif (ADV_HOR==3)
2082 
2083  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2084  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2085 
2086  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2087  -dtt_dxi* &
2088  ( min(vx_c_help, 0.0_dp) &
2089  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2090  *insq_g11_sgx(j,i) &
2091  +max(vx_c_help, 0.0_dp) &
2092  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2093  *insq_g11_sgx(j,i-1) ) &
2094  -dtt_deta* &
2095  ( min(vy_c_help, 0.0_dp) &
2096  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2097  *insq_g22_sgy(j,i) &
2098  +max(vy_c_help, 0.0_dp) &
2099  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2100  *insq_g22_sgy(j-1,i) )
2101 
2102 #endif
2103 
2104 end do
2105 
2106 kc=kcmax
2107 lgs_a0(kc) = 0.0_dp
2108 lgs_a1(kc) = 1.0_dp
2109 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2110  ! zero water content at the ice surface
2111 
2112 end subroutine calc_temp_enth_2_c
2113 
2114 !-------------------------------------------------------------------------------
2115 !> Computation of temperature and age for an ice column with a temperate base
2116 !! with the enthalpy method:
2117 !! Set-up of the equations for the age of ice.
2118 !<------------------------------------------------------------------------------
2119 subroutine calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
2120  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
2121  adv_vert_sg, abs_adv_vert_sg, &
2122  dtime_temp, dtt_dxi, dtt_deta, &
2123  dtt_2dxi, dtt_2deta, &
2124  i, j, &
2125  lgs_a0, lgs_a1, lgs_a2, lgs_b)
2126 
2127 implicit none
2128 
2129 integer(i4b), intent(in) :: i, j
2130 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
2131  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
2132 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
2133  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
2134  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2135 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
2136 
2137 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2138  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2139  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2140  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2141 
2142 integer(i4b) :: kc
2143 real(dp) :: vx_c_help, vy_c_help
2144 real(dp) :: adv_vert_help
2145 
2146 !-------- Initialisation --------
2147 
2148 lgs_a0 = 0.0_dp
2149 lgs_a1 = 0.0_dp
2150 lgs_a2 = 0.0_dp
2151 lgs_b = 0.0_dp
2152 
2153 !-------- Actual computation --------
2154 
2155 kc=0 ! adv_vert_sg(0) <= 0
2156 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2157 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2158 
2159 #if (ADV_HOR==2)
2160 
2161 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2162  -dtt_2dxi* &
2163  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2164  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2165  *insq_g11_sgx(j,i) &
2166  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2167  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2168  *insq_g11_sgx(j,i-1) ) &
2169  -dtt_2deta* &
2170  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2171  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2172  *insq_g22_sgy(j,i) &
2173  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2174  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2175  *insq_g22_sgy(j-1,i) )
2176 
2177 #elif (ADV_HOR==3)
2178 
2179 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2180 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2181 
2182 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2183  -dtt_dxi* &
2184  ( min(vx_c_help, 0.0_dp) &
2185  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2186  *insq_g11_sgx(j,i) &
2187  +max(vx_c_help, 0.0_dp) &
2188  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2189  *insq_g11_sgx(j,i-1) ) &
2190  -dtt_deta* &
2191  ( min(vy_c_help, 0.0_dp) &
2192  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2193  *insq_g22_sgy(j,i) &
2194  +max(vy_c_help, 0.0_dp) &
2195  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2196  *insq_g22_sgy(j-1,i) )
2197 
2198 #endif
2199 
2200 do kc=1, kcmax-1
2201 
2202 #if (ADV_VERT==1)
2203 
2204  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2205  -ci1(kc)*ci2(kc-1)
2206  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2207  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2208  -ci1(kc)*ci2(kc)
2209 
2210 #elif (ADV_VERT==2)
2211 
2212  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2213  lgs_a1(kc) = 1.0_dp &
2214  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2215  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2216  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2217 
2218 #elif (ADV_VERT==3)
2219 
2220  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2221 
2222  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2223  lgs_a1(kc) = 1.0_dp &
2224  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2225  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2226 
2227 #endif
2228 
2229 #if (ADV_HOR==2)
2230 
2231  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2232  -dtt_2dxi* &
2233  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2234  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2235  *insq_g11_sgx(j,i) &
2236  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2237  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2238  *insq_g11_sgx(j,i-1) ) &
2239  -dtt_2deta* &
2240  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2241  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2242  *insq_g22_sgy(j,i) &
2243  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2244  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2245  *insq_g22_sgy(j-1,i) )
2246 
2247 #elif (ADV_HOR==3)
2248 
2249  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2250  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2251 
2252  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2253  -dtt_dxi* &
2254  ( min(vx_c_help, 0.0_dp) &
2255  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2256  *insq_g11_sgx(j,i) &
2257  +max(vx_c_help, 0.0_dp) &
2258  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2259  *insq_g11_sgx(j,i-1) ) &
2260  -dtt_deta* &
2261  ( min(vy_c_help, 0.0_dp) &
2262  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2263  *insq_g22_sgy(j,i) &
2264  +max(vy_c_help, 0.0_dp) &
2265  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2266  *insq_g22_sgy(j-1,i) )
2267 
2268 #endif
2269 
2270 end do
2271 
2272 kc=kcmax
2273 if (as_perp(j,i) >= 0.0_dp) then
2274  lgs_a0(kc) = 0.0_dp
2275  lgs_a1(kc) = 1.0_dp
2276  lgs_b(kc) = 0.0_dp
2277 else
2278  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2279  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2280  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2281  ! assumed/enforced
2282 #if (ADV_HOR==2)
2283 
2284  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2285  -dtt_2dxi* &
2286  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2287  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2288  *insq_g11_sgx(j,i) &
2289  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2290  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2291  *insq_g11_sgx(j,i-1) ) &
2292  -dtt_2deta* &
2293  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2294  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2295  *insq_g22_sgy(j,i) &
2296  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2297  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2298  *insq_g22_sgy(j-1,i) )
2299 
2300 #elif (ADV_HOR==3)
2301 
2302  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2303  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2304 
2305  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2306  -dtt_dxi* &
2307  ( min(vx_c_help, 0.0_dp) &
2308  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2309  *insq_g11_sgx(j,i) &
2310  +max(vx_c_help, 0.0_dp) &
2311  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2312  *insq_g11_sgx(j,i-1) ) &
2313  -dtt_deta* &
2314  ( min(vy_c_help, 0.0_dp) &
2315  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2316  *insq_g22_sgy(j,i) &
2317  +max(vy_c_help, 0.0_dp) &
2318  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2319  *insq_g22_sgy(j-1,i) )
2320 
2321 #endif
2322 
2323 end if
2324 
2325 end subroutine calc_temp_enth_2_d
2326 
2327 !-------------------------------------------------------------------------------
2328 !> Computation of temperature, age, water content and enthalpy for an
2329 !! ice-free column.
2330 !<------------------------------------------------------------------------------
2331 subroutine calc_temp_enth_r(atr1, alb1, i, j)
2332 
2333 use sico_maths_m, only : tri_sle
2334 use enth_temp_omega_m, only : enth_fct_temp_omega
2335 
2336 implicit none
2337 
2338 integer(i4b), intent(in) :: i, j
2339 real(dp), intent(in) :: atr1, alb1
2340 
2341 integer(i4b) :: kc, kt, kr
2342 real(dp) :: ctr1, clb1
2343 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2344  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2345  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2346  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2347  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2348 real(dp) :: enth_val
2349 
2350 !-------- Abbreviations --------
2351 
2352 ctr1 = atr1
2353 clb1 = alb1*q_geo(j,i)
2354 
2355 !-------- Set up the equations for the bedrock temperature --------
2356 
2357 kr=0
2358 lgs_a1(kr) = 1.0_dp
2359 lgs_a2(kr) = -1.0_dp
2360 lgs_b(kr) = clb1
2361 
2362 #if (Q_LITHO==1)
2363 ! (coupled heat-conducting bedrock)
2364 
2365 do kr=1, krmax-1
2366  lgs_a0(kr) = - ctr1
2367  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2368  lgs_a2(kr) = - ctr1
2369  lgs_b(kr) = temp_r(kr,j,i)
2370 end do
2371 
2372 #elif (Q_LITHO==0)
2373 ! (no coupled heat-conducting bedrock)
2374 
2375 do kr=1, krmax-1
2376  lgs_a0(kr) = 1.0_dp
2377  lgs_a1(kr) = 0.0_dp
2378  lgs_a2(kr) = -1.0_dp
2379  lgs_b(kr) = 2.0_dp*clb1
2380 end do
2381 
2382 #endif
2383 
2384 kr=krmax
2385 lgs_a0(kr) = 0.0_dp
2386 lgs_a1(kr) = 1.0_dp
2387 lgs_b(kr) = temp_s(j,i)
2388 
2389 !-------- Solve system of linear equations --------
2390 
2391 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2392 
2393 !-------- Assign the result --------
2394 
2395 do kr=0, krmax
2396  temp_r_neu(kr,j,i) = lgs_x(kr)
2397 end do
2398 
2399 !-------- Water content, age and enthalpy
2400 ! in the non-existing lower (kt) ice layer --------
2401 
2402 enth_val = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2403 
2404 do kt=0, ktmax
2405  omega_t_neu(kt,j,i) = 0.0_dp
2406  age_t_neu(kt,j,i) = 0.0_dp
2407  enth_t_neu(kt,j,i) = enth_val
2408 end do
2409 
2410 !-------- Temperature, age, water content and enthalpy
2411 ! in the non-existing upper (kc) ice layer --------
2412 
2413 do kc=0, kcmax
2414  temp_c_neu(kc,j,i) = temp_s(j,i)
2415  age_c_neu(kc,j,i) = 0.0_dp
2416  omega_c_neu(kc,j,i) = 0.0_dp
2417  enth_c_neu(kc,j,i) = enth_val
2418 end do
2419 
2420 end subroutine calc_temp_enth_r
2421 
2422 !-------------------------------------------------------------------------------
2423 !> Computation of temperature and age for ice shelves (floating ice)
2424 !! with the enthalpy method.
2425 !<------------------------------------------------------------------------------
2426 subroutine calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
2427  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
2428  ai1, ai2, &
2429  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
2430 
2431 use ice_material_properties_m, only : kappa_val, c_val, viscosity
2432 use sico_maths_m, only : tri_sle
2433 use enth_temp_omega_m, only : enth_fct_temp_omega, &
2434  temp_fct_enth, omega_fct_enth
2435 
2436 implicit none
2437 
2438 integer(i4b), intent(in) :: i, j
2439 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
2440  at3_1(0:kcmax), at3_2(0:kcmax), &
2441  at4_1(0:kcmax), at4_2(0:kcmax), &
2442  at5(0:kcmax), at6(0:kcmax), at7, &
2443  ai1(0:kcmax), ai2(0:kcmax), &
2444  atr1, alb1
2445 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
2446 
2447 integer(i4b) :: kc, kt, kr
2448 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
2449  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
2450 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
2451  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2452 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
2453 real(dp) :: temp_c_help(0:kcmax)
2454 real(dp) :: vx_c_help, vy_c_help
2455 real(dp) :: adv_vert_help
2456 real(dp) :: dtt_dxi, dtt_deta
2457 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2458  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2459  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2460  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2461  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2462 real(dp), parameter :: zero=0.0_dp
2463 
2464 !-------- Check for boundary points --------
2465 
2466 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
2467  stop ' >>> calc_temp_enth_ssa: Boundary points not allowed.'
2468 
2469 !-------- Abbreviations --------
2470 
2471 ctr1 = atr1
2472 clb1 = alb1*q_geo(j,i)
2473 
2474 #if (ADV_VERT==1)
2475 
2476 do kc=1, kcmax-1
2477  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
2478 end do
2479 
2480 kc=0
2481 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2482  ! only needed for kc=0 ...
2483 kc=kcmax-1
2484 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2485  ! ... and kc=KCMAX-1
2486 
2487 #elif (ADV_VERT==2 || ADV_VERT==3)
2488 
2489 do kc=0, kcmax-1
2490  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2491 end do
2492 
2493 #endif
2494 
2495 do kc=0, kcmax
2496 
2497  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
2498  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
2499  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
2500  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
2501  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
2502  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
2503  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
2504  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
2505  ce5(kc) = at5(kc)/h_c(j,i)
2506  ce7(kc) = 2.0_dp*at7 &
2507  *viscosity(de_ssa(j,i), &
2508  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2509  enh_c(kc,j,i), 0_i2b) &
2510  *de_ssa(j,i)**2
2511  ci1(kc) = ai1(kc)/h_c(j,i)
2512 
2513 end do
2514 
2515 #if (ADV_VERT==1)
2516 
2517 kc=0
2518 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2519 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2520 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2521 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2522 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2523 kc=kcmax-1
2524 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2525 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2526 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2527 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2528 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2529 
2530 #elif (ADV_VERT==2 || ADV_VERT==3)
2531 
2532 do kc=0, kcmax-1
2533  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2534  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2535  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2536  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2537  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2538 end do
2539 
2540 #endif
2541 
2542 do kc=0, kcmax-1
2543  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2544  ce6(kc) = at6(kc) &
2545  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
2546  ci2(kc) = ai2(kc)/h_c(j,i)
2547 end do
2548 
2549 #if (ADV_HOR==3)
2550 dtt_dxi = 2.0_dp*dtt_2dxi
2551 dtt_deta = 2.0_dp*dtt_2deta
2552 #endif
2553 
2554 !-------- Set up the equations for the bedrock temperature --------
2555 
2556 kr=0
2557 lgs_a1(kr) = 1.0_dp
2558 lgs_a2(kr) = -1.0_dp
2559 lgs_b(kr) = clb1
2560 
2561 #if (Q_LITHO==1)
2562 ! (coupled heat-conducting bedrock)
2563 
2564 do kr=1, krmax-1
2565  lgs_a0(kr) = - ctr1
2566  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2567  lgs_a2(kr) = - ctr1
2568  lgs_b(kr) = temp_r(kr,j,i)
2569 end do
2570 
2571 #elif (Q_LITHO==0)
2572 ! (no coupled heat-conducting bedrock)
2573 
2574 do kr=1, krmax-1
2575  lgs_a0(kr) = 1.0_dp
2576  lgs_a1(kr) = 0.0_dp
2577  lgs_a2(kr) = -1.0_dp
2578  lgs_b(kr) = 2.0_dp*clb1
2579 end do
2580 
2581 #endif
2582 
2583 kr=krmax
2584 lgs_a0(kr) = 0.0_dp
2585 lgs_a1(kr) = 1.0_dp
2586 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
2587 
2588 !-------- Solve system of linear equations --------
2589 
2590 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2591 
2592 !-------- Assign the result --------
2593 
2594 do kr=0, krmax
2595  temp_r_neu(kr,j,i) = lgs_x(kr)
2596 end do
2597 
2598 !-------- Set up the equations for the ice temperature --------
2599 
2600 kc=0
2601 lgs_a1(kc) = 1.0_dp
2602 lgs_a2(kc) = 0.0_dp
2603 lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i)-delta_tm_sw, 0.0_dp)
2604  ! zero water content assumed
2605 
2606 do kc=1, kcmax-1
2607 
2608 #if (ADV_VERT==1)
2609 
2610  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2611  -ce5(kc)*ce6(kc-1)
2612  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
2613  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2614  -ce5(kc)*ce6(kc)
2615 
2616 #elif (ADV_VERT==2)
2617 
2618  lgs_a0(kc) &
2619  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2620  -ce5(kc)*ce6(kc-1)
2621  lgs_a1(kc) &
2622  = 1.0_dp &
2623  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2624  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2625  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2626  lgs_a2(kc) &
2627  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2628  -ce5(kc)*ce6(kc)
2629 
2630 #elif (ADV_VERT==3)
2631 
2632  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2633 
2634  lgs_a0(kc) &
2635  = -max(adv_vert_help, 0.0_dp) &
2636  -ce5(kc)*ce6(kc-1)
2637  lgs_a1(kc) &
2638  = 1.0_dp &
2639  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2640  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2641  lgs_a2(kc) &
2642  = min(adv_vert_help, 0.0_dp) &
2643  -ce5(kc)*ce6(kc)
2644 
2645 #endif
2646 
2647 #if (ADV_HOR==2)
2648 
2649  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2650  -dtt_2dxi* &
2651  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2652  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2653  *insq_g11_sgx(j,i) &
2654  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2655  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2656  *insq_g11_sgx(j,i-1) ) &
2657  -dtt_2deta* &
2658  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2659  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2660  *insq_g22_sgy(j,i) &
2661  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2662  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2663  *insq_g22_sgy(j-1,i) )
2664 
2665 #elif (ADV_HOR==3)
2666 
2667  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2668  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2669 
2670  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2671  -dtt_dxi* &
2672  ( min(vx_c_help, 0.0_dp) &
2673  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2674  *insq_g11_sgx(j,i) &
2675  +max(vx_c_help, 0.0_dp) &
2676  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2677  *insq_g11_sgx(j,i-1) ) &
2678  -dtt_deta* &
2679  ( min(vy_c_help, 0.0_dp) &
2680  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2681  *insq_g22_sgy(j,i) &
2682  +max(vy_c_help, 0.0_dp) &
2683  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2684  *insq_g22_sgy(j-1,i) )
2685 
2686 #endif
2687 
2688 end do
2689 
2690 kc=kcmax
2691 lgs_a0(kc) = 0.0_dp
2692 lgs_a1(kc) = 1.0_dp
2693 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2694  ! zero water content assumed
2695 
2696 !-------- Solve system of linear equations --------
2697 
2698 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2699 
2700 !-------- Assign the result --------
2701 
2702 do kc=0, kcmax
2703  enth_c_neu(kc,j,i) = lgs_x(kc)
2704  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
2705  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
2706 end do
2707 
2708 !-------- Set enthalpy and water content in the redundant,
2709 ! lower (kt) ice layer to the value at the ice base --------
2710 
2711 do kt=0, ktmax
2712  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
2713  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
2714 end do
2715 
2716 !-------- Water drainage from the non-existing temperate ice --------
2717 
2718 q_tld(j,i) = 0.0_dp
2719 
2720 !-------- Set up the equations for the age of cold ice --------
2721 
2722 kc=0 ! adv_vert_sg(0) <= 0
2723 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2724 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2725 
2726 #if (ADV_HOR==2)
2727 
2728 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2729  -dtt_2dxi* &
2730  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2731  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2732  *insq_g11_sgx(j,i) &
2733  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2734  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2735  *insq_g11_sgx(j,i-1) ) &
2736  -dtt_2deta* &
2737  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2738  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2739  *insq_g22_sgy(j,i) &
2740  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2741  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2742  *insq_g22_sgy(j-1,i) )
2743 
2744 #elif (ADV_HOR==3)
2745 
2746 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2747 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2748 
2749 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2750  -dtt_dxi* &
2751  ( min(vx_c_help, 0.0_dp) &
2752  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2753  *insq_g11_sgx(j,i) &
2754  +max(vx_c_help, 0.0_dp) &
2755  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2756  *insq_g11_sgx(j,i-1) ) &
2757  -dtt_deta* &
2758  ( min(vy_c_help, 0.0_dp) &
2759  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2760  *insq_g22_sgy(j,i) &
2761  +max(vy_c_help, 0.0_dp) &
2762  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2763  *insq_g22_sgy(j-1,i) )
2764 
2765 #endif
2766 
2767 do kc=1, kcmax-1
2768 
2769 #if (ADV_VERT==1)
2770 
2771  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2772  -ci1(kc)*ci2(kc-1)
2773  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2774  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2775  -ci1(kc)*ci2(kc)
2776 
2777 #elif (ADV_VERT==2)
2778 
2779  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2780  lgs_a1(kc) = 1.0_dp &
2781  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2782  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2783  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2784 
2785 #elif (ADV_VERT==3)
2786 
2787  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2788 
2789  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2790  lgs_a1(kc) = 1.0_dp &
2791  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2792  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2793 
2794 #endif
2795 
2796 #if (ADV_HOR==2)
2797 
2798  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2799  -dtt_2dxi* &
2800  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2801  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2802  *insq_g11_sgx(j,i) &
2803  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2804  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2805  *insq_g11_sgx(j,i-1) ) &
2806  -dtt_2deta* &
2807  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2808  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2809  *insq_g22_sgy(j,i) &
2810  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2811  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2812  *insq_g22_sgy(j-1,i) )
2813 
2814 #elif (ADV_HOR==3)
2815 
2816  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2817  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2818 
2819  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2820  -dtt_dxi* &
2821  ( min(vx_c_help, 0.0_dp) &
2822  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2823  *insq_g11_sgx(j,i) &
2824  +max(vx_c_help, 0.0_dp) &
2825  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2826  *insq_g11_sgx(j,i-1) ) &
2827  -dtt_deta* &
2828  ( min(vy_c_help, 0.0_dp) &
2829  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2830  *insq_g22_sgy(j,i) &
2831  +max(vy_c_help, 0.0_dp) &
2832  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2833  *insq_g22_sgy(j-1,i) )
2834 
2835 #endif
2836 
2837 end do
2838 
2839 kc=kcmax
2840 if (as_perp(j,i) >= zero) then
2841  lgs_a0(kc) = 0.0_dp
2842  lgs_a1(kc) = 1.0_dp
2843  lgs_b(kc) = 0.0_dp
2844 else
2845  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2846  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2847  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2848  ! assumed/enforced
2849 #if (ADV_HOR==2)
2850 
2851  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2852  -dtt_2dxi* &
2853  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2854  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2855  *insq_g11_sgx(j,i) &
2856  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2857  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2858  *insq_g11_sgx(j,i-1) ) &
2859  -dtt_2deta* &
2860  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2861  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2862  *insq_g22_sgy(j,i) &
2863  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2864  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2865  *insq_g22_sgy(j-1,i) )
2866 
2867 #elif (ADV_HOR==3)
2868 
2869  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2870  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2871 
2872  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2873  -dtt_dxi* &
2874  ( min(vx_c_help, 0.0_dp) &
2875  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2876  *insq_g11_sgx(j,i) &
2877  +max(vx_c_help, 0.0_dp) &
2878  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2879  *insq_g11_sgx(j,i-1) ) &
2880  -dtt_deta* &
2881  ( min(vy_c_help, 0.0_dp) &
2882  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2883  *insq_g22_sgy(j,i) &
2884  +max(vy_c_help, 0.0_dp) &
2885  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2886  *insq_g22_sgy(j-1,i) )
2887 
2888 #endif
2889 
2890 end if
2891 
2892 !-------- Solve system of linear equations --------
2893 
2894 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2895 
2896 !-------- Assign the result,
2897 ! restriction to interval [0, AGE_MAX yr] --------
2898 
2899 do kc=0, kcmax
2900 
2901  age_c_neu(kc,j,i) = lgs_x(kc)
2902 
2903  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
2904  age_c_neu(kc,j,i) = 0.0_dp
2905  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
2906  age_c_neu(kc,j,i) = age_max*year_sec
2907 
2908 end do
2909 
2910 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
2911 
2912 do kt=0, ktmax
2913  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
2914 end do
2915 
2916 end subroutine calc_temp_enth_ssa
2917 
2918 !-------------------------------------------------------------------------------
2919 
2920 end module calc_temp_enth_m
2921 !
calc_temp_enth_m::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_m.F90:833
sico_variables_m::vx_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
Definition: sico_variables_m.F90:402
sico_variables_m::temp_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
Definition: sico_variables_m.F90:410
sico_variables_m::kappa_r
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
Definition: sico_variables_m.F90:546
sico_variables_m::delta_tm_sw
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
Definition: sico_variables_m.F90:538
sico_variables_m::temp_r_neu
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:474
sico_variables_m::de_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) de_c
de_c(kc,j,i): Full effective strain rate in the upper (kc) ice domain
Definition: sico_variables_m.F90:501
sico_variables_m::age_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
Definition: sico_variables_m.F90:416
sico_variables_m::n_cts
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
Definition: sico_variables_m.F90:65
sico_variables_m::insq_g22_sgy
real(dp), dimension(0:jmax, 0:imax) insq_g22_sgy
insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
Definition: sico_variables_m.F90:165
sico_variables_m::enth_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:479
sico_variables_m::eaz_c
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
Definition: sico_variables_m.F90:134
sico_variables_m::q_geo
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
Definition: sico_variables_m.F90:299
sico_variables_m::vx_b_g
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:274
sico_variables_m::enh_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
Definition: sico_variables_m.F90:428
sico_variables_m::temp_r
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
Definition: sico_variables_m.F90:472
sico_variables_m::age_t_neu
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:459
sico_variables_m::h_r
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
Definition: sico_variables_m.F90:542
ice_material_properties_m
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
Definition: ice_material_properties_m.F90:39
sico_variables_m::eps
real(dp), parameter eps
eps: Small number
Definition: sico_variables_m.F90:610
sico_variables_m::omega_max
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
Definition: sico_variables_m.F90:540
sico_variables_m::insq_g11_sgx
real(dp), dimension(0:jmax, 0:imax) insq_g11_sgx
insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
Definition: sico_variables_m.F90:163
enth_temp_omega_m
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
Definition: enth_temp_omega_m.F90:37
sico_maths_m
Several mathematical tools used by SICOPOLIS.
Definition: sico_maths_m.F90:35
sico_variables_m::insq_g22_g
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
Definition: sico_variables_m.F90:153
sico_variables_m::dh_c_dtau
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
Definition: sico_variables_m.F90:242
ice_material_properties_m::kappa_val
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
Definition: ice_material_properties_m.F90:247
sico_variables_m::maske
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
Definition: sico_variables_m.F90:49
sico_variables_m::dh_c_dxi_g
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:220
enth_temp_omega_m::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_m.F90:269
sico_variables_m::as_perp
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
Definition: sico_variables_m.F90:320
sico_variables_m::dzs_deta_g
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:224
sico_variables_m::temp_t_m
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
Definition: sico_variables_m.F90:455
sico_variables_m::de_ssa
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
Definition: sico_variables_m.F90:431
sico_variables_m::h_t_neu
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:351
sico_maths_m::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_maths_m.F90:134
ice_material_properties_m::creep
real(dp) function, public creep(sigma_val)
Creep response function for ice.
Definition: ice_material_properties_m.F90:330
sico_variables_m::temp_s
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
Definition: sico_variables_m.F90:334
sico_variables_m::g
real(dp) g
G: Acceleration due to gravity.
Definition: sico_variables_m.F90:531
sico_variables_m::kc_cts
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
Definition: sico_variables_m.F90:69
sico_variables_m::q_tld
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Definition: sico_variables_m.F90:307
sico_vars_m
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
enth_temp_omega_m::omega_fct_enth
real(dp) function, public omega_fct_enth(enth_val, temp_m_val)
Water content as a function of enthalpy.
Definition: enth_temp_omega_m.F90:293
sico_variables_m::ea
real(dp) ea
ea: Abbreviation for exp(aa)
Definition: sico_variables_m.F90:132
calc_temp_enth_m::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_m.F90:1231
sico_variables_m::enth_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
Definition: sico_variables_m.F90:477
sico_variables_m::flag_aa_nonzero
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Definition: sico_variables_m.F90:130
sico_variables_m::age_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:418
sico_variables_m::aa
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Definition: sico_variables_m.F90:126
sico_variables_m::rho_c_r
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
Definition: sico_variables_m.F90:544
sico_types_m
Declarations of kind types for SICOPOLIS.
Definition: sico_types_m.F90:35
calc_temp_enth_m::calc_temp_enth_2_a1
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: ...
Definition: calc_temp_enth_m.F90:1680
sico_variables_m::dh_c_deta_g
real(dp), dimension(0:jmax, 0:imax) dh_c_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:230
calc_temp_enth_m::calc_temp_enth_2_b
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: ...
Definition: calc_temp_enth_m.F90:1876
sico_variables_m::kc_cts_neu
integer(i2b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:71
sico_variables_m::temp_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:412
sico_variables_m::h_c
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
Definition: sico_variables_m.F90:188
sico_variables_m::vy_b_g
real(dp), dimension(0:jmax, 0:imax) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:276
sico_variables_m::omega_t_neu
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:453
sico_variables_m::dzm_dtau
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
Definition: sico_variables_m.F90:236
sico_variables_m::omega_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
Definition: sico_variables_m.F90:481
sico_variables_m::temp_c_m
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_m
temp_c_m(kc,j,i): Melting temperature in the upper (kc) ice domain
Definition: sico_variables_m.F90:414
sico_variables_m::enth_t_neu
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:488
sico_variables_m::nue
real(dp) nue
NUE: Water diffusivity in ice.
Definition: sico_variables_m.F90:533
ice_material_properties_m::c_val
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
Definition: ice_material_properties_m.F90:289
sico_variables_m::p_weert_inv
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
Definition: sico_variables_m.F90:250
sico_variables_m::flag_shelfy_stream
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
Definition: sico_variables_m.F90:113
sico_variables_m::vz_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
Definition: sico_variables_m.F90:408
calc_temp_enth_m::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_m.F90:1019
sico_variables_m::h_c_neu
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:349
enth_temp_omega_m::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_m.F90:254
calc_temp_enth_m
Computation of temperature, water content and age with the enthalpy method.
Definition: calc_temp_enth_m.F90:35
calc_temp_enth_m::calc_temp_enth_ssa
subroutine calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, alb1, ai1, ai2, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature and age for ice shelves (floating ice) with the enthalpy method...
Definition: calc_temp_enth_m.F90:2430
sico_variables_m::dzm_dxi_g
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:216
sico_variables_m::dzs_dxi_g
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:214
sico_variables_m::vy_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
Definition: sico_variables_m.F90:446
sico_variables_m::beta
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
Definition: sico_variables_m.F90:535
calc_temp_enth_m::calc_temp_enth_2_d
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: ...
Definition: calc_temp_enth_m.F90:2126
sico_variables_m::dzm_deta_g
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:226
calc_temp_enth_m::calc_temp_enth
subroutine, public calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Main subroutine of calc_temp_enth_m: Computation of temperature, water content and age with the entha...
Definition: calc_temp_enth_m.F90:54
sico_variables_m::insq_g11_g
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
Definition: sico_variables_m.F90:151
sico_variables_m::c_drag
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
Definition: sico_variables_m.F90:256
sico_variables_m::vy_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
Definition: sico_variables_m.F90:405
ice_material_properties_m::viscosity
real(dp) function, public 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: ice_material_properties_m.F90:406
sico_variables_m::zm_neu
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:343
ice_material_properties_m::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_properties_m.F90:217
sico_variables_m::rho
real(dp) rho
RHO: Density of ice.
Definition: sico_variables_m.F90:523
sico_variables_m::n_cts_neu
integer(i2b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:67
sico_variables_m::sigma_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
Definition: sico_variables_m.F90:426
sico_variables_m::zb
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
Definition: sico_variables_m.F90:174
sico_variables_m::omega_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:483
sico_variables_m::zeta_c
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
Definition: sico_variables_m.F90:119
sico_variables_m::h_t
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Definition: sico_variables_m.F90:192
sico_variables_m
Declarations of global variables for SICOPOLIS.
Definition: sico_variables_m.F90:35
sico_variables_m::vx_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
Definition: sico_variables_m.F90:443
Generated on Thu Aug 3 2017 21:26:01 for SICOPOLIS V3.3 by   doxygen 1.8.11