SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
output2.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : o u t p u t 2
4 !
5 !> @file
6 !!
7 !! Writing of time-series data on file in ASCII format.
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 !> Writing of time-series data on file in ASCII format.
34 !<------------------------------------------------------------------------------
35 subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)
36 
37 use sico_types
38 use sico_variables
39 
40 implicit none
41 
42 real(dp), intent(in) :: time, dxi, deta, delta_ts, glac_index, z_sl
43 
44 integer(i4b) :: i, j
45 integer(i4b) :: n_base, n_tempbase
46 real(dp) :: time_val, &
47  v_tot, v_grounded, v_floating, &
48  a_tot, a_grounded, a_floating, &
49  h_max, zs_max, v_temp, v_fw, v_sle, &
50  a_temp, v_bm, v_tld, h_t_max, vs_max, tbh_max
51 real(dp) :: x_pos, y_pos
52 real(dp), dimension(0:JMAX,0:IMAX) :: h_ice, h_t_help
53 real(dp) :: v_gr_redu, a_surf, rhosw_rho_ratio
54 real(dp) :: vs_help, tbh_help
55 real(dp) :: h_ave_sed, tbh_ave_sed, atb_sed
56 real(dp) :: sum_area_sed
57 
58 !-------- Ice thickness --------
59 
60 h_ice = h_c + h_t
61 
62 !-------- Thickness of the temperate layer --------
63 
64 #if CALCMOD==1
65 do i=1, imax-1
66 do j=1, jmax-1
67  h_t_help(j,i) = h_t(j,i)
68 end do
69 end do
70 #elif CALCMOD==0
71 do i=1, imax-1
72 do j=1, jmax-1
73  h_t_help(j,i) = h_c(j,i)*eaz_c_quotient(kc_cts(j,i))
74 end do
75 end do
76 #endif
77 
78 !-------- Maximum ice elevation and thickness, ice volume,
79 ! sea-level equivalent, volume of the temperate ice,
80 ! area, area covered by temperate ice,
81 ! freshwater production due to melting and calving,
82 ! water drainage due to basal melting,
83 ! water drainage from the temperate layer,
84 ! maximum thickness of the temperate layer,
85 ! maximum surface velocity,
86 ! maximum basal temperature rel. to pmp --------
87 
88 #if ( defined(ANT) \
89  || defined(grl) \
90  || defined(scand) \
91  || defined(nhem) \
92  || defined(tibet) \
93  || defined(asf) \
94  || defined(emtp2sge) \
95  || defined(heino) ) /* terrestrial ice sheet */
96 
97 rhosw_rho_ratio = rho_sw/rho
98 
99 #elif ( defined(NMARS) \
100  || defined(smars) ) /* martian ice sheet */
101 
102 rhosw_rho_ratio = 0.0_dp ! dummy value
103 
104 #else
105 
106 stop ' output2: No valid domain specified!'
107 
108 #endif
109 
110 v_grounded = 0.0_dp
111 v_gr_redu = 0.0_dp
112 v_floating = 0.0_dp
113 a_grounded = 0.0_dp
114 a_floating = 0.0_dp
115 h_max = 0.0_dp
116 zs_max = -1.11e+11_dp
117 v_temp = 0.0_dp
118 v_fw = 0.0_dp
119 a_temp = 0.0_dp
120 v_bm = 0.0_dp
121 v_tld = 0.0_dp
122 h_t_max = 0.0_dp
123 vs_max = 0.0_dp
124 tbh_max = -273.15_dp
125 
126 do i=1, imax-1
127 do j=1, jmax-1
128 
129  if (maske(j,i) == 0) then ! grounded ice
130 
131  if (zs(j,i) > zs_max) zs_max = zs(j,i)
132  if (h_ice(j,i) > h_max ) h_max = h_ice(j,i)
133  if (h_t_help(j,i) > h_t_max) h_t_max = h_t_help(j,i)
134 
135  v_grounded = v_grounded + h_ice(j,i) *area(j,i)
136  v_temp = v_temp + h_t_help(j,i)*area(j,i)
137 
138 #if ( defined(ANT) \
139  || defined(grl) \
140  || defined(scand) \
141  || defined(nhem) \
142  || defined(tibet) \
143  || defined(asf) \
144  || defined(emtp2sge) \
145  || defined(heino) ) /* terrestrial ice sheet */
146 
147  if (zl0(j,i) < z_sl) &
148  v_gr_redu = v_gr_redu + rhosw_rho_ratio*(z_sl-zl0(j,i))*area(j,i)
149 
150 #endif
151 
152  a_grounded = a_grounded + area(j,i)
153 
154  if (n_cts(j,i) /= -1) a_temp = a_temp + area(j,i)
155 
156  v_bm = v_bm + q_bm(j,i) *area(j,i)
157  v_tld = v_tld + q_tld(j,i) *area(j,i)
158 
159  vs_help = sqrt( &
160  0.25_dp*(vx_c(kcmax,j,i)+vx_c(kcmax,j,i-1))**2 &
161  +0.25_dp*(vy_c(kcmax,j,i)+vy_c(kcmax,j-1,i))**2 )
162  if (vs_help > vs_max) vs_max = vs_help
163 
164  if (n_cts(j,i) >= 0) then ! temperate base
165  tbh_max = 0.0_dp
166  else ! cold base
167  tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), -eps)
168  if (tbh_help > tbh_max) tbh_max = tbh_help
169  end if
170 
171  else if (maske(j,i) == 3) then ! floating ice
172  ! (basal temperature assumed to be below
173  ! the pressure melting point for pure ice)
174 
175  if (zs(j,i) > zs_max) zs_max = zs(j,i)
176  if (h_ice(j,i) > h_max) h_max = h_ice(j,i)
177 
178  v_floating = v_floating + h_ice(j,i)*area(j,i)
179  a_floating = a_floating + area(j,i)
180 
181  v_bm = v_bm + q_bm(j,i) *area(j,i)
182 
183  vs_help = sqrt( &
184  0.25_dp*(vx_c(kcmax,j,i)+vx_c(kcmax,j,i-1))**2 &
185  +0.25_dp*(vy_c(kcmax,j,i)+vy_c(kcmax,j-1,i))**2 )
186  if (vs_help > vs_max) vs_max = vs_help
187 
188  tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), -eps)
189  if (tbh_help > tbh_max) tbh_max = tbh_help
190 
191  end if
192 
193  if ( (maske(j,i) == 0).or.(maske_help(j,i) == 0) &
194  .or.(maske(j,i) == 3).or.(maske_help(j,i) == 3) ) then
195  ! grounded or floating ice before or after the time step
196  v_fw = v_fw + accum(j,i)*area(j,i) &
197  - (dzs_dtau(j,i)-dzb_dtau(j,i))*area(j,i)
198  end if
199 
200 end do
201 end do
202 
203 !-------- Conversion --------
204 
205 #if ( defined(ANT) \
206  || defined(grl) \
207  || defined(scand) \
208  || defined(nhem) \
209  || defined(tibet) \
210  || defined(asf) \
211  || defined(emtp2sge) \
212  || defined(heino) ) /* terrestrial ice sheet */
213 
214 a_surf = 3.61132e+14_dp ! global ocean area, in m^2
215 v_sle = (v_grounded-v_gr_redu)*(rho/rho_sw)/a_surf
216  ! m^3 ice equiv./m^2 -> m sea water equiv.
217 
218 #elif ( defined(NMARS) \
219  || defined(smars) ) /* martian ice sheet */
220 
221 a_surf = 1.4441e+14_dp ! surface area of Mars, in m^2
222 v_sle = v_grounded*(rho_i/rho_w)*(1.0_dp-frac_dust)/a_surf
223  ! m^3 (ice+dust) equiv./m^2 -> m water equiv.
224 
225 #endif
226 
227 #if ( !defined(OUT_TIMES) || OUT_TIMES==1 )
228 time_val = time /year_sec ! s -> a
229 #elif OUT_TIMES == 2
230 time_val = (time+year_zero) /year_sec ! s -> a
231 #else
232 stop ' output2: OUT_TIMES must be either 1 or 2!'
233 #endif
234 
235 v_tot = v_grounded + v_floating ! in m^3
236 a_tot = a_grounded + a_floating ! in m^2
237 
238 if (zs_max < -1.1e+11_dp) zs_max = 0.0_dp
239 
240 vs_max = vs_max *year_sec ! m/s -> m/a
241 
242 #if ( defined(ANT) \
243  || defined(grl) \
244  || defined(scand) \
245  || defined(nhem) \
246  || defined(tibet) \
247  || defined(asf) \
248  || defined(emtp2sge) \
249  || defined(heino) ) /* terrestrial ice sheet */
250 
251 v_fw = v_fw *year_sec *(rho/rho_w)
252  ! m^3/s ice equiv. -> m^3/a water equiv.
253 v_bm = v_bm *year_sec *(rho/rho_w)
254  ! m^3/s ice equiv. -> m^3/a water equiv.
255 v_tld = v_tld *year_sec *(rho/rho_w)
256  ! m^3/s ice equiv. -> m^3/a water equiv.
257 
258 #elif ( defined(NMARS) \
259  || defined(smars) ) /* martian ice sheet */
260 
261 v_fw = v_fw *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
262  ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
263 v_bm = v_bm *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
264  ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
265 v_tld = v_tld *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
266  ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
267 #endif
268 
269 !-------- Writing of data on time-series file --------
270 
271 #if ( defined(ANT) \
272  || defined(grl) \
273  || defined(scand) \
274  || defined(nhem) \
275  || defined(tibet) \
276  || defined(asf) \
277  || defined(emtp2sge) \
278  || defined(heino) ) /* terrestrial ice sheet */
279 
280 if (forcing_flag == 1) then
281 
282 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
283 
284  write(12,1010) time_val, delta_ts, z_sl, &
285  h_max, zs_max, v_grounded, v_temp, v_fw, v_sle, &
286  a_grounded, a_temp, v_bm, v_tld, h_t_max, vs_max
287 
288 #elif OUTSER_V_A==2
289 
290  write(12,1010) time_val, delta_ts, z_sl, &
291  v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
292  h_max, zs_max, v_temp, v_fw, v_sle, &
293  a_temp, v_bm, v_tld, h_t_max, vs_max
294 
295 #else
296  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
297 #endif
298 
299 else if (forcing_flag == 2) then
300 
301 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
302 
303  write(12,1010) time_val, glac_index, z_sl, &
304  h_max, zs_max, v_grounded, v_temp, v_fw, v_sle, &
305  a_grounded, a_temp, v_bm, v_tld, h_t_max, vs_max
306 
307 #elif OUTSER_V_A==2
308 
309  write(12,1010) time_val, glac_index, z_sl, &
310  v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
311  h_max, zs_max, v_temp, v_fw, v_sle, &
312  a_temp, v_bm, v_tld, h_t_max, vs_max
313 
314 #else
315  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
316 #endif
317 
318 else if (forcing_flag == 3) then
319 
320 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
321 
322  write(12,1010) time_val, 1.11e+11_dp, z_sl, &
323  h_max, zs_max, v_grounded, v_temp, v_fw, v_sle, &
324  a_grounded, a_temp, v_bm, v_tld, h_t_max, vs_max
325 
326 #elif OUTSER_V_A==2
327 
328  write(12,1010) time_val, 1.11e+11_dp, z_sl, &
329  v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
330  h_max, zs_max, v_temp, v_fw, v_sle, &
331  a_temp, v_bm, v_tld, h_t_max, vs_max
332 
333 #else
334  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
335 #endif
336 
337 end if
338 
339 #elif ( defined(NMARS) \
340  || defined(smars) ) /* martian ice sheet */
341 
342 if (forcing_flag == 1) then
343 
344  write(12,1010) time_val, delta_ts, z_sl, &
345  h_max, zs_max, v_tot, v_temp, v_fw, tbh_max, &
346  a_tot, a_temp, v_bm, v_tld, h_t_max, vs_max
347 
348 else if (forcing_flag == 2) then
349 
350  write(12,1010) time_val, glac_index, z_sl, &
351  h_max, zs_max, v_tot, v_temp, v_fw, tbh_max, &
352  a_tot, a_temp, v_bm, v_tld, h_t_max, vs_max
353 
354 else if (forcing_flag == 3) then
355 
356  write(12,1010) time_val, 1.11e+11_dp, z_sl, &
357  h_max, zs_max, v_tot, v_temp, v_fw, tbh_max, &
358  a_tot, a_temp, v_bm, v_tld, h_t_max, vs_max
359 
360 end if
361 
362 #endif
363 
364 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
365 
366 1010 format(1pe13.6,2(1pe13.4),/,13x,6(1pe13.4),/,13x,6(1pe13.4),/)
367 
368 #elif OUTSER_V_A==2
369 
370 1010 format(1pe13.6,2(1pe13.4),/,13x,6(1pe13.4),/, &
371  26x,5(1pe13.4),/,26x,5(1pe13.4),/)
372 
373 #else
374  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
375 #endif
376 
377 !-------- Time-series data for the sediment region of HEINO --------
378 
379 #if defined(HEINO)
380 
381 ! ------ Average ice thickness, average basal temperature rel. to pmp,
382 ! temperate basal area
383 
384 h_ave_sed = 0.0_dp
385 tbh_ave_sed = 0.0_dp
386 atb_sed = 0.0_dp
387 sum_area_sed = 0.0_dp
388 
389 do i=1, imax-1
390 do j=1, jmax-1
391 
392  if (maske_sedi(j,i) == 7) then ! sediment region
393 
394  sum_area_sed = sum_area_sed + area(j,i)
395 
396  h_ave_sed = h_ave_sed + area(j,i)*h_ice(j,i)
397 
398  if (n_cts(j,i) /= -1) then ! temperate base
399  tbh_help = 0.0_dp
400  else ! cold base
401  tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), -eps)
402  end if
403  tbh_ave_sed = tbh_ave_sed + area(j,i)*tbh_help
404 
405  if (n_cts(j,i) /= -1) atb_sed = atb_sed + area(j,i)
406 
407  end if
408 
409 end do
410 end do
411 
412 if (sum_area_sed > eps) then
413  h_ave_sed = h_ave_sed / sum_area_sed
414  tbh_ave_sed = tbh_ave_sed / sum_area_sed
415 else
416  stop ' output2: No sediment area found!'
417 end if
418 
419 ! ------ Writing of data on time-series file
420 
421 if (forcing_flag == 1) then
422  write(15,1011) time_val, delta_ts, z_sl, h_ave_sed, tbh_ave_sed, atb_sed
423 else if (forcing_flag == 2) then
424  write(15,1011) time_val, glac_index, z_sl, h_ave_sed, tbh_ave_sed, atb_sed
425 else if (forcing_flag == 3) then
426  write(15,1011) time_val, 1.11e+11_dp, z_sl, h_ave_sed, tbh_ave_sed, atb_sed
427 end if
428 
429 1011 format(1pe13.6,2(1pe13.4),/,13x,3(1pe13.4),/)
430 
431 #endif
432 
433 end subroutine output2
434 !
Generated on Mon May 20 2013 10:18:42 for SICOPOLIS V3.0 by   doxygen 1.8.1