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