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