SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_temp_poly.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ p o l y
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age in polythermal mode.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of temperature, water content and age in polythermal mode.
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
36  dtime_temp, mean_accum_inv)
37 
38 use sico_types
39 use sico_variables
40 
41 implicit none
42 
43 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
44 real(dp), intent(in) :: dtime_temp
45 real(dp), intent(in) :: mean_accum_inv
46 
47 integer(i4b) :: i, j, kc, kt, kr, ii, jj
48 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
49  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
50  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
51  acb1, acb2, acb3, acb4, &
52  ai1(0:kcmax), ai2(0:kcmax), ai3, ai4, &
53  atr1, am1, am2, alb1
54 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
55 real(dp) :: dtime_temp_inv, dtt_2dxi, dtt_2deta
56 real(dp) :: temp_c_help(0:kcmax)
57 real(dp) :: fact_thick
58 real(dp) :: time_lag_cts
59 real(dp) :: vol_t, vol_t_smooth, korrfakt_t, &
60  dh_t_smooth(0:jmax,0:imax)
61 
62 !-------- Term abbreviations
63 
64 at7 = 2.0_dp/rho*dtime_temp
65 
66 aw1 = dtime_temp/dzeta_t
67 aw2 = dtime_temp/dzeta_t
68 aw3 = dtime_temp/dzeta_t
69 aw4 = dtime_temp/dzeta_t
70 aw5 = nue/rho*dtime_temp/(dzeta_t**2)
71 aw7 = 2.0_dp/(rho*l)*dtime_temp
72 aw8 = beta**2/(rho*l) &
73  *(kappa_val(0.0_dp)-kappa_val(-1.0_dp))*dtime_temp
74 aw9 = beta/l*dtime_temp
75 
76 ai3 = agediff*dtime_temp/(dzeta_t**2)
77 ai4 = deform*dzeta_c/(ea-1.0_dp)
78 
79 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
80 
81 am1 = deform*beta*dzeta_c/(ea-1.0_dp)
82 am2 = deform*l*rho*dzeta_c/(ea-1.0_dp)
83 
84 acb1 = (ea-1.0_dp)/deform/dzeta_c
85 acb2 = kappa_r/h_r/dzeta_r
86 acb3 = rho*g
87 acb4 = rho*g
88 
89 alb1 = h_r/kappa_r*dzeta_r
90 
91 aqtld = dzeta_t/dtime_temp
92 
93 dtt_2dxi = 0.5_dp*dtime_temp/dxi
94 dtt_2deta = 0.5_dp*dtime_temp/deta
95 
96 dtime_temp_inv = 1.0_dp/dtime_temp
97 
98 do kc=0, kcmax
99  at1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
100  at2_1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
101  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(deform*eaz_c(kc)) &
102  *dtime_temp/dzeta_c
103  at3_1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
104  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(deform*eaz_c(kc)) &
105  *dtime_temp/dzeta_c
106  at4_1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
107  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(deform*eaz_c(kc)) &
108  *dtime_temp/dzeta_c
109  at5(kc) = (ea-1.0_dp)/(rho*deform*eaz_c(kc)) &
110  *dtime_temp/dzeta_c
111  if (kc /= kcmax) then
112  at6(kc) = (ea-1.0_dp) &
113  /(deform*exp(deform*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
114  /dzeta_c
115  else
116  at6(kc) = 0.0_dp
117  end if
118  ai1(kc) = agediff*(ea-1.0_dp)/(deform*eaz_c(kc)) &
119  *dtime_temp/dzeta_c
120  if (kc /= kcmax) then
121  ai2(kc) = (ea-1.0_dp) &
122  /(deform*exp(deform*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
123  /dzeta_c
124  else
125  ai2(kc) = 0.0_dp
126  end if
127 end do
128 
129 !-------- Computation loop for temperature, water content and age --------
130 
131 do i=1, imax-1 ! skipping domain margins
132 do j=1, jmax-1 ! skipping domain margins
133 
134  if (maske(j,i)==0) then ! glaciated land
135 
136 ! ------ Old vertical column cold
137 
138  if (n_cts(j,i).eq.-1) then
139 
140  n_cts_neu(j,i) = n_cts(j,i)
141  zm_neu(j,i) = zb(j,i)
142  h_c_neu(j,i) = h_c(j,i)
143  h_t_neu(j,i) = h_t(j,i)
144 
145  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
146  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
147  acb3, acb4, alb1, ai1, ai2, ai4, &
148  mean_accum_inv, &
149  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
150 
151 ! ---- Check whether base has become temperate
152 
153  if (temp_c_neu(0,j,i).gt.temp_c_m(0,j,i)) then
154 
155  n_cts_neu(j,i) = 0
156 
157  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
158  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
159  ai1, ai2, ai4, mean_accum_inv, &
160  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
161 
162  end if
163 
164 ! ---- Check whether even temperate layer has formed
165 
166  if ( &
167  ( n_cts_neu(j,i).eq.0 ).and. &
168  ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
169  .gt.(am1*h_c_neu(j,i)) ) &
170  ) then
171 
172  n_cts_neu(j,i) = 1
173  zm_neu(j,i) = zb(j,i)+0.001_dp
174  h_c_neu(j,i) = h_c(j,i)-0.001_dp
175  h_t_neu(j,i) = h_t(j,i)+0.001_dp
176 ! ! CTS preliminarily positioned 1 mm above ice base --------
177 
178  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
179  at4_1, at4_2, at5, at6, at7, atr1, &
180  am1, am2, alb1, &
181  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
182  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
183  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
184 
185  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
186  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
187  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
188  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
189  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
190  i, j)
191 
192  end if
193 
194 ! ------ Old vertical column with temperate base
195 
196  else if (n_cts(j,i).eq.0) then
197 
198  n_cts_neu(j,i) = n_cts(j,i)
199  zm_neu(j,i) = zb(j,i)
200  h_c_neu(j,i) = h_c(j,i)
201  h_t_neu(j,i) = h_t(j,i)
202 
203  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
204  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
205  ai1, ai2, ai4, mean_accum_inv, &
206  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
207 
208 ! ---- Check whether temperate base becomes cold
209 
210  if ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
211  .lt. (am1*h_c(j,i)) ) then
212 
213  n_cts_neu(j,i) = -1
214 
215  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
216  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
217  acb3, acb4, alb1, ai1, ai2, ai4, &
218  mean_accum_inv, &
219  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
220 
221  if (temp_c_neu(0,j,i).ge.temp_c_m(0,j,i)) then
222 
223  n_cts_neu(j,i) = 0
224 
225  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
226  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
227  ai1, ai2, ai4, mean_accum_inv, &
228  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
229 
230  end if
231 
232  end if
233 
234 ! ---- Check whether temperate layer has formed
235 
236  if ( &
237  ( n_cts_neu(j,i).eq.0 ).and. &
238  ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
239  .gt.(am1*h_c_neu(j,i)) ) &
240  ) then
241 
242  n_cts_neu(j,i) = 1
243  zm_neu(j,i) = zb(j,i)+0.001_dp
244  h_c_neu(j,i) = h_c(j,i)-0.001_dp
245  h_t_neu(j,i) = h_t(j,i)+0.001_dp
246 ! ! CTS preliminarily positioned 1 mm above ice base --------
247 
248  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
249  at4_1, at4_2, at5, at6, at7, atr1, &
250  am1, am2, alb1, &
251  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
252  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
253  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
254 
255  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
256  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
257  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
258  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
259  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
260  i, j)
261 
262  end if
263 
264 ! ------ Old vertical column with temperate base and temperate layer
265 
266  else ! n_cts(j,i).eq.1
267 
268  n_cts_neu(j,i) = n_cts(j,i)
269  zm_neu(j,i) = zm(j,i)
270  h_c_neu(j,i) = h_c(j,i)
271  h_t_neu(j,i) = h_t(j,i)
272 
273  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
274  at4_1, at4_2, at5, at6, at7, atr1, &
275  am1, am2, alb1, &
276  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
277  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
278  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
279 
280  if ( (temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))).gt.0.0_dp ) &
281  then
282  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
283  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
284  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
285  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
286  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
287  i, j)
288  else
289  call shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
290  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
291  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
292  ai1, ai2, ai3, ai4, mean_accum_inv, dzeta_t, &
293  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
294  i, j)
295  end if
296 
297  end if
298 
299 #if MARGIN==3
300 
301  else if (maske(j,i)==3) then ! floating ice
302 
303  n_cts_neu(j,i) = -1
304  zm_neu(j,i) = zb(j,i)
305  h_c_neu(j,i) = h_c(j,i) + h_t(j,i)
306  h_t_neu(j,i) = 0.0_dp
307 
308  call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
309  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
310  ai1, ai2, ai4, mean_accum_inv, &
311  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
312 
313 ! ------ Reset temperatures above melting to the melting point
314 ! (should not occur, but just in case)
315 
316  do kc=0, kcmax
317  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
318  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
319  end do
320 
321 #endif
322 
323  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
324 
325  n_cts_neu(j,i) = -1
326  zm_neu(j,i) = zb(j,i)
327  h_c_neu(j,i) = h_c(j,i)
328  h_t_neu(j,i) = h_t(j,i)
329 
330  call calc_temp_r(atr1, alb1, i, j)
331 
332 endif
333 
334 end do
335 end do ! End of computation loop
336 
337 !-------- Extrapolate values on margins --------
338 
339 ! ------ Lower left corner
340 
341 i=0
342 j=0
343 
344 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
345  ! glaciated land or floating ice
346  ii=i+1
347  jj=j+1
348 
349  do kc=0,kcmax
350  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
351  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
352  end do
353 
354  do kt=0,ktmax
355  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
356  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
357  end do
358 
359  do kr=0,krmax
360  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
361  end do
362 
363  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
364  h_c_neu(j,i) = h_c(j,i)
365  h_t_neu(j,i) = h_t(j,i)
366  zm_neu(j,i) = zb(j,i)
367 
368 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
369 
370  n_cts_neu(j,i) = -1
371  zm_neu(j,i) = zb(j,i)
372  h_c_neu(j,i) = h_c(j,i)
373  h_t_neu(j,i) = h_t(j,i)
374 
375  call calc_temp_r(atr1, alb1, i, j)
376 
377 end if
378 
379 ! ------ Lower right corner
380 
381 i=imax
382 j=0
383 
384 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
385  ! glaciated land or floating ice
386  ii=i-1
387  jj=j+1
388 
389  do kc=0,kcmax
390  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
391  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
392  end do
393 
394  do kt=0,ktmax
395  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
396  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
397  end do
398 
399  do kr=0,krmax
400  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
401  end do
402 
403  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
404  h_c_neu(j,i) = h_c(j,i)
405  h_t_neu(j,i) = h_t(j,i)
406  zm_neu(j,i) = zb(j,i)
407 
408 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
409 
410  n_cts_neu(j,i) = -1
411  zm_neu(j,i) = zb(j,i)
412  h_c_neu(j,i) = h_c(j,i)
413  h_t_neu(j,i) = h_t(j,i)
414 
415  call calc_temp_r(atr1, alb1, i, j)
416 
417 end if
418 
419 ! ------ Upper left corner
420 
421 i=0
422 j=jmax
423 
424 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
425  ! glaciated land or floating ice
426  ii=i+1
427  jj=j-1
428 
429  do kc=0,kcmax
430  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
431  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
432  end do
433 
434  do kt=0,ktmax
435  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
436  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
437  end do
438 
439  do kr=0,krmax
440  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
441  end do
442 
443  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
444  h_c_neu(j,i) = h_c(j,i)
445  h_t_neu(j,i) = h_t(j,i)
446  zm_neu(j,i) = zb(j,i)
447 
448 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
449 
450  n_cts_neu(j,i) = -1
451  zm_neu(j,i) = zb(j,i)
452  h_c_neu(j,i) = h_c(j,i)
453  h_t_neu(j,i) = h_t(j,i)
454 
455  call calc_temp_r(atr1, alb1, i, j)
456 
457 end if
458 
459 ! ------ Upper right corner
460 
461 i=imax
462 j=jmax
463 
464 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
465  ! glaciated land or floating ice
466  ii=i-1
467  jj=j-1
468 
469  do kc=0,kcmax
470  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
471  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
472  end do
473 
474  do kt=0,ktmax
475  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
476  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
477  end do
478 
479  do kr=0,krmax
480  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
481  end do
482 
483  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
484  h_c_neu(j,i) = h_c(j,i)
485  h_t_neu(j,i) = h_t(j,i)
486  zm_neu(j,i) = zb(j,i)
487 
488 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
489 
490  n_cts_neu(j,i) = -1
491  zm_neu(j,i) = zb(j,i)
492  h_c_neu(j,i) = h_c(j,i)
493  h_t_neu(j,i) = h_t(j,i)
494 
495  call calc_temp_r(atr1, alb1, i, j)
496 
497 end if
498 
499 ! ------ Lower and upper margins
500 
501 do i=1, imax-1
502 
503 ! ---- Lower margin
504 
505  j=0
506 
507  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
508  ! glaciated land or floating ice
509  ii=i
510  jj=j+1
511 
512  do kc=0,kcmax
513  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
514  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
515  end do
516 
517  do kt=0,ktmax
518  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
519  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
520  end do
521 
522  do kr=0,krmax
523  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
524  end do
525 
526  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
527  h_c_neu(j,i) = h_c(j,i)
528  h_t_neu(j,i) = h_t(j,i)
529  zm_neu(j,i) = zb(j,i)
530 
531  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
532 
533  n_cts_neu(j,i) = -1
534  zm_neu(j,i) = zb(j,i)
535  h_c_neu(j,i) = h_c(j,i)
536  h_t_neu(j,i) = h_t(j,i)
537 
538  call calc_temp_r(atr1, alb1, i, j)
539 
540  end if
541 
542 ! ---- Upper margin
543 
544  j=jmax
545 
546  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
547  ! glaciated land or floating ice
548  ii=i
549  jj=j-1
550 
551  do kc=0,kcmax
552  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
553  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
554  end do
555 
556  do kt=0,ktmax
557  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
558  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
559  end do
560 
561  do kr=0,krmax
562  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
563  end do
564 
565  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
566  h_c_neu(j,i) = h_c(j,i)
567  h_t_neu(j,i) = h_t(j,i)
568  zm_neu(j,i) = zb(j,i)
569 
570  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
571 
572  n_cts_neu(j,i) = -1
573  zm_neu(j,i) = zb(j,i)
574  h_c_neu(j,i) = h_c(j,i)
575  h_t_neu(j,i) = h_t(j,i)
576 
577  call calc_temp_r(atr1, alb1, i, j)
578 
579  end if
580 
581 end do
582 
583 ! ------ Left and right margins
584 
585 do j=1, jmax-1
586 
587 ! ---- Left margin
588 
589  i=0
590 
591  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
592  ! glaciated land or floating ice
593  ii=i+1
594  jj=j
595 
596  do kc=0,kcmax
597  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
598  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
599  end do
600 
601  do kt=0,ktmax
602  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
603  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
604  end do
605 
606  do kr=0,krmax
607  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
608  end do
609 
610  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
611  h_c_neu(j,i) = h_c(j,i)
612  h_t_neu(j,i) = h_t(j,i)
613  zm_neu(j,i) = zb(j,i)
614 
615  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
616 
617  n_cts_neu(j,i) = -1
618  zm_neu(j,i) = zb(j,i)
619  h_c_neu(j,i) = h_c(j,i)
620  h_t_neu(j,i) = h_t(j,i)
621 
622  call calc_temp_r(atr1, alb1, i, j)
623 
624  end if
625 
626 ! ---- Right margin
627 
628  i=imax
629 
630  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
631  ! glaciated land or floating ice
632  ii=i-1
633  jj=j
634 
635  do kc=0,kcmax
636  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
637  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
638  end do
639 
640  do kt=0,ktmax
641  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
642  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
643  end do
644 
645  do kr=0,krmax
646  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
647  end do
648 
649  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
650  h_c_neu(j,i) = h_c(j,i)
651  h_t_neu(j,i) = h_t(j,i)
652  zm_neu(j,i) = zb(j,i)
653 
654  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
655 
656  n_cts_neu(j,i) = -1
657  zm_neu(j,i) = zb(j,i)
658  h_c_neu(j,i) = h_c(j,i)
659  h_t_neu(j,i) = h_t(j,i)
660 
661  call calc_temp_r(atr1, alb1, i, j)
662 
663  end if
664 
665 end do
666 
667 !-------- Smoothing of H_t_neu with numerical diffusion --------
668 
669 ! ------ Volume of temperate ice without smoothing
670 
671 vol_t = 0.0_dp
672 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
673 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
674  if (n_cts_neu(j,i).eq.1) then
675  vol_t = vol_t + h_t_neu(j,i)*area(j,i)
676  end if
677 end do
678 end do
679 
680 ! ------ Smoothing
681 
682 do i=1, imax-1
683 do j=1, jmax-1
684  if (n_cts_neu(j,i).ne.-1) then
685 
686  dh_t_smooth(j,i) = numdiff_h_t* ( -4.0_dp*h_t_neu(j,i) &
687  +h_t_neu(j,i+1)+h_t_neu(j,i-1) &
688  +h_t_neu(j+1,i)+h_t_neu(j-1,i) )
689  if (dh_t_smooth(j,i).gt.0.001_dp) n_cts_neu(j,i) = 1
690 
691  end if
692 end do
693 end do
694 
695 do i=1, imax-1
696 do j=1, jmax-1
697  if (n_cts_neu(j,i).eq.1) then
698  h_t_neu(j,i) = h_t_neu(j,i) + dh_t_smooth(j,i)
699  end if
700 end do
701 end do
702 
703 ! ------ Volume of temperate ice with smoothing
704 
705 vol_t_smooth = 0.0_dp
706 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
707 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
708  if (n_cts_neu(j,i).eq.1) then
709  vol_t_smooth = vol_t_smooth + h_t_neu(j,i)*area(j,i)
710  end if
711 end do
712 end do
713 
714 ! ------ Correction so that volume is not changed by the smoothing
715 
716 if (vol_t_smooth.gt.0.0_dp) then
717 
718  korrfakt_t = vol_t/vol_t_smooth
719  do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
720  do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
721  if (n_cts_neu(j,i).eq.1) then
722  h_t_neu(j,i) = h_t_neu(j,i)*korrfakt_t
723 ! zm_neu(j,i) = zb(j,i) + H_t_neu(j,i)
724 ! H_c_neu(j,i) = zs(j,i) - zm_neu(j,i)
725  end if
726  end do
727  end do
728 
729 end if
730 
731 !-------- Numerical time lag for evolution of H_t_neu --------
732 
733 time_lag_cts = tau_cts*year_sec ! yr --> s
734 
735 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
736 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
737 
738  if (n_cts_neu(j,i).eq.1) then
739 
740  h_t_neu(j,i) = ( time_lag_cts*h_t(j,i) &
741  + dtime_temp*h_t_neu(j,i) ) &
742  /(time_lag_cts+dtime_temp)
743 
744  zm_neu(j,i) = zb(j,i) + h_t_neu(j,i)
745  h_c_neu(j,i) = zs(j,i) - zm_neu(j,i)
746 
747  end if
748 
749 end do
750 end do
751 
752 end subroutine calc_temp_poly
753 !
Generated on Thu Sep 12 2013 09:45:11 for SICOPOLIS V3.1 by   doxygen 1.8.2