SICOPOLIS V3.3
sico_init_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ i n i t _ m
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 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 module sico_init_m
36 
37  use sico_types_m
39  use sico_vars_m
40 
41  implicit none
42 
43  public
44 
45 contains
46 
47 !-------------------------------------------------------------------------------
48 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
49 !<------------------------------------------------------------------------------
50 subroutine sico_init(delta_ts, glac_index, &
51  mean_accum, &
52  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
53  time, time_init, time_end, time_output, &
54  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
55  z_sl, dzsl_dtau, z_mar, &
56  ii, jj, nn, &
57  ndat2d, ndat3d, n_output, &
58  runname)
59 
60  use compare_float_m
64 
65  use read_m, only : read_phys_para, read_kei
66 
67  use boundary_m
69  use calc_enhance_m
70  use calc_vxy_m
71  use calc_vz_m
72  use calc_dxyz_m
74 
75  use output_m
76 
77 implicit none
78 
79 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
80  jj((imax+1)*(jmax+1)), &
81  nn(0:jmax,0:imax)
82 integer(i4b), intent(out) :: ndat2d, ndat3d
83 integer(i4b), intent(out) :: n_output
84 real(dp), intent(out) :: delta_ts, glac_index
85 real(dp), intent(out) :: mean_accum
86 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
87  dtime_out, dtime_ser
88 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
89 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
90 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
91 character(len=100), intent(out) :: runname
92 
93 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
94 integer(i4b) :: ios
95 integer(i4b) :: ierr
96 integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_ref
97 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
98 real(dp) :: time_init0, time_end0, time_output0(100)
99 real(dp) :: d_dummy
100 character(len=100) :: anfdatname
101 character(len=256) :: filename_with_path
102 character(len=256) :: shell_command
103 character :: ch_dummy
104 logical :: flag_init_output, flag_3d_output
105 
106 character(len=64), parameter :: fmt1 = '(a)', &
107  fmt2 = '(a,i4)', &
108  fmt2a = '(a,i0)', &
109  fmt3 = '(a,es12.4)'
110 
111 character(len= 8) :: ch_imax
112 character(len=128) :: fmt4
113 
114 write(ch_imax, fmt='(i8)') imax
115 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
116 
117 write(unit=6, fmt='(a)') ' '
118 write(unit=6, fmt='(a)') ' -------- sico_init --------'
119 write(unit=6, fmt='(a)') ' '
120 
121 !-------- Name of the computational domain --------
122 
123 #if (defined(ANT))
124 ch_domain_long = 'Antarctica'
125 ch_domain_short = 'ant'
126 
127 #elif (defined(ASF))
128 ch_domain_long = 'Austfonna'
129 ch_domain_short = 'asf'
130 
131 #elif (defined(EMTP2SGE))
132 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
133 ch_domain_short = 'emtp2sge'
134 
135 #elif (defined(GRL))
136 ch_domain_long = 'Greenland'
137 ch_domain_short = 'grl'
138 
139 #elif (defined(NHEM))
140 ch_domain_long = 'Northern hemisphere'
141 ch_domain_short = 'nhem'
142 
143 #elif (defined(SCAND))
144 ch_domain_long = 'Scandinavia and Eurasia'
145 ch_domain_short = 'scand'
146 
147 #elif (defined(TIBET))
148 ch_domain_long = 'Tibet'
149 ch_domain_short = 'tibet'
150 
151 #elif (defined(NMARS))
152 ch_domain_long = 'North polar cap of Mars'
153 ch_domain_short = 'nmars'
154 
155 #elif (defined(SMARS))
156 ch_domain_long = 'South polar cap of Mars'
157 ch_domain_short = 'smars'
158 
159 #elif (defined(XYZ))
160 ch_domain_long = 'XYZ'
161 ch_domain_short = 'xyz'
162 #if (defined(HEINO))
163 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
164 #endif
165 
166 #else
167 stop ' >>> sico_init: No valid domain specified!'
168 #endif
169 
170 !-------- Some initial values --------
171 
172 n_output = 0
173 
174 dtime = 0.0_dp
175 dtime_temp = 0.0_dp
176 dtime_wss = 0.0_dp
177 dtime_out = 0.0_dp
178 dtime_ser = 0.0_dp
179 
180 time = 0.0_dp
181 time_init = 0.0_dp
182 time_end = 0.0_dp
183 time_output = 0.0_dp
184 
185 !-------- Initialisation of the Library of Iterative Solvers Lis,
186 ! if required --------
187 
188 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
189  call lis_initialize(ierr)
190 #endif
191 
192 !-------- Read physical parameters --------
193 
194 call read_phys_para()
195 
196 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
197 
198 ! ------ Some auxiliary quantities required for the enthalpy method
199 
200 call calc_c_int_table(c, -190, 10, l)
202 
203 !-------- Compatibility check of the SICOPOLIS version with the header file
204 
205 if ( trim(version) /= trim(sico_version) ) &
206  stop ' >>> sico_init: SICOPOLIS version not compatible with header file!'
207 
208 !-------- Check whether the dynamics and thermodynamics modes are defined
209 
210 #if (!defined(DYNAMICS))
211 stop ' >>> sico_init: DYNAMICS not defined in the header file!'
212 #endif
213 
214 #if (!defined(CALCMOD))
215 stop ' >>> sico_init: CALCMOD not defined in the header file!'
216 #endif
217 
218 #if (defined(ENTHMOD))
219 write(6, fmt='(a)') ' >>> sico_init: ENTHMOD must not be defined any more.'
220 write(6, fmt='(a)') ' Please update your header file!'
221 stop
222 #endif
223 
224 !-------- Compatibility check of the horizontal resolution with the
225 ! number of grid points --------
226 
227 #if (GRID==0 || GRID==1)
228 
229  stop ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
230 
231 #elif (GRID==2)
232 
233 
234 if ( (approx_equal(dlambda, 1.0_dp/3.0_dp, eps_sp_dp)) &
235  .and.(approx_equal(dphi , 1.0_dp/3.0_dp, eps_sp_dp)) ) then
236 
237  if ((imax /= 135).or.(jmax /= 51)) then
238  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
239  end if
240 
241 else if ( (approx_equal(dlambda, 1.0_dp/6.0_dp, eps_sp_dp)) &
242  .and.(approx_equal(dphi , 1.0_dp/6.0_dp, eps_sp_dp)) ) then
243 
244  if ((imax /= 270).or.(jmax /= 102)) then
245  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
246  end if
247 
248 else if ( (approx_equal(dlambda, 1.0_dp/12.0_dp, eps_sp_dp)) &
249  .and.(approx_equal(dphi , 1.0_dp/12.0_dp, eps_sp_dp)) ) then
250 
251  if ((imax /= 540).or.(jmax /= 204)) then
252  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
253  end if
254 
255 else
256  stop ' >>> sico_init: DLAMBDA / DPHI wrong!'
257 end if
258 
259 #endif
260 
261 !-------- Compatibility check of the thermodynamics mode
262 ! (cold vs. polythermal vs. enthalpy method)
263 ! and the number of grid points in the lower (kt) ice domain --------
264 
265 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
266 
267 if (ktmax > 2) then
268  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
269  write(6, fmt='(a)') ' the separate kt domain is redundant.'
270  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
271  write(6, fmt='(a)') ' '
272 end if
273 
274 #endif
275 
276 !-------- Compatibility check of surface-temperature and precipitation
277 ! determination by interpolation between present and LGM values
278 ! with a glacial index --------
279 
280 #if (TSURFACE == 5 && ACCSURFACE != 5)
281 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
282 #endif
283 
284 #if (TSURFACE != 5 && ACCSURFACE == 5)
285 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
286 #endif
287 
288 !-------- Compatibility check of discretization schemes for the horizontal and
289 ! vertical advection terms in the temperature and age equations --------
290 
291 #if (ADV_HOR==1)
292 stop ' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!'
293 #endif
294 
295 !-------- Check whether for the shallow shelf
296 ! or shelfy stream approximation
297 ! the chosen grid is Cartesian coordinates
298 ! without distortion correction (GRID==0) --------
299 
300 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
301 #if (GRID != 0)
302 write(6, fmt='(a)') ' >>> sico_init: Distortion correction not implemented'
303 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
304 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
305 write(6, fmt='(a)') ' -> GRID==0 required!'
306 stop
307 #endif
308 #endif
309 
310 !-------- Setting of forcing flag --------
311 
312 #if (TSURFACE <= 4)
313 
314 forcing_flag = 1 ! forcing by delta_ts
315 
316 #elif (TSURFACE == 5)
317 
318 forcing_flag = 2 ! forcing by glac_index
319 
320 #endif
321 
322 !-------- Initialization of numerical time steps --------
323 
324 dtime0 = dtime0
325 dtime_temp0 = dtime_temp0
326 #if (REBOUND==2)
327 dtime_wss0 = dtime_wss0
328 #endif
329 
330 !-------- Further initializations --------
331 
332 dzeta_c = 1.0_dp/real(kcmax,dp)
333 dzeta_t = 1.0_dp/real(ktmax,dp)
334 dzeta_r = 1.0_dp/real(krmax,dp)
335 
336 ndat2d = 1
337 ndat3d = 1
338 
339 !-------- General abbreviations --------
340 
341 ! ------ kc domain
342 
343 if (deform >= eps) then
344 
345  flag_aa_nonzero = .true. ! non-equidistant grid
346 
347  aa = deform
348  ea = exp(aa)
349 
350  kc=0
351  zeta_c(kc) = 0.0_dp
352  eaz_c(kc) = 1.0_dp
353  eaz_c_quotient(kc) = 0.0_dp
354 
355  do kc=1, kcmax-1
356  zeta_c(kc) = kc*dzeta_c
357  eaz_c(kc) = exp(aa*zeta_c(kc))
358  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
359  end do
360 
361  kc=kcmax
362  zeta_c(kc) = 1.0_dp
363  eaz_c(kc) = exp(aa)
364  eaz_c_quotient(kc) = 1.0_dp
365 
366 else
367 
368  flag_aa_nonzero = .false. ! equidistant grid
369 
370  aa = 0.0_dp
371  ea = 1.0_dp
372 
373  kc=0
374  zeta_c(kc) = 0.0_dp
375  eaz_c(kc) = 1.0_dp
376  eaz_c_quotient(kc) = 0.0_dp
377 
378  do kc=1, kcmax-1
379  zeta_c(kc) = kc*dzeta_c
380  eaz_c(kc) = 1.0_dp
381  eaz_c_quotient(kc) = zeta_c(kc)
382  end do
383 
384  kc=kcmax
385  zeta_c(kc) = 1.0_dp
386  eaz_c(kc) = 1.0_dp
387  eaz_c_quotient(kc) = 1.0_dp
388 
389 end if
390 
391 ! ------ kt domain
392 
393 kt=0
394 zeta_t(kt) = 0.0_dp
395 
396 do kt=1, ktmax-1
397  zeta_t(kt) = kt*dzeta_t
398 end do
399 
400 kt=ktmax
401 zeta_t(kt) = 1.0_dp
402 
403 ! ------ kr domain
404 
405 kr=0
406 zeta_r(kr) = 0.0_dp
407 
408 do kr=1, krmax-1
409  zeta_r(kr) = kr*dzeta_r
410 end do
411 
412 kr=krmax
413 zeta_r(kr) = 1.0_dp
414 
415 !-------- Construction of a vector (with index n) from a 2-d array
416 ! (with indices i, j) by diagonal numbering --------
417 
418 n=1
419 
420 do m=0, imax+jmax
421  do i=m, 0, -1
422  j = m-i
423  if ((i <= imax).and.(j <= jmax)) then
424  ii(n) = i
425  jj(n) = j
426  nn(j,i) = n
427  n=n+1
428  end if
429  end do
430 end do
431 
432 !-------- Specification of current simulation --------
433 
434 runname = runname
435 anfdatname = anfdatname
436 
437 #if (defined(YEAR_ZERO))
439 #else
440 year_zero = 2000.0_dp ! default value 2000 CE
441 #endif
442 
443 time_init0 = time_init0
444 time_end0 = time_end0
445 dtime_ser0 = dtime_ser0
446 
447 #if (OUTPUT==1 || OUTPUT==3)
448 dtime_out0 = dtime_out0
449 #endif
450 
451 #if (OUTPUT==2 || OUTPUT==3)
452 n_output = n_output
453 time_output0( 1) = time_out0_01
454 time_output0( 2) = time_out0_02
455 time_output0( 3) = time_out0_03
456 time_output0( 4) = time_out0_04
457 time_output0( 5) = time_out0_05
458 time_output0( 6) = time_out0_06
459 time_output0( 7) = time_out0_07
460 time_output0( 8) = time_out0_08
461 time_output0( 9) = time_out0_09
462 time_output0(10) = time_out0_10
463 time_output0(11) = time_out0_11
464 time_output0(12) = time_out0_12
465 time_output0(13) = time_out0_13
466 time_output0(14) = time_out0_14
467 time_output0(15) = time_out0_15
468 time_output0(16) = time_out0_16
469 time_output0(17) = time_out0_17
470 time_output0(18) = time_out0_18
471 time_output0(19) = time_out0_19
472 time_output0(20) = time_out0_20
473 #endif
474 
475 !-------- Write log file --------
476 
477 shell_command = 'if [ ! -d'
478 shell_command = trim(shell_command)//' '//outpath
479 shell_command = trim(shell_command)//' '//'] ; then mkdir'
480 shell_command = trim(shell_command)//' '//outpath
481 shell_command = trim(shell_command)//' '//'; fi'
482 call system(trim(shell_command))
483  ! Check whether directory OUTPATH exists. If not, it is created.
484 
485 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
486 
487 open(10, iostat=ios, file=trim(filename_with_path), status='new')
488 
489 if (ios /= 0) stop ' >>> sico_init: Error when opening the log file!'
490 
491 write(10, fmt=trim(fmt1)) 'Computational domain:'
492 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
493 write(10, fmt=trim(fmt1)) ' '
494 
495 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
496 write(10, fmt=trim(fmt1)) ' '
497 
498 write(10, fmt=trim(fmt2)) 'imax =', imax
499 write(10, fmt=trim(fmt2)) 'jmax =', jmax
500 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
501 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
502 write(10, fmt=trim(fmt2)) 'krmax =', krmax
503 write(10, fmt=trim(fmt1)) ' '
504 
505 write(10, fmt=trim(fmt3)) 'a =', aa
506 write(10, fmt=trim(fmt1)) ' '
507 
508 #if (GRID==0 || GRID==1)
509 stop ' >>> sico_init: GRID==0 or 1 not allowed for this application!'
510 #elif (GRID==2)
511 write(10, fmt=trim(fmt3)) 'lambda0 =', lambda_0
512 write(10, fmt=trim(fmt3)) 'phi0 =', phi_0
513 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
514 write(10, fmt=trim(fmt3)) 'dphi =', dphi
515 #endif
516 write(10, fmt=trim(fmt1)) ' '
517 
518 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
519 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
520 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
521 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
522 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
523 #if (REBOUND==2)
524 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
525 #endif
526 write(10, fmt=trim(fmt1)) ' '
527 
528 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
529 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
530 #if (ANF_DAT==1)
531 #if (defined(ZB_PRESENT_FILE))
532 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
533 #endif
534 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
535 #endif
536 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
537 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
538 #if (ANF_DAT==1 && defined(TEMP_INIT))
539 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
540 #endif
541 #if (ANF_DAT==3)
542 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
543 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
544 #endif
545 write(10, fmt=trim(fmt1)) ' '
546 
547 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
548 write(10, fmt=trim(fmt1)) ' '
549 
550 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
551 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
552 #if (CALCTHK==2 || CALCTHK==5)
553 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
554 #if (ITER_MAX_SOR>0)
555 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
556 #endif
557 #endif
558 write(10, fmt=trim(fmt1)) ' '
559 #endif
560 
561 write(10, fmt=trim(fmt1)) 'temp_mm_present file = '//temp_mm_present_file
562 #if (TSURFACE==1)
563 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
564 #elif (TSURFACE==3)
565 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
566 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
567 #elif (TSURFACE==4)
568 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
569 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
570 #elif (TSURFACE==5)
571 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
572 write(10, fmt=trim(fmt1)) 'temp_mm_anom file = '//temp_mm_anom_file
573 write(10, fmt=trim(fmt3)) 'temp_mm_anom fact = ', temp_mm_anom_fact
574 #endif
575 
576 write(10, fmt=trim(fmt1)) 'precip_mm_present file = '//precip_mm_present_file
577 #if (ACCSURFACE==1)
578 write(10, fmt=trim(fmt3)) 'accfact =', accfact
579 #elif (ACCSURFACE==2 || ACCSURFACE==3)
580 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
581 #elif (ACCSURFACE==5)
582 write(10, fmt=trim(fmt1)) 'precip_mm_anom file = '//precip_mm_anom_file
583 write(10, fmt=trim(fmt3)) 'precip_mm_anom fact = ', precip_mm_anom_fact
584 #endif
585 #if (ACCSURFACE <= 3)
586 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
587 #if (ELEV_DESERT == 1)
588 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
589 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
590 #endif
591 #endif
592 
593 #if (ABLSURFACE==3)
594 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
595 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
596 #endif
597 
598 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
599 #if (SEA_LEVEL==1)
600 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
601 #elif (SEA_LEVEL==3)
602 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
603 #endif
604 write(10, fmt=trim(fmt1)) ' '
605 
606 #if (MARGIN==2)
607 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
608 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
609 write(10, fmt=trim(fmt1)) ' '
610 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
611  || marine_ice_calving==6 || marine_ice_calving==7)
612 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
613 write(10, fmt=trim(fmt1)) ' '
614 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
615 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
616 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
617 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
618 write(10, fmt=trim(fmt1)) ' '
619 #endif
620 #elif (MARGIN==3)
621 #if (ICE_SHELF_CALVING==2)
622 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
623 write(10, fmt=trim(fmt1)) ' '
624 #endif
625 #endif
626 
627 #if (defined(BASAL_HYDROLOGY))
628 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
629 #endif
630 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
631 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
632 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
633 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
634 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
635 #if (defined(TIME_RAMP_UP_SLIDE))
636 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
637 #endif
638 #if (SLIDE_LAW==2)
639 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
640 #endif
641 #if (BASAL_HYDROLOGY==1)
642 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
643 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
644 #endif
645 write(10, fmt=trim(fmt1)) ' '
646 
647 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
648 #if (Q_GEO_MOD==1)
649 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
650 #elif (Q_GEO_MOD==2)
651 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
652 #endif
653 write(10, fmt=trim(fmt1)) ' '
654 
655 #if (defined(MARINE_ICE_BASAL_MELTING))
656 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
657 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
658 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
659 #endif
660 write(10, fmt=trim(fmt1)) ' '
661 #endif
662 
663 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
664 #if (REBOUND==1)
665 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
666 #endif
667 #if (REBOUND==1 || REBOUND==2)
668 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
669 #if (TIME_LAG_MOD==1)
670 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
671 #elif (TIME_LAG_MOD==2)
672 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
673 #else
674 stop ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
675 #endif
676 #endif
677 #if (REBOUND==2)
678 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
679 #if (FLEX_RIG_MOD==1)
680 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
681 #elif (FLEX_RIG_MOD==2)
682 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
683 #else
684 stop ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
685 #endif
686 #endif
687 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
688 write(10, fmt=trim(fmt1)) ' '
689 
690 #if (FLOW_LAW==2)
691 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
692 write(10, fmt=trim(fmt1)) ' '
693 #endif
694 #if (FIN_VISC==2)
695 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
696 write(10, fmt=trim(fmt1)) ' '
697 #endif
698 
699 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
700 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
701 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
702 #endif
703 #if (ENHMOD==2 || ENHMOD==3)
704 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
705 #endif
706 #if (ENHMOD==2)
707 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
708 #endif
709 #if (ENHMOD==3)
710 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
711 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
712 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
713 #endif
714 #if (ENHMOD==4 || ENHMOD==5)
715 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
716 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
717 #endif
718 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
719 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
720 #endif
721 write(10, fmt=trim(fmt1)) ' '
722 
723 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
724 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
725 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
726 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
727 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
728 #if (defined(QBM_MIN))
729 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
730 #elif (defined(QB_MIN)) /* obsolete */
731 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
732 #endif
733 #if (defined(QBM_MAX))
734 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
735 #elif (defined(QB_MAX)) /* obsolete */
736 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
737 #endif
738 write(10, fmt=trim(fmt3)) 'age_min =', age_min
739 write(10, fmt=trim(fmt3)) 'age_max =', age_max
740 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
741 #if (ADV_VERT==1)
742 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
743 #endif
744 write(10, fmt=trim(fmt1)) ' '
745 
746 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
747 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
748 #if (defined(LIS_OPTS))
749 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
750 #endif
751 #if (defined(N_ITER_SSA))
752 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
753 #endif
754 #if (defined(RELAX_FACT_SSA))
755 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
756 #endif
757 #endif
758 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
759 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
760 #endif
761 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
762 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
763 #endif
764 write(10, fmt=trim(fmt1)) ' '
765 
766 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
767 #if (CALCMOD==-1 && defined(TEMP_CONST))
768 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
769 #endif
770 #if (CALCMOD==-1 && defined(AGE_CONST))
771 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
772 #endif
773 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
774 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
775 #endif
776 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
777 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
778 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
779 #if (MARGIN==2)
780 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
781 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
782 #elif (MARGIN==3)
783 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
784 #endif
785 #if (defined(THK_EVOL))
786 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
787 #else
788 stop ' >>> sico_init: Define THK_EVOL in header file!'
789 #endif
790 #if (defined(CALCTHK))
791 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
792 #else
793 stop ' >>> sico_init: Define CALCTHK in header file!'
794 #endif
795 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
796 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
797 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
798 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
799 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
800 #if (ACCSURFACE==5)
801 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
802 #endif
803 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
804 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
805 #if (defined(MB_ACCOUNT))
806 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
807 #endif
808 write(10, fmt=trim(fmt1)) ' '
809 
810 #if (defined(OUT_TIMES))
811 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
812 #endif
813 #if (defined(OUTPUT_INIT))
814 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
815 #endif
816 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
817 #if (OUTPUT==1 || OUTPUT==3)
818 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
819 #endif
820 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
821 #if (OUTPUT==1 || OUTPUT==2)
822 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
823 #endif
824 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
825 #if (OUTPUT==2 || OUTPUT==3)
826 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
827 do n=1, n_output
828  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
829 end do
830 #endif
831 write(10, fmt=trim(fmt1)) ' '
832 
833 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
834 
835 close(10, status='keep')
836 
837 !-------- Conversion of time quantities --------
838 
839 year_zero = year_zero*year_sec ! a --> s
840 time_init = time_init0*year_sec ! a --> s
841 time_end = time_end0*year_sec ! a --> s
842 dtime = dtime0*year_sec ! a --> s
843 dtime_temp = dtime_temp0*year_sec ! a --> s
844 #if (REBOUND==2)
845 dtime_wss = dtime_wss0*year_sec ! a --> s
846 #endif
847 dtime_ser = dtime_ser0*year_sec ! a --> s
848 #if (OUTPUT==1 || OUTPUT==3)
849 dtime_out = dtime_out0*year_sec ! a --> s
850 #endif
851 #if (OUTPUT==2 || OUTPUT==3)
852 do n=1, n_output
853  time_output(n) = time_output0(n)*year_sec ! a --> s
854 end do
855 #endif
856 
857 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) &
858  stop ' >>> sico_init: dtime_temp must be a multiple of dtime!'
859 
860 #if (REBOUND==2)
861 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) &
862  stop ' >>> sico_init: dtime_wss must be a multiple of dtime!'
863 #endif
864 
865 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) &
866  stop ' >>> sico_init: dtime_ser must be a multiple of dtime!'
867 
868 #if (OUTPUT==1 || OUTPUT==3)
869 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) &
870  stop ' >>> sico_init: dtime_out must be a multiple of dtime!'
871 #endif
872 
873 time = time_init
874 
875 !-------- Reading of measurements for present monthly-mean precipitation --------
876 
877 #if (GRID==0 || GRID==1)
878 
879 stop ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
880 
881 #elif (GRID==2)
882 
883 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
884  trim(precip_mm_present_file)
885 
886 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
887 
888 #endif
889 
890 if (ios /= 0) stop ' >>> sico_init: Error when opening the precipitation file!'
891 
892 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
893 
894 do n=1, 12 ! month counter
895  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
896  do j=jmax, 0, -1
897  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
898  end do
899 end do
900 
901 close(21, status='keep')
902 
903 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
904 
905 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
906  ! mm/a water equiv. -> m/s ice equiv.
907 
908 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
909 
910 #if (ACCSURFACE==5)
911 
912 #if (GRID==2)
913 
914 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
915  trim(precip_mm_anom_file)
916 
917 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
918 
919 #endif
920 
921 if (ios /= 0) stop ' >>> sico_init: Error when opening the precip anomaly file!'
922 
923 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
924 
925 do n=1, 12 ! month counter
926  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
927  do j=jmax, 0, -1
928  read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
929  end do
930 end do
931 
932 close(21, status='keep')
933 
934 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
935 
936 do i=0, imax
937 do j=0, jmax
938 
939 #if (PRECIP_ANOM_INTERPOL==1)
940  do n=1, 12 ! month counter
941  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
942  end do
943 #elif (PRECIP_ANOM_INTERPOL==2)
944  do n=1, 12 ! month counter
945  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
946  end do
947 #else
948  stop ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
949 #endif
950 
951 end do
952 end do
953 
954 #endif
955 
956 !-------- Mean accumulation --------
957 
958 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
959  ! mm/a water equiv. -> m/s ice equiv.
960 
961 !-------- Reading of present topography mask --------
962 
963 #if (GRID==0 || GRID==1)
964 
965 stop ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
966 
967 #elif (GRID==2)
968 
969 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
970  trim(mask_present_file)
971 
972 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
973 
974 #endif
975 
976 if (ios /= 0) stop ' >>> sico_init: Error when opening the mask file!'
977 
978 do m=1, 6; read(24, fmt='(a)') ch_dummy; end do
979 
980 do j=jmax, 0, -1
981  read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
982 end do
983 
984 close(24, status='keep')
985 
986 !-------- Reading of data for present monthly-mean surface temperature --------
987 
988 #if (GRID==2)
989 
990 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
991  trim(temp_mm_present_file)
992 
993 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
994 
995 #endif
996 
997 if (ios /= 0) stop ' >>> sico_init: Error when opening the temperature file!'
998 
999 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1000 
1001 do n=1, 12 ! month counter
1002  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1003  do j=jmax, 0, -1
1004  read(21, fmt=*) (temp_mm_present(j,i,n), i=0,imax)
1005  end do
1006 end do
1007 
1008 close(21, status='keep')
1009 
1010 !-------- Reading of LGM monthly-mean surface-temperature anomalies --------
1011 
1012 #if (TSURFACE==5)
1013 
1014 #if (GRID==2)
1015 
1016 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1017  trim(temp_mm_anom_file)
1018 
1019 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1020 
1021 #endif
1022 
1023 if (ios /= 0) stop ' >>> sico_init: Error when opening the temperature anomaly file!'
1024 
1025 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1026 
1027 do n=1, 12 ! month counter
1028  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1029  do j=jmax, 0, -1
1030  read(21, fmt=*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
1031  end do
1032 end do
1033 
1034 close(21, status='keep')
1035 
1036 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
1037 
1038 #endif
1039 
1040 !-------- Present reference elevation
1041 ! (for precipitation and surface-temperature data) --------
1042 
1043 #if (GRID==0 || GRID==1)
1044 
1045 stop ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
1046 
1047 #elif (GRID==2)
1048 
1049 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1050  trim(zs_present_file)
1051 
1052 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1053 
1054 #endif
1055 
1056 if (ios /= 0) stop ' >>> sico_init: Error when opening the zs file!'
1057 
1058 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1059 
1060 do j=jmax, 0, -1
1061  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1062 end do
1063 
1064 close(21, status='keep')
1065 
1066 ! ------ Reset bathymetry data (sea floor elevation) to the
1067 ! sea surface
1068 
1069 do i=0, imax
1070 do j=0, jmax
1071  if (maske_ref(j,i) >= 2) zs_ref(j,i) = 0.0_dp
1072 end do
1073 end do
1074 
1075 !-------- Read data for delta_ts --------
1076 
1077 #if (TSURFACE==4)
1078 
1079 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1080 
1081 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1082 
1083 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for delta_ts!'
1084 
1085 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1086 
1087 if (ch_dummy /= '#') then
1088  write(6, fmt=*) ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max'
1089  write(6, fmt=*) ' not defined in data file for delta_ts!'
1090  stop
1091 end if
1092 
1094 
1095 allocate(griptemp(0:ndata_grip))
1096 
1097 do n=0, ndata_grip
1098  read(21, fmt=*) d_dummy, griptemp(n)
1099 end do
1100 
1101 close(21, status='keep')
1102 
1103 #endif
1104 
1105 !-------- Read data for the glacial index --------
1106 
1107 #if (TSURFACE==5)
1108 
1109 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1110 
1111 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1112 
1113 if (ios /= 0) stop ' >>> sico_init: Error when opening the glacial-index file!'
1114 
1115 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1116 
1117 if (ch_dummy /= '#') then
1118  write(6, fmt=*) ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max'
1119  write(6, fmt=*) ' not defined in glacial-index file!'
1120  stop
1121 end if
1122 
1124 
1125 allocate(glacial_index(0:ndata_gi))
1126 
1127 do n=0, ndata_gi
1128  read(21, fmt=*) d_dummy, glacial_index(n)
1129 end do
1130 
1131 close(21, status='keep')
1132 
1133 #endif
1134 
1135 !-------- Read data for z_sl --------
1136 
1137 #if (SEA_LEVEL==3)
1138 
1139 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1140 
1141 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1142 
1143 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for z_sl!'
1144 
1145 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1146 
1147 if (ch_dummy /= '#') then
1148  write(6, fmt=*) ' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1149  write(6, fmt=*) ' not defined in data file for z_sl!'
1150  stop
1151 end if
1152 
1154 
1155 allocate(specmap_zsl(0:ndata_specmap))
1156 
1157 do n=0, ndata_specmap
1158  read(21, fmt=*) d_dummy, specmap_zsl(n)
1159 end do
1160 
1161 close(21, status='keep')
1162 
1163 #endif
1164 
1165 !-------- Determination of the geothermal heat flux --------
1166 
1167 #if (Q_GEO_MOD==1)
1168 
1169 ! ------ Constant value
1170 
1171 do i=0, imax
1172 do j=0, jmax
1173  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1174 end do
1175 end do
1176 
1177 #elif (Q_GEO_MOD==2)
1178 
1179 ! ------ Read data from file
1180 
1181 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1182  trim(q_geo_file)
1183 
1184 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1185 
1186 if (ios /= 0) stop ' >>> sico_init: Error when opening the qgeo file!'
1187 
1188 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1189 
1190 do j=jmax, 0, -1
1191  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1192 end do
1193 
1194 close(21, status='keep')
1195 
1196 do i=0, imax
1197 do j=0, jmax
1198  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1199 end do
1200 end do
1201 
1202 #endif
1203 
1204 !-------- Reading of tabulated kei function--------
1205 
1206 #if (REBOUND==0 || REBOUND==1)
1207 
1208 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1209  ! dummy values
1210 #elif (REBOUND==2)
1211 
1212 call read_kei()
1213 
1214 #endif
1215 
1216 !-------- Determination of the time lag
1217 ! of the relaxing asthenosphere --------
1218 
1219 #if (REBOUND==1 || REBOUND==2)
1220 
1221 #if (TIME_LAG_MOD==1)
1222 
1223 time_lag_asth = time_lag*year_sec ! a -> s
1224 
1225 #elif (TIME_LAG_MOD==2)
1226 
1227 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1228  trim(time_lag_file)
1229 
1230 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1231 
1232 if (ios /= 0) stop ' >>> sico_init: Error when opening the time-lag file!'
1233 
1234 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1235 
1236 do j=jmax, 0, -1
1237  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1238 end do
1239 
1240 close(29, status='keep')
1241 
1242 time_lag_asth = time_lag_asth*year_sec ! a -> s
1243 
1244 #endif
1245 
1246 #elif (REBOUND==0)
1247 
1248 time_lag_asth = 0.0_dp ! dummy values
1249 
1250 #endif
1251 
1252 !-------- Determination of the flexural rigidity of the lithosphere --------
1253 
1254 #if (REBOUND==2)
1255 
1256 #if (FLEX_RIG_MOD==1)
1257 
1258 flex_rig_lith = flex_rig
1259 
1260 #elif (FLEX_RIG_MOD==2)
1261 
1262 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1263  trim(flex_rig_file)
1264 
1265 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1266 
1267 if (ios /= 0) stop ' >>> sico_init: Error when opening the flex-rig file!'
1268 
1269 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1270 
1271 do j=jmax, 0, -1
1272  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1273 end do
1274 
1275 close(29, status='keep')
1276 
1277 #endif
1278 
1279 #elif (REBOUND==0 || REBOUND==1)
1280 
1281 flex_rig_lith = 0.0_dp ! dummy values
1282 
1283 #endif
1284 
1285 !-------- Definition of initial values --------
1286 
1287 ! ------ Present topography
1288 
1289 #if (ANF_DAT==1)
1290 
1291 stop ' >>> sico_init: ANF_DAT==1 not allowed for Tibet application!'
1292 
1293 ! ------ Ice-free, relaxed bedrock
1294 
1295 #elif (ANF_DAT==2)
1296 
1297 call topography2(dxi, deta)
1298 
1299 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1300 
1301 call boundary(time_init, dtime, dxi, deta, &
1302  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1303 
1304 as_perp_apl = 0.0_dp
1305 
1306 q_bm = 0.0_dp
1307 q_tld = 0.0_dp
1308 q_b_tot = 0.0_dp
1309 
1310 p_b_w = 0.0_dp
1311 h_w = 0.0_dp
1312 
1313 call init_temp_water_age_2()
1314 
1315 #if (ENHMOD==1)
1316  call calc_enhance_1()
1317 #elif (ENHMOD==2)
1318  call calc_enhance_2()
1319 #elif (ENHMOD==3)
1320  call calc_enhance_3(time_init)
1321 #elif (ENHMOD==4)
1322  call calc_enhance_4()
1323 #elif (ENHMOD==5)
1324  call calc_enhance_5()
1325 #else
1326  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1327 #endif
1328 
1329 ! ------ Read initial state from output of previous simulation
1330 
1331 #elif (ANF_DAT==3)
1332 
1333 call topography3(dxi, deta, z_sl, anfdatname)
1334 
1335 call boundary(time_init, dtime, dxi, deta, &
1336  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1337 
1338 where ((maske==0_i2b).or.(maske==3_i2b))
1339  ! grounded or floating ice
1341 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1342  as_perp_apl = 0.0_dp
1343 end where
1344 
1345 q_b_tot = q_bm + q_tld
1346 
1347 #endif
1348 
1349 !-------- Inner-point flag --------
1350 
1351 flag_inner_point = .true.
1352 
1353 flag_inner_point(0,:) = .false.
1354 flag_inner_point(jmax,:) = .false.
1355 
1356 flag_inner_point(:,0) = .false.
1357 flag_inner_point(:,imax) = .false.
1358 
1359 !-------- Grounding line and calving front flags --------
1360 
1361 flag_grounding_line_1 = .false.
1362 flag_grounding_line_2 = .false.
1363 
1364 flag_calving_front_1 = .false.
1365 flag_calving_front_2 = .false.
1366 
1367 #if (MARGIN==1 || MARGIN==2)
1368 
1369 continue
1370 
1371 #elif (MARGIN==3)
1372 
1373 do i=1, imax-1
1374 do j=1, jmax-1
1375 
1376  if ( (maske(j,i)==0_i2b) & ! grounded ice
1377  .and. &
1378  ( (maske(j,i+1)==3_i2b) & ! with
1379  .or.(maske(j,i-1)==3_i2b) & ! one
1380  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1381  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1382  ) &
1383  flag_grounding_line_1(j,i) = .true.
1384 
1385  if ( (maske(j,i)==3_i2b) & ! floating ice
1386  .and. &
1387  ( (maske(j,i+1)==0_i2b) & ! with
1388  .or.(maske(j,i-1)==0_i2b) & ! one
1389  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1390  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1391  ) &
1392  flag_grounding_line_2(j,i) = .true.
1393 
1394  if ( (maske(j,i)==3_i2b) & ! floating ice
1395  .and. &
1396  ( (maske(j,i+1)==2_i2b) & ! with
1397  .or.(maske(j,i-1)==2_i2b) & ! one
1398  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1399  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1400  ) &
1401  flag_calving_front_1(j,i) = .true.
1402 
1403  if ( (maske(j,i)==2_i2b) & ! ocean
1404  .and. &
1405  ( (maske(j,i+1)==3_i2b) & ! with
1406  .or.(maske(j,i-1)==3_i2b) & ! one
1407  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1408  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1409  ) &
1410  flag_calving_front_2(j,i) = .true.
1411 
1412 end do
1413 end do
1414 
1415 do i=1, imax-1
1416 
1417  j=0
1418  if ( (maske(j,i)==2_i2b) & ! ocean
1419  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1420  ) & ! floating ice point
1421  flag_calving_front_2(j,i) = .true.
1422 
1423  j=jmax
1424  if ( (maske(j,i)==2_i2b) & ! ocean
1425  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1426  ) & ! floating ice point
1427  flag_calving_front_2(j,i) = .true.
1428 
1429 end do
1430 
1431 do j=1, jmax-1
1432 
1433  i=0
1434  if ( (maske(j,i)==2_i2b) & ! ocean
1435  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1436  ) & ! floating ice point
1437  flag_calving_front_2(j,i) = .true.
1438 
1439  i=imax
1440  if ( (maske(j,i)==2_i2b) & ! ocean
1441  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1442  ) & ! floating ice point
1443  flag_calving_front_2(j,i) = .true.
1444 
1445 end do
1446 
1447 #else
1448 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1449 #endif
1450 
1451 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1452 
1453 #if (GRID==2)
1454 
1455 do ir=-imax, imax
1456 do jr=-jmax, jmax
1457 
1458  dist_dxdy(jr,ir) = sqrt( (sq_g11_g(jmax/2,imax/2)*real(ir,dp)*dxi)**2 &
1459  + (sq_g22_g(jmax/2,imax/2)*real(jr,dp)*deta)**2 )
1460 
1461  ! This uses the metric tensor in the center of the Tibet
1462  ! domain for the entire domain; DIRTY TRICK !!!
1463 
1464 end do
1465 end do
1466 
1467 #endif
1468 
1469 !-------- Initial velocities --------
1470 
1471 call calc_temp_melt()
1472 
1473 #if (DYNAMICS==1 || DYNAMICS==2)
1474 
1475 call calc_vxy_b_sia(time, z_sl)
1476 call calc_vxy_sia(dzeta_c, dzeta_t)
1477 
1478 #if (MARGIN==3 || DYNAMICS==2)
1479 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1480 #endif
1481 
1482 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1483 
1484 #if (MARGIN==3)
1485 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1486 #endif
1487 
1488 #elif (DYNAMICS==0)
1489 
1490 call calc_vxy_static()
1491 call calc_vz_static()
1492 
1493 #else
1494  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1495 #endif
1496 
1497 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1498 
1499 !-------- Initial enthalpies --------
1500 
1501 #if (CALCMOD==0 || CALCMOD==-1)
1502 
1503 do i=0, imax
1504 do j=0, jmax
1505 
1506  do kc=0, kcmax
1507  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1508  end do
1509 
1510  do kt=0, ktmax
1511  enth_t(kt,j,i) = enth_c(0,j,i)
1512  end do
1513 
1514 end do
1515 end do
1516 
1517 #elif (CALCMOD==1)
1518 
1519 do i=0, imax
1520 do j=0, jmax
1521 
1522  do kc=0, kcmax
1523  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1524  end do
1525 
1526  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1527  do kt=0, ktmax
1528  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1529  end do
1530  else
1531  do kt=0, ktmax
1532  enth_t(kt,j,i) = enth_c(0,j,i)
1533  end do
1534  end if
1535 
1536 end do
1537 end do
1538 
1539 #elif (CALCMOD==2 || CALCMOD==3)
1540 
1541 do i=0, imax
1542 do j=0, jmax
1543 
1544  do kc=0, kcmax
1545  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1546  end do
1547 
1548  do kt=0, ktmax
1549  enth_t(kt,j,i) = enth_c(0,j,i)
1550  end do
1551 
1552 end do
1553 end do
1554 
1555 #else
1556 
1557 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1558 
1559 #endif
1560 
1561 !-------- Initialize time-series files --------
1562 
1563 ! ------ Time-series file for the ice sheet on the whole
1564 
1565 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1566 
1567 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1568 
1569 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1570 
1571 if (forcing_flag == 1) then
1572 
1573  write(12,1102)
1574  write(12,1103)
1575 
1576  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1577  ' V(m^3) V_g(m^3) V_f(m^3)', &
1578  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1579  ' V_sle(m) V_t(m^3)', &
1580  ' A_t(m^2)',/, &
1581  ' H_max(m) H_t_max(m)', &
1582  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1583  1103 format('----------------------------------------------------', &
1584  '---------------------------------------')
1585 
1586 else if (forcing_flag == 2) then
1587 
1588  write(12,1112)
1589  write(12,1113)
1590 
1591  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1592  ' V(m^3) V_g(m^3) V_f(m^3)', &
1593  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1594  ' V_sle(m) V_t(m^3)', &
1595  ' A_t(m^2)',/, &
1596  ' H_max(m) H_t_max(m)', &
1597  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1598  1113 format('----------------------------------------------------', &
1599  '---------------------------------------')
1600 
1601 end if
1602 
1603 ! ------ Time-series file for deep boreholes
1604 
1605 n_core = 0 ! No boreholes defined
1606 
1607 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1608 
1609 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1610 
1611 write(14,'(1x,a)') '---------------------'
1612 write(14,'(1x,a)') 'No boreholes defined.'
1613 write(14,'(1x,a)') '---------------------'
1614 
1615 !-------- Output of the initial state --------
1616 
1617 #if (defined(OUTPUT_INIT))
1618 
1619 #if (OUTPUT_INIT==0)
1620  flag_init_output = .false.
1621 #elif (OUTPUT_INIT==1)
1622  flag_init_output = .true.
1623 #else
1624  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1625 #endif
1626 
1627 #else
1628  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1629 #endif
1630 
1631 #if (OUTPUT==1)
1632 
1633 #if (ERGDAT==0)
1634  flag_3d_output = .false.
1635 #elif (ERGDAT==1)
1636  flag_3d_output = .true.
1637 #else
1638  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1639 #endif
1640 
1641  if (flag_init_output) &
1642  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1643  flag_3d_output, ndat2d, ndat3d)
1644 
1645 #elif (OUTPUT==2)
1646 
1647 if (time_output(1) <= time_init+eps) then
1648 
1649 #if (ERGDAT==0)
1650  flag_3d_output = .false.
1651 #elif (ERGDAT==1)
1652  flag_3d_output = .true.
1653 #else
1654  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1655 #endif
1656 
1657  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1658  flag_3d_output, ndat2d, ndat3d)
1659 
1660 end if
1661 
1662 #elif (OUTPUT==3)
1663 
1664  flag_3d_output = .false.
1665 
1666  if (flag_init_output) &
1667  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1668  flag_3d_output, ndat2d, ndat3d)
1669 
1670 if (time_output(1) <= time_init+eps) then
1671 
1672  flag_3d_output = .true.
1673 
1674  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1675  flag_3d_output, ndat2d, ndat3d)
1676 
1677 end if
1678 
1679 #else
1680  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1681 #endif
1682 
1683 if (flag_init_output) then
1684  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1685  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1686 end if
1687 
1688 end subroutine sico_init
1689 
1690 !-------------------------------------------------------------------------------
1691 !> Definition of the initial surface and bedrock topography
1692 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1693 !! For ice-free initial topography with relaxed lithosphere
1694 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1695 !<------------------------------------------------------------------------------
1696 subroutine topography1(dxi, deta)
1697 
1698 implicit none
1699 
1700 real(dp), intent(out) :: dxi, deta
1701 
1702 ! integer(i4b) :: i, j
1703 ! integer(i4b) :: ios
1704 ! real(dp) :: xi0, eta0
1705 
1706 dxi=0.0_dp; deta=0.0_dp ! dummy values
1707 
1708 write(6,'(/1x,a)') '>>> topography1: ANF_DAT==1 not defined for this domain!'
1709 stop
1710 
1711 end subroutine topography1
1712 
1713 !-------------------------------------------------------------------------------
1714 !> Definition of the initial surface and bedrock topography
1715 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1716 !! For ice-free initial topography with relaxed lithosphere.
1717 !<------------------------------------------------------------------------------
1718 subroutine topography2(dxi, deta)
1719 
1720 #if (GRID==0 || GRID==1)
1722 #endif
1723 
1724  use metric_m
1725  use topograd_m
1726 
1727 implicit none
1728 
1729 real(dp), intent(out) :: dxi, deta
1730 
1731 integer(i4b) :: i, j, n
1732 integer(i4b) :: ios
1733 real(dp) :: xi0, eta0
1734 character :: ch_dummy
1735 
1736 character(len= 8) :: ch_imax
1737 character(len=128) :: fmt4
1738 character(len=256) :: filename_with_path
1739 
1740 write(ch_imax, fmt='(i8)') imax
1741 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1742 
1743 !-------- Read topography --------
1744 
1745 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1746  trim(zl0_file)
1747 
1748 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1749 
1750 if (ios /= 0) stop ' >>> topography2: Error when opening the zl0 file!'
1751 
1752 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1753  trim(mask_present_file)
1754 
1755 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1756 
1757 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
1758 
1759 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1760 
1761 do j=jmax, 0, -1
1762  read(23, fmt=*) (zl0(j,i), i=0,imax)
1763 end do
1764 
1765 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1766 
1767 do j=jmax, 0, -1
1768  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1769 end do
1770 
1771 close(23, status='keep')
1772 close(24, status='keep')
1773 
1774 !-------- Further stuff --------
1775 
1776 dxi = dlambda *pi_180 ! deg -> rad
1777 deta = dphi *pi_180 ! deg -> rad
1778 
1779 xi0 = lambda_0 *pi_180 ! deg -> rad
1780 eta0 = phi_0 *pi_180 ! deg -> rad
1781 
1782 do i=0, imax
1783 do j=0, jmax
1784 
1785  if (maske(j,i) <= 1) then
1786  maske(j,i) = 1
1787  zs(j,i) = zl0(j,i)
1788  zb(j,i) = zl0(j,i)
1789  zl(j,i) = zl0(j,i)
1790  else ! (maske(j,i) >= 2)
1791  maske(j,i) = 2
1792 #if (MARGIN==1 || MARGIN==2)
1793  zs(j,i) = zl0(j,i)
1794  zb(j,i) = zl0(j,i)
1795 #elif (MARGIN==3)
1796  zs(j,i) = 0.0_dp ! present-day
1797  zb(j,i) = 0.0_dp ! sea level
1798 #endif
1799  zl(j,i) = zl0(j,i)
1800  end if
1801 
1802  xi(i) = xi0 + real(i,dp)*dxi
1803  eta(j) = eta0 + real(j,dp)*deta
1804 
1805  zm(j,i) = zb(j,i)
1806  n_cts(j,i) = -1
1807  kc_cts(j,i) = 0
1808 
1809  h_c(j,i) = 0.0_dp
1810  h_t(j,i) = 0.0_dp
1811 
1812  dzs_dtau(j,i) = 0.0_dp
1813  dzm_dtau(j,i) = 0.0_dp
1814  dzb_dtau(j,i) = 0.0_dp
1815  dzl_dtau(j,i) = 0.0_dp
1816  dh_c_dtau(j,i) = 0.0_dp
1817  dh_t_dtau(j,i) = 0.0_dp
1818 
1819 end do
1820 end do
1821 
1822 maske_old = maske
1823 
1824 !-------- Geographic coordinates, metric tensor,
1825 ! gradients of the topography --------
1826 
1827 do i=0, imax
1828 do j=0, jmax
1829 
1830 #if (GRID==0 || GRID==1) /* Stereographic projection */
1831 
1832  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1833  lambda0, phi0, lambda(j,i), phi(j,i))
1834 
1835 #elif (GRID==2) /* Geographic coordinates */
1836 
1837  lambda(j,i) = xi(i)
1838  phi(j,i) = eta(j)
1839 
1840 #endif
1841 
1842 end do
1843 end do
1844 
1845 call metric()
1846 
1847 #if (TOPOGRAD==0)
1848 call topograd_1(dxi, deta, 1)
1849 #elif (TOPOGRAD==1)
1850 call topograd_2(dxi, deta, 1)
1851 #endif
1852 
1853 !-------- Corresponding area of grid points --------
1854 
1855 do i=0, imax
1856 do j=0, jmax
1857  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1858 end do
1859 end do
1860 
1861 end subroutine topography2
1862 
1863 !-------------------------------------------------------------------------------
1864 !> Definition of the initial surface and bedrock topography
1865 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1866 !! For initial topography from previous simulation.
1867 !<------------------------------------------------------------------------------
1868 subroutine topography3(dxi, deta, z_sl, anfdatname)
1869 
1870  use read_m, only : read_erg_nc
1871 
1872 #if (GRID==0 || GRID==1)
1874 #endif
1875 
1876  use metric_m
1877  use topograd_m
1878 
1879 implicit none
1880 
1881 character(len=100), intent(in) :: anfdatname
1882 
1883 real(dp), intent(out) :: dxi, deta, z_sl
1884 
1885 integer(i4b) :: i, j, n
1886 integer(i4b) :: ios
1887 character(len=256) :: filename_with_path
1888 character :: ch_dummy
1889 
1890 !-------- Read data from time-slice file of previous simulation --------
1891 
1892 call read_erg_nc(z_sl, anfdatname)
1893 
1894 !-------- Read topography of the relaxed bedrock --------
1895 
1896 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1897  trim(zl0_file)
1898 
1899 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1900 
1901 if (ios /= 0) stop ' >>> topography3: Error when opening the zl0 file!'
1902 
1903 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1904 
1905 do j=jmax, 0, -1
1906  read(23, fmt=*) (zl0(j,i), i=0,imax)
1907 end do
1908 
1909 close(23, status='keep')
1910 
1911 !-------- Further stuff --------
1912 
1913 dxi = dlambda *pi_180 ! deg -> rad
1914 deta = dphi *pi_180 ! deg -> rad
1915 
1916 !-------- Geographic coordinates, metric tensor,
1917 ! gradients of the topography --------
1918 
1919 do i=0, imax
1920 do j=0, jmax
1921 
1922 #if (GRID==0 || GRID==1) /* Stereographic projection */
1923 
1924  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1925  lambda0, phi0, lambda(j,i), phi(j,i))
1926 
1927 #elif (GRID==2) /* Geographic coordinates */
1928 
1929  lambda(j,i) = xi(i)
1930  phi(j,i) = eta(j)
1931 
1932 #endif
1933 
1934 end do
1935 end do
1936 
1937 call metric()
1938 
1939 #if (TOPOGRAD==0)
1940 call topograd_1(dxi, deta, 1)
1941 #elif (TOPOGRAD==1)
1942 call topograd_2(dxi, deta, 1)
1943 #endif
1944 
1945 !-------- Corresponding area of grid points --------
1946 
1947 do i=0, imax
1948 do j=0, jmax
1949  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1950 end do
1951 end do
1952 
1953 end subroutine topography3
1954 
1955 !-------------------------------------------------------------------------------
1956 
1957 end module sico_init_m
1958 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
subroutine, public read_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:50
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:765
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_present
temp_mm_present(j,i,n): Present-day mean monthly surface temperature
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:59
real(dp) rho_w
RHO_W: Density of pure water.
subroutine, public calc_enhance_2()
Computation of the flow enhancement factor. Case ENHMOD==2: two different values depending on age for...
real(dp) a
A: Semi-major axis of the planet.
subroutine, public calc_vz_floating(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz for floating ice.
Definition: calc_vz_m.F90:259
subroutine, public output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in native binary or NetCDF format.
Definition: output_m.F90:59
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:813
integer(i2b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz_m.F90:37
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
subroutine, public stereo_inv_ellipsoid(x_val, y_val, A, B, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for an ellipsoidal planet.
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_m.F90:56
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
Definition: output_m.F90:4721
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
Computation of the flow enhancement factor.
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
Computation of the horizontal velocity vx, vy.
Definition: calc_vxy_m.F90:35
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Initial temperature, water content and age.
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
Computation of the melting and basal temperatures.
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
Definition: topograd_m.F90:201
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Definition: metric_m.F90:37
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
subroutine, public calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_m.F90:374
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_tmp_min, n_tmp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
Comparison of single- or double-precision floating-point numbers.
subroutine, public calc_enhance_1()
Computation of the flow enhancement factor. Case ENHMOD==1: constant for grounded ice...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_m.F90:478
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:56
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
subroutine, public output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format (and optionally in NetCDF format).
Definition: output_m.F90:3377
real(dp), dimension(-jmax:jmax,-imax:imax) dist_dxdy
dist_dxdy(jr,ir): Distance between grid points with delta_i=ir, delta_j=jr
Calculation of topography gradients on the staggered grid and on the grid points (including length re...
Definition: topograd_m.F90:39
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
subroutine, public metric()
Main routine of module metric_m: Definition of the components g11 and g22 of the metric tensor of the...
Definition: metric_m.F90:54
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) r_t
R_T: Coefficient of the water-content dependence in the rate factor for temperate ice...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
subroutine, public calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Main subroutine of calc_dxyz_m: Computation of all components of the strain-rate tensor, the full effective strain rate and the shear fraction.
Definition: calc_dxyz_m.F90:56
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
real(dp) l
L: Latent heat of ice.
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
Writing of output data on files.
Definition: output_m.F90:36
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
subroutine, public calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
Computation of the temperature integral of the specific heat of ice as a table (c_int_table). Further, definition of the latent heat of ice.
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
real(dp), dimension(0:jmax, 0:imax, 12) precip_lgm_anom
precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly precipitation rate at the...
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_lgm_anom
temp_mm_lgm_anom(j,i,n): LGM anomaly (difference LGM - present) of the mean monthly surface temperatu...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
subroutine, public calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
Definition: calc_vxy_m.F90:53
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) flex_rig_lith
flex_rig_lith(j,i): Flexural rigidity of the lithosphere
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:jmax, 0:imax) sq_g22_g
sq_g22_g(j,i): Square root of the coefficient g22 of the metric tensor on grid point (i...