SICOPOLIS V5-dev  Revision 1173
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, 40.0_dp, eps_sp_dp)) then
230 
231  if ((imax /= 150).or.(jmax /= 70)) then
232  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
233  end if
234 
235 else if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
236 
237  if ((imax /= 300).or.(jmax /= 140)) then
238  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
239  end if
240 
241 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
242 
243  if ((imax /= 600).or.(jmax /= 280)) 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 Scandinavia 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 ' >>> 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 Scandinavia 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 Scandinavia 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 Scandinavia 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 Scandinavia 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 scand 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 smb_corr = 0.0_dp
1325 
1326 q_bm = 0.0_dp
1327 q_tld = 0.0_dp
1328 q_b_tot = 0.0_dp
1329 
1330 p_b_w = 0.0_dp
1331 h_w = 0.0_dp
1332 
1333 call init_temp_water_age_2()
1334 
1335 #if (ENHMOD==1)
1336  call calc_enhance_1()
1337 #elif (ENHMOD==2)
1338  call calc_enhance_2()
1339 #elif (ENHMOD==3)
1340  call calc_enhance_3(time_init)
1341 #elif (ENHMOD==4)
1342  call calc_enhance_4()
1343 #elif (ENHMOD==5)
1344  call calc_enhance_5()
1345 #else
1346  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1347 #endif
1348 
1349 ! ------ Read initial state from output of previous simulation
1350 
1351 #elif (ANF_DAT==3)
1352 
1353 call topography3(dxi, deta, z_sl, anfdatname)
1354 
1355 call boundary(time_init, dtime, dxi, deta, &
1356  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1357 
1358 where ((maske==0_i2b).or.(maske==3_i2b))
1359  ! grounded or floating ice
1361 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1362  as_perp_apl = 0.0_dp
1363 end where
1364 
1365 smb_corr = 0.0_dp
1366 
1367 q_b_tot = q_bm + q_tld
1368 
1369 #endif
1370 
1371 !-------- Inner-point flag --------
1372 
1373 flag_inner_point = .true.
1374 
1375 flag_inner_point(0,:) = .false.
1376 flag_inner_point(jmax,:) = .false.
1377 
1378 flag_inner_point(:,0) = .false.
1379 flag_inner_point(:,imax) = .false.
1380 
1381 !-------- Grounding line and calving front flags --------
1382 
1383 flag_grounding_line_1 = .false.
1384 flag_grounding_line_2 = .false.
1385 
1386 flag_calving_front_1 = .false.
1387 flag_calving_front_2 = .false.
1388 
1389 #if (MARGIN==1 || MARGIN==2)
1390 
1391 continue
1392 
1393 #elif (MARGIN==3)
1394 
1395 do i=1, imax-1
1396 do j=1, jmax-1
1397 
1398  if ( (maske(j,i)==0_i2b) & ! grounded ice
1399  .and. &
1400  ( (maske(j,i+1)==3_i2b) & ! with
1401  .or.(maske(j,i-1)==3_i2b) & ! one
1402  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1403  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1404  ) &
1405  flag_grounding_line_1(j,i) = .true.
1406 
1407  if ( (maske(j,i)==3_i2b) & ! floating ice
1408  .and. &
1409  ( (maske(j,i+1)==0_i2b) & ! with
1410  .or.(maske(j,i-1)==0_i2b) & ! one
1411  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1412  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1413  ) &
1414  flag_grounding_line_2(j,i) = .true.
1415 
1416  if ( (maske(j,i)==3_i2b) & ! floating ice
1417  .and. &
1418  ( (maske(j,i+1)==2_i2b) & ! with
1419  .or.(maske(j,i-1)==2_i2b) & ! one
1420  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1421  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1422  ) &
1423  flag_calving_front_1(j,i) = .true.
1424 
1425  if ( (maske(j,i)==2_i2b) & ! ocean
1426  .and. &
1427  ( (maske(j,i+1)==3_i2b) & ! with
1428  .or.(maske(j,i-1)==3_i2b) & ! one
1429  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1430  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1431  ) &
1432  flag_calving_front_2(j,i) = .true.
1433 
1434 end do
1435 end do
1436 
1437 do i=1, imax-1
1438 
1439  j=0
1440  if ( (maske(j,i)==2_i2b) & ! ocean
1441  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1442  ) & ! floating ice point
1443  flag_calving_front_2(j,i) = .true.
1444 
1445  j=jmax
1446  if ( (maske(j,i)==2_i2b) & ! ocean
1447  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1448  ) & ! floating ice point
1449  flag_calving_front_2(j,i) = .true.
1450 
1451 end do
1452 
1453 do j=1, jmax-1
1454 
1455  i=0
1456  if ( (maske(j,i)==2_i2b) & ! ocean
1457  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1458  ) & ! floating ice point
1459  flag_calving_front_2(j,i) = .true.
1460 
1461  i=imax
1462  if ( (maske(j,i)==2_i2b) & ! ocean
1463  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1464  ) & ! floating ice point
1465  flag_calving_front_2(j,i) = .true.
1466 
1467 end do
1468 
1469 #else
1470 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1471 #endif
1472 
1473 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1474 
1475 #if (GRID==0 || GRID==1)
1476 
1477 do ir=-imax, imax
1478 do jr=-jmax, jmax
1479  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1480  ! distortion due to stereographic projection not accounted for
1481 end do
1482 end do
1483 
1484 #endif
1485 
1486 !-------- Initial velocities --------
1487 
1488 call calc_temp_melt()
1489 
1490 #if (DYNAMICS==1 || DYNAMICS==2)
1491 
1492 call calc_vxy_b_sia(time, z_sl)
1493 call calc_vxy_sia(dzeta_c, dzeta_t)
1494 
1495 #if (MARGIN==3 || DYNAMICS==2)
1496 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1497 #endif
1498 
1499 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1500 
1501 #if (MARGIN==3)
1502 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1503 #endif
1504 
1505 #elif (DYNAMICS==0)
1506 
1507 call calc_vxy_static()
1508 call calc_vz_static()
1509 
1510 #else
1511  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1512 #endif
1513 
1514 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1515 
1516 !-------- Initial enthalpies --------
1517 
1518 #if (CALCMOD==0 || CALCMOD==-1)
1519 
1520 do i=0, imax
1521 do j=0, jmax
1522 
1523  do kc=0, kcmax
1524  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1525  end do
1526 
1527  do kt=0, ktmax
1528  enth_t(kt,j,i) = enth_c(0,j,i)
1529  end do
1530 
1531 end do
1532 end do
1533 
1534 #elif (CALCMOD==1)
1535 
1536 do i=0, imax
1537 do j=0, jmax
1538 
1539  do kc=0, kcmax
1540  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1541  end do
1542 
1543  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1544  do kt=0, ktmax
1545  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1546  end do
1547  else
1548  do kt=0, ktmax
1549  enth_t(kt,j,i) = enth_c(0,j,i)
1550  end do
1551  end if
1552 
1553 end do
1554 end do
1555 
1556 #elif (CALCMOD==2 || CALCMOD==3)
1557 
1558 do i=0, imax
1559 do j=0, jmax
1560 
1561  do kc=0, kcmax
1562  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1563  end do
1564 
1565  do kt=0, ktmax
1566  enth_t(kt,j,i) = enth_c(0,j,i)
1567  end do
1568 
1569 end do
1570 end do
1571 
1572 #else
1573 
1574 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1575 
1576 #endif
1577 
1578 !-------- Initialize time-series files --------
1579 
1580 ! ------ Time-series file for the ice sheet on the whole
1581 
1582 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1583 
1584 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1585 
1586 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1587 
1588 if (forcing_flag == 1) then
1589 
1590  write(12,1102)
1591  write(12,1103)
1592 
1593  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1594  ' V(m^3) V_g(m^3) V_f(m^3)', &
1595  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1596  ' V_sle(m) V_t(m^3)', &
1597  ' A_t(m^2)',/, &
1598  ' H_max(m) H_t_max(m)', &
1599  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1600  1103 format('----------------------------------------------------', &
1601  '---------------------------------------')
1602 
1603 else if (forcing_flag == 2) then
1604 
1605  write(12,1112)
1606  write(12,1113)
1607 
1608  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1609  ' V(m^3) V_g(m^3) V_f(m^3)', &
1610  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1611  ' V_sle(m) V_t(m^3)', &
1612  ' A_t(m^2)',/, &
1613  ' H_max(m) H_t_max(m)', &
1614  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1615  1113 format('----------------------------------------------------', &
1616  '---------------------------------------')
1617 
1618 end if
1619 
1620 ! ------ Time-series file for deep boreholes
1621 
1622 n_core = 0 ! No boreholes defined
1623 
1624 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1625 
1626 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1627 
1628 write(14,'(1x,a)') '---------------------'
1629 write(14,'(1x,a)') 'No boreholes defined.'
1630 write(14,'(1x,a)') '---------------------'
1631 
1632 !-------- Output of the initial state --------
1633 
1634 #if (defined(OUTPUT_INIT))
1635 
1636 #if (OUTPUT_INIT==0)
1637  flag_init_output = .false.
1638 #elif (OUTPUT_INIT==1)
1639  flag_init_output = .true.
1640 #else
1641  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1642 #endif
1643 
1644 #else
1645  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1646 #endif
1647 
1648 #if (OUTPUT==1)
1649 
1650 #if (ERGDAT==0)
1651  flag_3d_output = .false.
1652 #elif (ERGDAT==1)
1653  flag_3d_output = .true.
1654 #else
1655  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1656 #endif
1657 
1658  if (flag_init_output) &
1659  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1660  flag_3d_output, ndat2d, ndat3d)
1661 
1662 #elif (OUTPUT==2)
1663 
1664 if (time_output(1) <= time_init+eps) then
1665 
1666 #if (ERGDAT==0)
1667  flag_3d_output = .false.
1668 #elif (ERGDAT==1)
1669  flag_3d_output = .true.
1670 #else
1671  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1672 #endif
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 #elif (OUTPUT==3)
1680 
1681  flag_3d_output = .false.
1682 
1683  if (flag_init_output) &
1684  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1685  flag_3d_output, ndat2d, ndat3d)
1686 
1687 if (time_output(1) <= time_init+eps) then
1688 
1689  flag_3d_output = .true.
1690 
1691  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1692  flag_3d_output, ndat2d, ndat3d)
1693 
1694 end if
1695 
1696 #else
1697  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1698 #endif
1699 
1700 if (flag_init_output) then
1701  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1702  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1703 end if
1704 
1705 end subroutine sico_init
1706 
1707 !-------------------------------------------------------------------------------
1708 !> Definition of the initial surface and bedrock topography
1709 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1710 !! For ice-free initial topography with relaxed lithosphere
1711 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1712 !<------------------------------------------------------------------------------
1713 subroutine topography1(dxi, deta)
1714 
1715 implicit none
1716 
1717 real(dp), intent(out) :: dxi, deta
1718 
1719 ! integer(i4b) :: i, j
1720 ! integer(i4b) :: ios
1721 ! real(dp) :: xi0, eta0
1722 
1723 dxi=0.0_dp; deta=0.0_dp ! dummy values
1724 
1725 write(6,'(/1x,a)') '>>> topography1: ANF_DAT==1 not defined for this domain!'
1726 stop
1727 
1728 end subroutine topography1
1729 
1730 !-------------------------------------------------------------------------------
1731 !> Definition of the initial surface and bedrock topography
1732 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1733 !! For ice-free initial topography with relaxed lithosphere.
1734 !<------------------------------------------------------------------------------
1735 subroutine topography2(dxi, deta)
1736 
1737 #if (GRID==0 || GRID==1)
1739 #endif
1740 
1741  use metric_m
1742  use topograd_m
1743 
1744 implicit none
1745 
1746 real(dp), intent(out) :: dxi, deta
1747 
1748 integer(i4b) :: i, j, n
1749 integer(i4b) :: ios
1750 real(dp) :: xi0, eta0
1751 character :: ch_dummy
1752 
1753 character(len= 8) :: ch_imax
1754 character(len=128) :: fmt4
1755 character(len=256) :: filename_with_path
1756 
1757 write(ch_imax, fmt='(i8)') imax
1758 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1759 
1760 !-------- Read topography --------
1761 
1762 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1763  trim(zl0_file)
1764 
1765 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1766 
1767 if (ios /= 0) stop ' >>> topography2: Error when opening the zl0 file!'
1768 
1769 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1770 
1771 do j=jmax, 0, -1
1772  read(23, fmt=*) (zl0(j,i), i=0,imax)
1773 end do
1774 
1775 close(23, status='keep')
1776 
1777 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1778  trim(mask_present_file)
1779 
1780 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1781 
1782 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
1783 
1784 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1785 
1786 do j=jmax, 0, -1
1787  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1788 end do
1789 
1790 close(24, status='keep')
1791 
1792 !-------- Further stuff --------
1793 
1794 dxi = dx *1000.0_dp ! km -> m
1795 deta = dx *1000.0_dp ! km -> m
1796 
1797 xi0 = x0 *1000.0_dp ! km -> m
1798 eta0 = y0 *1000.0_dp ! km -> m
1799 
1800 do i=0, imax
1801 do j=0, jmax
1802 
1803  if (maske(j,i) <= 1) then
1804  maske(j,i) = 1
1805  zs(j,i) = zl0(j,i)
1806  zb(j,i) = zl0(j,i)
1807  zl(j,i) = zl0(j,i)
1808  else ! (maske(j,i) >= 2)
1809  maske(j,i) = 2
1810 #if (MARGIN==1 || MARGIN==2)
1811  zs(j,i) = zl0(j,i)
1812  zb(j,i) = zl0(j,i)
1813 #elif (MARGIN==3)
1814  zs(j,i) = 0.0_dp ! present-day
1815  zb(j,i) = 0.0_dp ! sea level
1816 #endif
1817  zl(j,i) = zl0(j,i)
1818  end if
1819 
1820  xi(i) = xi0 + real(i,dp)*dxi
1821  eta(j) = eta0 + real(j,dp)*deta
1822 
1823  zm(j,i) = zb(j,i)
1824  n_cts(j,i) = -1
1825  kc_cts(j,i) = 0
1826 
1827  h_c(j,i) = 0.0_dp
1828  h_t(j,i) = 0.0_dp
1829 
1830  dzs_dtau(j,i) = 0.0_dp
1831  dzm_dtau(j,i) = 0.0_dp
1832  dzb_dtau(j,i) = 0.0_dp
1833  dzl_dtau(j,i) = 0.0_dp
1834  dh_c_dtau(j,i) = 0.0_dp
1835  dh_t_dtau(j,i) = 0.0_dp
1836 
1837 end do
1838 end do
1839 
1840 maske_old = maske
1841 
1842 !-------- Geographic coordinates, metric tensor,
1843 ! gradients of the topography --------
1844 
1845 do i=0, imax
1846 do j=0, jmax
1847 
1848 #if (GRID==0 || GRID==1) /* Stereographic projection */
1849 
1850  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1851  lambda0, phi0, lambda(j,i), phi(j,i))
1852 
1853 #elif (GRID==2) /* Geographic coordinates */
1854 
1855  lambda(j,i) = xi(i)
1856  phi(j,i) = eta(j)
1857 
1858 #endif
1859 
1860 end do
1861 end do
1862 
1863 call metric()
1864 
1865 #if (TOPOGRAD==0)
1866 call topograd_1(dxi, deta, 1)
1867 #elif (TOPOGRAD==1)
1868 call topograd_2(dxi, deta, 1)
1869 #endif
1870 
1871 !-------- Corresponding area of grid points --------
1872 
1873 do i=0, imax
1874 do j=0, jmax
1875  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1876 end do
1877 end do
1878 
1879 end subroutine topography2
1880 
1881 !-------------------------------------------------------------------------------
1882 !> Definition of the initial surface and bedrock topography
1883 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1884 !! For initial topography from previous simulation.
1885 !<------------------------------------------------------------------------------
1886 subroutine topography3(dxi, deta, z_sl, anfdatname)
1887 
1888  use read_m, only : read_erg_nc
1889 
1890 #if (GRID==0 || GRID==1)
1892 #endif
1893 
1894  use metric_m
1895  use topograd_m
1896 
1897 implicit none
1898 
1899 character(len=100), intent(in) :: anfdatname
1900 
1901 real(dp), intent(out) :: dxi, deta, z_sl
1902 
1903 integer(i4b) :: i, j, n
1904 integer(i4b) :: ios
1905 character(len=256) :: filename_with_path
1906 character :: ch_dummy
1907 
1908 !-------- Read data from time-slice file of previous simulation --------
1909 
1910 call read_erg_nc(z_sl, anfdatname)
1911 
1912 !-------- Read topography of the relaxed bedrock --------
1913 
1914 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1915  trim(zl0_file)
1916 
1917 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1918 
1919 if (ios /= 0) stop ' >>> topography3: Error when opening the zl0 file!'
1920 
1921 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1922 
1923 do j=jmax, 0, -1
1924  read(23, fmt=*) (zl0(j,i), i=0,imax)
1925 end do
1926 
1927 close(23, status='keep')
1928 
1929 !-------- Further stuff --------
1930 
1931 dxi = dx *1000.0_dp ! km -> m
1932 deta = dx *1000.0_dp ! km -> m
1933 
1934 !-------- Geographic coordinates, metric tensor,
1935 ! gradients of the topography --------
1936 
1937 do i=0, imax
1938 do j=0, jmax
1939 
1940 #if (GRID==0 || GRID==1) /* Stereographic projection */
1941 
1942  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1943  lambda0, phi0, lambda(j,i), phi(j,i))
1944 
1945 #elif (GRID==2) /* Geographic coordinates */
1946 
1947  lambda(j,i) = xi(i)
1948  phi(j,i) = eta(j)
1949 
1950 #endif
1951 
1952 end do
1953 end do
1954 
1955 call metric()
1956 
1957 #if (TOPOGRAD==0)
1958 call topograd_1(dxi, deta, 1)
1959 #elif (TOPOGRAD==1)
1960 call topograd_2(dxi, deta, 1)
1961 #endif
1962 
1963 !-------- Corresponding area of grid points --------
1964 
1965 do i=0, imax
1966 do j=0, jmax
1967  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1968 end do
1969 end do
1970 
1971 end subroutine topography3
1972 
1973 !-------------------------------------------------------------------------------
1974 
1975 end module sico_init_m
1976 !
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:775
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:823
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:4751
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), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
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:3407
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...