SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_temp_cold.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ c o l d
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age in cold-ice 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 and age in cold-ice mode.
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp_cold(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, 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 
59 !-------- Term abbreviations
60 
61 at7 = 2.0_dp/rho*dtime_temp
62 
63 aw1 = dtime_temp/dzeta_t
64 aw2 = dtime_temp/dzeta_t
65 aw3 = dtime_temp/dzeta_t
66 aw4 = dtime_temp/dzeta_t
67 aw5 = nue/rho*dtime_temp/(dzeta_t**2)
68 aw7 = 2.0_dp/(rho*l)*dtime_temp
69 aw8 = beta**2/(rho*l) &
70  *(kappa_val(0.0_dp)-kappa_val(-1.0_dp))*dtime_temp
71 aw9 = beta/l*dtime_temp
72 
73 ai3 = agediff*dtime_temp/(dzeta_t**2)
74 ai4 = deform*dzeta_c/(ea-1.0_dp)
75 
76 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
77 
78 am1 = deform*beta*dzeta_c/(ea-1.0_dp)
79 am2 = deform*l*rho*dzeta_c/(ea-1.0_dp)
80 
81 acb1 = (ea-1.0_dp)/deform/dzeta_c
82 acb2 = kappa_r/h_r/dzeta_r
83 acb3 = rho*g
84 acb4 = rho*g
85 
86 alb1 = h_r/kappa_r*dzeta_r
87 
88 aqtld = dzeta_t/dtime_temp
89 
90 dtt_2dxi = 0.5_dp*dtime_temp/dxi
91 dtt_2deta = 0.5_dp*dtime_temp/deta
92 
93 dtime_temp_inv = 1.0_dp/dtime_temp
94 
95 do kc=0, kcmax
96  at1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
97  at2_1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
98  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(deform*eaz_c(kc)) &
99  *dtime_temp/dzeta_c
100  at3_1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
101  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(deform*eaz_c(kc)) &
102  *dtime_temp/dzeta_c
103  at4_1(kc) = (ea-1.0_dp)/(deform*eaz_c(kc))*dtime_temp/dzeta_c
104  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(deform*eaz_c(kc)) &
105  *dtime_temp/dzeta_c
106  at5(kc) = (ea-1.0_dp)/(rho*deform*eaz_c(kc)) &
107  *dtime_temp/dzeta_c
108  if (kc /= kcmax) then
109  at6(kc) = (ea-1.0_dp) &
110  /(deform*exp(deform*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
111  /dzeta_c
112  else
113  at6(kc) = 0.0_dp
114  end if
115  ai1(kc) = agediff*(ea-1.0_dp)/(deform*eaz_c(kc)) &
116  *dtime_temp/dzeta_c
117  if (kc /= kcmax) then
118  ai2(kc) = (ea-1.0_dp) &
119  /(deform*exp(deform*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
120  /dzeta_c
121  else
122  ai2(kc) = 0.0_dp
123  end if
124 end do
125 
126 !-------- Computation loop for temperature and age --------
127 
128 do i=1, imax-1 ! skipping domain margins
129 do j=1, jmax-1 ! skipping domain margins
130 
131  if (maske(j,i)==0) then ! glaciated land
132 
133  n_cts_neu(j,i) = -1
134  zm_neu(j,i) = zb(j,i)
135  h_c_neu(j,i) = h_c(j,i)
136  h_t_neu(j,i) = h_t(j,i)
137 
138  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
139  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
140  acb3, acb4, alb1, ai1, ai2, ai4, &
141  mean_accum_inv, &
142  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
143 
144 ! ------ Reset temperatures above melting to the melting point,
145 ! look for the CTS
146 
147  kc_cts(j,i) = 0
148 
149  if (temp_c_neu(0,j,i).gt.temp_c_m(0,j,i)) then
150  n_cts_neu(j,i) = 0
151  kc_cts(j,i) = 0
152  temp_c_neu(0,j,i) = temp_c_m(0,j,i)
153  temp_r_neu(krmax,j,i) = temp_c_m(0,j,i)
154  end if
155 
156  do kc=1, kcmax
157  if (temp_c_neu(kc,j,i).gt.temp_c_m(kc,j,i)) then
158  kc_cts(j,i) = kc
159  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
160  end if
161  end do
162 
163 #if MARGIN==3
164 
165  else if (maske(j,i)==3) then ! floating ice
166 
167  n_cts_neu(j,i) = -1
168  kc_cts(j,i) = 0
169  zm_neu(j,i) = zb(j,i)
170  h_c_neu(j,i) = h_c(j,i)
171  h_t_neu(j,i) = 0.0_dp
172 
173  call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
174  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
175  ai1, ai2, ai4, mean_accum_inv, &
176  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
177 
178 ! ------ Reset temperatures above melting to the melting point
179 ! (should not occur, but just in case)
180 
181  do kc=0, kcmax
182  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
183  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
184  end do
185 
186 #endif
187 
188  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
189 
190  n_cts_neu(j,i) = -1
191  kc_cts(j,i) = 0
192  zm_neu(j,i) = zb(j,i)
193  h_c_neu(j,i) = h_c(j,i)
194  h_t_neu(j,i) = h_t(j,i)
195 
196  call calc_temp_r(atr1, alb1, i, j)
197 
198  end if
199 
200 end do
201 end do
202 
203 !-------- Extrapolate values on margins --------
204 
205 ! ------ Lower left corner
206 
207 i=0
208 j=0
209 
210 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
211  ! glaciated land or floating ice
212  ii=i+1
213  jj=j+1
214 
215  do kc=0,kcmax
216  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
217  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
218  end do
219 
220  do kr=0,krmax
221  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
222  end do
223 
224  n_cts_neu(j,i) = n_cts_neu(jj,ii)
225  kc_cts(j,i) = kc_cts(jj,ii)
226  zm_neu(j,i) = zb(j,i)
227  h_c_neu(j,i) = h_c(j,i)
228  h_t_neu(j,i) = h_t(j,i)
229 
230 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
231 
232  n_cts_neu(j,i) = -1
233  kc_cts(j,i) = 0
234  zm_neu(j,i) = zb(j,i)
235  h_c_neu(j,i) = h_c(j,i)
236  h_t_neu(j,i) = h_t(j,i)
237 
238  call calc_temp_r(atr1, alb1, i, j)
239 
240 end if
241 
242 ! ------ Lower right corner
243 
244 i=imax
245 j=0
246 
247 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
248  ! glaciated land or floating ice
249  ii=i-1
250  jj=j+1
251 
252  do kc=0,kcmax
253  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
254  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
255  end do
256 
257  do kr=0,krmax
258  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
259  end do
260 
261  n_cts_neu(j,i) = n_cts_neu(jj,ii)
262  kc_cts(j,i) = kc_cts(jj,ii)
263  zm_neu(j,i) = zb(j,i)
264  h_c_neu(j,i) = h_c(j,i)
265  h_t_neu(j,i) = h_t(j,i)
266 
267 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
268 
269  n_cts_neu(j,i) = -1
270  kc_cts(j,i) = 0
271  zm_neu(j,i) = zb(j,i)
272  h_c_neu(j,i) = h_c(j,i)
273  h_t_neu(j,i) = h_t(j,i)
274 
275  call calc_temp_r(atr1, alb1, i, j)
276 
277 end if
278 
279 ! ------ Upper left corner
280 
281 i=0
282 j=jmax
283 
284 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
285  ! glaciated land or floating ice
286  ii=i+1
287  jj=j-1
288 
289  do kc=0,kcmax
290  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
291  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
292  end do
293 
294  do kr=0,krmax
295  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
296  end do
297 
298  n_cts_neu(j,i) = n_cts_neu(jj,ii)
299  kc_cts(j,i) = kc_cts(jj,ii)
300  zm_neu(j,i) = zb(j,i)
301  h_c_neu(j,i) = h_c(j,i)
302  h_t_neu(j,i) = h_t(j,i)
303 
304 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
305 
306  n_cts_neu(j,i) = -1
307  kc_cts(j,i) = 0
308  zm_neu(j,i) = zb(j,i)
309  h_c_neu(j,i) = h_c(j,i)
310  h_t_neu(j,i) = h_t(j,i)
311 
312  call calc_temp_r(atr1, alb1, i, j)
313 
314 end if
315 
316 ! ------ Upper right corner
317 
318 i=imax
319 j=jmax
320 
321 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
322  ! glaciated land or floating ice
323  ii=i-1
324  jj=j-1
325 
326  do kc=0,kcmax
327  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
328  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
329  end do
330 
331  do kr=0,krmax
332  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
333  end do
334 
335  n_cts_neu(j,i) = n_cts_neu(jj,ii)
336  kc_cts(j,i) = kc_cts(jj,ii)
337  zm_neu(j,i) = zb(j,i)
338  h_c_neu(j,i) = h_c(j,i)
339  h_t_neu(j,i) = h_t(j,i)
340 
341 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
342 
343  n_cts_neu(j,i) = -1
344  kc_cts(j,i) = 0
345  zm_neu(j,i) = zb(j,i)
346  h_c_neu(j,i) = h_c(j,i)
347  h_t_neu(j,i) = h_t(j,i)
348 
349  call calc_temp_r(atr1, alb1, i, j)
350 
351 end if
352 
353 ! ------ Lower and upper margins
354 
355 do i=1, imax-1
356 
357 ! ---- Lower margin
358 
359  j=0
360 
361  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
362  ! glaciated land or floating ice
363  ii=i
364  jj=j+1
365 
366  do kc=0,kcmax
367  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
368  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
369  end do
370 
371  do kr=0,krmax
372  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
373  end do
374 
375  n_cts_neu(j,i) = n_cts_neu(jj,ii)
376  kc_cts(j,i) = kc_cts(jj,ii)
377  zm_neu(j,i) = zb(j,i)
378  h_c_neu(j,i) = h_c(j,i)
379  h_t_neu(j,i) = h_t(j,i)
380 
381  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
382 
383  n_cts_neu(j,i) = -1
384  kc_cts(j,i) = 0
385  zm_neu(j,i) = zb(j,i)
386  h_c_neu(j,i) = h_c(j,i)
387  h_t_neu(j,i) = h_t(j,i)
388 
389  call calc_temp_r(atr1, alb1, i, j)
390 
391  end if
392 
393 ! ---- Upper margin
394 
395  j=jmax
396 
397  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
398  ! glaciated land or floating ice
399  ii=i
400  jj=j-1
401 
402  do kc=0,kcmax
403  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
404  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
405  end do
406 
407  do kr=0,krmax
408  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
409  end do
410 
411  n_cts_neu(j,i) = n_cts_neu(jj,ii)
412  kc_cts(j,i) = kc_cts(jj,ii)
413  zm_neu(j,i) = zb(j,i)
414  h_c_neu(j,i) = h_c(j,i)
415  h_t_neu(j,i) = h_t(j,i)
416 
417  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
418 
419  n_cts_neu(j,i) = -1
420  kc_cts(j,i) = 0
421  zm_neu(j,i) = zb(j,i)
422  h_c_neu(j,i) = h_c(j,i)
423  h_t_neu(j,i) = h_t(j,i)
424 
425  call calc_temp_r(atr1, alb1, i, j)
426 
427  end if
428 
429 end do
430 
431 ! ------ Left and right margins
432 
433 do j=1, jmax-1
434 
435 ! ---- Left margin
436 
437  i=0
438 
439  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
440  ! glaciated land or floating ice
441  ii=i+1
442  jj=j
443 
444  do kc=0,kcmax
445  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
446  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
447  end do
448 
449  do kr=0,krmax
450  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
451  end do
452 
453  n_cts_neu(j,i) = n_cts_neu(jj,ii)
454  kc_cts(j,i) = kc_cts(jj,ii)
455  zm_neu(j,i) = zb(j,i)
456  h_c_neu(j,i) = h_c(j,i)
457  h_t_neu(j,i) = h_t(j,i)
458 
459  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
460 
461  n_cts_neu(j,i) = -1
462  kc_cts(j,i) = 0
463  zm_neu(j,i) = zb(j,i)
464  h_c_neu(j,i) = h_c(j,i)
465  h_t_neu(j,i) = h_t(j,i)
466 
467  call calc_temp_r(atr1, alb1, i, j)
468 
469  end if
470 
471 ! ---- Right margin
472 
473  i=imax
474 
475  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
476  ! glaciated land or floating ice
477  ii=i-1
478  jj=j
479 
480  do kc=0,kcmax
481  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
482  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
483  end do
484 
485  do kr=0,krmax
486  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
487  end do
488 
489  n_cts_neu(j,i) = n_cts_neu(jj,ii)
490  kc_cts(j,i) = kc_cts(jj,ii)
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  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
496 
497  n_cts_neu(j,i) = -1
498  kc_cts(j,i) = 0
499  zm_neu(j,i) = zb(j,i)
500  h_c_neu(j,i) = h_c(j,i)
501  h_t_neu(j,i) = h_t(j,i)
502 
503  call calc_temp_r(atr1, alb1, i, j)
504 
505  end if
506 
507 end do
508 
509 end subroutine calc_temp_cold
510 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1