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 if (approx_equal(dx, 80.0_dp, eps_sp_dp)) then
230 
231  if ((imax /= 156).or.(jmax /= 156)) then
232  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
233  end if
234 
235 else if (approx_equal(dx, 40.0_dp, eps_sp_dp)) then
236 
237  if ((imax /= 312).or.(jmax /= 312)) then
238  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
239  end if
240 
241 else if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
242 
243  if ((imax /= 624).or.(jmax /= 624)) then
244  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
245  end if
246 
247 else
248  stop ' >>> sico_init: DX wrong!'
249 end if
250 
251 #elif (GRID==2)
252 
253  stop ' >>> sico_init: GRID==2 not allowed for the northern hemisphere application!'
254 
255 #endif
256 
257 !-------- Compatibility check of the thermodynamics mode
258 ! (cold vs. polythermal vs. enthalpy method)
259 ! and the number of grid points in the lower (kt) ice domain --------
260 
261 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
262 
263 if (ktmax > 2) then
264  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
265  write(6, fmt='(a)') ' the separate kt domain is redundant.'
266  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
267  write(6, fmt='(a)') ' '
268 end if
269 
270 #endif
271 
272 !-------- Compatibility check of surface-temperature and precipitation
273 ! determination by interpolation between present and LGM values
274 ! with a glacial index --------
275 
276 #if (TSURFACE == 5 && ACCSURFACE != 5)
277 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
278 #endif
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 !-------- Compatibility check of discretization schemes for the horizontal and
285 ! vertical advection terms in the temperature and age equations --------
286 
287 #if (ADV_HOR==1)
288 stop ' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!'
289 #endif
290 
291 !-------- Check whether for the shallow shelf
292 ! or shelfy stream approximation
293 ! the chosen grid is Cartesian coordinates
294 ! without distortion correction (GRID==0) --------
295 
296 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
297 #if (GRID != 0)
298 write(6, fmt='(a)') ' >>> sico_init: Distortion correction not implemented'
299 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
300 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
301 write(6, fmt='(a)') ' -> GRID==0 required!'
302 stop
303 #endif
304 #endif
305 
306 !-------- Setting of forcing flag --------
307 
308 #if (TSURFACE <= 4)
309 
310 forcing_flag = 1 ! forcing by delta_ts
311 
312 #elif (TSURFACE == 5)
313 
314 forcing_flag = 2 ! forcing by glac_index
315 
316 #endif
317 
318 !-------- Initialization of numerical time steps --------
319 
320 dtime0 = dtime0
321 dtime_temp0 = dtime_temp0
322 #if (REBOUND==2)
323 dtime_wss0 = dtime_wss0
324 #endif
325 
326 !-------- Further initializations --------
327 
328 dzeta_c = 1.0_dp/real(kcmax,dp)
329 dzeta_t = 1.0_dp/real(ktmax,dp)
330 dzeta_r = 1.0_dp/real(krmax,dp)
331 
332 ndat2d = 1
333 ndat3d = 1
334 
335 !-------- General abbreviations --------
336 
337 ! ------ kc domain
338 
339 if (deform >= eps) then
340 
341  flag_aa_nonzero = .true. ! non-equidistant grid
342 
343  aa = deform
344  ea = exp(aa)
345 
346  kc=0
347  zeta_c(kc) = 0.0_dp
348  eaz_c(kc) = 1.0_dp
349  eaz_c_quotient(kc) = 0.0_dp
350 
351  do kc=1, kcmax-1
352  zeta_c(kc) = kc*dzeta_c
353  eaz_c(kc) = exp(aa*zeta_c(kc))
354  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
355  end do
356 
357  kc=kcmax
358  zeta_c(kc) = 1.0_dp
359  eaz_c(kc) = exp(aa)
360  eaz_c_quotient(kc) = 1.0_dp
361 
362 else
363 
364  flag_aa_nonzero = .false. ! equidistant grid
365 
366  aa = 0.0_dp
367  ea = 1.0_dp
368 
369  kc=0
370  zeta_c(kc) = 0.0_dp
371  eaz_c(kc) = 1.0_dp
372  eaz_c_quotient(kc) = 0.0_dp
373 
374  do kc=1, kcmax-1
375  zeta_c(kc) = kc*dzeta_c
376  eaz_c(kc) = 1.0_dp
377  eaz_c_quotient(kc) = zeta_c(kc)
378  end do
379 
380  kc=kcmax
381  zeta_c(kc) = 1.0_dp
382  eaz_c(kc) = 1.0_dp
383  eaz_c_quotient(kc) = 1.0_dp
384 
385 end if
386 
387 ! ------ kt domain
388 
389 kt=0
390 zeta_t(kt) = 0.0_dp
391 
392 do kt=1, ktmax-1
393  zeta_t(kt) = kt*dzeta_t
394 end do
395 
396 kt=ktmax
397 zeta_t(kt) = 1.0_dp
398 
399 ! ------ kr domain
400 
401 kr=0
402 zeta_r(kr) = 0.0_dp
403 
404 do kr=1, krmax-1
405  zeta_r(kr) = kr*dzeta_r
406 end do
407 
408 kr=krmax
409 zeta_r(kr) = 1.0_dp
410 
411 !-------- Construction of a vector (with index n) from a 2-d array
412 ! (with indices i, j) by diagonal numbering --------
413 
414 n=1
415 
416 do m=0, imax+jmax
417  do i=m, 0, -1
418  j = m-i
419  if ((i <= imax).and.(j <= jmax)) then
420  ii(n) = i
421  jj(n) = j
422  nn(j,i) = n
423  n=n+1
424  end if
425  end do
426 end do
427 
428 !-------- Specification of current simulation --------
429 
430 runname = runname
431 anfdatname = anfdatname
432 
433 #if (defined(YEAR_ZERO))
435 #else
436 year_zero = 2000.0_dp ! default value 2000 CE
437 #endif
438 
439 time_init0 = time_init0
440 time_end0 = time_end0
441 dtime_ser0 = dtime_ser0
442 
443 #if (OUTPUT==1 || OUTPUT==3)
444 dtime_out0 = dtime_out0
445 #endif
446 
447 #if (OUTPUT==2 || OUTPUT==3)
448 n_output = n_output
449 time_output0( 1) = time_out0_01
450 time_output0( 2) = time_out0_02
451 time_output0( 3) = time_out0_03
452 time_output0( 4) = time_out0_04
453 time_output0( 5) = time_out0_05
454 time_output0( 6) = time_out0_06
455 time_output0( 7) = time_out0_07
456 time_output0( 8) = time_out0_08
457 time_output0( 9) = time_out0_09
458 time_output0(10) = time_out0_10
459 time_output0(11) = time_out0_11
460 time_output0(12) = time_out0_12
461 time_output0(13) = time_out0_13
462 time_output0(14) = time_out0_14
463 time_output0(15) = time_out0_15
464 time_output0(16) = time_out0_16
465 time_output0(17) = time_out0_17
466 time_output0(18) = time_out0_18
467 time_output0(19) = time_out0_19
468 time_output0(20) = time_out0_20
469 #endif
470 
471 !-------- Write log file --------
472 
473 shell_command = 'if [ ! -d'
474 shell_command = trim(shell_command)//' '//outpath
475 shell_command = trim(shell_command)//' '//'] ; then mkdir'
476 shell_command = trim(shell_command)//' '//outpath
477 shell_command = trim(shell_command)//' '//'; fi'
478 call system(trim(shell_command))
479  ! Check whether directory OUTPATH exists. If not, it is created.
480 
481 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
482 
483 open(10, iostat=ios, file=trim(filename_with_path), status='new')
484 
485 if (ios /= 0) stop ' >>> sico_init: Error when opening the log file!'
486 
487 write(10, fmt=trim(fmt1)) 'Computational domain:'
488 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
489 write(10, fmt=trim(fmt1)) ' '
490 
491 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
492 write(10, fmt=trim(fmt1)) ' '
493 
494 write(10, fmt=trim(fmt2)) 'imax =', imax
495 write(10, fmt=trim(fmt2)) 'jmax =', jmax
496 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
497 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
498 write(10, fmt=trim(fmt2)) 'krmax =', krmax
499 write(10, fmt=trim(fmt1)) ' '
500 
501 write(10, fmt=trim(fmt3)) 'a =', aa
502 write(10, fmt=trim(fmt1)) ' '
503 
504 #if (GRID==0 || GRID==1)
505 write(10, fmt=trim(fmt3)) 'x0 =', x0
506 write(10, fmt=trim(fmt3)) 'y0 =', y0
507 write(10, fmt=trim(fmt3)) 'dx =', dx
508 #elif (GRID==2)
509 stop ' >>> sico_init: GRID==2 not allowed for this application!'
510 #endif
511 write(10, fmt=trim(fmt1)) ' '
512 
513 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
514 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
515 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
516 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
517 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
518 #if (REBOUND==2)
519 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
520 #endif
521 write(10, fmt=trim(fmt1)) ' '
522 
523 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
524 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
525 #if (ANF_DAT==1)
526 #if (defined(ZB_PRESENT_FILE))
527 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
528 #endif
529 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
530 #endif
531 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
532 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
533 #if (ANF_DAT==1 && defined(TEMP_INIT))
534 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
535 #endif
536 #if (ANF_DAT==3)
537 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
538 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
539 #endif
540 write(10, fmt=trim(fmt1)) ' '
541 
542 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
543 write(10, fmt=trim(fmt1)) ' '
544 
545 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
546 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
547 #if (CALCTHK==2 || CALCTHK==5)
548 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
549 #if (ITER_MAX_SOR>0)
550 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
551 #endif
552 #endif
553 write(10, fmt=trim(fmt1)) ' '
554 #endif
555 
556 write(10, fmt=trim(fmt1)) 'temp_mm_present file = '//temp_mm_present_file
557 #if (TSURFACE==1)
558 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
559 #elif (TSURFACE==3)
560 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
561 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
562 #elif (TSURFACE==4)
563 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
564 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
565 #elif (TSURFACE==5)
566 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
567 write(10, fmt=trim(fmt1)) 'temp_mm_anom file = '//temp_mm_anom_file
568 write(10, fmt=trim(fmt3)) 'temp_mm_anom fact = ', temp_mm_anom_fact
569 #endif
570 
571 write(10, fmt=trim(fmt1)) 'precip_mm_present file = '//precip_mm_present_file
572 #if (ACCSURFACE==1)
573 write(10, fmt=trim(fmt3)) 'accfact =', accfact
574 #elif (ACCSURFACE==2 || ACCSURFACE==3)
575 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
576 #elif (ACCSURFACE==5)
577 write(10, fmt=trim(fmt1)) 'precip_mm_anom file = '//precip_mm_anom_file
578 write(10, fmt=trim(fmt3)) 'precip_mm_anom fact = ', precip_mm_anom_fact
579 #endif
580 #if (ACCSURFACE <= 3)
581 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
582 #if (ELEV_DESERT == 1)
583 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
584 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
585 #endif
586 #endif
587 
588 #if (ABLSURFACE==3)
589 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
590 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
591 #endif
592 
593 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
594 #if (SEA_LEVEL==1)
595 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
596 #elif (SEA_LEVEL==3)
597 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
598 #endif
599 write(10, fmt=trim(fmt1)) ' '
600 
601 #if (MARGIN==2)
602 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
603 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
604 write(10, fmt=trim(fmt1)) ' '
605 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
606  || marine_ice_calving==6 || marine_ice_calving==7)
607 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
608 write(10, fmt=trim(fmt1)) ' '
609 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
610 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
611 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
612 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
613 write(10, fmt=trim(fmt1)) ' '
614 #endif
615 #elif (MARGIN==3)
616 #if (ICE_SHELF_CALVING==2)
617 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
618 write(10, fmt=trim(fmt1)) ' '
619 #endif
620 #endif
621 
622 #if (defined(BASAL_HYDROLOGY))
623 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
624 #endif
625 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
626 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
627 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
628 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
629 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
630 #if (defined(TIME_RAMP_UP_SLIDE))
631 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
632 #endif
633 #if (SLIDE_LAW==2)
634 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
635 #endif
636 #if (BASAL_HYDROLOGY==1)
637 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
638 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
639 #endif
640 write(10, fmt=trim(fmt1)) ' '
641 
642 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
643 #if (Q_GEO_MOD==1)
644 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
645 #elif (Q_GEO_MOD==2)
646 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
647 #endif
648 write(10, fmt=trim(fmt1)) ' '
649 
650 #if (defined(MARINE_ICE_BASAL_MELTING))
651 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
652 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
653 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
654 #endif
655 write(10, fmt=trim(fmt1)) ' '
656 #endif
657 
658 #if (MARGIN==3)
659 write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
660 #if (FLOATING_ICE_BASAL_MELTING<=3)
661 write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
662 write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
663 #endif
664 write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
665 write(10, fmt=trim(fmt3)) 'z_abyss =', z_abyss
666 #if (FLOATING_ICE_BASAL_MELTING==4)
667 write(10, fmt=trim(fmt3)) 'temp_ocean =', temp_ocean
668 write(10, fmt=trim(fmt3)) 'Omega_qbm =', omega_qbm
669 write(10, fmt=trim(fmt3)) 'alpha_qbm =', alpha_qbm
670 #endif
671 #if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
672 write(10, fmt=trim(fmt3)) 'H_w_0 =', h_w_0
673 #endif
674 write(10, fmt=trim(fmt1)) ' '
675 #endif
676 
677 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
678 #if (REBOUND==1)
679 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
680 #endif
681 #if (REBOUND==1 || REBOUND==2)
682 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
683 #if (TIME_LAG_MOD==1)
684 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
685 #elif (TIME_LAG_MOD==2)
686 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
687 #else
688 stop ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
689 #endif
690 #endif
691 #if (REBOUND==2)
692 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
693 #if (FLEX_RIG_MOD==1)
694 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
695 #elif (FLEX_RIG_MOD==2)
696 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
697 #else
698 stop ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
699 #endif
700 #endif
701 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
702 write(10, fmt=trim(fmt1)) ' '
703 
704 #if (FLOW_LAW==2)
705 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
706 write(10, fmt=trim(fmt1)) ' '
707 #endif
708 #if (FIN_VISC==2)
709 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
710 write(10, fmt=trim(fmt1)) ' '
711 #endif
712 
713 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
714 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
715 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
716 #endif
717 #if (ENHMOD==2 || ENHMOD==3)
718 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
719 #endif
720 #if (ENHMOD==2)
721 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
722 #endif
723 #if (ENHMOD==3)
724 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
725 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
726 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
727 #endif
728 #if (ENHMOD==4 || ENHMOD==5)
729 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
730 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
731 #endif
732 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
733 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
734 #endif
735 write(10, fmt=trim(fmt1)) ' '
736 
737 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
738 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
739 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
740 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
741 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
742 #if (defined(QBM_MIN))
743 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
744 #elif (defined(QB_MIN)) /* obsolete */
745 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
746 #endif
747 #if (defined(QBM_MAX))
748 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
749 #elif (defined(QB_MAX)) /* obsolete */
750 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
751 #endif
752 write(10, fmt=trim(fmt3)) 'age_min =', age_min
753 write(10, fmt=trim(fmt3)) 'age_max =', age_max
754 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
755 #if (ADV_VERT==1)
756 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
757 #endif
758 write(10, fmt=trim(fmt1)) ' '
759 
760 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
761 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
762 #if (defined(LIS_OPTS))
763 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
764 #endif
765 #if (defined(N_ITER_SSA))
766 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
767 #endif
768 #if (defined(RELAX_FACT_SSA))
769 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
770 #endif
771 #endif
772 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
773 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
774 #endif
775 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
776 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
777 #endif
778 write(10, fmt=trim(fmt1)) ' '
779 
780 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
781 #if (CALCMOD==-1 && defined(TEMP_CONST))
782 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
783 #endif
784 #if (CALCMOD==-1 && defined(AGE_CONST))
785 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
786 #endif
787 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
788 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
789 #endif
790 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
791 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
792 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
793 #if (MARGIN==2)
794 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
795 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
796 #elif (MARGIN==3)
797 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
798 #endif
799 #if (defined(THK_EVOL))
800 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
801 #else
802 stop ' >>> sico_init: Define THK_EVOL in header file!'
803 #endif
804 #if (defined(CALCTHK))
805 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
806 #else
807 stop ' >>> sico_init: Define CALCTHK in header file!'
808 #endif
809 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
810 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
811 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
812 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
813 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
814 #if (ACCSURFACE==5)
815 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
816 #endif
817 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
818 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
819 #if (defined(MB_ACCOUNT))
820 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
821 #endif
822 write(10, fmt=trim(fmt1)) ' '
823 
824 #if (defined(OUT_TIMES))
825 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
826 #endif
827 #if (defined(OUTPUT_INIT))
828 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
829 #endif
830 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
831 #if (OUTPUT==1 || OUTPUT==3)
832 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
833 #endif
834 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
835 #if (OUTPUT==1 || OUTPUT==2)
836 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
837 #endif
838 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
839 #if (OUTPUT==2 || OUTPUT==3)
840 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
841 do n=1, n_output
842  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
843 end do
844 #endif
845 write(10, fmt=trim(fmt1)) ' '
846 
847 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
848 
849 close(10, status='keep')
850 
851 !-------- Conversion of time quantities --------
852 
853 year_zero = year_zero*year_sec ! a --> s
854 time_init = time_init0*year_sec ! a --> s
855 time_end = time_end0*year_sec ! a --> s
856 dtime = dtime0*year_sec ! a --> s
857 dtime_temp = dtime_temp0*year_sec ! a --> s
858 #if (REBOUND==2)
859 dtime_wss = dtime_wss0*year_sec ! a --> s
860 #endif
861 dtime_ser = dtime_ser0*year_sec ! a --> s
862 #if (OUTPUT==1 || OUTPUT==3)
863 dtime_out = dtime_out0*year_sec ! a --> s
864 #endif
865 #if (OUTPUT==2 || OUTPUT==3)
866 do n=1, n_output
867  time_output(n) = time_output0(n)*year_sec ! a --> s
868 end do
869 #endif
870 
871 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) &
872  stop ' >>> sico_init: dtime_temp must be a multiple of dtime!'
873 
874 #if (REBOUND==2)
875 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) &
876  stop ' >>> sico_init: dtime_wss must be a multiple of dtime!'
877 #endif
878 
879 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) &
880  stop ' >>> sico_init: dtime_ser must be a multiple of dtime!'
881 
882 #if (OUTPUT==1 || OUTPUT==3)
883 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) &
884  stop ' >>> sico_init: dtime_out must be a multiple of dtime!'
885 #endif
886 
887 time = time_init
888 
889 !-------- Reading of measurements for present monthly-mean precipitation --------
890 
891 #if (GRID==0 || GRID==1)
892 
893 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
894  trim(precip_mm_present_file)
895 
896 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
897 
898 #elif (GRID==2)
899 
900 stop ' >>> sico_init: GRID==2 not allowed for the northern hemisphere application!'
901 
902 #endif
903 
904 if (ios /= 0) stop ' >>> sico_init: Error when opening the precipitation file!'
905 
906 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
907 
908 do n=1, 12 ! month counter
909  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
910  do j=jmax, 0, -1
911  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
912  end do
913 end do
914 
915 close(21, status='keep')
916 
917 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
918 
919 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
920  ! mm/a water equiv. -> m/s ice equiv.
921 
922 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
923 
924 #if (ACCSURFACE==5)
925 
926 #if (GRID==0 || GRID==1)
927 
928 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
929  trim(precip_mm_anom_file)
930 
931 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
932 
933 #endif
934 
935 if (ios /= 0) stop ' >>> sico_init: Error when opening the precip anomaly file!'
936 
937 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
938 
939 do n=1, 12 ! month counter
940  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
941  do j=jmax, 0, -1
942  read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
943  end do
944 end do
945 
946 close(21, status='keep')
947 
948 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
949 
950 do i=0, imax
951 do j=0, jmax
952 
953 #if (PRECIP_ANOM_INTERPOL==1)
954  do n=1, 12 ! month counter
955  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
956  end do
957 #elif (PRECIP_ANOM_INTERPOL==2)
958  do n=1, 12 ! month counter
959  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
960  end do
961 #else
962  stop ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
963 #endif
964 
965 end do
966 end do
967 
968 #endif
969 
970 !-------- Mean accumulation --------
971 
972 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
973  ! mm/a water equiv. -> m/s ice equiv.
974 
975 !-------- Reading of present topography mask --------
976 
977 #if (GRID==0 || GRID==1)
978 
979 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
980  trim(mask_present_file)
981 
982 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
983 
984 #elif (GRID==2)
985 
986 stop ' >>> sico_init: GRID==2 not allowed for the northern hemisphere application!'
987 
988 #endif
989 
990 if (ios /= 0) stop ' >>> sico_init: Error when opening the mask file!'
991 
992 do m=1, 6; read(24, fmt='(a)') ch_dummy; end do
993 
994 do j=jmax, 0, -1
995  read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
996 end do
997 
998 close(24, status='keep')
999 
1000 !-------- Reading of data for present monthly-mean surface temperature --------
1001 
1002 #if (GRID==0 || GRID==1)
1003 
1004 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1005  trim(temp_mm_present_file)
1006 
1007 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1008 
1009 #elif (GRID==2)
1010 
1011 stop ' >>> sico_init: GRID==2 not allowed for the northern hemisphere application!'
1012 
1013 #endif
1014 
1015 if (ios /= 0) stop ' >>> sico_init: Error when opening the temperature file!'
1016 
1017 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1018 
1019 do n=1, 12 ! month counter
1020  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1021  do j=jmax, 0, -1
1022  read(21, fmt=*) (temp_mm_present(j,i,n), i=0,imax)
1023  end do
1024 end do
1025 
1026 close(21, status='keep')
1027 
1028 !-------- Reading of LGM monthly-mean surface-temperature anomalies --------
1029 
1030 #if (TSURFACE==5)
1031 
1032 #if (GRID==0 || GRID==1)
1033 
1034 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1035  trim(temp_mm_anom_file)
1036 
1037 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1038 
1039 #endif
1040 
1041 if (ios /= 0) stop ' >>> sico_init: Error when opening the temperature anomaly file!'
1042 
1043 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1044 
1045 do n=1, 12 ! month counter
1046  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1047  do j=jmax, 0, -1
1048  read(21, fmt=*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
1049  end do
1050 end do
1051 
1052 close(21, status='keep')
1053 
1054 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
1055 
1056 #endif
1057 
1058 !-------- Present reference elevation
1059 ! (for precipitation and surface-temperature data) --------
1060 
1061 #if (GRID==0 || GRID==1)
1062 
1063 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1064  trim(zs_present_file)
1065 
1066 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1067 
1068 #elif (GRID==2)
1069 
1070 stop ' >>> sico_init: GRID==2 not allowed for the northern hemisphere application!'
1071 
1072 #endif
1073 
1074 if (ios /= 0) stop ' >>> sico_init: Error when opening the zs file!'
1075 
1076 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1077 
1078 do j=jmax, 0, -1
1079  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1080 end do
1081 
1082 close(21, status='keep')
1083 
1084 ! ------ Reset bathymetry data (sea floor elevation) to the
1085 ! sea surface
1086 
1087 do i=0, imax
1088 do j=0, jmax
1089  if (maske_ref(j,i) >= 2) zs_ref(j,i) = 0.0_dp
1090 end do
1091 end do
1092 
1093 !-------- Read data for delta_ts --------
1094 
1095 #if (TSURFACE==4)
1096 
1097 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1098 
1099 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1100 
1101 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for delta_ts!'
1102 
1103 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1104 
1105 if (ch_dummy /= '#') then
1106  write(6, fmt=*) ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max'
1107  write(6, fmt=*) ' not defined in data file for delta_ts!'
1108  stop
1109 end if
1110 
1112 
1113 allocate(griptemp(0:ndata_grip))
1114 
1115 do n=0, ndata_grip
1116  read(21, fmt=*) d_dummy, griptemp(n)
1117 end do
1118 
1119 close(21, status='keep')
1120 
1121 #endif
1122 
1123 !-------- Read data for the glacial index --------
1124 
1125 #if (TSURFACE==5)
1126 
1127 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1128 
1129 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1130 
1131 if (ios /= 0) stop ' >>> sico_init: Error when opening the glacial-index file!'
1132 
1133 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1134 
1135 if (ch_dummy /= '#') then
1136  write(6, fmt=*) ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max'
1137  write(6, fmt=*) ' not defined in glacial-index file!'
1138  stop
1139 end if
1140 
1142 
1143 allocate(glacial_index(0:ndata_gi))
1144 
1145 do n=0, ndata_gi
1146  read(21, fmt=*) d_dummy, glacial_index(n)
1147 end do
1148 
1149 close(21, status='keep')
1150 
1151 #endif
1152 
1153 !-------- Read data for z_sl --------
1154 
1155 #if (SEA_LEVEL==3)
1156 
1157 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1158 
1159 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1160 
1161 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for z_sl!'
1162 
1163 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1164 
1165 if (ch_dummy /= '#') then
1166  write(6, fmt=*) ' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1167  write(6, fmt=*) ' not defined in data file for z_sl!'
1168  stop
1169 end if
1170 
1172 
1173 allocate(specmap_zsl(0:ndata_specmap))
1174 
1175 do n=0, ndata_specmap
1176  read(21, fmt=*) d_dummy, specmap_zsl(n)
1177 end do
1178 
1179 close(21, status='keep')
1180 
1181 #endif
1182 
1183 !-------- Determination of the geothermal heat flux --------
1184 
1185 #if (Q_GEO_MOD==1)
1186 
1187 ! ------ Constant value
1188 
1189 do i=0, imax
1190 do j=0, jmax
1191  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1192 end do
1193 end do
1194 
1195 #elif (Q_GEO_MOD==2)
1196 
1197 ! ------ Read data from file
1198 
1199 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1200  trim(q_geo_file)
1201 
1202 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1203 
1204 if (ios /= 0) stop ' >>> sico_init: Error when opening the qgeo file!'
1205 
1206 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1207 
1208 do j=jmax, 0, -1
1209  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1210 end do
1211 
1212 close(21, status='keep')
1213 
1214 do i=0, imax
1215 do j=0, jmax
1216  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1217 end do
1218 end do
1219 
1220 #endif
1221 
1222 !-------- Reading of tabulated kei function--------
1223 
1224 #if (REBOUND==0 || REBOUND==1)
1225 
1226 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1227  ! dummy values
1228 #elif (REBOUND==2)
1229 
1230 call read_kei()
1231 
1232 #endif
1233 
1234 !-------- Determination of the time lag
1235 ! of the relaxing asthenosphere --------
1236 
1237 #if (REBOUND==1 || REBOUND==2)
1238 
1239 #if (TIME_LAG_MOD==1)
1240 
1241 time_lag_asth = time_lag*year_sec ! a -> s
1242 
1243 #elif (TIME_LAG_MOD==2)
1244 
1245 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1246  trim(time_lag_file)
1247 
1248 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1249 
1250 if (ios /= 0) stop ' >>> sico_init: Error when opening the time-lag file!'
1251 
1252 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1253 
1254 do j=jmax, 0, -1
1255  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1256 end do
1257 
1258 close(29, status='keep')
1259 
1260 time_lag_asth = time_lag_asth*year_sec ! a -> s
1261 
1262 #endif
1263 
1264 #elif (REBOUND==0)
1265 
1266 time_lag_asth = 0.0_dp ! dummy values
1267 
1268 #endif
1269 
1270 !-------- Determination of the flexural rigidity of the lithosphere --------
1271 
1272 #if (REBOUND==2)
1273 
1274 #if (FLEX_RIG_MOD==1)
1275 
1276 flex_rig_lith = flex_rig
1277 
1278 #elif (FLEX_RIG_MOD==2)
1279 
1280 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1281  trim(flex_rig_file)
1282 
1283 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1284 
1285 if (ios /= 0) stop ' >>> sico_init: Error when opening the flex-rig file!'
1286 
1287 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1288 
1289 do j=jmax, 0, -1
1290  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1291 end do
1292 
1293 close(29, status='keep')
1294 
1295 #endif
1296 
1297 #elif (REBOUND==0 || REBOUND==1)
1298 
1299 flex_rig_lith = 0.0_dp ! dummy values
1300 
1301 #endif
1302 
1303 !-------- Definition of initial values --------
1304 
1305 ! ------ Present topography
1306 
1307 #if (ANF_DAT==1)
1308 
1309 stop ' >>> sico_init: ANF_DAT==1 not allowed for northern hemisphere application!'
1310 
1311 ! ------ Ice-free, relaxed bedrock
1312 
1313 #elif (ANF_DAT==2)
1314 
1315 call topography2(dxi, deta)
1316 
1317 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1318 
1319 call boundary(time_init, dtime, dxi, deta, &
1320  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1321 
1322 as_perp_apl = 0.0_dp
1323 
1324 q_bm = 0.0_dp
1325 q_tld = 0.0_dp
1326 q_b_tot = 0.0_dp
1327 
1328 p_b_w = 0.0_dp
1329 h_w = 0.0_dp
1330 
1331 call init_temp_water_age_2()
1332 
1333 #if (ENHMOD==1)
1334  call calc_enhance_1()
1335 #elif (ENHMOD==2)
1336  call calc_enhance_2()
1337 #elif (ENHMOD==3)
1338  call calc_enhance_3(time_init)
1339 #elif (ENHMOD==4)
1340  call calc_enhance_4()
1341 #elif (ENHMOD==5)
1342  call calc_enhance_5()
1343 #else
1344  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1345 #endif
1346 
1347 ! ------ Read initial state from output of previous simulation
1348 
1349 #elif (ANF_DAT==3)
1350 
1351 call topography3(dxi, deta, z_sl, anfdatname)
1352 
1353 call boundary(time_init, dtime, dxi, deta, &
1354  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1355 
1356 where ((maske==0_i2b).or.(maske==3_i2b))
1357  ! grounded or floating ice
1359 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1360  as_perp_apl = 0.0_dp
1361 end where
1362 
1363 q_b_tot = q_bm + q_tld
1364 
1365 #endif
1366 
1367 !-------- Inner-point flag --------
1368 
1369 flag_inner_point = .true.
1370 
1371 flag_inner_point(0,:) = .false.
1372 flag_inner_point(jmax,:) = .false.
1373 
1374 flag_inner_point(:,0) = .false.
1375 flag_inner_point(:,imax) = .false.
1376 
1377 !-------- Grounding line and calving front flags --------
1378 
1379 flag_grounding_line_1 = .false.
1380 flag_grounding_line_2 = .false.
1381 
1382 flag_calving_front_1 = .false.
1383 flag_calving_front_2 = .false.
1384 
1385 #if (MARGIN==1 || MARGIN==2)
1386 
1387 continue
1388 
1389 #elif (MARGIN==3)
1390 
1391 do i=1, imax-1
1392 do j=1, jmax-1
1393 
1394  if ( (maske(j,i)==0_i2b) & ! grounded ice
1395  .and. &
1396  ( (maske(j,i+1)==3_i2b) & ! with
1397  .or.(maske(j,i-1)==3_i2b) & ! one
1398  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1399  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1400  ) &
1401  flag_grounding_line_1(j,i) = .true.
1402 
1403  if ( (maske(j,i)==3_i2b) & ! floating ice
1404  .and. &
1405  ( (maske(j,i+1)==0_i2b) & ! with
1406  .or.(maske(j,i-1)==0_i2b) & ! one
1407  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1408  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1409  ) &
1410  flag_grounding_line_2(j,i) = .true.
1411 
1412  if ( (maske(j,i)==3_i2b) & ! floating ice
1413  .and. &
1414  ( (maske(j,i+1)==2_i2b) & ! with
1415  .or.(maske(j,i-1)==2_i2b) & ! one
1416  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1417  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1418  ) &
1419  flag_calving_front_1(j,i) = .true.
1420 
1421  if ( (maske(j,i)==2_i2b) & ! ocean
1422  .and. &
1423  ( (maske(j,i+1)==3_i2b) & ! with
1424  .or.(maske(j,i-1)==3_i2b) & ! one
1425  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1426  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1427  ) &
1428  flag_calving_front_2(j,i) = .true.
1429 
1430 end do
1431 end do
1432 
1433 do i=1, imax-1
1434 
1435  j=0
1436  if ( (maske(j,i)==2_i2b) & ! ocean
1437  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1438  ) & ! floating ice point
1439  flag_calving_front_2(j,i) = .true.
1440 
1441  j=jmax
1442  if ( (maske(j,i)==2_i2b) & ! ocean
1443  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1444  ) & ! floating ice point
1445  flag_calving_front_2(j,i) = .true.
1446 
1447 end do
1448 
1449 do j=1, jmax-1
1450 
1451  i=0
1452  if ( (maske(j,i)==2_i2b) & ! ocean
1453  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1454  ) & ! floating ice point
1455  flag_calving_front_2(j,i) = .true.
1456 
1457  i=imax
1458  if ( (maske(j,i)==2_i2b) & ! ocean
1459  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1460  ) & ! floating ice point
1461  flag_calving_front_2(j,i) = .true.
1462 
1463 end do
1464 
1465 #else
1466 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1467 #endif
1468 
1469 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1470 
1471 #if (GRID==0 || GRID==1)
1472 
1473 do ir=-imax, imax
1474 do jr=-jmax, jmax
1475  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1476  ! distortion due to stereographic projection not accounted for
1477 end do
1478 end do
1479 
1480 #endif
1481 
1482 !-------- Initial velocities --------
1483 
1484 call calc_temp_melt()
1485 
1486 #if (DYNAMICS==1 || DYNAMICS==2)
1487 
1488 call calc_vxy_b_sia(time, z_sl)
1489 call calc_vxy_sia(dzeta_c, dzeta_t)
1490 
1491 #if (MARGIN==3 || DYNAMICS==2)
1492 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1493 #endif
1494 
1495 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1496 
1497 #if (MARGIN==3)
1498 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1499 #endif
1500 
1501 #elif (DYNAMICS==0)
1502 
1503 call calc_vxy_static()
1504 call calc_vz_static()
1505 
1506 #else
1507  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1508 #endif
1509 
1510 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1511 
1512 !-------- Initial enthalpies --------
1513 
1514 #if (CALCMOD==0 || CALCMOD==-1)
1515 
1516 do i=0, imax
1517 do j=0, jmax
1518 
1519  do kc=0, kcmax
1520  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1521  end do
1522 
1523  do kt=0, ktmax
1524  enth_t(kt,j,i) = enth_c(0,j,i)
1525  end do
1526 
1527 end do
1528 end do
1529 
1530 #elif (CALCMOD==1)
1531 
1532 do i=0, imax
1533 do j=0, jmax
1534 
1535  do kc=0, kcmax
1536  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1537  end do
1538 
1539  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1540  do kt=0, ktmax
1541  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1542  end do
1543  else
1544  do kt=0, ktmax
1545  enth_t(kt,j,i) = enth_c(0,j,i)
1546  end do
1547  end if
1548 
1549 end do
1550 end do
1551 
1552 #elif (CALCMOD==2 || CALCMOD==3)
1553 
1554 do i=0, imax
1555 do j=0, jmax
1556 
1557  do kc=0, kcmax
1558  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1559  end do
1560 
1561  do kt=0, ktmax
1562  enth_t(kt,j,i) = enth_c(0,j,i)
1563  end do
1564 
1565 end do
1566 end do
1567 
1568 #else
1569 
1570 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1571 
1572 #endif
1573 
1574 !-------- Initialize time-series files --------
1575 
1576 ! ------ Time-series file for the ice sheet on the whole
1577 
1578 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1579 
1580 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1581 
1582 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1583 
1584 if (forcing_flag == 1) then
1585 
1586  write(12,1102)
1587  write(12,1103)
1588 
1589  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1590  ' V(m^3) V_g(m^3) V_f(m^3)', &
1591  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1592  ' V_sle(m) V_t(m^3)', &
1593  ' A_t(m^2)',/, &
1594  ' H_max(m) H_t_max(m)', &
1595  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1596  1103 format('----------------------------------------------------', &
1597  '---------------------------------------')
1598 
1599 else if (forcing_flag == 2) then
1600 
1601  write(12,1112)
1602  write(12,1113)
1603 
1604  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1605  ' V(m^3) V_g(m^3) V_f(m^3)', &
1606  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1607  ' V_sle(m) V_t(m^3)', &
1608  ' A_t(m^2)',/, &
1609  ' H_max(m) H_t_max(m)', &
1610  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1611  1113 format('----------------------------------------------------', &
1612  '---------------------------------------')
1613 
1614 end if
1615 
1616 ! ------ Time-series file for deep boreholes
1617 
1618 n_core = 0 ! No boreholes defined
1619 
1620 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1621 
1622 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1623 
1624 write(14,'(1x,a)') '---------------------'
1625 write(14,'(1x,a)') 'No boreholes defined.'
1626 write(14,'(1x,a)') '---------------------'
1627 
1628 !-------- Output of the initial state --------
1629 
1630 #if (defined(OUTPUT_INIT))
1631 
1632 #if (OUTPUT_INIT==0)
1633  flag_init_output = .false.
1634 #elif (OUTPUT_INIT==1)
1635  flag_init_output = .true.
1636 #else
1637  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1638 #endif
1639 
1640 #else
1641  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1642 #endif
1643 
1644 #if (OUTPUT==1)
1645 
1646 #if (ERGDAT==0)
1647  flag_3d_output = .false.
1648 #elif (ERGDAT==1)
1649  flag_3d_output = .true.
1650 #else
1651  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1652 #endif
1653 
1654  if (flag_init_output) &
1655  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1656  flag_3d_output, ndat2d, ndat3d)
1657 
1658 #elif (OUTPUT==2)
1659 
1660 if (time_output(1) <= time_init+eps) then
1661 
1662 #if (ERGDAT==0)
1663  flag_3d_output = .false.
1664 #elif (ERGDAT==1)
1665  flag_3d_output = .true.
1666 #else
1667  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1668 #endif
1669 
1670  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1671  flag_3d_output, ndat2d, ndat3d)
1672 
1673 end if
1674 
1675 #elif (OUTPUT==3)
1676 
1677  flag_3d_output = .false.
1678 
1679  if (flag_init_output) &
1680  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1681  flag_3d_output, ndat2d, ndat3d)
1682 
1683 if (time_output(1) <= time_init+eps) then
1684 
1685  flag_3d_output = .true.
1686 
1687  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1688  flag_3d_output, ndat2d, ndat3d)
1689 
1690 end if
1691 
1692 #else
1693  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1694 #endif
1695 
1696 if (flag_init_output) then
1697  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1698  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1699 end if
1700 
1701 end subroutine sico_init
1702 
1703 !-------------------------------------------------------------------------------
1704 !> Definition of the initial surface and bedrock topography
1705 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1706 !! For ice-free initial topography with relaxed lithosphere
1707 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1708 !<------------------------------------------------------------------------------
1709 subroutine topography1(dxi, deta)
1710 
1711 implicit none
1712 
1713 real(dp), intent(out) :: dxi, deta
1714 
1715 ! integer(i4b) :: i, j
1716 ! integer(i4b) :: ios
1717 ! real(dp) :: xi0, eta0
1718 
1719 dxi=0.0_dp; deta=0.0_dp ! dummy values
1720 
1721 write(6,'(/1x,a)') '>>> topography1: ANF_DAT==1 not defined for this domain!'
1722 stop
1723 
1724 end subroutine topography1
1725 
1726 !-------------------------------------------------------------------------------
1727 !> Definition of the initial surface and bedrock topography
1728 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1729 !! For ice-free initial topography with relaxed lithosphere.
1730 !<------------------------------------------------------------------------------
1731 subroutine topography2(dxi, deta)
1732 
1733 #if (GRID==0 || GRID==1)
1735 #endif
1736 
1737  use metric_m
1738  use topograd_m
1739 
1740 implicit none
1741 
1742 real(dp), intent(out) :: dxi, deta
1743 
1744 integer(i4b) :: i, j, n
1745 integer(i4b) :: ios
1746 real(dp) :: xi0, eta0
1747 character :: ch_dummy
1748 
1749 character(len= 8) :: ch_imax
1750 character(len=128) :: fmt4
1751 character(len=256) :: filename_with_path
1752 
1753 write(ch_imax, fmt='(i8)') imax
1754 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1755 
1756 !-------- Read topography --------
1757 
1758 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1759  trim(zl0_file)
1760 
1761 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1762 
1763 if (ios /= 0) stop ' >>> topography2: Error when opening the zl0 file!'
1764 
1765 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1766  trim(mask_present_file)
1767 
1768 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1769 
1770 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
1771 
1772 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1773 
1774 do j=jmax, 0, -1
1775  read(23, fmt=*) (zl0(j,i), i=0,imax)
1776 end do
1777 
1778 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1779 
1780 do j=jmax, 0, -1
1781  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1782 end do
1783 
1784 close(23, status='keep')
1785 close(24, status='keep')
1786 
1787 !-------- Further stuff --------
1788 
1789 dxi = dx *1000.0_dp ! km -> m
1790 deta = dx *1000.0_dp ! km -> m
1791 
1792 xi0 = x0 *1000.0_dp ! km -> m
1793 eta0 = y0 *1000.0_dp ! km -> m
1794 
1795 do i=0, imax
1796 do j=0, jmax
1797 
1798  if (maske(j,i) <= 1) then
1799  maske(j,i) = 1
1800  zs(j,i) = zl0(j,i)
1801  zb(j,i) = zl0(j,i)
1802  zl(j,i) = zl0(j,i)
1803  else ! (maske(j,i) >= 2)
1804  maske(j,i) = 2
1805 #if (MARGIN==1 || MARGIN==2)
1806  zs(j,i) = zl0(j,i)
1807  zb(j,i) = zl0(j,i)
1808 #elif (MARGIN==3)
1809  zs(j,i) = 0.0_dp ! present-day
1810  zb(j,i) = 0.0_dp ! sea level
1811 #endif
1812  zl(j,i) = zl0(j,i)
1813  end if
1814 
1815  xi(i) = xi0 + real(i,dp)*dxi
1816  eta(j) = eta0 + real(j,dp)*deta
1817 
1818  zm(j,i) = zb(j,i)
1819  n_cts(j,i) = -1
1820  kc_cts(j,i) = 0
1821 
1822  h_c(j,i) = 0.0_dp
1823  h_t(j,i) = 0.0_dp
1824 
1825  dzs_dtau(j,i) = 0.0_dp
1826  dzm_dtau(j,i) = 0.0_dp
1827  dzb_dtau(j,i) = 0.0_dp
1828  dzl_dtau(j,i) = 0.0_dp
1829  dh_c_dtau(j,i) = 0.0_dp
1830  dh_t_dtau(j,i) = 0.0_dp
1831 
1832 end do
1833 end do
1834 
1835 maske_old = maske
1836 
1837 !-------- Geographic coordinates, metric tensor,
1838 ! gradients of the topography --------
1839 
1840 do i=0, imax
1841 do j=0, jmax
1842 
1843 #if (GRID==0 || GRID==1) /* Stereographic projection */
1844 
1845  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1846  lambda0, phi0, lambda(j,i), phi(j,i))
1847 
1848 #elif (GRID==2) /* Geographic coordinates */
1849 
1850  lambda(j,i) = xi(i)
1851  phi(j,i) = eta(j)
1852 
1853 #endif
1854 
1855 end do
1856 end do
1857 
1858 call metric()
1859 
1860 #if (TOPOGRAD==0)
1861 call topograd_1(dxi, deta, 1)
1862 #elif (TOPOGRAD==1)
1863 call topograd_2(dxi, deta, 1)
1864 #endif
1865 
1866 !-------- Corresponding area of grid points --------
1867 
1868 do i=0, imax
1869 do j=0, jmax
1870  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1871 end do
1872 end do
1873 
1874 end subroutine topography2
1875 
1876 !-------------------------------------------------------------------------------
1877 !> Definition of the initial surface and bedrock topography
1878 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1879 !! For initial topography from previous simulation.
1880 !<------------------------------------------------------------------------------
1881 subroutine topography3(dxi, deta, z_sl, anfdatname)
1882 
1883  use read_m, only : read_erg_nc
1884 
1885 #if (GRID==0 || GRID==1)
1887 #endif
1888 
1889  use metric_m
1890  use topograd_m
1891 
1892 implicit none
1893 
1894 character(len=100), intent(in) :: anfdatname
1895 
1896 real(dp), intent(out) :: dxi, deta, z_sl
1897 
1898 integer(i4b) :: i, j, n
1899 integer(i4b) :: ios
1900 character(len=256) :: filename_with_path
1901 character :: ch_dummy
1902 
1903 !-------- Read data from time-slice file of previous simulation --------
1904 
1905 call read_erg_nc(z_sl, anfdatname)
1906 
1907 !-------- Read topography of the relaxed bedrock --------
1908 
1909 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1910  trim(zl0_file)
1911 
1912 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1913 
1914 if (ios /= 0) stop ' >>> topography3: Error when opening the zl0 file!'
1915 
1916 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1917 
1918 do j=jmax, 0, -1
1919  read(23, fmt=*) (zl0(j,i), i=0,imax)
1920 end do
1921 
1922 close(23, status='keep')
1923 
1924 !-------- Further stuff --------
1925 
1926 dxi = dx *1000.0_dp ! km -> m
1927 deta = dx *1000.0_dp ! km -> m
1928 
1929 !-------- Geographic coordinates, metric tensor,
1930 ! gradients of the topography --------
1931 
1932 do i=0, imax
1933 do j=0, jmax
1934 
1935 #if (GRID==0 || GRID==1) /* Stereographic projection */
1936 
1937  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1938  lambda0, phi0, lambda(j,i), phi(j,i))
1939 
1940 #elif (GRID==2) /* Geographic coordinates */
1941 
1942  lambda(j,i) = xi(i)
1943  phi(j,i) = eta(j)
1944 
1945 #endif
1946 
1947 end do
1948 end do
1949 
1950 call metric()
1951 
1952 #if (TOPOGRAD==0)
1953 call topograd_1(dxi, deta, 1)
1954 #elif (TOPOGRAD==1)
1955 call topograd_2(dxi, deta, 1)
1956 #endif
1957 
1958 !-------- Corresponding area of grid points --------
1959 
1960 do i=0, imax
1961 do j=0, jmax
1962  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1963 end do
1964 end do
1965 
1966 end subroutine topography3
1967 
1968 !-------------------------------------------------------------------------------
1969 
1970 end module sico_init_m
1971 !
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), 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...