SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
sico_init.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : s i c o _ i n i t
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
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 !> Initialisations for SICOPOLIS.
34 !<------------------------------------------------------------------------------
35 subroutine sico_init(delta_ts, glac_index, &
36  mean_accum, mean_accum_inv, &
37  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38  time, time_init, time_end, time_output, &
39  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40  z_sl, dzsl_dtau, z_mar, &
41  ii, jj, nn, &
42  ndat2d, ndat3d, n_output, &
43  runname)
44 
45 use sico_types
46 use sico_variables
47 
48 implicit none
49 
50 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
51  jj((imax+1)*(jmax+1)), &
52  nn(0:jmax,0:imax)
53 integer(i4b), intent(out) :: ndat2d, ndat3d
54 integer(i4b), intent(out) :: n_output
55 real(dp), intent(out) :: delta_ts, glac_index
56 real(dp), intent(out) :: mean_accum, mean_accum_inv
57 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
58  dtime_out, dtime_ser
59 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
60 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
61 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
62 character(len=100), intent(out) :: runname
63 
64 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
65 integer(i4b) :: ios
66 integer(i4b) :: ierr
67 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
68 real(dp) :: time_init0, time_end0, time_output0(100)
69 real(dp) :: d_dummy
70 character (len=100) :: anfdatname
71 character (len=256) :: shell_command
72 character :: ch_dummy
73 logical :: flag_3d_output
74 
75 character(len=64), parameter :: fmt1 = '(a)', &
76  fmt2 = '(a,i4)', &
77  fmt3 = '(a,es12.4)'
78 
79 character(len= 8) :: ch_imax
80 character(len=128) :: fmt4
81 
82 write(ch_imax, fmt='(i8)') imax
83 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
84 
85 !-------- Some initial values --------
86 
87 n_output = 0
88 
89 dtime = 0.0_dp
90 dtime_temp = 0.0_dp
91 dtime_wss = 0.0_dp
92 dtime_out = 0.0_dp
93 dtime_ser = 0.0_dp
94 
95 time = 0.0_dp
96 time_init = 0.0_dp
97 time_end = 0.0_dp
98 time_output = 0.0_dp
99 
100 !-------- Initialisation of the Library of Iterative Solvers Lis,
101 ! if required --------
102 
103 #if (CALCZS==4 || MARGIN==3)
104 
105 #include "lisf.h"
106 call lis_initialize(ierr)
107 
108 #endif
109 
110 !-------- Read physical parameters --------
111 
112 call phys_para()
113 
114 !-------- Compatibility check of the SICOPOLIS version with the header file
115 
116 if ( trim(version) /= trim(sico_version) ) &
117  stop ' sico_init: SICOPOLIS version not compatible with header file!'
118 
119 !-------- Compatibility check of the horizontal resolution with the
120 ! number of grid points --------
121 
122 #if (GRID==0 || GRID==1)
123 
124 if (abs(dx-40.0_dp) < eps) then
125 
126  if ((imax /= 150).or.(jmax /= 70)) then
127  stop ' sico_init: IMAX and/or JMAX wrong!'
128  end if
129 
130 else if (abs(dx-20.0_dp) < eps) then
131 
132  if ((imax /= 300).or.(jmax /= 140)) then
133  stop ' sico_init: IMAX and/or JMAX wrong!'
134  end if
135 
136 else if (abs(dx-10.0_dp) < eps) then
137 
138  if ((imax /= 600).or.(jmax /= 280)) then
139  stop ' sico_init: IMAX and/or JMAX wrong!'
140  end if
141 
142 else
143  stop ' sico_init: DX wrong!'
144 end if
145 
146 #elif GRID==2
147 
148  stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
149 
150 #endif
151 
152 !-------- Compatibility check of surface-temperature and precipitation
153 ! determination by interpolation between present and LGM values
154 ! with a glacial index --------
155 
156 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
157 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
158 #endif
159 
160 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
161 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
162 #endif
163 
164 !-------- Compatibility check of discretization schemes for the horizontal and
165 ! vertical advection terms in the temperature and age equations --------
166 
167 #if ADV_HOR==1
168 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
169 #endif
170 
171 !-------- Setting of forcing flag --------
172 
173 #if (TSURFACE <= 4)
174 
175 forcing_flag = 1 ! forcing by delta_ts
176 
177 #elif (TSURFACE == 5)
178 
179 forcing_flag = 2 ! forcing by glac_index
180 
181 #endif
182 
183 !-------- Initialization of numerical time steps --------
184 
185 dtime0 = dtime0
186 dtime_temp0 = dtime_temp0
187 #if REBOUND==2
188 dtime_wss0 = dtime_wss0
189 #endif
190 
191 !-------- Further initializations --------
192 
193 dzeta_c = 1.0_dp/real(kcmax,dp)
194 dzeta_t = 1.0_dp/real(ktmax,dp)
195 dzeta_r = 1.0_dp/real(krmax,dp)
196 
197 ndat2d = 1
198 ndat3d = 1
199 
200 !-------- Construction of a vector (with index n) from a 2-d array
201 ! (with indices i, j) by diagonal numbering --------
202 
203 n=1
204 
205 do m=0, imax+jmax
206  do i=m, 0, -1
207  j = m-i
208  if ((i <= imax).and.(j <= jmax)) then
209  ii(n) = i
210  jj(n) = j
211  nn(j,i) = n
212  n=n+1
213  end if
214  end do
215 end do
216 
217 !-------- Specification of current simulation --------
218 
219 runname = runname
220 anfdatname = anfdatname
221 
222 #if defined(YEAR_ZERO)
223 year_zero = year_zero
224 #else
225 year_zero = 2000.0_dp ! default value 2000 CE
226 #endif
227 
228 time_init0 = time_init0
229 time_end0 = time_end0
230 dtime_ser0 = dtime_ser0
231 
232 #if ( OUTPUT==1 || OUTPUT==3 )
233 dtime_out0 = dtime_out0
234 #endif
235 
236 #if ( OUTPUT==2 || OUTPUT==3 )
237 n_output = n_output
238 time_output0( 1) = time_out0_01
239 time_output0( 2) = time_out0_02
240 time_output0( 3) = time_out0_03
241 time_output0( 4) = time_out0_04
242 time_output0( 5) = time_out0_05
243 time_output0( 6) = time_out0_06
244 time_output0( 7) = time_out0_07
245 time_output0( 8) = time_out0_08
246 time_output0( 9) = time_out0_09
247 time_output0(10) = time_out0_10
248 time_output0(11) = time_out0_11
249 time_output0(12) = time_out0_12
250 time_output0(13) = time_out0_13
251 time_output0(14) = time_out0_14
252 time_output0(15) = time_out0_15
253 time_output0(16) = time_out0_16
254 time_output0(17) = time_out0_17
255 time_output0(18) = time_out0_18
256 time_output0(19) = time_out0_19
257 time_output0(20) = time_out0_20
258 #endif
259 
260 !-------- Write log file --------
261 
262 shell_command = 'if [ ! -d'
263 shell_command = trim(shell_command)//' '//outpath
264 shell_command = trim(shell_command)//' '//'] ; then mkdir'
265 shell_command = trim(shell_command)//' '//outpath
266 shell_command = trim(shell_command)//' '//'; fi'
267 call system(trim(shell_command))
268  ! Check whether directory OUTPATH exists. If not, it is created.
269 
270 open(10, iostat=ios, &
271  file=outpath//'/'//trim(runname)//'.log', &
272  status='new')
273 if (ios /= 0) &
274  stop 'Error when opening the log file!'
275 
276 write(10, fmt=trim(fmt1)) 'Computational domain:'
277 #if defined(ANT)
278 write(10, fmt=trim(fmt1)) 'Antarctica'
279 #elif defined(GRL)
280 write(10, fmt=trim(fmt1)) 'Greenland'
281 #elif defined(SCAND)
282 write(10, fmt=trim(fmt1)) 'Scandinavia and Eurasia'
283 #elif defined(NHEM)
284 write(10, fmt=trim(fmt1)) 'Northern hemisphere'
285 #elif defined(EMTP2SGE)
286 write(10, fmt=trim(fmt1)) 'EISMINT Phase 2 Simplified Geometry Experiment'
287 #elif defined(NMARS)
288 write(10, fmt=trim(fmt1)) 'North polar cap of Mars'
289 #elif defined(SMARS)
290 write(10, fmt=trim(fmt1)) 'South polar cap of Mars'
291 #else
292 stop ' sico_init: Subroutine sico_init: No valid domain specified!'
293 #endif
294 write(10, fmt=trim(fmt1)) ' '
295 
296 write(10, fmt=trim(fmt2)) 'imax =', imax
297 write(10, fmt=trim(fmt2)) 'jmax =', jmax
298 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
299 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
300 write(10, fmt=trim(fmt2)) 'krmax =', krmax
301 write(10, fmt=trim(fmt1)) ' '
302 
303 write(10, fmt=trim(fmt3)) 'a =', deform
304 write(10, fmt=trim(fmt1)) ' '
305 
306 #if (GRID==0 || GRID==1)
307 write(10, fmt=trim(fmt3)) 'x0 =', x0
308 write(10, fmt=trim(fmt3)) 'y0 =', y0
309 write(10, fmt=trim(fmt3)) 'dx =', dx
310 #elif GRID==2
311 stop ' sico_init: GRID==2 not allowed for this application!'
312 #endif
313 write(10, fmt=trim(fmt1)) ' '
314 
315 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
316 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
317 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
318 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
319 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
320 #if REBOUND==2
321 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
322 #endif
323 #if ( OUTPUT==1 || OUTPUT==3 )
324 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
325 #endif
326 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
327 write(10, fmt=trim(fmt1)) ' '
328 
329 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
330 #if ANF_DAT==1
331 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
332 #endif
333 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
334 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
335 write(10, fmt=trim(fmt1)) ' '
336 
337 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
338 write(10, fmt=trim(fmt1)) ' '
339 
340 #if (CALCZS==3 || CALCZS==4)
341 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
342 #if CALCZS==3
343 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
344 #endif
345 write(10, fmt=trim(fmt1)) ' '
346 #endif
347 
348 write(10, fmt=trim(fmt1)) 'temp_mm_present file = '//temp_mm_present_file
349 #if TSURFACE==1
350 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
351 #elif TSURFACE==3
352 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
353 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
354 #elif TSURFACE==4
355 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
356 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
357 #elif TSURFACE==5
358 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
359 write(10, fmt=trim(fmt1)) 'temp_mm_anom file = '//temp_mm_anom_file
360 write(10, fmt=trim(fmt3)) 'temp_mm_anom fact = ', temp_mm_anom_fact
361 #endif
362 
363 write(10, fmt=trim(fmt1)) 'precip_mm_present file = '//precip_mm_present_file
364 #if ACCSURFACE==1
365 write(10, fmt=trim(fmt3)) 'accfact =', accfact
366 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
367 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
368 #elif ( ACCSURFACE==5 )
369 write(10, fmt=trim(fmt1)) 'precip_mm_anom file = '//precip_mm_anom_file
370 write(10, fmt=trim(fmt3)) 'precip_mm_anom fact = ', precip_mm_anom_fact
371 #endif
372 #if (ACCSURFACE <= 3)
373 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
374 #if (ELEV_DESERT == 1)
375 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
376 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
377 #endif
378 #endif
379 
380 #if ABLSURFACE==3
381 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
382 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
383 #endif
384 
385 #if SEA_LEVEL==1
386 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
387 #elif SEA_LEVEL==3
388 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
389 #endif
390 write(10, fmt=trim(fmt1)) ' '
391 
392 #if MARGIN==2
393 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
394 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
395 write(10, fmt=trim(fmt1)) ' '
396 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
397  || marine_ice_calving==6 || marine_ice_calving==7 )
398 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
399 write(10, fmt=trim(fmt1)) ' '
400 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
401 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
402 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
403 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
404 write(10, fmt=trim(fmt1)) ' '
405 #endif
406 #elif MARGIN==3
407 #if ICE_SHELF_CALVING==2
408 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
409 write(10, fmt=trim(fmt1)) ' '
410 #endif
411 #endif
412 
413 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
414 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
415 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
416 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
417 #if SLIDE_LAW==2
418 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
419 #endif
420 write(10, fmt=trim(fmt1)) ' '
421 
422 #if Q_GEO_MOD==1
423 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
424 #elif Q_GEO_MOD==2
425 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
426 #endif
427 write(10, fmt=trim(fmt1)) ' '
428 
429 #if defined(MARINE_ICE_BASAL_MELTING)
430 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
431 #if ( MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3 )
432 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
433 #endif
434 write(10, fmt=trim(fmt1)) ' '
435 #endif
436 
437 #if REBOUND==1
438 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
439 write(10, fmt=trim(fmt1)) ' '
440 #endif
441 
442 #if FLOW_LAW==2
443 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
444 write(10, fmt=trim(fmt1)) ' '
445 #endif
446 #if FIN_VISC==2
447 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
448 write(10, fmt=trim(fmt1)) ' '
449 #endif
450 
451 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
452 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
453 #if ( ENHMOD==2 || ENHMOD==3 )
454 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
455 #endif
456 #if ENHMOD==2
457 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
458 #endif
459 #if ENHMOD==3
460 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
461 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
462 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
463 #endif
464 #if MARGIN==3
465 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
466 #endif
467 write(10, fmt=trim(fmt1)) ' '
468 
469 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
470 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
471 write(10, fmt=trim(fmt3)) 'H_min =', h_min
472 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
473 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
474 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
475 #if defined(QBM_MIN)
476 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
477 #elif defined(QB_MIN) /* obsolete */
478 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
479 #endif
480 #if defined(QBM_MAX)
481 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
482 #elif defined(QB_MAX) /* obsolete */
483 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
484 #endif
485 write(10, fmt=trim(fmt3)) 'age_min =', age_min
486 write(10, fmt=trim(fmt3)) 'age_max =', age_max
487 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
488 #if ADV_VERT==1
489 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
490 #endif
491 write(10, fmt=trim(fmt1)) ' '
492 
493 write(10, fmt=trim(fmt2)) 'GRID =', grid
494 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
495 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
496 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
497 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
498 #if MARGIN==2
499 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
500 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
501 #elif MARGIN==3
502 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
503 #endif
504 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
505 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
506 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
507 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
508 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
509 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
510 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
511 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
512 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
513 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
514 #if ACCSURFACE==5
515 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
516 #endif
517 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
518 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
519 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
520 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
521 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
522 write(10, fmt=trim(fmt1)) ' '
523 
524 #if defined(OUT_TIMES)
525 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
526 #endif
527 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
528 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
529 #if defined(OUTSER_V_A)
530 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
531 #endif
532 #if ( OUTPUT==1 || OUTPUT==2 )
533 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
534 #endif
535 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
536 write(10, fmt=trim(fmt1)) ' '
537 #if ( OUTPUT==2 || OUTPUT==3 )
538 write(10, fmt=trim(fmt2)) 'n_output =', n_output
539 do n=1, n_output
540  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
541 end do
542 write(10, fmt=trim(fmt1)) ' '
543 #endif
544 
545 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
546 
547 close(10, status='keep')
548 
549 !-------- Conversion of time quantities --------
550 
551 year_zero = year_zero*year_sec ! a --> s
552 time_init = time_init0*year_sec ! a --> s
553 time_end = time_end0*year_sec ! a --> s
554 dtime = dtime0*year_sec ! a --> s
555 dtime_temp = dtime_temp0*year_sec ! a --> s
556 #if REBOUND==2
557 dtime_wss = dtime_wss0*year_sec ! a --> s
558 #endif
559 dtime_ser = dtime_ser0*year_sec ! a --> s
560 #if ( OUTPUT==1 || OUTPUT==3 )
561 dtime_out = dtime_out0*year_sec ! a --> s
562 #endif
563 #if ( OUTPUT==2 || OUTPUT==3 )
564 do n=1, n_output
565  time_output(n) = time_output0(n)*year_sec ! a --> s
566 end do
567 #endif
568 
569 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
570  stop ' sico_init: dtime_temp must be a multiple of dtime!'
571 
572 #if REBOUND==2
573 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
574  stop ' sico_init: dtime_wss must be a multiple of dtime!'
575 #endif
576 
577 time = time_init
578 
579 !-------- Reading of measurements for present monthly-mean precipitation --------
580 
581 #if (GRID==0 || GRID==1)
582 
583 open(21, iostat=ios, &
584  file=inpath//'/scand/'//precip_mm_present_file, &
585  recl=8192, status='old')
586 
587 #elif GRID==2
588 
589 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
590 
591 #endif
592 
593 if (ios /= 0) stop ' sico_init: Error when opening the precipitation file!'
594 
595 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
596 
597 do n=1, 12 ! month counter
598  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
599  do j=jmax, 0, -1
600  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
601  end do
602 end do
603 
604 close(21, status='keep')
605 
606 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
607 
608 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
609  ! mm/a water equiv. -> m/s ice equiv.
610 
611 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
612 
613 #if ACCSURFACE==5
614 
615 #if (GRID==0 || GRID==1)
616 
617 open(21, iostat=ios, &
618  file=inpath//'/scand/'//precip_mm_anom_file, &
619  recl=8192, status='old')
620 
621 #endif
622 
623 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
624 
625 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
626 
627 do n=1, 12 ! month counter
628  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
629  do j=jmax, 0, -1
630  read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
631  end do
632 end do
633 
634 close(21, status='keep')
635 
636 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
637 
638 do i=0, imax
639 do j=0, jmax
640 
641 #if (PRECIP_ANOM_INTERPOL==1)
642  do n=1, 12 ! month counter
643  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
644  end do
645 #elif (PRECIP_ANOM_INTERPOL==2)
646  do n=1, 12 ! month counter
647  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
648  end do
649 #else
650  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
651 #endif
652 
653 end do
654 end do
655 
656 #endif
657 
658 !-------- Mean accumulation --------
659 
660 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
661  ! mm/a water equiv. -> m/s ice equiv.
662 
663 mean_accum_inv = 1.0_dp/mean_accum
664 
665 !-------- Reading of present topography mask --------
666 
667 #if (GRID==0 || GRID==1)
668 
669 open(24, iostat=ios, &
670  file=inpath//'/scand/'//mask_present_file, &
671  status='old')
672 
673 #elif GRID==2
674 
675 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
676 
677 #endif
678 
679 if (ios /= 0) &
680  stop ' sico_init: Error when opening the mask file!'
681 
682 do m=1, 6; read(24, fmt='(a)') ch_dummy; end do
683 
684 do j=jmax, 0, -1
685  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
686 end do
687 
688 close(24, status='keep')
689 
690 !-------- Reading of data for present monthly-mean surface temperature --------
691 
692 #if (GRID==0 || GRID==1)
693 
694 open(21, iostat=ios, &
695  file=inpath//'/scand/'//temp_mm_present_file, &
696  recl=8192, status='old')
697 
698 #elif GRID==2
699 
700 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
701 
702 #endif
703 
704 if (ios /= 0) stop ' sico_init: Error when opening the temperature file!'
705 
706 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
707 
708 do n=1, 12 ! month counter
709  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
710  do j=jmax, 0, -1
711  read(21, fmt=*) (temp_mm_present(j,i,n), i=0,imax)
712  end do
713 end do
714 
715 close(21, status='keep')
716 
717 !-------- Reading of LGM monthly-mean surface-temperature anomalies --------
718 
719 #if TSURFACE==5
720 
721 #if (GRID==0 || GRID==1)
722 
723 open(21, iostat=ios, &
724  file=inpath//'/scand/'//temp_mm_anom_file, &
725  recl=8192, status='old')
726 
727 #endif
728 
729 if (ios /= 0) stop ' sico_init: Error when opening the temperature anomaly file!'
730 
731 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
732 
733 do n=1, 12 ! month counter
734  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
735  do j=jmax, 0, -1
736  read(21, fmt=*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
737  end do
738 end do
739 
740 close(21, status='keep')
741 
742 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
743 
744 #endif
745 
746 !-------- Present reference elevation
747 ! (for precipitation and surface-temperature data) --------
748 
749 #if (GRID==0 || GRID==1)
750 
751 open(21, iostat=ios, &
752  file=inpath//'/scand/'//zs_present_file, &
753  recl=8192, status='old')
754 
755 #elif GRID==2
756 
757 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
758 
759 #endif
760 
761 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
762 
763 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
764 
765 do j=jmax, 0, -1
766  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
767 end do
768 
769 close(21, status='keep')
770 
771 ! ------ Reset bathymetry data (sea floor elevation) to the
772 ! sea surface
773 
774 do i=0, imax
775 do j=0, jmax
776  if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
777 end do
778 end do
779 
780 !-------- Read data for delta_ts --------
781 
782 #if TSURFACE==4
783 
784 open(21, iostat=ios, &
785  file=inpath//'/general/'//grip_temp_file, &
786  status='old')
787 if (ios /= 0) &
788  stop ' sico_init: Error when opening the data file for delta_ts!'
789 
790 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
791 
792 if (ch_dummy /= '#') then
793  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
794  write(6, fmt=*) ' not defined in data file for delta_ts!'
795  stop
796 end if
797 
798 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
799 
800 allocate(griptemp(0:ndata_grip))
801 
802 do n=0, ndata_grip
803  read(21, fmt=*) d_dummy, griptemp(n)
804 end do
805 
806 close(21, status='keep')
807 
808 #endif
809 
810 !-------- Read data for the glacial index --------
811 
812 #if TSURFACE==5
813 
814 open(21, iostat=ios, &
815  file=inpath//'/general/'//glac_ind_file, status='old')
816 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
817 
818 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
819 
820 if (ch_dummy /= '#') then
821  write(6, fmt=*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
822  write(6, fmt=*) ' not defined in glacial-index file!'
823  stop
824 end if
825 
826 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
827 
828 allocate(glacial_index(0:ndata_gi))
829 
830 do n=0, ndata_gi
831  read(21, fmt=*) d_dummy, glacial_index(n)
832 end do
833 
834 close(21, status='keep')
835 
836 #endif
837 
838 !-------- Read data for z_sl --------
839 
840 #if SEA_LEVEL==3
841 
842 open(21, iostat=ios, &
843  file=inpath//'/general/'//sea_level_file, &
844  status='old')
845 if (ios /= 0) &
846  stop ' sico_init: Error when opening the data file for z_sl!'
847 
848 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
849 
850 if (ch_dummy /= '#') then
851  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
852  write(6, fmt=*) ' not defined in data file for z_sl!'
853  stop
854 end if
855 
856 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
857 
858 allocate(specmap_zsl(0:ndata_specmap))
859 
860 do n=0, ndata_specmap
861  read(21, fmt=*) d_dummy, specmap_zsl(n)
862 end do
863 
864 close(21, status='keep')
865 
866 #endif
867 
868 !-------- Determination of the geothermal heat flux --------
869 
870 #if Q_GEO_MOD==1
871 
872 ! ------ Constant value
873 
874 do i=0, imax
875 do j=0, jmax
876  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
877 end do
878 end do
879 
880 #elif Q_GEO_MOD==2
881 
882 ! ------ Read data from file
883 
884 open(21, iostat=ios, &
885  file=inpath//'/scand/'//q_geo_file, &
886  recl=8192, status='old')
887 
888 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
889 
890 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
891 
892 do j=jmax, 0, -1
893  read(21, fmt=*) (q_geo(j,i), i=0,imax)
894 end do
895 
896 close(21, status='keep')
897 
898 do i=0, imax
899 do j=0, jmax
900  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
901 end do
902 end do
903 
904 #endif
905 
906 !-------- Reading of tabulated kei function--------
907 
908 #if (REBOUND==0 || REBOUND==1)
909 
910 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
911  ! dummy values
912 #elif REBOUND==2
913 
914 call read_kei()
915 
916 #endif
917 
918 !-------- Definition of initial values --------
919 
920 ! ------ Present topography
921 
922 #if ANF_DAT==1
923 
924 stop ' sico_init: ANF_DAT==1 not allowed for scand application!'
925 
926 ! ------ Ice-free, relaxed bedrock
927 
928 #elif ANF_DAT==2
929 
930 call topography2(dxi, deta)
931 
932 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
933 
934 call boundary(time_init, dtime, dxi, deta, &
935  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
936 
937 do i=0, imax
938 do j=0, jmax
939 
940  q_bm(j,i) = 0.0_dp
941  q_tld(j,i) = 0.0_dp
942  q_b_tot(j,i) = 0.0_dp
943 
944  p_b_w(j,i) = 0.0_dp
945  h_w(j,i) = 0.0_dp
946 
947  do kc=0, kcmax
948  temp_c(kc,j,i) = temp_s(j,i)
949  age_c(kc,j,i) = 15000.0_dp*year_sec
950  end do
951 
952  do kt=0, ktmax
953  omega_t(kt,j,i) = 0.0_dp
954  age_t(kt,j,i) = 15000.0_dp*year_sec
955  end do
956 
957  do kr=0, krmax
958  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
959  *(1.0_dp-real(kr,dp)/real(krmax,dp))
960 ! (linear temperature distribution according to the
961 ! geothermal heat flux)
962  end do
963 
964 end do
965 end do
966 
967 call calc_enhance(time_init)
968 
969 ! ------ Read initial state from output of previous simulation
970 
971 #elif ANF_DAT==3
972 
973 #if NETCDF==1
974 call topography3(dxi, deta, z_sl, anfdatname)
975 #elif NETCDF==2
976 call topography3_nc(dxi, deta, z_sl, anfdatname)
977 #else
978 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
979 #endif
980 
981 call boundary(time_init, dtime, dxi, deta, &
982  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
983 
984 q_b_tot = q_bm + q_tld
985 
986 call calc_enhance(time_init)
987 
988 #endif
989 
990 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
991 
992 #if (GRID==0 || GRID==1)
993 
994 do ir=-imax, imax
995 do jr=-jmax, jmax
996  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
997  ! distortion due to stereographic projection not accounted for
998 end do
999 end do
1000 
1001 #endif
1002 
1003 !-------- General abbreviations --------
1004 
1005 ea = exp(deform)
1006 
1007 do kc=0, kcmax
1008  zeta_c(kc) = kc*dzeta_c
1009  eaz_c(kc) = exp(deform*zeta_c(kc))
1010  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
1011 end do
1012 
1013 do kt=0, ktmax
1014  zeta_t(kt) = kt*dzeta_t
1015 end do
1016 
1017 !-------- Initial velocities --------
1018 
1019 call calc_temp_melt()
1020 
1021 call calc_vxy_b_sia(time, z_sl)
1022 call calc_vxy_sia(dzeta_c, dzeta_t)
1023 
1024 #if MARGIN==3
1025 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
1026 #endif
1027 
1028 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1029 
1030 #if MARGIN==3
1031 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1032 #endif
1033 
1034 !-------- Initialize time-series files --------
1035 
1036 ! ------ Time-series file for the ice sheet on the whole
1037 
1038 open(12, iostat=ios, &
1039  file=outpath//'/'//trim(runname)//'.ser', &
1040  status='new')
1041 if (ios /= 0) &
1042  stop ' sico_init: Error when opening the ser file!'
1043 
1044 if (forcing_flag == 1) then
1045 
1046  write(12,1102)
1047  write(12,1103)
1048 
1049 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1050 
1051  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1052  ' H_max(m) zs_max(m) V_g(m^3)', &
1053  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1054  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1055  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1056  1103 format('----------------------------------------------------', &
1057  '---------------------------------------')
1058 
1059 #elif OUTSER_V_A==2
1060 
1061  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1062  ' V(m^3) V_g(m^3) V_f(m^3)', &
1063  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1064  ' H_max(m) zs_max(m)', &
1065  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1066  ' A_t(m^2) V_bm(m^3/a)', &
1067  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1068  1103 format('----------------------------------------------------', &
1069  '---------------------------------------')
1070 
1071 #else
1072  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1073 #endif
1074 
1075 else if (forcing_flag == 2) then
1076 
1077  write(12,1112)
1078  write(12,1113)
1079 
1080 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1081 
1082  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1083  ' H_max(m) zs_max(m) V_g(m^3)', &
1084  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1085  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1086  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1087  1113 format('----------------------------------------------------', &
1088  '---------------------------------------')
1089 
1090 #elif OUTSER_V_A==2
1091 
1092  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1093  ' V(m^3) V_g(m^3) V_f(m^3)', &
1094  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1095  ' H_max(m) zs_max(m)', &
1096  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1097  ' A_t(m^2) V_bm(m^3/a)', &
1098  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1099  1113 format('----------------------------------------------------', &
1100  '---------------------------------------')
1101 
1102 #else
1103  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1104 #endif
1105 
1106 end if
1107 
1108 ! ------ Time-series file for the ice-thickness field
1109 
1110 #if OUTSER==2
1111 
1112 open(13, iostat=ios, &
1113  file=outpath//'/'//trim(runname)//'.thk', &
1114  status='new')
1115 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1116 
1117 if (forcing_flag == 1) then
1118 
1119  write(13,1104)
1120  write(13,1105)
1121 
1122  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1123  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1124  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1125  1105 format('----------------------------------------------------')
1126 
1127 else if (forcing_flag == 2) then
1128 
1129  write(13,1114)
1130  write(13,1115)
1131 
1132  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1133  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1134  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1135  1115 format('----------------------------------------------------')
1136 
1137 end if
1138 
1139 #endif
1140 
1141 !-------- Output of the initial state --------
1142 
1143 #if OUTPUT==1
1144 
1145 #if ERGDAT==0
1146  flag_3d_output = .false.
1147 #elif ERGDAT==1
1148  flag_3d_output = .true.
1149 #else
1150  stop ' sico_init: ERGDAT must be either 0 or 1!'
1151 #endif
1152 
1153 #if NETCDF==1
1154  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1155  flag_3d_output, ndat2d, ndat3d)
1156 #elif NETCDF==2
1157  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1158  flag_3d_output, ndat2d, ndat3d)
1159 #else
1160  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1161 #endif
1162 
1163 #elif OUTPUT==2
1164 
1165 if (time_output(1) <= time_init+eps) then
1166 
1167 #if ERGDAT==0
1168  flag_3d_output = .false.
1169 #elif ERGDAT==1
1170  flag_3d_output = .true.
1171 #else
1172  stop ' sico_init: ERGDAT must be either 0 or 1!'
1173 #endif
1174 
1175 #if NETCDF==1
1176  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1177  flag_3d_output, ndat2d, ndat3d)
1178 #elif NETCDF==2
1179  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1180  flag_3d_output, ndat2d, ndat3d)
1181 #else
1182  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1183 #endif
1184 
1185 end if
1186 
1187 #elif OUTPUT==3
1188 
1189  flag_3d_output = .false.
1190 
1191 #if NETCDF==1
1192  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1193  flag_3d_output, ndat2d, ndat3d)
1194 #elif NETCDF==2
1195  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1196  flag_3d_output, ndat2d, ndat3d)
1197 #else
1198  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1199 #endif
1200 
1201 if (time_output(1) <= time_init+eps) then
1202 
1203  flag_3d_output = .true.
1204 
1205 #if NETCDF==1
1206  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1207  flag_3d_output, ndat2d, ndat3d)
1208 #elif NETCDF==2
1209  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1210  flag_3d_output, ndat2d, ndat3d)
1211 #else
1212  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1213 #endif
1214 
1215 end if
1216 
1217 #else
1218  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1219 #endif
1220 
1221 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1222 #if OUTSER==2
1223 call output3(time_init, delta_ts, glac_index, z_sl)
1224 #endif
1225 
1226 end subroutine sico_init
1227 !
Generated on Thu Sep 12 2013 09:44:59 for SICOPOLIS V3.1 by   doxygen 1.8.2