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