SICOPOLIS V3.0
 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  sum_accum, 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 #if (NETCDF > 1)
48 use netcdf
49 #endif
50 
51 implicit none
52 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
53 integer(i4b) :: ios, ios1, ios2, ios3, ios4
54 integer(i4b) :: ndat2d, ndat3d
55 integer(i4b) :: n_output
56 integer(i4b) :: ierr
57 integer(i4b) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1)), nn(0:jmax,0:imax)
58 real(dp) :: delta_ts, glac_index
59 real(dp) :: sum_accum, mean_accum, mean_accum_inv
60 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
61 real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
62 real(dp) :: time_init0, time_end0, time_output0(100)
63 real(dp) :: time, time_init, time_end, time_output(100)
64 real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
65 real(dp) :: z_sl, dzsl_dtau, z_mar
66 real(dp) :: d_dummy
67 character (len=100) :: runname, anfdatname, target_topo_dat_name
68 character (len=256) :: shell_command
69 character :: ch_dummy
70 logical :: flag_3d_output
71 
72 #if (NETCDF > 1)
73 integer(i4b) :: dimid
74 ! dimid: Dimension ID
75 integer(i4b) :: ncv
76 ! ncv: Variable ID
77 #endif
78 
79 !-------- Initialisation of the Library of Iterative Solvers Lis,
80 ! if required --------
81 
82 #if (CALCZS==4 || MARGIN==3)
83 
84 #include "lisf.h"
85 call lis_initialize(ierr)
86 
87 #endif
88 
89 !-------- Read physical parameters --------
90 
91 call phys_para
92 
93 !-------- Compatibility check of the SICOPOLIS version with the header file
94 
95 if ( trim(version) /= trim(sico_version) ) &
96  stop ' sico_init: SICOPOLIS version not compatible with header file!'
97 
98 !-------- Compatibility check of the horizontal resolution with the
99 ! number of grid points --------
100 
101 #if (GRID==0 || GRID==1)
102 
103 if (abs(dx-80.0_dp) < eps) then
104 
105  if ( (imax /= 70).or.(jmax /= 60) & ! Regular domain
106  ) then
107  stop ' sico_init: IMAX and/or JMAX wrong!'
108  end if
109 
110 else if (abs(dx-40.0_dp) < eps) then
111 
112  if ( ( (imax /= 140).or.(jmax /= 120) ) & ! Regular domain
113  .and. &
114  ( (imax /= 150).or.(jmax /= 150) ) & ! SeaRISE domain
115  ) then
116  stop ' sico_init: IMAX and/or JMAX wrong!'
117  end if
118 
119 else if (abs(dx-20.0_dp) < eps) then
120 
121  if ( ( (imax /= 280).or.(jmax /= 240) ) & ! Regular domain
122  .and. &
123  ( (imax /= 300).or.(jmax /= 300) ) & ! SeaRISE domain
124  ) then
125  stop ' sico_init: IMAX and/or JMAX wrong!'
126  end if
127 
128 else if (abs(dx-10.0_dp) < eps) then
129 
130  if ( (imax /= 600).or.(jmax /= 600) & ! SeaRISE domain
131  ) then
132  stop ' sico_init: IMAX and/or JMAX wrong!'
133  end if
134 
135 else
136  stop ' sico_init: DX wrong!'
137 end if
138 
139 #elif GRID==2
140 
141 stop ' sico_init: GRID==2 not allowed for this application!'
142 
143 #endif
144 
145 !-------- Compatibility check of surface-temperature
146 ! and precipitation switches --------
147 
148 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
149 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
150 #endif
151 
152 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
153 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
154 #endif
155 
156 #if (TSURFACE == 6)
157 #if (ACCSURFACE != 6)
158 write(6,'(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
159 write(6,'(a)') ' and NETCDF==2 must be used together!'
160 stop
161 #elif (NETCDF != 2)
162 write(6,'(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
163 write(6,'(a)') ' and NETCDF==2 must be used together!'
164 stop
165 #endif
166 #endif
167 
168 #if (ACCSURFACE == 6)
169 #if (TSURFACE != 6)
170 write(6,'(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
171 write(6,'(a)') ' and NETCDF==2 must be used together!'
172 stop
173 #elif (NETCDF != 2)
174 write(6,'(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
175 write(6,'(a)') ' and NETCDF==2 must be used together!'
176 stop
177 #endif
178 #endif
179 
180 !-------- Compatibility check of discretization schemes for the horizontal and
181 ! vertical advection terms in the temperature and age equations --------
182 
183 #if ADV_HOR==1
184 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
185 #endif
186 
187 #if ((ADV_HOR == 3) && (ADV_VERT != 3))
188 stop ' sico_init: Options ADV_HOR==3 and ADV_VERT==3 must be used together!'
189 #endif
190 
191 #if ((ADV_HOR != 3) && (ADV_VERT == 3))
192 stop ' sico_init: Options ADV_HOR==3 and ADV_VERT==3 must be used together!'
193 #endif
194 
195 !-------- Check whether for coupled ice-sheet/ice-shelf dynamics
196 ! (MARGIN==3) the chosen grid is in Cartesian coordinates
197 ! without distortion correction (GRID==0) --------
198 
199 #if ( MARGIN == 3 && GRID != 0 )
200 write(6,'(a)') ' sico_init: Option MARGIN==3 requires GRID==0'
201 write(6,'(a)') ' because distortion correction is not'
202 write(6,'(a)') ' implemented for the shallow shelf approximation!'
203 stop
204 #endif
205 
206 !-------- Setting of forcing flag --------
207 
208 #if (TSURFACE <= 4)
209 
210 forcing_flag = 1 ! forcing by delta_ts
211 
212 #elif (TSURFACE == 5)
213 
214 forcing_flag = 2 ! forcing by glac_index
215 
216 #elif (TSURFACE == 6)
217 
218 forcing_flag = 3 ! forcing by time-dependent surface temperature
219  ! and precipitation data
220 
221 #endif
222 
223 !-------- Initialization of numerical time steps --------
224 
225 dtime0 = dtime0
226 dtime_temp0 = dtime_temp0
227 #if REBOUND==2
228 dtime_wss0 = dtime_wss0
229 #endif
230 
231 !-------- Further initializations --------
232 
233 dzeta_c = 1.0_dp/real(kcmax,dp)
234 dzeta_t = 1.0_dp/real(ktmax,dp)
235 dzeta_r = 1.0_dp/real(krmax,dp)
236 
237 ndat2d = 1
238 ndat3d = 1
239 
240 !-------- Construction of a vector (with index n) from a 2-d array
241 ! (with indices i, j) by diagonal numbering --------
242 
243 n=1
244 
245 do m=0, imax+jmax
246  do i=m, 0, -1
247  j = m-i
248  if ((i.le.imax).and.(j.le.jmax)) then
249  ii(n) = i
250  jj(n) = j
251  nn(j,i) = n
252  n=n+1
253  end if
254  end do
255 end do
256 
257 !-------- Specification of current simulation --------
258 
259 runname = runname
260 anfdatname = anfdatname
261 
262 #if defined(YEAR_ZERO)
263 year_zero = year_zero
264 #else
265 year_zero = 2000.0_dp ! default value 2000 CE
266 #endif
267 
268 time_init0 = time_init0
269 time_end0 = time_end0
270 dtime_ser0 = dtime_ser0
271 
272 #if ( OUTPUT==1 || OUTPUT==3 )
273 dtime_out0 = dtime_out0
274 #endif
275 
276 #if ( OUTPUT==2 || OUTPUT==3 )
277 n_output = n_output
278 time_output0( 1) = time_out0_01
279 time_output0( 2) = time_out0_02
280 time_output0( 3) = time_out0_03
281 time_output0( 4) = time_out0_04
282 time_output0( 5) = time_out0_05
283 time_output0( 6) = time_out0_06
284 time_output0( 7) = time_out0_07
285 time_output0( 8) = time_out0_08
286 time_output0( 9) = time_out0_09
287 time_output0(10) = time_out0_10
288 time_output0(11) = time_out0_11
289 time_output0(12) = time_out0_12
290 time_output0(13) = time_out0_13
291 time_output0(14) = time_out0_14
292 time_output0(15) = time_out0_15
293 time_output0(16) = time_out0_16
294 time_output0(17) = time_out0_17
295 time_output0(18) = time_out0_18
296 time_output0(19) = time_out0_19
297 time_output0(20) = time_out0_20
298 #endif
299 
300 !-------- Write log file --------
301 
302 shell_command = 'if [ ! -d'
303 shell_command = trim(shell_command)//' '//outpath
304 shell_command = trim(shell_command)//' '//'] ; then mkdir'
305 shell_command = trim(shell_command)//' '//outpath
306 shell_command = trim(shell_command)//' '//'; fi'
307 call system(trim(shell_command))
308  ! Check whether directory OUTPATH exists. If not, it is created.
309 
310 open(10, iostat=ios, &
311  file=outpath//'/'//trim(runname)//'.log', &
312  status='new')
313 if (ios /= 0) &
314  stop ' sico_init: Error when opening the log file!'
315 
316 write(10,1000) 'Computational domain:'
317 #if defined(ANT)
318 write(10,1000) 'Antarctica'
319 #elif defined(GRL)
320 write(10,1000) 'Greenland'
321 #elif defined(SCAND)
322 write(10,1000) 'Scandinavia and Eurasia'
323 #elif defined(NHEM)
324 write(10,1000) 'Northern hemisphere'
325 #elif defined(EMTP2SGE)
326 write(10,1000) 'EISMINT Phase 2 Simplified Geometry Experiment'
327 #elif defined(NMARS)
328 write(10,1000) 'North polar cap of Mars'
329 #elif defined(SMARS)
330 write(10,1000) 'South polar cap of Mars'
331 #else
332 stop ' sico_init: No valid domain specified!'
333 #endif
334 write(10,1000) ' '
335 
336 write(10,1001) 'imax =', imax
337 write(10,1001) 'jmax =', jmax
338 write(10,1001) 'kcmax =', kcmax
339 write(10,1001) 'ktmax =', ktmax
340 write(10,1001) 'krmax =', krmax
341 write(10,1000) ' '
342 
343 write(10,1002) 'a =', deform
344 write(10,1000) ' '
345 
346 #if (GRID==0 || GRID==1)
347 write(10,1002) 'x0 =', x0
348 write(10,1002) 'y0 =', y0
349 write(10,1002) 'dx =', dx
350 #elif GRID==2
351 stop ' sico_init: GRID==2 not allowed for this application!'
352 #endif
353 write(10,1000) ' '
354 
355 write(10,1002) 'year_zero =', year_zero
356 write(10,1002) 'time_init =', time_init0
357 write(10,1002) 'time_end =', time_end0
358 write(10,1002) 'dtime =', dtime0
359 write(10,1002) 'dtime_temp =', dtime_temp0
360 #if REBOUND==2
361 write(10,1002) 'dtime_wss =', dtime_wss0
362 #endif
363 #if ( OUTPUT==1 || OUTPUT==3 )
364 write(10,1002) 'dtime_out =', dtime_out0
365 #endif
366 write(10,1002) 'dtime_ser =', dtime_ser0
367 write(10,1000) ' '
368 
369 write(10,1000) 'zs_present file = '//zs_present_file
370 #if ANF_DAT==1
371 write(10,1000) 'zl_present file = '//zl_present_file
372 #endif
373 write(10,1000) 'zl0 file = '//zl0_file
374 write(10,1000) 'mask_present file = '//mask_present_file
375 #if ANF_DAT==3
376 write(10,1000) 'Initial-value file = '//anfdatname
377 write(10,1000) 'Path to initial-value file = '//anfdatpath
378 #endif
379 write(10,1000) ' '
380 
381 #if ZS_EVOL==2
382 write(10,1002) 'time_target_topo_init =', time_target_topo_init0
383 write(10,1002) 'time_target_topo_final =', time_target_topo_final0
384 write(10,1000) 'Target-topography file = '//target_topo_dat_name
385 write(10,1000) 'Path to target-topography file = '//target_topo_dat_path
386 write(10,1000) ' '
387 #endif
388 
389 #if ZS_EVOL==3
390 write(10,1000) 'Maximum ice extent mask file = '//mask_maxextent_file
391 write(10,1000) ' '
392 #endif
393 
394 write(10,1000) 'Physical-parameter file = '//phys_para_file
395 write(10,1000) ' '
396 
397 #if (CALCZS==3 || CALCZS==4)
398 write(10,1002) 'ovi_weight =', ovi_weight
399 #if CALCZS==3
400 write(10,1002) 'omega_sor =', omega_sor
401 #endif
402 write(10,1000) ' '
403 #endif
404 
405 #if TSURFACE==1
406 write(10,1002) 'delta_ts0 =', delta_ts0
407 #elif TSURFACE==3
408 write(10,1002) 'sine_amplit =', sine_amplit
409 write(10,1002) 'sine_period =', sine_period
410 #elif TSURFACE==4
411 write(10,1000) 'GRIP file = '//grip_temp_file
412 write(10,1002) 'grip_temp_fact =', grip_temp_fact
413 #elif TSURFACE==5
414 write(10,1000) 'Glacial-index file = '//glac_ind_file
415 write(10,1000) 'temp_ma_anom file = '//temp_ma_anom_file
416 write(10,1002) 'temp_ma_anom fact = ', temp_ma_anom_fact
417 write(10,1000) 'temp_mj_anom file = '//temp_mj_anom_file
418 write(10,1002) 'temp_mj_anom fact = ', temp_mj_anom_fact
419 #elif TSURFACE==6
420 write(10,1000) 'temp_precip_anom file = '//temp_precip_anom_file
421 write(10,1002) 'temp_ma_anom fact = ', temp_ma_anom_fact
422 write(10,1002) 'temp_mj_anom fact = ', temp_mj_anom_fact
423 #endif
424 
425 write(10,1000) 'precip_present file = '//precip_present_file
426 #if ACCSURFACE==1
427 write(10,1002) 'accfact =', accfact
428 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
429 write(10,1002) 'gamma_s =', gamma_s
430 #elif ( ACCSURFACE==5 )
431 write(10,1000) 'precip_anom file = '//precip_anom_file
432 write(10,1002) 'precip_anom fact = ', precip_anom_fact
433 #elif ( ACCSURFACE==6 )
434 write(10,1000) 'temp_precip_anom file = '//temp_precip_anom_file
435 write(10,1002) 'precip_anom fact = ', precip_anom_fact
436 #endif
437 #if (ACCSURFACE <= 3)
438 write(10,1001) 'ELEV_DESERT =', elev_desert
439 #if (ELEV_DESERT == 1)
440 write(10,1002) 'gamma_p =', gamma_p
441 write(10,1002) 'zs_thresh =', zs_thresh
442 #endif
443 #endif
444 
445 #if ABLSURFACE==3
446 write(10,1002) 'lambda_lti =', lambda_lti
447 write(10,1002) 'temp_lti =', temp_lti
448 #endif
449 
450 #if SEA_LEVEL==1
451 write(10,1002) 'z_sl0 =', z_sl0
452 #elif SEA_LEVEL==3
453 write(10,1000) 'sea-level file = '//sea_level_file
454 #endif
455 write(10,1000) ' '
456 
457 #if MARGIN==2
458 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
459 write(10,1002) 'z_mar =', z_mar
460 write(10,1000) ' '
461 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
462  || marine_ice_calving==6 || marine_ice_calving==7 )
463 write(10,1002) 'fact_z_mar =', fact_z_mar
464 write(10,1000) ' '
465 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
466 write(10,1002) 'calv_uw_coeff =', calv_uw_coeff
467 write(10,1002) 'r1_calv_uw =', r1_calv_uw
468 write(10,1002) 'r2_calv_uw =', r2_calv_uw
469 write(10,1000) ' '
470 #endif
471 #elif MARGIN==3
472 #if ICE_SHELF_CALVING==2
473 write(10,1002) 'H_calv =', h_calv
474 write(10,1000) ' '
475 #endif
476 #endif
477 
478 write(10,1002) 'c_slide =', c_slide
479 write(10,1002) 'gamma_slide =', gamma_slide
480 write(10,1001) 'p_weert =', p_weert
481 write(10,1001) 'q_weert =', q_weert
482 #if SLIDE_LAW==2
483 write(10,1002) 'red_pres_limit_fact =', red_pres_limit_fact
484 #endif
485 write(10,1000) ' '
486 
487 #if ( defined(SEDI_SLIDE) && SEDI_SLIDE==2 )
488 write(10,1000) 'Sediment-mask file = '//mask_sedi_file
489 write(10,1000) ' '
490 
491 write(10,1002) 'c_slide_sedi =', c_slide_sedi
492 write(10,1002) 'gamma_slide_sedi =', gamma_slide_sedi
493 write(10,1001) 'p_weert_sedi =', p_weert_sedi
494 write(10,1001) 'q_weert_sedi =', q_weert_sedi
495 write(10,1000) ' '
496 #endif
497 
498 #if Q_GEO_MOD==1
499 write(10,1002) 'q_geo =', q_geo
500 #elif Q_GEO_MOD==2
501 write(10,1000) 'q_geo file = '//q_geo_file
502 #endif
503 write(10,1000) ' '
504 
505 #if defined(MARINE_ICE_BASAL_MELTING)
506 write(10,1001) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
507 #if ( MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3 )
508 write(10,1002) 'qbm_marine =', qbm_marine
509 #endif
510 write(10,1000) ' '
511 #endif
512 
513 #if MARGIN==3
514 write(10,1001) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
515 #if FLOATING_ICE_BASAL_MELTING<=3
516 write(10,1002) 'qbm_float_1 =', qbm_float_1
517 write(10,1002) 'qbm_float_2 =', qbm_float_2
518 write(10,1002) 'qbm_float_3 =', qbm_float_3
519 write(10,1002) 'z_abyss =', z_abyss
520 #endif
521 write(10,1000) ' '
522 #endif
523 
524 #if REBOUND==1
525 write(10,1002) 'frac_llra =', frac_llra
526 write(10,1000) ' '
527 #endif
528 
529 #if FLOW_LAW==2
530 write(10,1002) 'gr_size =', gr_size
531 write(10,1000) ' '
532 #endif
533 #if FIN_VISC==2
534 write(10,1002) 'sigma_res =', sigma_res
535 write(10,1000) ' '
536 #endif
537 
538 write(10,1001) 'ENHMOD =', enhmod
539 write(10,1002) 'enh_fact =', enh_fact
540 #if ( ENHMOD==2 || ENHMOD==3 )
541 write(10,1002) 'enh_intg =', enh_intg
542 #endif
543 #if ENHMOD==2
544 write(10,1002) 'age_trans =', age_trans_0
545 #endif
546 #if ENHMOD==3
547 write(10,1002) 'date_trans1 =', date_trans1_0
548 write(10,1002) 'date_trans2 =', date_trans2_0
549 write(10,1002) 'date_trans3 =', date_trans3_0
550 #endif
551 #if MARGIN==3
552 write(10,1002) 'enh_shelf =', enh_shelf
553 #endif
554 write(10,1000) ' '
555 
556 write(10,1002) 'numdiff_H_t =', numdiff_h_t
557 write(10,1002) 'tau_cts =', tau_cts
558 write(10,1002) 'H_min =', h_min
559 write(10,1002) 'vh_max =', vh_max
560 write(10,1002) 'hd_min =', hd_min
561 write(10,1002) 'hd_max =', hd_max
562 #if defined(QBM_MIN)
563 write(10,1002) 'qbm_min =', qbm_min
564 #elif defined(QB_MIN) /* obsolete */
565 write(10,1002) 'qb_min =', qb_min
566 #endif
567 #if defined(QBM_MAX)
568 write(10,1002) 'qbm_max =', qbm_max
569 #elif defined(QB_MAX) /* obsolete */
570 write(10,1002) 'qb_max =', qb_max
571 #endif
572 write(10,1002) 'age_min =', age_min
573 write(10,1002) 'age_max =', age_max
574 write(10,1002) 'mean_accum =', mean_accum
575 #if ADV_VERT==1
576 write(10,1002) 'age_diff =', agediff
577 write(10,1002) 'm_age =', m_age
578 #endif
579 write(10,1000) ' '
580 
581 write(10,1001) 'GRID =', grid
582 write(10,1001) 'CALCMOD =', calcmod
583 write(10,1001) 'FLOW_LAW =', flow_law
584 write(10,1001) 'FIN_VISC =', fin_visc
585 write(10,1001) 'MARGIN =', margin
586 #if MARGIN==2
587 write(10,1001) 'MARINE_ICE_FORMATION =', marine_ice_formation
588 write(10,1001) 'MARINE_ICE_CALVING =', marine_ice_calving
589 #elif MARGIN==3
590 write(10,1001) 'ICE_SHELF_CALVING =', ice_shelf_calving
591 #endif
592 write(10,1001) 'ANF_DAT =', anf_dat
593 write(10,1001) 'REBOUND =', rebound
594 write(10,1001) 'Q_LITHO =', q_litho
595 write(10,1001) 'ZS_EVOL =', zs_evol
596 write(10,1001) 'CALCZS =', calczs
597 write(10,1001) 'ADV_HOR =', adv_hor
598 write(10,1001) 'ADV_VERT =', adv_vert
599 write(10,1001) 'TOPOGRAD =', topograd
600 write(10,1001) 'TEMP_PRESENT_PARA =', temp_present_para
601 write(10,1001) 'TSURFACE =', tsurface
602 write(10,1001) 'ACCSURFACE =', accsurface
603 #if ACCSURFACE==5
604 write(10,1001) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
605 #endif
606 write(10,1001) 'SOLID_PRECIP =', solid_precip
607 write(10,1001) 'ABLSURFACE =', ablsurface
608 write(10,1001) 'SEA_LEVEL =', sea_level
609 write(10,1001) 'SLIDE_LAW =', slide_law
610 #if defined(SEDI_SLIDE)
611 write(10,1001) 'SEDI_SLIDE =', sedi_slide
612 #endif
613 write(10,1001) 'Q_GEO_MOD =', q_geo_mod
614 write(10,1000) ' '
615 
616 #if defined(OUT_TIMES)
617 write(10,1001) 'OUT_TIMES =', out_times
618 #endif
619 write(10,1001) 'OUTPUT =', output
620 write(10,1001) 'OUTSER =', outser
621 #if defined(OUTSER_V_A)
622 write(10,1001) 'OUTSER_V_A =', outser_v_a
623 #endif
624 #if ( OUTPUT==1 || OUTPUT==2 )
625 write(10,1001) 'ERGDAT =', ergdat
626 #endif
627 write(10,1001) 'NETCDF =', netcdf
628 write(10,1000) ' '
629 #if ( OUTPUT==2 || OUTPUT==3 )
630 write(10,1001) 'n_output =', n_output
631 do n=1, n_output
632  write(10,1002) 'time_output =', time_output0(n)
633 end do
634 write(10,1000) ' '
635 #endif
636 
637 write(10,1000) 'Program version and date: '//version//' / '//date
638 
639  1000 format(a)
640  1001 format(a,i4)
641  1002 format(a,es12.4)
642 
643 close(10, status='keep')
644 
645 !-------- Conversion of time quantities --------
646 
647 year_zero = year_zero*year_sec ! a --> s
648 time_init = time_init0*year_sec ! a --> s
649 time_end = time_end0*year_sec ! a --> s
650 dtime = dtime0*year_sec ! a --> s
651 dtime_temp = dtime_temp0*year_sec ! a --> s
652 #if REBOUND==2
653 dtime_wss = dtime_wss0*year_sec ! a --> s
654 #endif
655 dtime_ser = dtime_ser0*year_sec ! a --> s
656 #if ( OUTPUT==1 || OUTPUT==3 )
657 dtime_out = dtime_out0*year_sec ! a --> s
658 #endif
659 #if ( OUTPUT==2 || OUTPUT==3 )
660 do n=1, n_output
661  time_output(n) = time_output0(n)*year_sec ! a --> s
662 end do
663 #endif
664 
665 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
666  stop ' sico_init: dtime_temp must be a multiple of dtime!'
667 
668 #if REBOUND==2
669 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
670  stop ' sico_init: dtime_wss must be a multiple of dtime!'
671 #endif
672 
673 #if ZS_EVOL==2
674 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
675 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
676 #endif
677 
678 !-------- Reading of present mean-annual precipitation rate --------
679 
680 #if (GRID==0 || GRID==1)
681 
682 open(21, iostat=ios, &
683  file=inpath//'/ant/'//precip_present_file, &
684  recl=8192, status='old')
685 
686 #elif GRID==2
687 
688 stop ' sico_init: GRID==2 not allowed for this application!'
689 
690 #endif
691 
692 if (ios /= 0) stop ' sico_init: Error when opening the precip file!'
693 
694 do n=1, 6; read(21,'(a)') ch_dummy; end do
695 
696 do j=jmax, 0, -1
697  read(21,*) (precip_ma_present(j,i), i=0,imax)
698 end do
699 
700 close(21, status='keep')
701 
702 precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
703  ! mm/a water equiv. -> m/s ice equiv.
704 
705 !-------- Present monthly precipitation rates
706 ! (assumed to be equal to the mean annual precipitation rate) --------
707 
708 do n=1, 12 ! month counter
709 do i=0, imax
710 do j=0, jmax
711  precip_present(j,i,n) = precip_ma_present(j,i)
712 end do
713 end do
714 end do
715 
716 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
717 
718 #if ACCSURFACE==5
719 
720 #if (GRID==0 || GRID==1)
721 
722 open(21, iostat=ios, &
723  file=inpath//'/ant/'//precip_anom_file, &
724  recl=8192, status='old')
725 
726 #endif
727 
728 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
729 
730 do n=1, 6; read(21,'(a)') ch_dummy; end do
731 
732 do j=jmax, 0, -1
733  read(21,*) (precip_ma_lgm_anom(j,i), i=0,imax)
734 end do
735 
736 close(21, status='keep')
737 
738 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
739 
740 !-------- LGM monthly precipitation-rate anomalies (assumed to be
741 ! equal to the mean annual precipitation-rate anomaly) --------
742 
743 do n=1, 12 ! month counter
744 do i=0, imax
745 do j=0, jmax
746 
747  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
748  ! positive values ensured
749 
750 #if (PRECIP_ANOM_INTERPOL==1)
751  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
752 #elif (PRECIP_ANOM_INTERPOL==2)
753  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
754 #else
755  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
756 #endif
757 
758 end do
759 end do
760 end do
761 
762 #endif
763 
764 !-------- Mean accumulation --------
765 
766 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
767  ! mm/a water equiv. -> m/s ice equiv.
768 
769 mean_accum_inv = 1.0_dp/mean_accum
770 
771 !-------- Read present topography mask --------
772 
773 #if (GRID==0 || GRID==1)
774 
775 open(24, iostat=ios, &
776  file=inpath//'/ant/'//mask_present_file, &
777  status='old')
778 
779 #elif GRID==2
780 
781 stop ' sico_init: GRID==2 not allowed for this application!'
782 
783 #endif
784 
785 if (ios /= 0) stop ' sico_init: Error when opening the mask file!'
786 
787 do n=1, 6; read(24,'(a)') ch_dummy; end do
788 
789 do j=jmax, 0, -1
790  read(24,2300) (maske_help(j,i), i=0,imax)
791 end do
792 
793 close(24, status='keep')
794 
795 2300 format(imax(i1),i1)
796 
797 !-------- Read rock/sediment mask --------
798 
799 #if ( defined(SEDI_SLIDE) && SEDI_SLIDE==2 )
800 
801 open(24, iostat=ios, &
802  file=inpath//'/ant/'//mask_sedi_file, &
803  recl=8192, status='old')
804 
805 if (ios /= 0) stop ' sico_init: Error when opening the rock/sediment mask file!'
806 
807 do n=1, 6; read(24,'(a)') ch_dummy; end do
808 
809 do j=jmax, 0, -1
810  read(24,2300) (maske_sedi(j,i), i=0,imax)
811 end do
812 
813 close(24, status='keep')
814 
815 #endif
816 
817 !-------- Present reference elevation (for precipitation data) --------
818 
819 #if (GRID==0 || GRID==1)
820 
821 open(21, iostat=ios, &
822  file=inpath//'/ant/'//zs_present_file, &
823  recl=8192, status='old')
824 
825 #elif GRID==2
826 
827 stop ' sico_init: GRID==2 not allowed for this application!'
828 
829 #endif
830 
831 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
832 
833 do n=1, 6; read(21,'(a)') ch_dummy; end do
834 
835 do j=jmax, 0, -1
836  read(21,*) (zs_ref(j,i), i=0,imax)
837 end do
838 
839 close(21, status='keep')
840 
841 ! ------ Reset bathymetry data (sea floor elevation) to the
842 ! sea surface
843 
844 do i=0, imax
845 do j=0, jmax
846  if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
847 end do
848 end do
849 
850 !-------- Reading of the prescribed target topography --------
851 
852 #if ZS_EVOL==2
853 
854 target_topo_dat_name = trim(target_topo_dat_name)
855 
856 #if NETCDF==1
857 write(6,'(a)') ' sico_init: Reading of target topography'
858 write(6,'(a)') ' only implemented for NetCDF files (NETCDF==2)!'
859 stop
860 #elif NETCDF==2
861 call read_target_topo_nc(target_topo_dat_name)
862 #else
863 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
864 #endif
865 
866 #endif
867 
868 !-------- Reading of the maximum ice extent mask --------
869 
870 #if ZS_EVOL==3
871 
872 #if (GRID==0 || GRID==1)
873 
874 open(24, iostat=ios, &
875  file=inpath//'/ant/'//mask_maxextent_file, &
876  recl=8192, status='old')
877 
878 #elif GRID==2
879 
880 stop ' sico_init: GRID==2 not allowed for this application!'
881 
882 #endif
883 
884 if (ios /= 0) &
885  stop ' sico_init: Error when opening the maximum ice extent mask file!'
886 
887 do n=1, 6; read(24,'(a)') ch_dummy; end do
888 
889 do j=jmax, 0, -1
890  read(24,2300) (maske_maxextent(j,i), i=0,imax)
891 end do
892 
893 close(24, status='keep')
894 
895 #endif
896 
897 !-------- Reading of LGM mean-annual and mean-January (summer)
898 ! surface-temperature anomalies --------
899 
900 #if TSURFACE==5
901 
902 #if (GRID==0 || GRID==1)
903 
904 open(21, iostat=ios, &
905  file=inpath//'/ant/'//temp_ma_anom_file, &
906  recl=8192, status='old')
907 
908 #endif
909 
910 if (ios /= 0) stop ' sico_init: Error when opening the temp_ma anomaly file!'
911 
912 #if (GRID==0 || GRID==1)
913 
914 open(22, iostat=ios, &
915  file=inpath//'/ant/'//temp_mj_anom_file, &
916  recl=8192, status='old')
917 
918 #endif
919 
920 if (ios /= 0) stop ' sico_init: Error when opening the temp_mj anomaly file!'
921 
922 do n=1, 6; read(21,'(a)') ch_dummy; end do
923 do n=1, 6; read(22,'(a)') ch_dummy; end do
924 
925 do j=jmax, 0, -1
926  read(21,*) (temp_ma_lgm_anom(j,i), i=0,imax)
927  read(22,*) (temp_mj_lgm_anom(j,i), i=0,imax)
928 end do
929 
930 close(21, status='keep')
931 close(22, status='keep')
932 
933 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
934 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
935 
936 #endif
937 
938 !-------- Read data for delta_ts --------
939 
940 #if TSURFACE==4
941 
942 open(21, iostat=ios, &
943  file=inpath//'/general/'//grip_temp_file, &
944  status='old')
945 if (ios /= 0) &
946  stop ' sico_init: Error when opening the data file for delta_ts!'
947 
948 read(21,*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
949 
950 if (ch_dummy /= '#') then
951  write(6,*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
952  write(6,*) ' not defined in data file for delta_ts!'
953  stop
954 end if
955 
956 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
957 
958 allocate(griptemp(0:ndata_grip))
959 
960 do n=0, ndata_grip
961  read(21,*) d_dummy, griptemp(n)
962 end do
963 
964 close(21, status='keep')
965 
966 #endif
967 
968 !-------- Read data for the glacial index --------
969 
970 #if TSURFACE==5
971 
972 open(21, iostat=ios, &
973  file=inpath//'/general/'//glac_ind_file, status='old')
974 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
975 
976 read(21,*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
977 
978 if (ch_dummy /= '#') then
979  write(6,*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
980  write(6,*) ' not defined in glacial-index file!'
981  stop
982 end if
983 
984 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
985 
986 allocate(glacial_index(0:ndata_gi))
987 
988 do n=0, ndata_gi
989  read(21,*) d_dummy, glacial_index(n)
990 end do
991 
992 close(21, status='keep')
993 
994 #endif
995 
996 !-------- Reading of time-dependent surface-temperature
997 ! and precipitation anomaly data --------
998 
999 #if ( (TSURFACE==6) && (ACCSURFACE==6) )
1000 
1001 call check( nf90_open(inpath//'/ant/'//temp_precip_anom_file, &
1002  nf90_nowrite, ncid_temp_precip) )
1003 
1004 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid) )
1005 
1006 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1007  len = ndata_temp_precip) )
1008 
1009 ndata_temp_precip = ndata_temp_precip-1
1010 
1011 allocate(temp_precip_time(0:ndata_temp_precip))
1012 
1013 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv) )
1014 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time) ) ! in a
1015 
1016 temp_precip_time_min = temp_precip_time(0) ! in a
1017 temp_precip_time_stp = temp_precip_time(1)-temp_precip_time(0)
1018  ! in a; constant time step assumed
1019 temp_precip_time_max = temp_precip_time(ndata_temp_precip) ! in a
1020 
1021 #endif
1022 
1023 !-------- Read data for z_sl --------
1024 
1025 #if SEA_LEVEL==3
1026 
1027 open(21, iostat=ios, &
1028  file=inpath//'/general/'//sea_level_file, &
1029  status='old')
1030 if (ios /= 0) &
1031  stop ' sico_init: Error when opening the data file for z_sl!'
1032 
1033 read(21,*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1034 
1035 if (ch_dummy /= '#') then
1036  write(6,*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1037  write(6,*) ' not defined in data file for z_sl!'
1038  stop
1039 end if
1040 
1041 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
1042 
1043 allocate(specmap_zsl(0:ndata_specmap))
1044 
1045 do n=0, ndata_specmap
1046  read(21,*) d_dummy, specmap_zsl(n)
1047 end do
1048 
1049 close(21, status='keep')
1050 
1051 #endif
1052 
1053 !-------- Determination of the geothermal heat flux --------
1054 
1055 #if Q_GEO_MOD==1
1056 
1057 ! ------ Constant value
1058 
1059 do i=0, imax
1060 do j=0, jmax
1061  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1062 end do
1063 end do
1064 
1065 #elif Q_GEO_MOD==2
1066 
1067 ! ------ Read data from file
1068 
1069 open(21, iostat=ios, &
1070  file=inpath//'/ant/'//q_geo_file, &
1071  recl=8192, status='old')
1072 
1073 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
1074 
1075 do n=1, 6; read(21,'(a)') ch_dummy; end do
1076 
1077 do j=jmax, 0, -1
1078  read(21,*) (q_geo(j,i), i=0,imax)
1079 end do
1080 
1081 close(21, status='keep')
1082 
1083 do i=0, imax
1084 do j=0, jmax
1085  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1086 end do
1087 end do
1088 
1089 #endif
1090 
1091 !-------- Reading of tabulated kei function--------
1092 
1093 #if (REBOUND==0 || REBOUND==1)
1094 
1095 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1096  ! dummy values
1097 #elif REBOUND==2
1098 
1099 call read_kei
1100 
1101 #endif
1102 
1103 !-------- Definition of initial values --------
1104 
1105 ! ------ Present topography
1106 
1107 #if ANF_DAT==1
1108 
1109 call topography1(dxi, deta)
1110 
1111 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1112 
1113 call boundary(time_init, dtime, dxi, deta, &
1114  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1115 
1116 do i=0, imax
1117 do j=0, jmax
1118 
1119  q_bm(j,i) = 0.0_dp
1120  q_tld(j,i) = 0.0_dp
1121  q_b_tot(j,i) = 0.0_dp
1122 
1123  p_b_w(j,i) = 0.0_dp
1124  h_w(j,i) = 0.0_dp
1125 
1126  do kc=0, kcmax
1127  temp_c(kc,j,i) = -10.0_dp
1128  age_c(kc,j,i) = 15000.0_dp*year_sec
1129  end do
1130 
1131  do kt=0, ktmax
1132  omega_t(kt,j,i) = 0.0_dp
1133  age_t(kt,j,i) = 15000.0_dp*year_sec
1134  end do
1135 
1136  do kr=0, krmax
1137  temp_r(kr,j,i) = -10.0_dp+q_geo(j,i)*h_r/kappa_r &
1138  *(1.0_dp-real(kr,dp)/real(krmax,dp))
1139 ! (linear temperature distribution according to the
1140 ! geothermal heat flux)
1141  end do
1142 
1143 end do
1144 end do
1145 
1146 call calc_enhance(time_init)
1147 
1148 ! ------ Ice-free, relaxed bedrock
1149 
1150 #elif ANF_DAT==2
1151 
1152 call topography2(dxi, deta)
1153 
1154 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1155 
1156 call boundary(time_init, dtime, dxi, deta, &
1157  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1158 
1159 do i=0, imax
1160 do j=0, jmax
1161 
1162  q_bm(j,i) = 0.0_dp
1163  q_tld(j,i) = 0.0_dp
1164  q_b_tot(j,i) = 0.0_dp
1165 
1166  p_b_w(j,i) = 0.0_dp
1167  h_w(j,i) = 0.0_dp
1168 
1169  do kc=0, kcmax
1170  temp_c(kc,j,i) = temp_s(j,i)
1171  age_c(kc,j,i) = 15000.0_dp*year_sec
1172  end do
1173 
1174  do kt=0, ktmax
1175  omega_t(kt,j,i) = 0.0_dp
1176  age_t(kt,j,i) = 15000.0_dp*year_sec
1177  end do
1178 
1179  do kr=0, krmax
1180  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
1181  *(1.0_dp-real(kr,dp)/real(krmax,dp))
1182 ! (linear temperature distribution according to the
1183 ! geothermal heat flux)
1184  end do
1185 
1186 end do
1187 end do
1188 
1189 call calc_enhance(time_init)
1190 
1191 ! ------ Read initial state from output of previous simulation
1192 
1193 #elif ANF_DAT==3
1194 
1195 #if NETCDF==1
1196 call topography3(dxi, deta, z_sl, anfdatname)
1197 #elif NETCDF==2
1198 call topography3_nc(dxi, deta, z_sl, anfdatname)
1199 #else
1200 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1201 #endif
1202 
1203 call boundary(time_init, dtime, dxi, deta, &
1204  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1205 
1206 q_b_tot = q_bm + q_tld
1207 
1208 call calc_enhance(time_init)
1209 
1210 #endif
1211 
1212 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1213 
1214 #if (GRID==0 || GRID==1)
1215 
1216 do ir=-imax, imax
1217 do jr=-jmax, jmax
1218  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1219  ! distortion due to stereographic projection not accounted for
1220 end do
1221 end do
1222 
1223 #endif
1224 
1225 !-------- General abbreviations --------
1226 
1227 ea = exp(deform)
1228 
1229 do kc=0, kcmax
1230  zeta_c(kc) = kc*dzeta_c
1231  eaz_c(kc) = exp(deform*zeta_c(kc))
1232  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
1233 end do
1234 
1235 do kt=0, ktmax
1236  zeta_t(kt) = kt*dzeta_t
1237 end do
1238 
1239 !-------- Initial velocities --------
1240 
1241 call calc_temp_melt
1242 
1243 call calc_vxy_b_sia(time, z_sl)
1244 call calc_vxy_sia(dzeta_c, dzeta_t)
1245 
1246 #if MARGIN==3
1247 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
1248 #endif
1249 
1250 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1251 
1252 #if MARGIN==3
1253 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1254 #endif
1255 
1256 !-------- Initialize time-series files --------
1257 
1258 ! ------ Time-series file for the ice sheet on the whole
1259 
1260 open(12, iostat=ios, &
1261  file=outpath//'/'//trim(runname)//'.ser', &
1262  status='new')
1263 if (ios /= 0) &
1264  stop ' sico_init: Error when opening the ser file!'
1265 
1266 if (forcing_flag == 1) then
1267 
1268  write(12,1102)
1269  write(12,1103)
1270 
1271 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1272 
1273  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1274  ' H_max(m) zs_max(m) V_g(m^3)', &
1275  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1276  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1277  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1278  1103 format('----------------------------------------------------', &
1279  '---------------------------------------')
1280 
1281 #elif OUTSER_V_A==2
1282 
1283  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1284  ' V(m^3) V_g(m^3) V_f(m^3)', &
1285  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1286  ' H_max(m) zs_max(m)', &
1287  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1288  ' A_t(m^2) V_bm(m^3/a)', &
1289  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1290  1103 format('----------------------------------------------------', &
1291  '---------------------------------------')
1292 
1293 #else
1294  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1295 #endif
1296 
1297 else if (forcing_flag == 2) then
1298 
1299  write(12,1112)
1300  write(12,1113)
1301 
1302 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1303 
1304  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1305  ' H_max(m) zs_max(m) V_g(m^3)', &
1306  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1307  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1308  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1309  1113 format('----------------------------------------------------', &
1310  '---------------------------------------')
1311 
1312 #elif OUTSER_V_A==2
1313 
1314  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1315  ' V(m^3) V_g(m^3) V_f(m^3)', &
1316  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1317  ' H_max(m) zs_max(m)', &
1318  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1319  ' A_t(m^2) V_bm(m^3/a)', &
1320  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1321  1113 format('----------------------------------------------------', &
1322  '---------------------------------------')
1323 
1324 #else
1325  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1326 #endif
1327 
1328 else if (forcing_flag == 3) then
1329 
1330  write(12,1122)
1331  write(12,1123)
1332 
1333 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1334 
1335  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1336  ' H_max(m) zs_max(m) V_g(m^3)', &
1337  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1338  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1339  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1340  1123 format('----------------------------------------------------', &
1341  '---------------------------------------')
1342 
1343 #elif OUTSER_V_A==2
1344 
1345  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1346  ' V(m^3) V_g(m^3) V_f(m^3)', &
1347  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1348  ' H_max(m) zs_max(m)', &
1349  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1350  ' A_t(m^2) V_bm(m^3/a)', &
1351  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1352  1123 format('----------------------------------------------------', &
1353  '---------------------------------------')
1354 
1355 #else
1356  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1357 #endif
1358 
1359 end if
1360 
1361 ! ------ Time-series file for the ice-thickness field
1362 
1363 #if OUTSER==2
1364 
1365 open(13, iostat=ios, &
1366  file=outpath//'/'//trim(runname)//'.thk', &
1367  status='new')
1368 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1369 
1370 if (forcing_flag == 1) then
1371 
1372  write(13,1104)
1373  write(13,1105)
1374 
1375  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1376  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1377  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1378  1105 format('----------------------------------------------------')
1379 
1380 else if (forcing_flag == 2) then
1381 
1382  write(13,1114)
1383  write(13,1115)
1384 
1385  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1386  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1387  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1388  1115 format('----------------------------------------------------')
1389 
1390 else if (forcing_flag == 3) then
1391 
1392  write(13,1124)
1393  write(13,1125)
1394 
1395  1124 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1396  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1397  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1398  1125 format('----------------------------------------------------')
1399 
1400 end if
1401 
1402 #endif
1403 
1404 ! ------ Time-series file for deep boreholes
1405 
1406 #if OUTSER==3
1407 
1408 n_core = 6 ! Vostok, Dome A, Dome C, Dome F, Kohnen, Byrd
1409 
1410 allocate(lambda_core(n_core), phi_core(n_core), &
1411  x_core(n_core), y_core(n_core))
1412 
1413 phi_core(1) = -78.467_dp *pi_180 ! Geographical position of Vostok,
1414 lambda_core(1) = 106.8_dp *pi_180 ! conversion deg -> rad
1415 call num_coord(lambda_core(1), phi_core(1), x_core(1), y_core(1))
1416 
1417 phi_core(2) = -80.367_dp *pi_180 ! Geographical position of Dome A,
1418 lambda_core(2) = 77.35_dp *pi_180 ! conversion deg -> rad
1419 call num_coord(lambda_core(2), phi_core(2), x_core(2), y_core(2))
1420 
1421 phi_core(3) = -75.1_dp *pi_180 ! Geographical position of Dome C,
1422 lambda_core(3) = 123.4_dp *pi_180 ! conversion deg -> rad
1423 call num_coord(lambda_core(3), phi_core(3), x_core(3), y_core(3))
1424 
1425 phi_core(4) = -77.317_dp *pi_180 ! Geographical position of Dome F,
1426 lambda_core(4) = 39.7_dp *pi_180 ! conversion deg -> rad
1427 call num_coord(lambda_core(4), phi_core(4), x_core(4), y_core(4))
1428 
1429 phi_core(5) = -75.0_dp *pi_180 ! Geographical position of Kohnen,
1430 lambda_core(5) = 0.067_dp *pi_180 ! conversion deg -> rad
1431 call num_coord(lambda_core(5), phi_core(5), x_core(5), y_core(5))
1432 
1433 phi_core(6) = -80.017_dp *pi_180 ! Geographical position of Byrd,
1434 lambda_core(6) = -119.517_dp *pi_180 ! conversion deg -> rad
1435 call num_coord(lambda_core(6), phi_core(6), x_core(6), y_core(6))
1436 
1437 open(14, iostat=ios, &
1438  file=outpath//'/'//trim(runname)//'.core', &
1439  status='new')
1440 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
1441 
1442 if (forcing_flag == 1) then
1443 
1444  write(14,1106)
1445  write(14,1107)
1446 
1447  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1448  ' H_Vo(m) H_DA(m) H_DC(m)', &
1449  ' H_DF(m) H_Ko(m) H_By(m)',/, &
1450  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
1451  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
1452  ' T_Vo(C) T_DA(C) T_DC(C)', &
1453  ' T_DF(C) T_Ko(C) T_By(C)',/, &
1454  ' Rb_Vo(W/m2) Rb_DA(W/m2) Rb_DC(W/m2)', &
1455  ' Rb_DF(W/m2) Rb_Ko(W/m2) Rb_By(W/m2)')
1456  1107 format('----------------------------------------------------', &
1457  '---------------------------------------')
1458 
1459 else if (forcing_flag == 2) then
1460 
1461  write(14,1116)
1462  write(14,1117)
1463 
1464  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
1465  ' H_Vo(m) H_DA(m) H_DC(m)', &
1466  ' H_DF(m) H_Ko(m) H_By(m)',/, &
1467  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
1468  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
1469  ' T_Vo(C) T_DA(C) T_DC(C)', &
1470  ' T_DF(C) T_Ko(C) T_By(C)',/, &
1471  ' Rb_Vo(W/m2) Rb_DA(W/m2) Rb_DC(W/m2)', &
1472  ' Rb_DF(W/m2) Rb_Ko(W/m2) Rb_By(W/m2)')
1473  1117 format('----------------------------------------------------', &
1474  '---------------------------------------')
1475 
1476 else if (forcing_flag == 3) then
1477 
1478  write(14,1126)
1479  write(14,1127)
1480 
1481  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1482  ' H_Vo(m) H_DA(m) H_DC(m)', &
1483  ' H_DF(m) H_Ko(m) H_By(m)',/, &
1484  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
1485  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
1486  ' T_Vo(C) T_DA(C) T_DC(C)', &
1487  ' T_DF(C) T_Ko(C) T_By(C)',/, &
1488  ' Rb_Vo(W/m2) Rb_DA(W/m2) Rb_DC(W/m2)', &
1489  ' Rb_DF(W/m2) Rb_Ko(W/m2) Rb_By(W/m2)')
1490  1127 format('----------------------------------------------------', &
1491  '---------------------------------------')
1492 
1493 end if
1494 
1495 #endif
1496 
1497 !-------- Output of the initial state --------
1498 
1499 #if OUTPUT==1
1500 
1501 #if ERGDAT==0
1502  flag_3d_output = .false.
1503 #elif ERGDAT==1
1504  flag_3d_output = .true.
1505 #else
1506  stop ' sico_init: ERGDAT must be either 0 or 1!'
1507 #endif
1508 
1509 #if NETCDF==1
1510  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1511  flag_3d_output, ndat2d, ndat3d)
1512 #elif NETCDF==2
1513  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1514  flag_3d_output, ndat2d, ndat3d)
1515 #else
1516  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1517 #endif
1518 
1519 #elif OUTPUT==2
1520 
1521 if (time_output(1) <= time_init+eps) then
1522 
1523 #if ERGDAT==0
1524  flag_3d_output = .false.
1525 #elif ERGDAT==1
1526  flag_3d_output = .true.
1527 #else
1528  stop ' sico_init: ERGDAT must be either 0 or 1!'
1529 #endif
1530 
1531 #if NETCDF==1
1532  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1533  flag_3d_output, ndat2d, ndat3d)
1534 #elif NETCDF==2
1535  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1536  flag_3d_output, ndat2d, ndat3d)
1537 #else
1538  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1539 #endif
1540 
1541 end if
1542 
1543 #elif OUTPUT==3
1544 
1545  flag_3d_output = .false.
1546 
1547 #if NETCDF==1
1548  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1549  flag_3d_output, ndat2d, ndat3d)
1550 #elif NETCDF==2
1551  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1552  flag_3d_output, ndat2d, ndat3d)
1553 #else
1554  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1555 #endif
1556 
1557 if (time_output(1) <= time_init+eps) then
1558 
1559  flag_3d_output = .true.
1560 
1561 #if NETCDF==1
1562  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1563  flag_3d_output, ndat2d, ndat3d)
1564 #elif NETCDF==2
1565  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1566  flag_3d_output, ndat2d, ndat3d)
1567 #else
1568  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1569 #endif
1570 
1571 end if
1572 
1573 #else
1574  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1575 #endif
1576 
1577 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1578 #if OUTSER==2
1579 call output3(time_init, delta_ts, glac_index, z_sl)
1580 #endif
1581 #if OUTSER==3
1582 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1583 #endif
1584 
1585 end subroutine sico_init
1586 !
Generated on Mon May 20 2013 10:18:36 for SICOPOLIS V3.0 by   doxygen 1.8.1