SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_age_c.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ a g e _ c
4 !
5 !> @file
6 !!
7 !! Computation of age for a cold 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 ice column.
34 !<------------------------------------------------------------------------------
35 subroutine calc_age_c(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) :: dtime_temp, dxi, deta, dzeta_c, dtau
89 
90 real(dp), parameter :: zero=0.0_dp
91 
92 !-------- Check for boundary points --------
93 
94 if ((i.eq.0).or.(i.eq.imax).or.(j.eq.0).or.(j.eq.jmax)) &
95  stop ' calc_age_c: Boundary points not allowed.'
96 
97 !-------- Some definitions --------
98 
99 dxi = dx *1000.0_dp
100 deta = dx *1000.0_dp
101 dzeta_c = 1.0_dp/dble(kcmax)
102 dtau = dtime_temp
103 
104 g11_g22_g_inv(j,i) = 1.0_dp / ( sq_g11_g(j,i)**2 * sq_g22_g(j,i)**2 )
105 
106 d_sq_g11_g22_g_dxi(j,i) = ( sq_g11_g(j,i+1) * sq_g22_g(j,i+1) &
107  - sq_g11_g(j,i-1) * sq_g22_g(j,i-1) ) &
108  / ( 2.0_dp * dxi )
109 
110 d_sq_g11_g22_g_deta(j,i) = ( sq_g11_g(j+1,i) * sq_g22_g(j+1,i) &
111  - sq_g11_g(j-1,i) * sq_g22_g(j-1,i) ) &
112  / ( 2.0_dp * deta )
113 
114 eaz_c_quarter = exp(deform*(0.25_dp)*dzeta_c)
115 
116 !-------- Further definitions --------
117 
118 do kc=0, kcmax-1
119 
120 vx_c_half(kc,j,i) = ( vz_c(kc,j,i) + vz_c(kc,j,i-1) &
121  + vz_c(kc+1,j,i) + vz_c(kc+1,j,i-1) ) * 0.25_dp
122 vy_c_half(kc,j,i) = ( vz_c(kc,j,i) + vz_c(kc,j-1,i) &
123  + vz_c(kc+1,j,i) + vz_c(kc+1,j-1,i) ) * 0.25_dp
124 
125 age_c_half(kc,j,i) = ( age_c(kc,j,i) + age_c(kc+1,j,i) ) * 0.5_dp
126 eaz_c_half(kc) = exp(deform*(kc+0.5_dp)*dzeta_c)
127 
128 t_c(kc,j,i) = ( (ea-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc))*dzm_dtau(j,i) &
129  + (eaz_c_half(kc)-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc)) &
130  * dh_c_dtau(j,i) )
131 x_c(kc,j,i) = insq_g11_g(j,i) &
132  * ( (ea-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc))*dzm_dxi(j,i) &
133  + (eaz_c_half(kc)-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc)) &
134  * dh_c_dxi(j,i) )
135 y_c(kc,j,i) = insq_g22_g(j,i) &
136  * ( (ea-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc))*dzm_deta(j,i) &
137  + (eaz_c_half(kc)-1.0_dp)/(h_c(j,i)*deform*eaz_c_half(kc)) &
138  * dh_c_deta(j,i))
139 
140 end do
141 
142 !-------- vz at the grid point --------
143 
144  vz_c_g(0,j,i) = vz_c(0,j,i)
145 do kc=1, kcmax
146  vz_c_g(kc,j,i) = ( vz_c(kc,j,i) + vz_c(kc-1,j,i) ) * 0.5_dp
147 end do
148 
149 do kc=0, kcmax
150 
151 !-------- vx and vy at grid points --------
152 
153 vx_c_g(kc,j,i) = ( vx_c(kc,j,i) + vx_c(kc,j,i-1) ) * 0.5_dp
154 vy_c_g(kc,j,i) = ( vy_c(kc,j,i) + vy_c(kc,j-1,i) ) * 0.5_dp
155 
156 !-------- Abbreviations --------
157 
158 t_c_mean(kc,j,i) = ( (ea-1.0_dp)/(h_c(j,i)*eaz_c(kc))*dzm_dtau(j,i) &
159  - 1.0_dp/(h_c(j,i)*eaz_c(kc))*dh_c_dtau(j,i) )
160 x_c_mean(kc,j,i) = insq_g11_g(j,i) &
161  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c(kc))*dzm_dxi(j,i) &
162  - 1.0_dp/(h_c(j,i)*eaz_c(kc))*dh_c_dxi(j,i) )
163 y_c_mean(kc,j,i) = insq_g22_g(j,i) &
164  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c(kc))*dzm_deta(j,i) &
165  - 1.0_dp/(h_c(j,i)*eaz_c(kc))*dh_c_deta(j,i) )
166 
167 !-------- The terms of the mean values --------
168 
169 mean_c_t(kc,j,i) = dtau * t_c_mean(kc,j,i) * age_c(kc,j,i)
170 
171 mean_c_f(kc,j,i) = dtau * ( x_c_mean(kc,j,i) &
172  - d_sq_g11_g22_g_dxi(j,i) &
173  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
174  * age_c(kc,j,i) * vx_c_g(kc,j,i)
175 
176 mean_c_g(kc,j,i) = dtau * ( y_c_mean(kc,j,i) &
177  - d_sq_g11_g22_g_deta(j,i) &
178  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
179  * age_c(kc,j,i) * vy_c_g(kc,j,i)
180 
181 mean_c_h(kc,j,i) = dtau * ( ea-1.0_dp ) / ( h_c(j,i) * eaz_c(kc) ) &
182  * age_c(kc,j,i) * vz_c_g(kc,j,i)
183 
184 !-------- The advection terms --------
185 
186 adv_c_x(kc,j,i) = dtau/(2.0_dp*dxi) &
187  * ( age_c(kc,j,i)*insq_g11_g(j,i) &
188  * (abs(vx_c(kc,j,i))+vx_c(kc,j,i)) &
189  - age_c(kc,j,i+1)*insq_g11_g(j,i+1) &
190  * (abs(vx_c(kc,j,i))-vx_c(kc,j,i)) &
191  - age_c(kc,j,i-1)*insq_g11_g(j,i-1) &
192  * (abs(vx_c(kc,j,i-1))+vx_c(kc,j,i-1)) &
193  + age_c(kc,j,i)*insq_g11_g(j,i) &
194  * (abs(vx_c(kc,j,i-1))-vx_c(kc,j,i-1)) )
195 
196 adv_c_y(kc,j,i) = dtau/(2.0_dp*deta) &
197  * ( age_c(kc,j,i)*insq_g22_g(j,i) &
198  * (abs(vy_c(kc,j,i))+vy_c(kc,j,i)) &
199  - age_c(kc,j+1,i)*insq_g22_g(j+1,i) &
200  * (abs(vy_c(kc,j,i))-vy_c(kc,j,i)) &
201  - age_c(kc,j-1,i)*insq_g22_g(j-1,i) &
202  * (abs(vy_c(kc,j-1,i))+vy_c(kc,j-1,i)) &
203  + age_c(kc,j,i)*insq_g22_g(j,i) &
204  * (abs(vy_c(kc,j-1,i))-vy_c(kc,j-1,i)) )
205 
206 end do
207 
208 
209 do kc=1, kcmax-1
210 
211 !-------- The vertical advection term --------
212 
213 adv_c_z(kc,j,i) = dtau/(2.0_dp*dzeta_c) * (ea-1.0_dp)/(h_c(j,i)*deform) &
214  * ( age_c(kc,j,i)/eaz_c(kc) &
215  * (abs(vz_c(kc,j,i))+vz_c(kc,j,i)) &
216  - age_c(kc+1,j,i)/eaz_c(kc+1) &
217  * (abs(vz_c(kc,j,i))-vz_c(kc,j,i)) &
218  - age_c(kc-1,j,i)/eaz_c(kc-1) &
219  * (abs(vz_c(kc-1,j,i))+vz_c(kc-1,j,i)) &
220  + age_c(kc,j,i)/eaz_c(kc) &
221  * (abs(vz_c(kc-1,j,i))-vz_c(kc-1,j,i)) )
222 
223 !-------- Additional terms --------
224 
225 add_t_c_1(kc,j,i) = dtau/dzeta_c &
226  * ( t_c(kc,j,i) * age_c_half(kc,j,i) &
227  - t_c(kc-1,j,i) * age_c_half(kc-1,j,i) )
228 add_t_c_2(kc,j,i) = dtau/dzeta_c &
229  * ( x_c(kc,j,i)*age_c_half(kc,j,i)*vx_c_half(kc,j,i) &
230  - x_c(kc-1,j,i)*age_c_half(kc-1,j,i) &
231  *vx_c_half(kc-1,j,i) )
232 add_t_c_3(kc,j,i) = dtau/dzeta_c &
233  * ( y_c(kc,j,i)*age_c_half(kc,j,i)*vy_c_half(kc,j,i) &
234  - y_c(kc-1,j,i)*age_c_half(kc-1,j,i) &
235  *vy_c_half(kc-1,j,i) )
236 
237 end do
238 
239 
240 !-------- Abbreviations at kc=0 --------
241 
242 t_c_0(j,i) = ( (ea-1.0_dp)/(h_c(j,i)*deform)*dzm_dtau(j,i) )
243 x_c_0(j,i) = insq_g11_g(j,i) &
244  * ( (ea-1.0_dp)/(h_c(j,i)*deform)*dzm_dxi(j,i) )
245 y_c_0(j,i) = insq_g22_g(j,i) &
246  * ( (ea-1.0_dp)/(h_c(j,i)*deform)*dzm_deta(j,i) )
247 
248 age_c_quarter(j,i) = (3.0_dp*age_c(0,j,i)+age_c(1,j,i))/4.0_dp
249 
250 vx_c_quarter(j,i) = (3.0_dp*vx_c_g(0,j,i)+vx_c_g(1,j,i))/4.0_dp
251 vy_c_quarter(j,i) = (3.0_dp*vy_c_g(0,j,i)+vy_c_g(1,j,i))/4.0_dp
252 
253 t_c_mean_quarter(j,i) = ( (ea-1.0_dp)/(h_c(j,i)*eaz_c_quarter)*dzm_dtau(j,i) &
254  - 1.0_dp/(h_c(j,i)*eaz_c_quarter)*dh_c_dtau(j,i) )
255 x_c_mean_quarter(j,i) = insq_g11_g(j,i) &
256  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c_quarter)*dzm_dxi(j,i) &
257  - 1.0_dp/(h_c(j,i)*eaz_c_quarter)*dh_c_dxi(j,i) )
258 y_c_mean_quarter(j,i) = insq_g22_g(j,i) &
259  * ( (ea-1.0_dp)/(h_c(j,i)*eaz_c_quarter)*dzm_deta(j,i) &
260  - 1.0_dp/(h_c(j,i)*eaz_c_quarter)*dh_c_deta(j,i) )
261 
262 !-------- The vertical advection term at kc= --------
263 
264 adv_c_z_0(j,i) = (2.0_dp*dtau)/dzeta_c * (ea - 1.0_dp)/(h_c(j,i)*deform) &
265  * ( age_c(0,j,i)/2.0_dp &
266  * (abs(vz_c(0,j,i))+vz_c(0,j,i)) &
267  - age_c(1,j,i)/(2.0_dp*eaz_c(1)) &
268  * (abs(vz_c(0,j,i))-vz_c(0,j,i)) &
269  - age_c(0,j,i) &
270  * vz_c(0,j,i))
271 
272 !-------- The terms of the mean values at kc=0 --------
273 
274 mean_c_t_quarter(j,i) = dtau * t_c_mean_quarter(j,i) * age_c_quarter(j,i)
275 
276 mean_c_f_quarter(j,i) = dtau * ( x_c_mean_quarter(j,i) &
277  - d_sq_g11_g22_g_dxi(j,i) &
278  * sq_g22_g(j,i) * g11_g22_g_inv(j,i) ) &
279  * age_c_quarter(j,i) * vx_c_quarter(j,i)
280 
281 mean_c_g_quarter(j,i) = dtau * ( y_c_mean_quarter(j,i) &
282  - d_sq_g11_g22_g_deta(j,i) &
283  * sq_g11_g(j,i) * g11_g22_g_inv(j,i) ) &
284  * age_c_quarter(j,i) * vy_c_quarter(j,i)
285 
286 mean_c_h_quarter(j,i) = dtau * ( ea-1.0_dp)/( h_c(j,i)*eaz_c_quarter ) &
287  * age_c_quarter(j,i) * vz_c_g(0,j,i)
288 
289 !-------- Additional terms at kc=0 --------
290 
291 add_t_c_1_0(j,i) = 2.0_dp*dtau/dzeta_c &
292  * ( t_c(0,j,i)*age_c_half(0,j,i) &
293  - t_c_0(j,i) * age_c(0,j,i) )
294 add_t_c_2_0(j,i) = 2.0_dp*dtau/dzeta_c &
295  * ( x_c(0,j,i)*age_c_half(0,j,i)*vx_c_half(0,j,i) &
296  - x_c_0(j,i)*age_c(0,j,i)*vx_c_g(0,j,i) )
297 add_t_c_3_0(j,i) = 2.0_dp*dtau/dzeta_c &
298  * ( y_c(0,j,i)*age_c_half(0,j,i)*vy_c_half(0,j,i) &
299  - y_c_0(j,i)*age_c(0,j,i)*vy_c_g(0,j,i) )
300 
301 !-------- Solve the age equation at the free surface --------
302 
303 kc=kcmax
304 
305  age_c_neu(kc,j,i) = 0.0_dp
306 
307 !--------
308 
309 do kc=1, kcmax-1
310 
311  age_c_neu(kc,j,i) = age_c(kc,j,i) + dtau &
312  - adv_c_x(kc,j,i) &
313  - adv_c_y(kc,j,i) &
314  - adv_c_z(kc,j,i) &
315  + mean_c_t(kc,j,i) &
316  + mean_c_f(kc,j,i) &
317  + mean_c_g(kc,j,i) &
318  - mean_c_h(kc,j,i) &
319  + add_t_c_1(kc,j,i) &
320  + add_t_c_2(kc,j,i) &
321  + add_t_c_3(kc,j,i)
322 
323 end do
324 
325 !-------- Solve the age equation at bottom layer --------
326 
327 kc=0
328 
329  age_c_neu(0,j,i) = age_c(0,j,i) + dtau &
330  - adv_c_x(0,j,i) &
331  - adv_c_y(0,j,i) &
332  - adv_c_z_0(j,i) &
333  + mean_c_t_quarter(j,i) &
334  + mean_c_f_quarter(j,i) &
335  + mean_c_g_quarter(j,i) &
336  - mean_c_h_quarter(j,i) &
337  + add_t_c_1_0(j,i) &
338  + add_t_c_2_0(j,i) &
339  + add_t_c_3_0(j,i)
340 
341 !-------- Assign the result,
342 ! restriction to interval [0, AGE_MAX yr] --------
343 
344 do kc=0, kcmax
345 
346  if ( (maske_neu(j,i).eq.1).or.(maske_neu(j,i).eq.2) ) &
347  age_c_neu(kc,j,i) = 0.0_dp
348  if (age_c_neu(kc,j,i).lt.zero) &
349  age_c_neu(kc,j,i) = 0.0_dp
350  if (age_c_neu(kc,j,i).gt.(age_max*year_sec)) &
351  age_c_neu(kc,j,i) = age_max*year_sec
352 
353 end do
354 
355 !-------- Age of the ice in the non-existing temperate layer --------
356 
357 do kt=0, ktmax
358  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
359 end do
360 
361 end subroutine calc_age_c
362 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1