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