SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_age_t.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ a g e _ t
4 !
5 !> @file
6 !!
7 !! Computation of age for a cold / temperate ice column.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Bernd Muegge, 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 age for a cold / temperate ice column.
34 !<------------------------------------------------------------------------------
35 subroutine calc_age_t(dtime_temp, i, j)
36 
37 use sico_types
38 use sico_variables
39 
40 implicit none
41 integer(i4b) :: i, j, kc, kt
42 real(dp) :: g11_g22_g_inv(0:jmax,0:imax)
43 real(dp) :: d_sq_g11_g22_g_dxi(0:jmax,0:imax), &
44  d_sq_g11_g22_g_deta(0:jmax,0:imax)
45 
46 real(dp) :: eaz_c_half(0:kcmax)
47 real(dp) :: vx_c_g(0:kcmax,0:jmax,0:imax), &
48  vy_c_g(0:kcmax,0:jmax,0:imax), &
49  vz_c_g(0:kcmax,0:jmax,0:imax)
50 real(dp) :: vx_c_half(0:kcmax,0:jmax,0:imax), &
51  vy_c_half(0:kcmax,0:jmax,0:imax)
52 real(dp) :: age_c_half(0:kcmax,0:jmax,0:imax)
53 real(dp) :: t_c_mean(0:kcmax,0:jmax,0:imax), &
54  x_c_mean(0:kcmax,0:jmax,0:imax), &
55  y_c_mean(0:kcmax,0:jmax,0:imax)
56 real(dp) :: t_c_mean_quarter(0:jmax,0:imax), &
57  x_c_mean_quarter(0:jmax,0:imax), &
58  y_c_mean_quarter(0:jmax,0:imax)
59 real(dp) :: t_c(0:kcmax,0:jmax,0:imax), &
60  x_c(0:kcmax,0:jmax,0:imax), &
61  y_c(0:kcmax,0:jmax,0:imax)
62 real(dp) :: t_c_0(0:jmax,0:imax), &
63  x_c_0(0:jmax,0:imax), &
64  y_c_0(0:jmax,0:imax)
65 real(dp) :: vx_c_quarter(0:jmax,0:imax), &
66  vy_c_quarter(0:jmax,0:imax)
67 real(dp) :: age_c_quarter(0:jmax,0:imax)
68 real(dp) :: adv_c_x(0:kcmax,0:jmax,0:imax), &
69  adv_c_y(0:kcmax,0:jmax,0:imax), &
70  adv_c_z(0:kcmax,0:jmax,0:imax)
71 real(dp) :: adv_c_z_0(0:jmax,0:imax)
72 real(dp) :: mean_c_t(0:kcmax,0:jmax,0:imax), &
73  mean_c_f(0:kcmax,0:jmax,0:imax), &
74  mean_c_g(0:kcmax,0:jmax,0:imax), &
75  mean_c_h(0:kcmax,0:jmax,0:imax)
76 real(dp) :: mean_c_t_quarter(0:jmax,0:imax), &
77  mean_c_f_quarter(0:jmax,0:imax), &
78  mean_c_g_quarter(0:jmax,0:imax), &
79  mean_c_h_quarter(0:jmax,0:imax)
80 real(dp) :: add_t_c_1(0:kcmax,0:jmax,0:imax), &
81  add_t_c_2(0:kcmax,0:jmax,0:imax), &
82  add_t_c_3(0:kcmax,0:jmax,0:imax)
83 real(dp) :: add_t_c_1_0(0:jmax,0:imax), &
84  add_t_c_2_0(0:jmax,0:imax), &
85  add_t_c_3_0(0:jmax,0:imax)
86 real(dp) :: eaz_c_quarter
87 
88 real(dp) :: vx_t_g(0:ktmax,0:jmax,0:imax), &
89  vy_t_g(0:ktmax,0:jmax,0:imax), &
90  vz_t_g(0:ktmax,0:jmax,0:imax)
91 real(dp) :: vx_t_half(0:ktmax,0:jmax,0:imax), &
92  vy_t_half(0:ktmax,0:jmax,0:imax)
93 real(dp) :: age_t_half(0:ktmax,0:jmax,0:imax)
94 real(dp) :: t_t_mean(0:jmax,0:imax), &
95  x_t_mean(0:jmax,0:imax), &
96  y_t_mean(0:jmax,0:imax)
97 real(dp) :: t_t(0:ktmax,0:jmax,0:imax), &
98  x_t(0:ktmax,0:jmax,0:imax), &
99  y_t(0:ktmax,0:jmax,0:imax)
100 real(dp) :: t_t_0(0:jmax,0:imax), &
101  x_t_0(0:jmax,0:imax), &
102  y_t_0(0:jmax,0:imax)
103 real(dp) :: t_t_max(0:jmax,0:imax), &
104  x_t_max(0:jmax,0:imax), &
105  y_t_max(0:jmax,0:imax)
106 real(dp) :: vx_t_quarter(0:jmax,0:imax), &
107  vy_t_quarter(0:jmax,0:imax)
108 real(dp) :: vx_t_threequarter(0:jmax,0:imax), &
109  vy_t_threequarter(0:jmax,0:imax)
110 real(dp) :: age_t_quarter(0:jmax,0:imax)
111 real(dp) :: age_t_threequarter(0:jmax,0:imax)
112 real(dp) :: adv_t_x(0:ktmax,0:jmax,0:imax), &
113  adv_t_y(0:ktmax,0:jmax,0:imax), &
114  adv_t_z(0:ktmax,0:jmax,0:imax)
115 real(dp) :: adv_t_z_0(0:jmax,0:imax)
116 real(dp) :: adv_t_z_max(0:jmax,0:imax)
117 real(dp) :: mean_t_t(0:ktmax,0:jmax,0:imax), &
118  mean_t_f(0:ktmax,0:jmax,0:imax), &
119  mean_t_g(0:ktmax,0:jmax,0:imax)
120 real(dp) :: mean_t_t_quarter(0:jmax,0:imax), &
121  mean_t_f_quarter(0:jmax,0:imax), &
122  mean_t_g_quarter(0:jmax,0:imax)
123 real(dp) :: mean_t_t_threequarter(0:jmax,0:imax), &
124  mean_t_f_threequarter(0:jmax,0:imax), &
125  mean_t_g_threequarter(0:jmax,0:imax)
126 real(dp) :: add_t_t_1(0:ktmax,0:jmax,0:imax), &
127  add_t_t_2(0:ktmax,0:jmax,0:imax), &
128  add_t_t_3(0:ktmax,0:jmax,0:imax)
129 real(dp) :: add_t_t_1_0(0:jmax,0:imax), &
130  add_t_t_2_0(0:jmax,0:imax), &
131  add_t_t_3_0(0:jmax,0:imax)
132 real(dp) :: add_t_t_1_max(0:jmax,0:imax), &
133  add_t_t_2_max(0:jmax,0:imax), &
134  add_t_t_3_max(0:jmax,0:imax)
135 
136 real(dp) :: dtime_temp, dxi, deta, dzeta_c, dzeta_t, dtau
137 real(dp), parameter :: zero=0.0_dp
138 
139 !-------- Check for boundary points --------
140 
141 if ((i.eq.0).or.(i.eq.imax).or.(j.eq.0).or.(j.eq.jmax)) &
142  stop ' calc_age_t: Boundary points not allowed.'
143 
144 !######## TEMPERATE ICE ##############################################
145 
146 !-------- Some definitions --------
147 
148 dxi = dx *1000.0_dp
149 deta = dx *1000.0_dp
150 dzeta_c = 1.0_dp/real(kcmax,dp)
151 dzeta_t = 1.0_dp/real(ktmax,dp)
152 dtau = dtime_temp
153 
154 g11_g22_g_inv(j,i) = 1.0_dp / ( sq_g11_g(j,i)**2 * sq_g22_g(j,i)**2 )
155 
156 d_sq_g11_g22_g_dxi(j,i) = ( sq_g11_g(j,i+1) * sq_g22_g(j,i+1) &
157  - sq_g11_g(j,i-1) * sq_g22_g(j,i-1) ) &
158  / ( 2.0_dp * dxi )
159 
160 d_sq_g11_g22_g_deta(j,i) = ( sq_g11_g(j+1,i) * sq_g22_g(j+1,i) &
161  - sq_g11_g(j-1,i) * sq_g22_g(j-1,i) ) &
162  / ( 2.0_dp * deta )
163 
164 !-------- Further definitions --------
165 
166 do kt=0, ktmax-1
167 
168 vx_t_half(kt,j,i) = ( vz_t(kt,j,i) + vz_t(kt,j,i-1) &
169  + vz_t(kt+1,j,i) + vz_t(kt+1,j,i-1) ) * 0.25_dp
170 vy_t_half(kt,j,i) = ( vz_t(kt,j,i) + vz_t(kt,j-1,i) &
171  + vz_t(kt+1,j,i) + vz_t(kt+1,j-1,i) ) * 0.25_dp
172 
173 age_t_half(kt,j,i) = ( age_t(kt,j,i) + age_t(kt+1,j,i) ) * 0.5_dp
174 
175 t_t(kt,j,i) = ( dzb_dtau(j,i) + (kt+0.5_dp)*dzeta_t * dh_t_dtau(j,i) )
176 x_t(kt,j,i) = insq_g11_g(j,i) &
177  * ( dzb_dxi(j,i) + (kt+0.5_dp)*dzeta_t * dh_t_dxi(j,i) )
178 y_t(kt,j,i) = insq_g22_g(j,i) &
179  * ( dzb_deta(j,i) + (kt+0.5_dp)*dzeta_t * dh_t_deta(j,i))
180 
181 end do
182 
183 !-------- vz at the grid point --------
184 
185  vz_t_g(0,j,i) = vz_t(0,j,i)
186 do kt=1, ktmax
187  vz_t_g(kt,j,i) = ( vz_t(kt,j,i) + vz_t(kt-1,j,i) ) * 0.5_dp
188 end do
189 
190 do kt=0, ktmax
191 
192 !-------- vx and vy at grid points --------
193 
194 vx_t_g(kt,j,i) = ( vx_t(kt,j,i) + vx_t(kt,j,i-1) ) * 0.5_dp
195 vy_t_g(kt,j,i) = ( vy_t(kt,j,i) + vy_t(kt,j-1,i) ) * 0.5_dp
196 
197 !-------- Abbreviations --------
198 
199 t_t_mean(j,i) = ( 1.0_dp/h_t(j,i)*dh_t_dtau(j,i) )
200 x_t_mean(j,i) = insq_g11_g(j,i) &
201  * ( 1.0_dp/h_t(j,i)*dh_t_dxi(j,i) )
202 y_t_mean(j,i) = insq_g22_g(j,i) &
203  * ( 1.0_dp/h_t(j,i)*dh_t_deta(j,i) )
204 
205 !-------- The terms of the mean values --------
206 
207 mean_t_t(kt,j,i) = dtau * t_t_mean(j,i) * age_t(kt,j,i)
208 
209 mean_t_f(kt,j,i) = dtau * ( x_t_mean(j,i) &
210  - d_sq_g11_g22_g_dxi(j,i) &
211  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
212  * age_t(kt,j,i) * vx_t_g(kt,j,i)
213 
214 mean_t_g(kt,j,i) = dtau * ( y_t_mean(j,i) &
215  - d_sq_g11_g22_g_deta(j,i) &
216  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
217  * age_t(kt,j,i) * vy_t_g(kt,j,i)
218 
219 !-------- The advection terms --------
220 
221 adv_t_x(kt,j,i) = dtau/(2.0_dp*dxi) &
222  * ( age_t(kt,j,i)*insq_g11_g(j,i) &
223  * (abs(vx_t(kt,j,i))+vx_t(kt,j,i)) &
224  - age_t(kt,j,i+1)*insq_g11_g(j,i+1) &
225  * (abs(vx_t(kt,j,i))-vx_t(kt,j,i)) &
226  - age_t(kt,j,i-1)*insq_g11_g(j,i-1) &
227  * (abs(vx_t(kt,j,i-1))+vx_t(kt,j,i-1)) &
228  + age_t(kt,j,i)*insq_g11_g(j,i) &
229  * (abs(vx_t(kt,j,i-1))-vx_t(kt,j,i-1)) )
230 
231 adv_t_y(kt,j,i) = dtau/(2.0_dp*deta) &
232  * ( age_t(kt,j,i)*insq_g22_g(j,i) &
233  * (abs(vy_t(kt,j,i))+vy_t(kt,j,i)) &
234  - age_t(kt,j+1,i)*insq_g22_g(j+1,i) &
235  * (abs(vy_t(kt,j,i))-vy_t(kt,j,i)) &
236  - age_t(kt,j-1,i)*insq_g22_g(j-1,i) &
237  * (abs(vy_t(kt,j-1,i))+vy_t(kt,j-1,i)) &
238  + age_t(kt,j,i)*insq_g22_g(j,i) &
239  * (abs(vy_t(kt,j-1,i))-vy_t(kt,j-1,i)) )
240 
241 end do
242 
243 
244 do kt=1, ktmax-1
245 
246 !-------- The vertical advection term --------
247 
248 adv_t_z(kt,j,i) = dtau/(2.0_dp*dzeta_t) * 1.0_dp/h_t(j,i) &
249  * ( age_t(kt,j,i) * (abs(vz_t(kt,j,i))+vz_t(kt,j,i)) &
250  - age_t(kt+1,j,i) * (abs(vz_t(kt,j,i))-vz_t(kt,j,i)) &
251  - age_t(kt-1,j,i) * (abs(vz_t(kt-1,j,i))+vz_t(kt-1,j,i)) &
252  + age_t(kt,j,i) * (abs(vz_t(kt-1,j,i))-vz_t(kt-1,j,i)) )
253 
254 !-------- Additional terms --------
255 
256 add_t_t_1(kt,j,i) = dtau/dzeta_t &
257  * ( t_t(kt,j,i) * age_t_half(kt,j,i) &
258  - t_t(kt-1,j,i) * age_t_half(kt-1,j,i) )
259 add_t_t_2(kt,j,i) = dtau/dzeta_t &
260  * ( x_t(kt,j,i)*age_t_half(kt,j,i)*vx_t_half(kt,j,i) &
261  - x_t(kt-1,j,i)*age_t_half(kt-1,j,i) &
262  *vx_t_half(kt-1,j,i) )
263 add_t_t_3(kt,j,i) = dtau/dzeta_t &
264  * ( y_t(kt,j,i)*age_t_half(kt,j,i)*vy_t_half(kt,j,i) &
265  - y_t(kt-1,j,i)*age_t_half(kt-1,j,i) &
266  *vy_t_half(kt-1,j,i) )
267 
268 end do
269 
270 
271 !-------- Abbreviations at kt=0 --------
272 
273 t_t_0(j,i) = dzb_dtau(j,i)
274 x_t_0(j,i) = insq_g11_g(j,i) * dzb_dxi(j,i)
275 y_t_0(j,i) = insq_g22_g(j,i) * dzb_deta(j,i)
276 
277 age_t_quarter(j,i) = (3.0_dp*age_t(0,j,i)+age_t(1,j,i))/4.0_dp
278 
279 vx_t_quarter(j,i) = (3.0_dp*vx_t_g(0,j,i)+vx_t_g(1,j,i))/4.0_dp
280 vy_t_quarter(j,i) = (3.0_dp*vy_t_g(0,j,i)+vy_t_g(1,j,i))/4.0_dp
281 
282 !-------- The vertical advection term at kt=0 --------
283 
284 adv_t_z_0(j,i) = (2.0_dp*dtau)/dzeta_t * 1.0_dp/h_t(j,i) &
285  * ( age_t(0,j,i)/2.0_dp * (abs(vz_t(0,j,i))+vz_t(0,j,i)) &
286  - age_t(1,j,i)/2.0_dp * (abs(vz_t(0,j,i))-vz_t(0,j,i)) &
287  - age_t(0,j,i) * vz_t(0,j,i))
288 
289 !-------- The terms of the mean values at kt=0 --------
290 
291 mean_t_t_quarter(j,i) = dtau * t_t_mean(j,i) * age_t_quarter(j,i)
292 
293 mean_t_f_quarter(j,i) = dtau * ( x_t_mean(j,i) &
294  - d_sq_g11_g22_g_dxi(j,i) &
295  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
296  * age_t_quarter(j,i) * vx_t_quarter(j,i)
297 
298 mean_t_g_quarter(j,i) = dtau * ( y_t_mean(j,i) &
299  - d_sq_g11_g22_g_deta(j,i) &
300  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
301  * age_t_quarter(j,i) * vy_t_quarter(j,i)
302 
303 !-------- Additional terms at kt=0 --------
304 
305 add_t_t_1_0(j,i) = 2.0_dp*dtau/dzeta_t &
306  * ( t_t(0,j,i)*age_t_half(0,j,i) &
307  - t_t_0(j,i) * age_t(0,j,i) )
308 add_t_t_2_0(j,i) = 2.0_dp*dtau/dzeta_t &
309  * ( x_t(0,j,i)*age_t_half(0,j,i)*vx_t_half(0,j,i) &
310  - x_t_0(j,i)*age_t(0,j,i)*vx_t_g(0,j,i) )
311 add_t_t_3_0(j,i) = 2.0_dp*dtau/dzeta_t &
312  * ( y_t(0,j,i)*age_t_half(0,j,i)*vy_t_half(0,j,i) &
313  - y_t_0(j,i)*age_t(0,j,i)*vy_t_g(0,j,i) )
314 
315 !-------- Abbreviations at kt=KTMAX --------
316 
317 t_t_max(j,i) = dzb_dtau(j,i) + dh_t_dtau(j,i)
318 x_t_max(j,i) = insq_g11_g(j,i) * ( dzb_dxi(j,i) + dh_t_dxi(j,i) )
319 y_t_max(j,i) = insq_g22_g(j,i) * ( dzb_deta(j,i) + dh_t_deta(j,i) )
320 
321 
322 age_t_threequarter(j,i) = (3.0_dp*age_t(ktmax,j,i)+age_t(ktmax-1,j,i))/4.0_dp
323 
324 vx_t_threequarter(j,i) = (3.0_dp*vx_t_g(ktmax,j,i)+vx_t_g(ktmax-1,j,i))/4.0_dp
325 vy_t_threequarter(j,i) = (3.0_dp*vy_t_g(ktmax,j,i)+vy_t_g(ktmax-1,j,i))/4.0_dp
326 
327 !-------- The vertical advection term at kt=KTMAX --------
328 
329 adv_t_z_max(j,i) = (2.0_dp*dtau)/dzeta_t * 1.0_dp/h_t(j,i) &
330  * ( age_t(ktmax,j,i) * vz_t(ktmax-1,j,i) &
331  - age_t(ktmax-1,j,i)/2.0_dp &
332  * (abs(vz_t(ktmax-1,j,i))+vz_t(ktmax-1,j,i)) &
333  + age_t(ktmax,j,i)/2.0_dp &
334  * (abs(vz_t(ktmax-1,j,i))-vz_t(ktmax-1,j,i)) )
335 
336 
337 !-------- The terms of the mean values at kt=KTMAX --------
338 
339 mean_t_t_threequarter(j,i) = dtau * t_t_mean(j,i) * age_t_threequarter(j,i)
340 
341 mean_t_f_threequarter(j,i) = dtau * ( x_t_mean(j,i) &
342  - d_sq_g11_g22_g_dxi(j,i) &
343  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
344  * age_t_threequarter(j,i) * vx_t_threequarter(j,i)
345 
346 mean_t_g_threequarter(j,i) = dtau * ( y_t_mean(j,i) &
347  - d_sq_g11_g22_g_deta(j,i) &
348  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
349  * age_t_threequarter(j,i) * vy_t_threequarter(j,i)
350 
351 !-------- Additional terms at kt=KTMAX --------
352 
353 add_t_t_1_max(j,i) = 2.0_dp*dtau/dzeta_t &
354  * ( t_t_max(j,i)*age_t(ktmax,j,i) &
355  - t_t(ktmax-1,j,i)*age_t_half(ktmax-1,j,i) )
356 add_t_t_2_max(j,i) = 2.0_dp*dtau/dzeta_t &
357  * ( x_t_max(j,i)*age_t(ktmax,j,i)*vx_t_g(0,j,i) &
358  - x_t(ktmax-1,j,i)*age_t_half(ktmax-1,j,i) &
359  *vx_t_half(ktmax-1,j,i) )
360 add_t_t_3_max(j,i) = 2.0_dp*dtau/dzeta_t &
361  * ( y_t_max(j,i)*age_t(ktmax,j,i)*vy_t_g(0,j,i) &
362  - y_t(ktmax-1,j,i)*age_t_half(ktmax-1,j,i) &
363  *vy_t_half(ktmax-1,j,i) )
364 
365 !######## COLD ICE ###################################################
366 
367 eaz_c_quarter = exp(deform*(0.25_dp)*dzeta_c)
368 
369 !-------- Further definitions --------
370 
371 do kc=0, kcmax-1
372 
373 vx_c_half(kc,j,i) = ( vz_c(kc,j,i) + vz_c(kc,j,i-1) &
374  + vz_c(kc+1,j,i) + vz_c(kc+1,j,i-1) ) * 0.25_dp
375 vy_c_half(kc,j,i) = ( vz_c(kc,j,i) + vz_c(kc,j-1,i) &
376  + vz_c(kc+1,j,i) + vz_c(kc+1,j-1,i) ) * 0.25_dp
377 
378 age_c_half(kc,j,i) = ( age_c(kc,j,i) + age_c(kc+1,j,i) ) * 0.5_dp
379 eaz_c_half(kc) = exp(deform*(kc+0.5_dp)*dzeta_c)
380 
381 t_c(kc,j,i) = ( (ea-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc))*dzm_dtau(j,i) &
382  + (eaz_c_half(kc)-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc)) &
383  * dh_c_dtau(j,i) )
384 x_c(kc,j,i) = insq_g11_g(j,i) &
385  * ( (ea-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc))*dzm_dxi(j,i) &
386  + (eaz_c_half(kc)-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc)) &
387  * dh_c_dxi(j,i) )
388 y_c(kc,j,i) = insq_g22_g(j,i) &
389  * ( (ea-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc))*dzm_deta(j,i) &
390  + (eaz_c_half(kc)-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc)) &
391  * dh_c_deta(j,i))
392 
393 end do
394 
395 !-------- vz at the grid point --------
396 
397  vz_c_g(0,j,i) = vz_c(0,j,i)
398 do kc=1, kcmax
399  vz_c_g(kc,j,i) = ( vz_c(kc,j,i) + vz_c(kc-1,j,i) ) * 0.5_dp
400 end do
401 
402 do kc=0, kcmax
403 
404 !-------- vx and vy at grid points --------
405 
406 vx_c_g(kc,j,i) = ( vx_c(kc,j,i) + vx_c(kc,j,i-1) ) * 0.5_dp
407 vy_c_g(kc,j,i) = ( vy_c(kc,j,i) + vy_c(kc,j-1,i) ) * 0.5_dp
408 
409 !-------- Abbreviations --------
410 
411 t_c_mean(kc,j,i) = ( (ea-1.0_dp)/(h_c(j,i)*eaz_c(kc))*dzm_dtau(j,i) &
412  - 1.0_dp/(h_c(j,i)*eaz_c(kc))*dh_c_dtau(j,i) )
413 x_c_mean(kc,j,i) = insq_g11_g(j,i) &
414  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c(kc))*dzm_dxi(j,i) &
415  - 1.0_dp/(h_c(j,i)*eaz_c(kc))*dh_c_dxi(j,i) )
416 y_c_mean(kc,j,i) = insq_g22_g(j,i) &
417  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c(kc))*dzm_deta(j,i) &
418  - 1.0_dp/(h_c(j,i)*eaz_c(kc))*dh_c_deta(j,i) )
419 
420 !-------- The terms of the mean values --------
421 
422 mean_c_t(kc,j,i) = dtau * t_c_mean(kc,j,i) * age_c(kc,j,i)
423 
424 mean_c_f(kc,j,i) = dtau * ( x_c_mean(kc,j,i) &
425  - d_sq_g11_g22_g_dxi(j,i) &
426  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
427  * age_c(kc,j,i) * vx_c_g(kc,j,i)
428 
429 mean_c_g(kc,j,i) = dtau * ( y_c_mean(kc,j,i) &
430  - d_sq_g11_g22_g_deta(j,i) &
431  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
432  * age_c(kc,j,i) * vy_c_g(kc,j,i)
433 
434 mean_c_h(kc,j,i) = dtau * ( ea-1.0_dp)/(h_c(j,i) *eaz_c(kc) ) &
435  * age_c(kc,j,i) * vz_c_g(kc,j,i)
436 
437 !-------- The advection terms --------
438 
439 adv_c_x(kc,j,i) = dtau/(2.0_dp*dxi) &
440  * ( age_c(kc,j,i)*insq_g11_g(j,i) &
441  * (abs(vx_c(kc,j,i))+vx_c(kc,j,i)) &
442  - age_c(kc,j,i+1)*insq_g11_g(j,i+1) &
443  * (abs(vx_c(kc,j,i))-vx_c(kc,j,i)) &
444  - age_c(kc,j,i-1)*insq_g11_g(j,i-1) &
445  * (abs(vx_c(kc,j,i-1))+vx_c(kc,j,i-1)) &
446  + age_c(kc,j,i)*insq_g11_g(j,i) &
447  * (abs(vx_c(kc,j,i-1))-vx_c(kc,j,i-1)) )
448 
449 adv_c_y(kc,j,i) = dtau/(2.0_dp*deta) &
450  * ( age_c(kc,j,i)*insq_g22_g(j,i) &
451  * (abs(vy_c(kc,j,i))+vy_c(kc,j,i)) &
452  - age_c(kc,j+1,i)*insq_g22_g(j+1,i) &
453  * (abs(vy_c(kc,j,i))-vy_c(kc,j,i)) &
454  - age_c(kc,j-1,i)*insq_g22_g(j-1,i) &
455  * (abs(vy_c(kc,j-1,i))+vy_c(kc,j-1,i)) &
456  + age_c(kc,j,i)*insq_g22_g(j,i) &
457  * (abs(vy_c(kc,j-1,i))-vy_c(kc,j-1,i)) )
458 
459 end do
460 
461 
462 do kc=1, kcmax-1
463 
464 !-------- The vertical advection term --------
465 
466 adv_c_z(kc,j,i) = dtau/(2.0_dp*dzeta_c) * (ea - 1.0_dp)/(h_c(j,i)*deform) &
467  * ( age_c(kc,j,i)/eaz_c(kc) &
468  * (abs(vz_c(kc,j,i))+vz_c(kc,j,i)) &
469  - age_c(kc+1,j,i)/eaz_c(kc+1) &
470  * (abs(vz_c(kc,j,i))-vz_c(kc,j,i)) &
471  - age_c(kc-1,j,i)/eaz_c(kc-1) &
472  * (abs(vz_c(kc-1,j,i))+vz_c(kc-1,j,i)) &
473  + age_c(kc,j,i)/eaz_c(kc) &
474  * (abs(vz_c(kc-1,j,i))-vz_c(kc-1,j,i)) )
475 
476 !-------- Additional terms --------
477 
478 add_t_c_1(kc,j,i) = dtau/dzeta_c &
479  * ( t_c(kc,j,i) * age_c_half(kc,j,i) &
480  - t_c(kc-1,j,i) * age_c_half(kc-1,j,i) )
481 add_t_c_2(kc,j,i) = dtau/dzeta_c &
482  * ( x_c(kc,j,i)*age_c_half(kc,j,i)*vx_c_half(kc,j,i) &
483  - x_c(kc-1,j,i)*age_c_half(kc-1,j,i) &
484  *vx_c_half(kc-1,j,i) )
485 add_t_c_3(kc,j,i) = dtau/dzeta_c &
486  * ( y_c(kc,j,i)*age_c_half(kc,j,i)*vy_c_half(kc,j,i) &
487  - y_c(kc-1,j,i)*age_c_half(kc-1,j,i) &
488  *vy_c_half(kc-1,j,i) )
489 
490 end do
491 
492 
493 !-------- Abbreviations at kc=0 --------
494 
495 t_c_0(j,i) = ( (ea-1.0_dp)/(h_c(j,i)*deform)*dzm_dtau(j,i) )
496 x_c_0(j,i) = insq_g11_g(j,i) &
497  * ( (ea-1.0_dp)/(h_c(j,i)*deform)*dzm_dxi(j,i) )
498 y_c_0(j,i) = insq_g22_g(j,i) &
499  * ( (ea-1.0_dp)/(h_c(j,i)*deform)*dzm_deta(j,i) )
500 
501 age_c_quarter(j,i) = (3.0_dp*age_c(0,j,i)+age_c(1,j,i))/4.0_dp
502 
503 vx_c_quarter(j,i) = (3.0_dp*vx_c_g(0,j,i)+vx_c_g(1,j,i))/4.0_dp
504 vy_c_quarter(j,i) = (3.0_dp*vy_c_g(0,j,i)+vy_c_g(1,j,i))/4.0_dp
505 
506 t_c_mean_quarter(j,i) = ( (ea-1.0_dp)/(h_c(j,i)*eaz_c_quarter)*dzm_dtau(j,i) &
507  - 1.0_dp/(h_c(j,i)*eaz_c_quarter)*dh_c_dtau(j,i) )
508 x_c_mean_quarter(j,i) = insq_g11_g(j,i) &
509  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c_quarter)*dzm_dxi(j,i) &
510  - 1.0_dp/(h_c(j,i)*eaz_c_quarter)*dh_c_dxi(j,i) )
511 y_c_mean_quarter(j,i) = insq_g22_g(j,i) &
512  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c_quarter)*dzm_deta(j,i) &
513  - 1.0_dp/(h_c(j,i)*eaz_c_quarter)*dh_c_deta(j,i) )
514 
515 !-------- The vertical advection term at kc=0 --------
516 
517 adv_c_z_0(j,i) = (2.0_dp*dtau)/dzeta_c * (ea - 1.0_dp)/(h_c(j,i)*deform) &
518  * ( age_c(0,j,i)/2.0_dp &
519  * (abs(vz_c(0,j,i))+vz_c(0,j,i)) &
520  - age_c(1,j,i)/(2.0_dp*eaz_c(1)) &
521  * (abs(vz_c(0,j,i))-vz_c(0,j,i)) &
522  - age_c(0,j,i) &
523  * vz_c(0,j,i))
524 
525 !-------- The terms of the mean values at kc=0 --------
526 
527 mean_c_t_quarter(j,i) = dtau * t_c_mean_quarter(j,i) * age_c_quarter(j,i)
528 
529 mean_c_f_quarter(j,i) = dtau * ( x_c_mean_quarter(j,i) &
530  - d_sq_g11_g22_g_dxi(j,i) &
531  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
532  * age_c_quarter(j,i) * vx_c_quarter(j,i)
533 
534 mean_c_g_quarter(j,i) = dtau * ( y_c_mean_quarter(j,i) &
535  - d_sq_g11_g22_g_deta(j,i) &
536  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
537  * age_c_quarter(j,i) * vy_c_quarter(j,i)
538 
539 mean_c_h_quarter(j,i) = dtau * ( ea-1)/(h_c(j,i) *eaz_c_quarter ) &
540  * age_c_quarter(j,i) * vz_c_g(0,j,i)
541 
542 !-------- Additional terms at kc=0 --------
543 
544 add_t_c_1_0(j,i) = 2.0_dp*dtau/dzeta_c &
545  * ( t_c(0,j,i)*age_c_half(0,j,i) &
546  - t_c_0(j,i) * age_c(0,j,i) )
547 add_t_c_2_0(j,i) = 2.0_dp*dtau/dzeta_c &
548  * ( x_c(0,j,i)*age_c_half(0,j,i)*vx_c_half(0,j,i) &
549  - x_c_0(j,i)*age_c(0,j,i)*vx_c_g(0,j,i) )
550 add_t_c_3_0(j,i) = 2.0_dp*dtau/dzeta_c &
551  * ( y_c(0,j,i)*age_c_half(0,j,i)*vy_c_half(0,j,i) &
552  - y_c_0(j,i)*age_c(0,j,i)*vy_c_g(0,j,i) )
553 
554 !-------- AGE COMPUTATION --------------------------------------------
555 
556 !-------- Solve the age equation at the free surface --------
557 
558 kc=kcmax
559 
560  age_c_neu(kc,j,i) = 0.0_dp
561 
562 !--------
563 
564 do kc=1, kcmax-1
565 
566  age_c_neu(kc,j,i) = age_c(kc,j,i) + dtau &
567  - adv_c_x(kc,j,i) &
568  - adv_c_y(kc,j,i) &
569  - adv_c_z(kc,j,i) &
570  + mean_c_t(kc,j,i) &
571  + mean_c_f(kc,j,i) &
572  + mean_c_g(kc,j,i) &
573  - mean_c_h(kc,j,i) &
574  + add_t_c_1(kc,j,i) &
575  + add_t_c_2(kc,j,i) &
576  + add_t_c_3(kc,j,i)
577 
578 end do
579 
580 !-------- Solve the age equation at the CTS --------
581 
582 kc=0
583 kt=ktmax
584 
585 if (vz_c(0,j,i).le.zero) then
586 
587  age_c_neu(0,j,i) = age_c(0,j,i) + dtau &
588  - adv_c_x(0,j,i) &
589  - adv_c_y(0,j,i) &
590  - adv_c_z_0(j,i) &
591  + mean_c_t_quarter(j,i) &
592  + mean_c_f_quarter(j,i) &
593  + mean_c_g_quarter(j,i) &
594  - mean_c_h_quarter(j,i) &
595  + add_t_c_1_0(j,i) &
596  + add_t_c_2_0(j,i) &
597  + add_t_c_3_0(j,i)
598 
599  age_t_neu(ktmax,j,i) = age_c(0,j,i)
600 
601 else if (vz_t(ktmax-1,j,i).ge.zero) then
602 
603  age_t_neu(ktmax,j,i) = age_t(ktmax,j,i) + dtau &
604  - adv_t_x(ktmax,j,i) &
605  - adv_t_y(ktmax,j,i) &
606  - adv_t_z_max(j,i) &
607  + mean_t_t_threequarter(j,i) &
608  + mean_t_f_threequarter(j,i) &
609  + mean_t_g_threequarter(j,i) &
610  + add_t_t_1_max(j,i) &
611  + add_t_t_2_max(j,i) &
612  + add_t_t_3_max(j,i)
613 
614  age_c(0,j,i) = age_t_neu(ktmax,j,i)
615 
616 end if
617 
618 !--------
619 
620 do kt=1, ktmax-1
621 
622  age_t_neu(kt,j,i) = age_t(kt,j,i) + dtau &
623  - adv_t_x(kt,j,i) &
624  - adv_t_y(kt,j,i) &
625  - adv_t_z(kt,j,i) &
626  + mean_t_t(kt,j,i) &
627  + mean_t_f(kt,j,i) &
628  + mean_t_g(kt,j,i) &
629  + add_t_t_1(kt,j,i) &
630  + add_t_t_2(kt,j,i) &
631  + add_t_t_3(kt,j,i)
632 
633 end do
634 
635 !-------- Solve the age equation at bottom layer --------
636 
637 kt=0
638 
639  age_t_neu(0,j,i) = age_t(0,j,i) + dtau &
640  - adv_t_x(0,j,i) &
641  - adv_t_y(0,j,i) &
642  - adv_t_z_0(j,i) &
643  + mean_t_t_quarter(j,i) &
644  + mean_t_f_quarter(j,i) &
645  + mean_t_g_quarter(j,i) &
646  + add_t_t_1_0(j,i) &
647  + add_t_t_2_0(j,i) &
648  + add_t_t_3_0(j,i)
649 
650 !-------- Assign the result,
651 ! restriction to interval [0, AGE_MAX yr] --------
652 
653 do kc=0, kcmax
654 
655  if ( (maske_neu(j,i).eq.1).or.(maske_neu(j,i).eq.2) ) &
656  age_c_neu(kc,j,i) = 0.0_dp
657  if (age_c_neu(kc,j,i).lt.zero) &
658  age_c_neu(kc,j,i) = 0.0_dp
659  if (age_c_neu(kc,j,i).gt.(age_max*year_sec)) &
660  age_c_neu(kc,j,i) = age_max*year_sec
661 
662 end do
663 
664 do kt=0, ktmax
665 
666  if ( (maske_neu(j,i).eq.1).or.(maske_neu(j,i).eq.2) ) &
667  age_t_neu(kt,j,i) = 0.0_dp
668  if (age_t_neu(kt,j,i).lt.zero) &
669  age_t_neu(kt,j,i) = 0.0_dp
670  if (age_t_neu(kt,j,i).gt.(age_max*year_sec)) &
671  age_t_neu(kt,j,i) = age_max*year_sec
672 
673 end do
674 
675 end subroutine calc_age_t
676 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1