SICOPOLIS V5-dev  Revision 1368
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-2018 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  use error_m
41 
42  implicit none
43 
44  public
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
50 !<------------------------------------------------------------------------------
51 subroutine sico_init(delta_ts, glac_index, &
52  mean_accum, &
53  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
54  time, time_init, time_end, time_output, &
55  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
56  z_sl, dzsl_dtau, z_mar, &
57  ii, jj, nn, &
58  ndat2d, ndat3d, n_output, &
59  runname)
60 
61 #ifndef ALLOW_OPENAD /* Normal */
62  use compare_float_m
63 #endif /* Normal */
64 
68 
69 #if (NETCDF > 1)
70  use netcdf
71  use nc_check_m
72 #endif
73 
74 #if (DISC>0 && !defined(EXEC_MAKE_C_DIS_0))
76 #endif
77 #if (DISC>0 && defined(EXEC_MAKE_C_DIS_0))
79 #endif
80 
81 #if (GRID==0 || GRID==1)
83 #endif
84 
86 
87 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
88  use read_m, only : read_ad_data
89 #endif
90 
91  use boundary_m
93  use calc_enhance_m
95  use calc_vxy_m
96  use calc_vz_m
97  use calc_dxyz_m
99 
100 #ifndef ALLOW_OPENAD /* Normal */
101  use output_m
102 #endif /* Normal */
103 
104 implicit none
105 
106 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
107  jj((imax+1)*(jmax+1)), &
108  nn(0:jmax,0:imax)
109 integer(i4b), intent(out) :: ndat2d, ndat3d
110 integer(i4b), intent(out) :: n_output
111 real(dp), intent(out) :: delta_ts, glac_index
112 real(dp), intent(out) :: mean_accum
113 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
114  dtime_out, dtime_ser
115 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
116 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
117 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
118 
119 character(len=100), intent(out) :: runname
120 
121 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
122 integer(i4b) :: ios, ios1, ios2, ios3, ios4
123 integer(i4b) :: ierr
124 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
125 real(dp) :: time_init0, time_end0, time_output0(100)
126 real(dp) :: d_dummy
127 character(len=100) :: anfdatname, target_topo_dat_name
128 character(len=256) :: filename_with_path
129 character(len=256) :: shell_command
130 character(len=3) :: ch_nc_test
131 character :: ch_dummy
132 logical :: flag_precip_monthly_mean
133 logical :: flag_init_output, flag_3d_output
134 
135 #if (defined(SMB_CORR_FILE))
136 real(sp), dimension(0:IMAX,0:JMAX) :: smb_corr_in_conv
137 #endif
138 
139 #if (defined(INITMIP_SMB_ANOM_FILE))
140 real(dp), dimension(0:IMAX,0:JMAX) :: smb_anom_initmip_conv
141 #endif
142 
143 #if (NETCDF > 1)
144 integer(i4b) :: dimid, ncid, ncv
145 ! dimid: Dimension ID
146 ! ncid: File ID
147 ! ncv: Variable ID
148 #endif
149 
150 real(dp), parameter :: inv_twelve = 1.0_dp/12.0_dp
151 
152 character(len=64), parameter :: thisroutine = 'sico_init'
153 
154 character(len=64), parameter :: fmt1 = '(a)', &
155  fmt2 = '(a,i4)', &
156  fmt2a = '(a,i0)', &
157  fmt3 = '(a,es12.4)'
158 
159 character(len= 8) :: ch_imax
160 character(len=128) :: fmt4
161 
162 write(ch_imax, fmt='(i8)') imax
163 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
164 
165 write(unit=6, fmt='(a)') ' '
166 write(unit=6, fmt='(a)') ' -------- sico_init --------'
167 write(unit=6, fmt='(a)') ' '
168 
169 !-------- Name of the computational domain --------
170 
171 #if (defined(ANT))
172 ch_domain_long = 'Antarctica'
173 ch_domain_short = 'ant'
174 
175 #elif (defined(ASF))
176 ch_domain_long = 'Austfonna'
177 ch_domain_short = 'asf'
178 
179 #elif (defined(EMTP2SGE))
180 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
181 ch_domain_short = 'emtp2sge'
182 
183 #elif (defined(GRL))
184 ch_domain_long = 'Greenland'
185 ch_domain_short = 'grl'
186 
187 #elif (defined(NHEM))
188 ch_domain_long = 'Northern hemisphere'
189 ch_domain_short = 'nhem'
190 
191 #elif (defined(SCAND))
192 ch_domain_long = 'Scandinavia and Eurasia'
193 ch_domain_short = 'scand'
194 
195 #elif (defined(TIBET))
196 ch_domain_long = 'Tibet'
197 ch_domain_short = 'tibet'
198 
199 #elif (defined(NMARS))
200 ch_domain_long = 'North polar cap of Mars'
201 ch_domain_short = 'nmars'
202 
203 #elif (defined(SMARS))
204 ch_domain_long = 'South polar cap of Mars'
205 ch_domain_short = 'smars'
206 
207 #elif (defined(XYZ))
208 ch_domain_long = 'XYZ'
209 ch_domain_short = 'xyz'
210 #if (defined(HEINO))
211 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
212 #endif
213 
214 #else
215 errormsg = ' >>> sico_init: No valid domain specified!'
216 call error(errormsg)
217 #endif
218 
219 !-------- Some initial values --------
220 
221 n_output = 0
222 
223 dtime = 0.0_dp
224 dtime_temp = 0.0_dp
225 dtime_wss = 0.0_dp
226 dtime_out = 0.0_dp
227 dtime_ser = 0.0_dp
228 
229 time = 0.0_dp
230 time_init = 0.0_dp
231 time_end = 0.0_dp
232 time_output = 0.0_dp
233 
234 !-------- Initialisation of the Library of Iterative Solvers Lis,
235 ! if required --------
236 
237 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
238 #ifndef ALLOW_OPENAD /* Normal */
239  call lis_initialize(ierr)
240 #else /* OpenAD */
241  call lis_init_f(ierr)
242 #endif /* Normal vs. OpenAD */
243 #endif
244 
245 !-------- Read physical parameters --------
246 
247 call read_phys_para()
248 
249 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
250 
251 ! ------ Some auxiliary quantities required for the enthalpy method
252 
253 call calc_c_int_table(c, -190, 10, l)
255 
256 !-------- Compatibility check of the SICOPOLIS version with the header file
257 
258 if ( trim(version) /= trim(sico_version) ) then
259  errormsg = ' >>> sico_init: ' &
260  //'SICOPOLIS version not compatible with header file!'
261  call error(errormsg)
262 end if
263 
264 !-------- Check whether the dynamics and thermodynamics modes are defined
265 
266 #if (!defined(DYNAMICS))
267 errormsg = ' >>> sico_init: DYNAMICS not defined in the header file!'
268 call error(errormsg)
269 #endif
270 
271 #if (!defined(CALCMOD))
272 errormsg = ' >>> sico_init: CALCMOD not defined in the header file!'
273 call error(errormsg)
274 #endif
275 
276 #if (defined(ENTHMOD))
277 errormsg = ' >>> sico_init: ENTHMOD must not be defined any more.' &
278  // end_of_line &
279  //' Please update your header file!'
280 call error(errormsg)
281 #endif
282 
283 !-------- Compatibility check of the horizontal resolution with the
284 ! number of grid points --------
285 
286 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
287 
288 #if (GRID==0 || GRID==1)
289 
290 if (approx_equal(dx, 40.0_dp, eps_sp_dp)) then
291 
292  if ( (imax /= 42).or.(jmax /= 72) ) then
293  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
294  call error(errormsg)
295  end if
296 
297 else if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
298 
299  if ( (imax /= 84).or.(jmax /= 144) ) then
300  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
301  call error(errormsg)
302  end if
303 
304 else if (approx_equal(dx, 16.0_dp, eps_sp_dp)) then
305 
306  if ( (imax /= 105).or.(jmax /= 180) ) then
307  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
308  call error(errormsg)
309  end if
310 
311 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
312 
313  if ( (imax /= 168).or.(jmax /= 288) ) then
314  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
315  call error(errormsg)
316  end if
317 
318 else if (approx_equal(dx, 8.0_dp, eps_sp_dp)) then
319 
320  if ( (imax /= 210).or.(jmax /= 360) ) then
321  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
322  call error(errormsg)
323  end if
324 
325 else if (approx_equal(dx, 5.0_dp, eps_sp_dp)) then
326 
327  if ( (imax /= 336).or.(jmax /= 576) ) then
328  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
329  call error(errormsg)
330  end if
331 
332 else if (approx_equal(dx, 4.0_dp, eps_sp_dp)) then
333 
334  if ( (imax /= 420).or.(jmax /= 720) ) then
335  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
336  call error(errormsg)
337  end if
338 
339 else
340  errormsg = ' >>> sico_init: DX wrong!'
341  call error(errormsg)
342 end if
343 
344 #elif (GRID==2)
345 
346 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
347 call error(errormsg)
348 
349 #endif
350 
351 #else /* OpenAD */
352 
353 print *, ' >>> sico_init: grid compatibility check not performed'
354 print *, ' in adjoint applications; check manually.'
355 
356 #endif /* Normal vs. OpenAD */
357 
358 !-------- Compatibility check of the thermodynamics mode
359 ! (cold vs. polythermal vs. enthalpy method)
360 ! and the number of grid points in the lower (kt) ice domain --------
361 
362 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
363 
364 if (ktmax > 2) then
365  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
366  write(6, fmt='(a)') ' the separate kt domain is redundant.'
367  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
368  write(6, fmt='(a)') ' '
369 end if
370 
371 #endif
372 
373 !-------- Compatibility check of surface-temperature
374 ! and precipitation switches --------
375 
376 #if (TSURFACE == 5 && ACCSURFACE != 5)
377 errormsg = ' >>> sico_init: ' &
378  //'Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
379 call error(errormsg)
380 #endif
381 
382 #if (TSURFACE != 5 && ACCSURFACE == 5)
383 errormsg = ' >>> sico_init: ' &
384  //'Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
385 call error(errormsg)
386 #endif
387 
388 #if (TSURFACE == 6)
389 #if (ACCSURFACE != 6)
390 errormsg = ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6' &
391  // end_of_line &
392  //' and NETCDF==2 must be used together!'
393 call error(errormsg)
394 #elif (NETCDF != 2)
395 errormsg = ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6' &
396  // end_of_line &
397  //' and NETCDF==2 must be used together!'
398 call error(errormsg)
399 #endif
400 #endif
401 
402 #if (ACCSURFACE == 6)
403 #if (TSURFACE != 6)
404 errormsg = ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6' &
405  // end_of_line &
406  //' and NETCDF==2 must be used together!'
407 call error(errormsg)
408 #elif (NETCDF != 2)
409 errormsg = ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6' &
410  // end_of_line &
411  //' and NETCDF==2 must be used together!'
412 call error(errormsg)
413 #endif
414 #endif
415 
416 #if (defined(INITMIP_SMB_ANOM_FILE) && NETCDF != 2)
417 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
418  errormsg = ' >>> sico_init: Defining INITMIP_SMB_ANOM_FILE' &
419  // end_of_line &
420  //' must be used together with NETCDF==2!'
421  call error(errormsg)
422 end if
423 #endif
424 
425 !-------- Compatibility check of discretization schemes for the horizontal and
426 ! vertical advection terms in the temperature and age equations --------
427 
428 #if (ADV_HOR==1)
429 errormsg = ' >>> sico_init: ' &
430  //'Option ADV_HOR==1 (central differences) not defined!'
431 call error(errormsg)
432 #endif
433 
434 !-------- Check whether for the shallow shelf
435 ! or shelfy stream approximation
436 ! the chosen grid is Cartesian coordinates
437 ! without distortion correction (GRID==0) --------
438 
439 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
440 #if (GRID != 0)
441 errormsg = ' >>> sico_init: Distortion correction not implemented' &
442  // end_of_line &
443  //' for the shallow shelf approximation (SSA)' &
444  // end_of_line &
445  //' or the shelfy stream approximation (SStA)' &
446  // end_of_line &
447  //' -> GRID==0 required!'
448 call error(errormsg)
449 #endif
450 #endif
451 
452 !-------- Check whether for the hidden ablation scheme
453 ! can be applied --------
454 
455 #if (MARGIN==3 && MB_ACCOUNT==1)
456 errormsg = ' >>> sico_init: MB_ACCOUNT==1 does not yet work for ice shelves!'
457 call error(errormsg)
458 #endif
459 
460 !-------- Setting of forcing flag --------
461 
462 #if (TSURFACE <= 4)
463 
464 forcing_flag = 1 ! forcing by delta_ts
465 
466 #elif (TSURFACE == 5)
467 
468 forcing_flag = 2 ! forcing by glac_index
469 
470 #elif (TSURFACE == 6)
471 
472 forcing_flag = 3 ! forcing by time-dependent surface temperature
473  ! and precipitation data
474 
475 #endif
476 
477 !-------- Initialization of numerical time steps --------
478 
479 dtime0 = dtime0
480 dtime_temp0 = dtime_temp0
481 #if (REBOUND==2)
482 dtime_wss0 = dtime_wss0
483 #endif
484 
485 !-------- Further initializations --------
486 
487 dzeta_c = 1.0_dp/real(kcmax,dp)
488 dzeta_t = 1.0_dp/real(ktmax,dp)
489 dzeta_r = 1.0_dp/real(krmax,dp)
490 
491 ndat2d = 1
492 ndat3d = 1
493 
494 !-------- General abbreviations --------
495 
496 ! ------ kc domain
497 
498 if (deform >= eps) then
499 
500  flag_aa_nonzero = .true. ! non-equidistant grid
501 
502  aa = deform
503  ea = exp(aa)
504 
505  kc=0
506  zeta_c(kc) = 0.0_dp
507  eaz_c(kc) = 1.0_dp
508  eaz_c_quotient(kc) = 0.0_dp
509 
510  do kc=1, kcmax-1
511  zeta_c(kc) = kc*dzeta_c
512  eaz_c(kc) = exp(aa*zeta_c(kc))
513  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
514  end do
515 
516  kc=kcmax
517  zeta_c(kc) = 1.0_dp
518  eaz_c(kc) = exp(aa)
519  eaz_c_quotient(kc) = 1.0_dp
520 
521 else
522 
523  flag_aa_nonzero = .false. ! equidistant grid
524 
525  aa = 0.0_dp
526  ea = 1.0_dp
527 
528  kc=0
529  zeta_c(kc) = 0.0_dp
530  eaz_c(kc) = 1.0_dp
531  eaz_c_quotient(kc) = 0.0_dp
532 
533  do kc=1, kcmax-1
534  zeta_c(kc) = kc*dzeta_c
535  eaz_c(kc) = 1.0_dp
536  eaz_c_quotient(kc) = zeta_c(kc)
537  end do
538 
539  kc=kcmax
540  zeta_c(kc) = 1.0_dp
541  eaz_c(kc) = 1.0_dp
542  eaz_c_quotient(kc) = 1.0_dp
543 
544 end if
545 
546 ! ------ kt domain
547 
548 kt=0
549 zeta_t(kt) = 0.0_dp
550 
551 do kt=1, ktmax-1
552  zeta_t(kt) = kt*dzeta_t
553 end do
554 
555 kt=ktmax
556 zeta_t(kt) = 1.0_dp
557 
558 ! ------ kr domain
559 
560 kr=0
561 zeta_r(kr) = 0.0_dp
562 
563 do kr=1, krmax-1
564  zeta_r(kr) = kr*dzeta_r
565 end do
566 
567 kr=krmax
568 zeta_r(kr) = 1.0_dp
569 
570 !-------- Construction of a vector (with index n) from a 2-d array
571 ! (with indices i, j) by diagonal numbering --------
572 
573 n=1
574 
575 do m=0, imax+jmax
576  do i=m, 0, -1
577  j = m-i
578  if ((i <= imax).and.(j <= jmax)) then
579  ii(n) = i
580  jj(n) = j
581  nn(j,i) = n
582  n=n+1
583  end if
584  end do
585 end do
586 
587 !-------- Specification of current simulation --------
588 
589 runname = runname
590 anfdatname = anfdatname
591 
592 #if (defined(YEAR_ZERO))
594 #else
595 year_zero = 2000.0_dp ! default value 2000 CE
596 #endif
597 
598 time_init0 = time_init0
599 time_end0 = time_end0
600 dtime_ser0 = dtime_ser0
601 
602 #if (OUTPUT==1 || OUTPUT==3)
603 dtime_out0 = dtime_out0
604 #endif
605 
606 #if (OUTPUT==2 || OUTPUT==3)
607 n_output = n_output
608 time_output0( 1) = time_out0_01
609 time_output0( 2) = time_out0_02
610 time_output0( 3) = time_out0_03
611 time_output0( 4) = time_out0_04
612 time_output0( 5) = time_out0_05
613 time_output0( 6) = time_out0_06
614 time_output0( 7) = time_out0_07
615 time_output0( 8) = time_out0_08
616 time_output0( 9) = time_out0_09
617 time_output0(10) = time_out0_10
618 time_output0(11) = time_out0_11
619 time_output0(12) = time_out0_12
620 time_output0(13) = time_out0_13
621 time_output0(14) = time_out0_14
622 time_output0(15) = time_out0_15
623 time_output0(16) = time_out0_16
624 time_output0(17) = time_out0_17
625 time_output0(18) = time_out0_18
626 time_output0(19) = time_out0_19
627 time_output0(20) = time_out0_20
628 #endif
629 
630 !-------- Write log file --------
631 
632 shell_command = 'if [ ! -d'
633 shell_command = trim(shell_command)//' '//outpath
634 shell_command = trim(shell_command)//' '//'] ; then mkdir'
635 shell_command = trim(shell_command)//' '//outpath
636 shell_command = trim(shell_command)//' '//'; fi'
637 call system(trim(shell_command))
638  ! Check whether directory OUTPATH exists. If not, it is created.
639 
640 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
641 
642 #if ( !defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD) ) /* Normal */
643 open(10, iostat=ios, file=trim(filename_with_path), status='new')
644 #else /* OpenAD */
645 open(10, iostat=ios, file=trim(filename_with_path))
646 #endif /* Normal vs. OpenAD */
647 
648 if (ios /= 0) then
649  errormsg = ' >>> sico_init: Error when opening the log file!'
650  call error(errormsg)
651 end if
652 
653 write(10, fmt=trim(fmt1)) 'Computational domain:'
654 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
655 write(10, fmt=trim(fmt1)) ' '
656 
657 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
658 write(10, fmt=trim(fmt1)) ' '
659 
660 write(10, fmt=trim(fmt2)) 'imax =', imax
661 write(10, fmt=trim(fmt2)) 'jmax =', jmax
662 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
663 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
664 write(10, fmt=trim(fmt2)) 'krmax =', krmax
665 write(10, fmt=trim(fmt1)) ' '
666 
667 write(10, fmt=trim(fmt3)) 'a =', aa
668 write(10, fmt=trim(fmt1)) ' '
669 
670 #if (GRID==0 || GRID==1)
671 write(10, fmt=trim(fmt3)) 'x0 =', x0
672 write(10, fmt=trim(fmt3)) 'y0 =', y0
673 write(10, fmt=trim(fmt3)) 'dx =', dx
674 #elif GRID==2
675 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
676 call error(errormsg)
677 #endif
678 write(10, fmt=trim(fmt1)) ' '
679 
680 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
681 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
682 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
683 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
684 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
685 #if (REBOUND==2)
686 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
687 #endif
688 #if (DISC>0)
689 write(10, fmt=trim(fmt3)) 'dtime_mar_coa =', dtime_mar_coa0
690 #endif
691 write(10, fmt=trim(fmt1)) ' '
692 
693 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
694 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
695 #if (ANF_DAT==1)
696 #if (defined(ZB_PRESENT_FILE))
697 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
698 #endif
699 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
700 #endif
701 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
702 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
703 #if (ANF_DAT==1 && defined(TEMP_INIT))
704 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
705 #endif
706 #if (ANF_DAT==3)
707 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
708 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
709 #endif
710 write(10, fmt=trim(fmt1)) ' '
711 
712 #if (THK_EVOL==2)
713 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
714 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
715 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
716 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
717 write(10, fmt=trim(fmt1)) ' '
718 #endif
719 
720 #if (THK_EVOL==3)
721 write(10, fmt=trim(fmt3)) 'target_topo_tau =', target_topo_tau0
722 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
723 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
724 write(10, fmt=trim(fmt1)) ' '
725 #endif
726 
727 #if (THK_EVOL==4)
728 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
729 write(10, fmt=trim(fmt1)) ' '
730 #endif
731 
732 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
733 write(10, fmt=trim(fmt1)) ' '
734 
735 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
736 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
737 #if (CALCTHK==2 || CALCTHK==5)
738 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
739 #if (ITER_MAX_SOR>0)
740 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
741 #endif
742 #endif
743 write(10, fmt=trim(fmt1)) ' '
744 #endif
745 
746 write(10, fmt=trim(fmt2a)) 'TEMP_PRESENT_PARA = ', temp_present_para
747 #if (defined(TEMP_PRESENT_OFFSET))
748 write(10, fmt=trim(fmt3)) 'TEMP_PRESENT_OFFSET =', temp_present_offset
749 #endif
750 write(10, fmt=trim(fmt1)) ' '
751 
752 write(10, fmt=trim(fmt2a)) 'TSURFACE = ', tsurface
753 #if (TSURFACE==1)
754 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
755 #elif (TSURFACE==3)
756 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
757 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
758 #elif (TSURFACE==4)
759 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
760 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
761 #elif (TSURFACE==5)
762 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
763 write(10, fmt=trim(fmt1)) 'temp_ma_anom file = '//temp_ma_anom_file
764 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
765 write(10, fmt=trim(fmt1)) 'temp_mj_anom file = '//temp_mj_anom_file
766 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
767 #elif (TSURFACE==6)
768 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
769 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
770 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
771 #endif
772 
773 #if (!defined(PRECIP_PRESENT_FILE))
774 errormsg = ' >>> sico_init: PRECIP_PRESENT_FILE not defined in the header file!'
775 call error(errormsg)
776 #endif
777 #if (!defined(PRECIP_MA_PRESENT_FILE))
778 errormsg = ' >>> sico_init: ' &
779  //'PRECIP_MA_PRESENT_FILE not defined in the header file!'
780 call error(errormsg)
781 #endif
782 #if (defined(PRECIP_PRESENT_FILE) && defined(PRECIP_MA_PRESENT_FILE))
783 if ( trim(adjustl(precip_present_file)) /= 'none' ) then
784  write(10, fmt=trim(fmt1)) 'precip_present_file = '//precip_present_file
785  flag_precip_monthly_mean = .true.
786 else if ( trim(adjustl(precip_ma_present_file)) /= 'none' ) then
787  write(10, fmt=trim(fmt1)) 'precip_ma_present_file = '//precip_ma_present_file
788  flag_precip_monthly_mean = .false.
789 else
790  errormsg = ' >>> sico_init: Neither PRECIP_PRESENT_FILE' &
791  // end_of_line &
792  //' nor PRECIP_MA_PRESENT_FILE' &
793  // end_of_line &
794  //' specified in the header file!'
795  call error(errormsg)
796 end if
797 #endif
798 #if (!defined(PRECIP_ZS_REF_FILE))
799 errormsg = ' >>> sico_init: PRECIP_ZS_REF_FILE not defined in the header file!'
800 call error(errormsg)
801 #else
802 write(10, fmt=trim(fmt1)) 'precip_zs_ref_file = '//precip_zs_ref_file
803 #endif
804 #if (ACCSURFACE==1)
805 write(10, fmt=trim(fmt3)) 'accfact =', accfact
806 #elif (ACCSURFACE==2 || ACCSURFACE==3)
807 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
808 #elif (ACCSURFACE==5)
809 write(10, fmt=trim(fmt1)) 'precip_anom file = '//precip_anom_file
810 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
811 #elif (ACCSURFACE==6)
812 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
813 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
814 #endif
815 #if (ACCSURFACE <= 3)
816 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
817 #if (ELEV_DESERT == 1)
818 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
819 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
820 #endif
821 #endif
822 
823 #if (ABLSURFACE==3)
824 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
825 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
826 #endif
827 
828 write(10, fmt=trim(fmt1)) ' '
829 
830 #if (defined(SMB_CORR_FILE))
831 if ( trim(adjustl(smb_corr_file)) /= 'none' ) then
832  write(10, fmt=trim(fmt1)) 'smb_corr_file = '//smb_corr_file
833  write(10, fmt=trim(fmt1)) ' '
834 end if
835 #endif
836 
837 #if (defined(INITMIP_CONST_SMB))
838 write(10, fmt=trim(fmt2a)) 'INITMIP_CONST_SMB = ', initmip_const_smb
839 #endif
840 #if (defined(INITMIP_SMB_ANOM_FILE))
841 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
842  write(10, fmt=trim(fmt1)) 'initmip_smb_anom file = '//initmip_smb_anom_file
843 end if
844 #endif
845 #if (defined(INITMIP_CONST_SMB) || defined(INITMIP_SMB_ANOM_FILE))
846 write(10, fmt=trim(fmt1)) ' '
847 #endif
848 
849 write(10, fmt=trim(fmt2a)) 'DISC = ', disc
850 #if (DISC>0)
851 write(10, fmt=trim(fmt3)) 'c_dis_0 =', c_dis_0
852 write(10, fmt=trim(fmt3)) 'c_dis_fac =', c_dis_fac
853 write(10, fmt=trim(fmt3)) 'm_H =', m_h
854 write(10, fmt=trim(fmt3)) 'm_D =', m_d
855 write(10, fmt=trim(fmt3)) 'r_mar_eff =', r_mar_eff
856 #if (defined(S_DIS))
857 write(10, fmt=trim(fmt3)) 's_dis =', s_dis
858 #endif
859 #if (defined(ALPHA_SUB))
860 write(10, fmt=trim(fmt3)) 'alpha_sub =', alpha_sub
861 #endif
862 #if (defined(ALPHA_O))
863 write(10, fmt=trim(fmt3)) 'alpha_o =', alpha_o
864 #endif
865 #endif
866 write(10, fmt=trim(fmt1)) ' '
867 
868 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
869 #if (SEA_LEVEL==1)
870 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
871 #elif (SEA_LEVEL==3)
872 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
873 #endif
874 write(10, fmt=trim(fmt1)) ' '
875 
876 #if (MARGIN==2)
877 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
878 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
879 write(10, fmt=trim(fmt1)) ' '
880 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
881  || marine_ice_calving==6 || marine_ice_calving==7)
882 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
883 write(10, fmt=trim(fmt1)) ' '
884 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
885 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
886 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
887 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
888 write(10, fmt=trim(fmt1)) ' '
889 #endif
890 #elif (MARGIN==3)
891 #if (ICE_SHELF_CALVING==2)
892 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
893 write(10, fmt=trim(fmt1)) ' '
894 #endif
895 #endif
896 
897 #if (defined(BASAL_HYDROLOGY))
898 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
899 #endif
900 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
901 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
902 write(10, fmt=trim(fmt3)) 'smw_coeff =', smw_coeff
903 write(10, fmt=trim(fmt2a)) 'r_smw = ', r_smw
904 write(10, fmt=trim(fmt2a)) 's_smw = ', s_smw
905 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
906 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
907 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
908 #if (defined(TIME_RAMP_UP_SLIDE))
909 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
910 #endif
911 #if (SLIDE_LAW==2 || SLIDE_LAW==3)
912 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
913 #endif
914 #if (BASAL_HYDROLOGY==1 \
915  && defined(hydro_slide_sat_fct) \
916  && defined(c_hw_slide) && defined(hw0_slide))
917 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
918 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
919 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
920 #endif
921 write(10, fmt=trim(fmt2a)) 'ICE_STREAM = ', ice_stream
922 #if (ICE_STREAM==2)
923 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
924 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
925 write(10, fmt=trim(fmt3)) 'smw_coeff_sedi =', smw_coeff_sedi
926 write(10, fmt=trim(fmt2a)) 'r_smw_sedi = ', r_smw_sedi
927 write(10, fmt=trim(fmt2a)) 's_smw_sedi = ', s_smw_sedi
928 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
929 write(10, fmt=trim(fmt2a)) 'p_weert_sedi = ', p_weert_sedi
930 write(10, fmt=trim(fmt2a)) 'q_weert_sedi = ', q_weert_sedi
931 #endif
932 write(10, fmt=trim(fmt1)) ' '
933 
934 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
935 #if (Q_GEO_MOD==1)
936 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
937 #elif (Q_GEO_MOD==2)
938 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
939 #endif
940 write(10, fmt=trim(fmt1)) ' '
941 
942 #if (defined(MARINE_ICE_BASAL_MELTING))
943 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
944 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
945 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
946 #endif
947 write(10, fmt=trim(fmt1)) ' '
948 #endif
949 
950 #if (MARGIN==3)
951 write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
952 #if (FLOATING_ICE_BASAL_MELTING<=3)
953 write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
954 write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
955 #endif
956 write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
957 write(10, fmt=trim(fmt3)) 'z_abyss =', z_abyss
958 #if (FLOATING_ICE_BASAL_MELTING==4)
959 write(10, fmt=trim(fmt3)) 'temp_ocean =', temp_ocean
960 write(10, fmt=trim(fmt3)) 'Omega_qbm =', omega_qbm
961 write(10, fmt=trim(fmt3)) 'alpha_qbm =', alpha_qbm
962 #endif
963 #if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
964 write(10, fmt=trim(fmt3)) 'H_w_0 =', h_w_0
965 #endif
966 write(10, fmt=trim(fmt1)) ' '
967 #endif
968 
969 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
970 #if (REBOUND==1)
971 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
972 #endif
973 #if (REBOUND==1 || REBOUND==2)
974 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
975 #if (TIME_LAG_MOD==1)
976 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
977 #elif (TIME_LAG_MOD==2)
978 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
979 #else
980 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
981 call error(errormsg)
982 #endif
983 #endif
984 #if (REBOUND==2)
985 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
986 #if (FLEX_RIG_MOD==1)
987 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
988 #elif (FLEX_RIG_MOD==2)
989 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
990 #else
991 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
992 call error(errormsg)
993 #endif
994 #endif
995 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
996 write(10, fmt=trim(fmt1)) ' '
997 
998 #if (FLOW_LAW==2)
999 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
1000 write(10, fmt=trim(fmt1)) ' '
1001 #endif
1002 #if (FIN_VISC==2)
1003 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
1004 write(10, fmt=trim(fmt1)) ' '
1005 #endif
1006 
1007 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
1008 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
1009 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
1010 #endif
1011 #if (ENHMOD==2 || ENHMOD==3)
1012 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
1013 #endif
1014 #if (ENHMOD==2)
1015 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
1016 #endif
1017 #if (ENHMOD==3)
1018 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
1019 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
1020 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
1021 #endif
1022 #if (ENHMOD==4 || ENHMOD==5)
1023 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
1024 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
1025 #endif
1026 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
1027 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
1028 #endif
1029 write(10, fmt=trim(fmt1)) ' '
1030 
1031 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
1032 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
1033 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
1034 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
1035 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
1036 #if (defined(QBM_MIN))
1037 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
1038 #elif (defined(QB_MIN)) /* obsolete */
1039 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
1040 #endif
1041 #if (defined(QBM_MAX))
1042 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
1043 #elif (defined(QB_MAX)) /* obsolete */
1044 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
1045 #endif
1046 write(10, fmt=trim(fmt3)) 'age_min =', age_min
1047 write(10, fmt=trim(fmt3)) 'age_max =', age_max
1048 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
1049 #if (ADV_VERT==1)
1050 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
1051 #endif
1052 write(10, fmt=trim(fmt1)) ' '
1053 
1054 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
1055 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
1056 #if (defined(LIS_OPTS))
1057 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
1058 #endif
1059 #if (defined(N_ITER_SSA))
1060 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
1061 #endif
1062 #if (defined(RELAX_FACT_SSA))
1063 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
1064 #endif
1065 #endif
1066 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
1067 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
1068 #endif
1069 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
1070 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
1071 #endif
1072 write(10, fmt=trim(fmt1)) ' '
1073 
1074 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
1075 #if (CALCMOD==-1 && defined(TEMP_CONST))
1076 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
1077 #endif
1078 #if (CALCMOD==-1 && defined(AGE_CONST))
1079 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
1080 #endif
1081 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
1082 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
1083 #endif
1084 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
1085 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
1086 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
1087 #if (MARGIN==2)
1088 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
1089 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
1090 #elif (MARGIN==3)
1091 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
1092 #endif
1093 #if (defined(THK_EVOL))
1094 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
1095 #else
1096 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
1097 call error(errormsg)
1098 #endif
1099 #if (defined(CALCTHK))
1100 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
1101 #else
1102 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
1103 call error(errormsg)
1104 #endif
1105 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
1106 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
1107 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
1108 write(10, fmt=trim(fmt1)) ' '
1109 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
1110 #if (ACCSURFACE==5)
1111 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
1112 #endif
1113 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
1114 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
1115 #if (defined(MB_ACCOUNT))
1116 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
1117 #endif
1118 write(10, fmt=trim(fmt1)) ' '
1119 
1120 #if (defined(OUT_TIMES))
1121 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
1122 #endif
1123 #if (defined(OUTPUT_INIT))
1124 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
1125 #endif
1126 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
1127 #if (OUTPUT==1 || OUTPUT==3)
1128 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
1129 #endif
1130 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
1131 #if (OUTPUT==1 || OUTPUT==2)
1132 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
1133 #endif
1134 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
1135 #if (OUTPUT==2 || OUTPUT==3)
1136 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
1137 do n=1, n_output
1138  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
1139 end do
1140 #endif
1141 write(10, fmt=trim(fmt1)) ' '
1142 
1143 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
1144 
1145 close(10, status='keep')
1146 
1147 !-------- Conversion of time quantities --------
1148 
1149 year_zero = year_zero*year_sec ! a --> s
1150 time_init = time_init0*year_sec ! a --> s
1151 time_end = time_end0*year_sec ! a --> s
1152 dtime = dtime0*year_sec ! a --> s
1153 dtime_temp = dtime_temp0*year_sec ! a --> s
1154 #if (REBOUND==2)
1155 dtime_wss = dtime_wss0*year_sec ! a --> s
1156 #endif
1157 dtime_ser = dtime_ser0*year_sec ! a --> s
1158 #if (OUTPUT==1 || OUTPUT==3)
1159 dtime_out = dtime_out0*year_sec ! a --> s
1160 #endif
1161 #if (OUTPUT==2 || OUTPUT==3)
1162 do n=1, n_output
1163  time_output(n) = time_output0(n)*year_sec ! a --> s
1164 end do
1165 #endif
1166 
1167 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
1168 
1169 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
1170  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
1171  call error(errormsg)
1172 end if
1173 
1174 #if (REBOUND==2)
1175 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) then
1176  errormsg = ' >>> sico_init: dtime_wss must be a multiple of dtime!'
1177  call error(errormsg)
1178 end if
1179 #endif
1180 
1181 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
1182  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
1183  call error(errormsg)
1184 end if
1185 
1186 #if (OUTPUT==1 || OUTPUT==3)
1187 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
1188  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
1189  call error(errormsg)
1190 end if
1191 #endif
1192 
1193 #else /* OpenAD */
1194 
1195 print *, ' >>> sico_init: not checking that time steps are '
1196 print *, ' multiples of each other in adjoint mode;'
1197 print *, ' check manually.'
1198 
1199 #endif /* Normal vs. OpenAD */
1200 
1201 #if (THK_EVOL==2)
1202 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
1203 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
1204 #endif
1205 
1206 #if (THK_EVOL==3)
1207 target_topo_tau = target_topo_tau0 *year_sec ! a --> s
1208 #endif
1209 
1210 time = time_init
1211 
1212 !-------- Reading of present-day
1213 ! monthly mean or mean annual precipitation rate --------
1214 
1215 #if (GRID==0 || GRID==1)
1216 
1217 #if (defined(PRECIP_PRESENT_FILE) && defined(PRECIP_MA_PRESENT_FILE))
1218 
1219 if (flag_precip_monthly_mean) then
1220  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1221  trim(precip_present_file)
1222 else
1223  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1224  trim(precip_ma_present_file)
1225 end if
1226 
1227 #endif
1228 
1229 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1230 
1231 #elif (GRID==2)
1232 
1233 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1234 call error(errormsg)
1235 
1236 #endif
1237 
1238 if (ios /= 0) then
1239  errormsg = ' >>> sico_init: Error when opening the precip file!'
1240  call error(errormsg)
1241 end if
1242 
1243 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1244 
1245 if (flag_precip_monthly_mean) then
1246 
1247  do n=1, 12 ! month counter
1248  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1249  do j=jmax, 0, -1
1250  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
1251  end do
1252  end do
1253 
1254  precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
1255  ! mm/a water equiv. -> m/s ice equiv.
1256 else
1257 
1258  do j=jmax, 0, -1
1259  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
1260  end do
1261 
1262  precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
1263  ! mm/a water equiv. -> m/s ice equiv.
1264 end if
1265 
1266 close(21, status='keep')
1267 
1268 ! ------ Computation of the still undefined present-day
1269 ! mean annual or monthly mean precipitation rate
1270 
1271 if (flag_precip_monthly_mean) then
1272 
1273  precip_ma_present = 0.0_dp ! initialisation
1274 
1275  do n=1, 12 ! month counter
1276  do i=0, imax
1277  do j=0, jmax
1278  precip_ma_present(j,i) = precip_ma_present(j,i) &
1279  + precip_present(j,i,n)*inv_twelve
1280  end do
1281  end do
1282  end do
1283 
1284 else
1285 
1286  do n=1, 12 ! month counter
1287  do i=0, imax
1288  do j=0, jmax
1289  precip_present(j,i,n) = precip_ma_present(j,i)
1290  ! monthly means assumed to be equal
1291  ! to the mean annual precipitation rate
1292  end do
1293  end do
1294  end do
1295 
1296 end if
1297 
1298 ! ------ Reference topography for present-day precipitation rate
1299 
1300 #if (GRID==0 || GRID==1)
1301 
1302 #if (defined(PRECIP_ZS_REF_FILE))
1303 
1304 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1305  trim(precip_zs_ref_file)
1306 
1307 #endif
1308 
1309 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1310 
1311 #elif (GRID==2)
1312 
1313 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1314 call error(errormsg)
1315 
1316 #endif
1317 
1318 if (ios /= 0) then
1319  errormsg = ' >>> sico_init: Error when opening the zs_ref file!'
1320  call error(errormsg)
1321 end if
1322 
1323 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1324 
1325 do j=jmax, 0, -1
1326  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1327 end do
1328 
1329 close(21, status='keep')
1330 
1331 do i=0, imax
1332 do j=0, jmax
1333  zs_ref(j,i) = max(zs_ref(j,i), 0.0_dp)
1334  ! resetting negative elevations (bathymetry data)
1335  ! to the present-day sea surface
1336 end do
1337 end do
1338 
1339 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
1340 
1341 #if (ACCSURFACE==5)
1342 
1343 #if (GRID==0 || GRID==1)
1344 
1345 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1346  trim(precip_anom_file)
1347 
1348 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1349 
1350 #endif
1351 
1352 if (ios /= 0) then
1353  errormsg = ' >>> sico_init: Error when opening the precip anomaly file!'
1354  call error(errormsg)
1355 end if
1356 
1357 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1358 
1359 do j=jmax, 0, -1
1360  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
1361 end do
1362 
1363 close(21, status='keep')
1364 
1365 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
1366 
1367 !-------- LGM monthly precipitation-rate anomalies (assumed to be
1368 ! equal to the mean annual precipitation-rate anomaly) --------
1369 
1370 do n=1, 12 ! month counter
1371 do i=0, imax
1372 do j=0, jmax
1373 
1374  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
1375  ! positive values ensured
1376 
1377 #if (PRECIP_ANOM_INTERPOL==1)
1378  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
1379 #elif (PRECIP_ANOM_INTERPOL==2)
1380  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1381 #else
1382  errormsg = ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1383  call error(errormsg)
1384 #endif
1385 
1386 end do
1387 end do
1388 end do
1389 
1390 #endif
1391 
1392 !-------- Mean accumulation --------
1393 
1394 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1395  ! mm/a water equiv. -> m/s ice equiv.
1396 
1397 !-------- Read rock/sediment mask --------
1398 
1399 #if (ICE_STREAM==2)
1400 
1401 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1402  trim(mask_sedi_file)
1403 
1404 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1405 
1406 if (ios /= 0) then
1407  errormsg = ' >>> sico_init: Error when opening the rock/sediment mask file!'
1408  call error(errormsg)
1409 end if
1410 
1411 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1412 
1413 do j=jmax, 0, -1
1414  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
1415 end do
1416 
1417 close(24, status='keep')
1418 
1419 #endif
1420 
1421 !-------- Reading of the prescribed target topography --------
1422 
1423 #if ((THK_EVOL==2) || (THK_EVOL==3))
1424 
1425 target_topo_dat_name = trim(target_topo_dat_name)
1426 
1427 #if (NETCDF==1)
1428 errormsg = ' >>> sico_init: Reading of target topography' &
1429  // end_of_line &
1430  //' only implemented for NetCDF files (NETCDF==2)!'
1431 call error(errormsg)
1432 #elif (NETCDF==2)
1433 call read_target_topo_nc(target_topo_dat_name)
1434 #else
1435 errormsg = ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1436 call error(errormsg)
1437 #endif
1438 
1439 #endif
1440 
1441 !-------- Reading of the maximum ice extent mask --------
1442 
1443 #if (THK_EVOL==4)
1444 
1445 #if (GRID==0 || GRID==1)
1446 
1447 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1448  trim(mask_maxextent_file)
1449 
1450 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1451 
1452 #elif (GRID==2)
1453 
1454 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1455 call error(errormsg)
1456 
1457 #endif
1458 
1459 if (ios /= 0) then
1460  errormsg = ' >>> sico_init: ' &
1461  //'Error when opening the maximum ice extent mask file!'
1462  call error(errormsg)
1463 end if
1464 
1465 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1466 
1467 do j=jmax, 0, -1
1468  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1469 end do
1470 
1471 close(24, status='keep')
1472 
1473 #endif
1474 
1475 !-------- Reading of LGM mean-annual and mean-July surface-temperature
1476 ! anomalies --------
1477 
1478 #if (TSURFACE==5)
1479 
1480 #if (GRID==0 || GRID==1)
1481 
1482 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1483  trim(temp_ma_anom_file)
1484 
1485 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1486 
1487 #endif
1488 
1489 if (ios /= 0) then
1490  errormsg = ' >>> sico_init: Error when opening the temp_ma anomaly file!'
1491  call error(errormsg)
1492 end if
1493 
1494 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1495 
1496 do j=jmax, 0, -1
1497  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
1498 end do
1499 
1500 close(21, status='keep')
1501 
1502 #if (GRID==0 || GRID==1)
1503 
1504 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1505  trim(temp_mj_anom_file)
1506 
1507 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1508 
1509 #endif
1510 
1511 if (ios /= 0) then
1512  errormsg = ' >>> sico_init: Error when opening the temp_mj anomaly file!'
1513  call error(errormsg)
1514 end if
1515 
1516 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1517 
1518 do j=jmax, 0, -1
1519  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
1520 end do
1521 
1522 close(22, status='keep')
1523 
1524 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
1525 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
1526 
1527 #endif
1528 
1529 !-------- Read data for delta_ts --------
1530 
1531 #if (TSURFACE==4)
1532 
1533 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1534 
1535 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1536 
1537 if (ios /= 0) then
1538  errormsg = ' >>> sico_init: Error when opening the data file for delta_ts!'
1539  call error(errormsg)
1540 end if
1541 
1542 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1543 
1544 if (ch_dummy /= '#') then
1545  errormsg = ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' &
1546  // end_of_line &
1547  //' not defined in data file for delta_ts!'
1548  call error(errormsg)
1549 end if
1550 
1552 
1553 if (allocated(griptemp)) then
1554  deallocate(griptemp)
1555 end if
1556 
1557 allocate(griptemp(0:ndata_grip))
1558 
1559 do n=0, ndata_grip
1560  read(21, fmt=*) d_dummy, griptemp(n)
1561 end do
1562 
1563 close(21, status='keep')
1564 
1565 #endif
1566 
1567 !-------- Read data for the glacial index --------
1568 
1569 #if (TSURFACE==5)
1570 
1571 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1572 
1573 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1574 
1575 if (ios /= 0) then
1576  errormsg = ' >>> sico_init: Error when opening the glacial-index file!'
1577  call error(errormsg)
1578 end if
1579 
1580 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1581 
1582 if (ch_dummy /= '#') then
1583  errormsg = ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max' &
1584  // end_of_line &
1585  //' not defined in glacial-index file!'
1586  call error(errormsg)
1587 end if
1588 
1590 
1591 allocate(glacial_index(0:ndata_gi))
1592 
1593 do n=0, ndata_gi
1594  read(21, fmt=*) d_dummy, glacial_index(n)
1595 end do
1596 
1597 close(21, status='keep')
1598 
1599 #endif
1600 
1601 !-------- Reading of time-dependent surface-temperature
1602 ! and precipitation anomaly data --------
1603 
1604 #if (TSURFACE==6 && ACCSURFACE==6)
1605 
1606 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1607  '/'//trim(temp_precip_anom_file)
1608 
1609 call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
1610  ncid_temp_precip), thisroutine )
1611 
1612 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1613  thisroutine )
1614 
1615 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1616  len = ndata_temp_precip), thisroutine )
1617 
1619 
1621 
1622 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1623 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1624  ! in a
1625 
1628  ! in a; constant time step assumed
1630 
1631 #endif
1632 
1633 !-------- Prescribed surface mass balance correction --------
1634 
1635 smb_corr_in = 0.0_dp
1636 
1637 #if (defined(SMB_CORR_FILE))
1638 
1639 if (trim(adjustl(smb_corr_file)) /= 'none') then
1640 
1641  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1642  trim(smb_corr_file)
1643 
1644  filename_with_path = adjustr(filename_with_path)
1645  n = len(filename_with_path)
1646  ch_nc_test = filename_with_path(n-2:n)
1647  filename_with_path = adjustl(filename_with_path)
1648 
1649  if (ch_nc_test == '.nc') then ! NetCDF file
1650 
1651 #if (NETCDF==1)
1652 
1653  errormsg = ' >>> sico_init: Reading of NetCDF SMB_CORR_FILE' &
1654  // end_of_line &
1655  //' requires NETCDF==2!'
1656  call error(errormsg)
1657 
1658 #elif (NETCDF==2)
1659 
1660  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1661  thisroutine )
1662 
1663  call check( nf90_inq_varid(ncid, 'DSMB', ncv) )
1664  call check( nf90_get_var(ncid, ncv, smb_corr_in_conv) )
1665 
1666  call check( nf90_close(ncid) )
1667 
1668  do i=0, imax
1669  do j=0, jmax
1670  smb_corr_in(j,i) = smb_corr_in_conv(i,j) /year_sec
1671  ! m/a ice equiv. -> m/s ice equiv.
1672  end do
1673  end do
1674 
1675 #else
1676  errormsg = ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1677  call error(errormsg)
1678 #endif
1679 
1680  else ! ASCII file
1681 
1682  open(21, iostat=ios, &
1683  file=trim(filename_with_path), recl=rcl1, status='old')
1684 
1685  if (ios /= 0) then
1686  errormsg = ' >>> sico_init: ' &
1687  //'Error when opening the SMB correction file!'
1688  call error(errormsg)
1689  end if
1690 
1691  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1692 
1693  do j=jmax, 0, -1
1694  read(21, fmt=*) (smb_corr_in(j,i), i=0,imax)
1695  end do
1696 
1697  close(21, status='keep')
1698 
1699  smb_corr_in = smb_corr_in /year_sec
1700  ! m/a ice equiv. -> m/s ice equiv.
1701 
1702  end if
1703 
1704 end if
1705 
1706 #endif
1707 
1708 !-------- Reading of ISMIP6 InitMIP SMB anomaly data --------
1709 
1710 #if (defined(INITMIP_SMB_ANOM_FILE))
1711 
1712 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
1713 
1714  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1715  '/'//trim(initmip_smb_anom_file)
1716 
1717  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1718  thisroutine )
1719 
1720  call check( nf90_inq_varid(ncid, 'DSMB', ncv) )
1721  call check( nf90_get_var(ncid, ncv, smb_anom_initmip_conv) )
1722 
1723  call check( nf90_close(ncid) )
1724 
1725  do i=0, imax
1726  do j=0, jmax
1727  smb_anom_initmip(j,i) = smb_anom_initmip_conv(i,j) /year_sec
1728  ! m/a ice equiv. -> m/s ice equiv.
1729  end do
1730  end do
1731 
1732 else
1733  smb_anom_initmip = 0.0_dp
1734 end if
1735 
1736 #endif
1737 
1738 !-------- Reading of the global annual temperature anomaly
1739 ! (for parameterising the sub-ocean temperature anomaly
1740 ! for the ice discharge parameterisation)
1741 
1742 #if (DISC==2)
1743 
1744 filename_with_path = trim(inpath)//'/general/dTg_paleo.dat'
1745 
1746 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1747 
1748 if (ios /= 0) then
1749  errormsg = ' >>> sico_init: ' &
1750  //'Error when opening the data file for dT_glann_CLIMBER!'
1751  call error(errormsg)
1752 end if
1753 
1754 read(21, fmt=*) ch_dummy, glann_time_min, glann_time_stp, glann_time_max
1755 
1756 if (ch_dummy /= '#') then
1757  errormsg = ' >>> sico_init: glann_time_min, glann_time_stp, glann_time_max' &
1758  // end_of_line &
1759  //' not defined in data file for dT_glann_CLIMBER!'
1760  call error(errormsg)
1761 end if
1762 
1763 ndata_glann = (glann_time_max-glann_time_min)/glann_time_stp
1764 
1765 allocate(dt_glann_climber(0:ndata_glann))
1766 
1767 do n=0, ndata_glann
1768  read(21, fmt=*) d_dummy, dt_glann_climber(n)
1769 end do
1770 
1771 close(21, status='keep')
1772 
1773 #endif
1774 
1775 !-------- Read data for z_sl --------
1776 
1777 #if (SEA_LEVEL==3)
1778 
1779 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1780 
1781 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1782 
1783 if (ios /= 0) then
1784  errormsg = ' >>> sico_init: Error when opening the data file for z_sl!'
1785  call error(errormsg)
1786 end if
1787 
1788 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1789 
1790 if (ch_dummy /= '#') then
1791  errormsg = ' >>> sico_init:' &
1792  // end_of_line &
1793  //' specmap_time_min, specmap_time_stp, specmap_time_max' &
1794  // end_of_line &
1795  //' not defined in data file for z_sl!'
1796  call error(errormsg)
1797 end if
1798 
1800 
1801 allocate(specmap_zsl(0:ndata_specmap))
1802 
1803 do n=0, ndata_specmap
1804  read(21, fmt=*) d_dummy, specmap_zsl(n)
1805 end do
1806 
1807 close(21, status='keep')
1808 
1809 #endif
1810 
1811 !-------- Determination of the geothermal heat flux --------
1812 
1813 #if (Q_GEO_MOD==1)
1814 
1815 ! ------ Constant value
1816 
1817 do i=0, imax
1818 do j=0, jmax
1819  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1820 end do
1821 end do
1822 
1823 #elif (Q_GEO_MOD==2)
1824 
1825 ! ------ Read data from file
1826 
1827 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1828  trim(q_geo_file)
1829 
1830 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1831 
1832 if (ios /= 0) then
1833  errormsg = ' >>> sico_init: Error when opening the qgeo file!'
1834  call error(errormsg)
1835 end if
1836 
1837 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1838 
1839 do j=jmax, 0, -1
1840  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1841 end do
1842 
1843 close(21, status='keep')
1844 
1845 do i=0, imax
1846 do j=0, jmax
1847  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1848 end do
1849 end do
1850 
1851 #endif
1852 
1853 !-------- Reading of tabulated kei function--------
1854 
1855 #if (REBOUND==0 || REBOUND==1)
1856 
1857 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1858  ! dummy values
1859 #elif (REBOUND==2)
1860 
1861 call read_kei()
1862 
1863 #endif
1864 
1865 !-------- Determination of the time lag
1866 ! of the relaxing asthenosphere --------
1867 
1868 #if (REBOUND==1 || REBOUND==2)
1869 
1870 #if (TIME_LAG_MOD==1)
1871 
1872 time_lag_asth = time_lag*year_sec ! a -> s
1873 
1874 #elif (TIME_LAG_MOD==2)
1875 
1876 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1877  trim(time_lag_file)
1878 
1879 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1880 
1881 if (ios /= 0) then
1882  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
1883  call error(errormsg)
1884 end if
1885 
1886 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1887 
1888 do j=jmax, 0, -1
1889  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1890 end do
1891 
1892 close(29, status='keep')
1893 
1894 time_lag_asth = time_lag_asth*year_sec ! a -> s
1895 
1896 #endif
1897 
1898 #elif (REBOUND==0)
1899 
1900 time_lag_asth = 0.0_dp ! dummy values
1901 
1902 #endif
1903 
1904 !-------- Determination of the flexural rigidity of the lithosphere --------
1905 
1906 #if (REBOUND==2)
1907 
1908 #if (FLEX_RIG_MOD==1)
1909 
1910 flex_rig_lith = flex_rig
1911 
1912 #elif (FLEX_RIG_MOD==2)
1913 
1914 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1915  trim(flex_rig_file)
1916 
1917 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1918 
1919 if (ios /= 0) then
1920  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1921  call error(errormsg)
1922 end if
1923 
1924 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1925 
1926 do j=jmax, 0, -1
1927  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1928 end do
1929 
1930 close(29, status='keep')
1931 
1932 #endif
1933 
1934 #elif (REBOUND==0 || REBOUND==1)
1935 
1936 flex_rig_lith = 0.0_dp ! dummy values
1937 
1938 #endif
1939 
1940 !-------- Definition of initial values --------
1941 
1942 ! ------ Present topography
1943 
1944 #if (ANF_DAT==1)
1945 
1946 call topography1(dxi, deta)
1947 
1948 #if (DISC>0) /* Ice discharge parameterization */
1949 call disc_param(dtime)
1950 call disc_fields()
1951 #endif
1952 
1953 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1954 
1955 call boundary(time_init, dtime, dxi, deta, &
1956  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1957 
1958 where ((maske==0_i1b).or.(maske==3_i1b))
1959  ! grounded or floating ice
1961 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1962  as_perp_apl = 0.0_dp
1963 end where
1964 
1965 smb_corr = 0.0_dp
1966 
1967 q_bm = 0.0_dp
1968 q_tld = 0.0_dp
1969 q_b_tot = 0.0_dp
1970 
1971 p_b_w = 0.0_dp
1972 h_w = 0.0_dp
1973 
1974 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1976 #elif (TEMP_INIT==2)
1978 #elif (TEMP_INIT==3)
1980 #elif (TEMP_INIT==4)
1982 #else
1983  errormsg = ' >>> sico_init: TEMP_INIT must be set to either 1, 2, 3 or 4!'
1984  call error(errormsg)
1985 #endif
1986 
1987 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
1988  call read_ad_data()
1989 #endif
1990 
1991 #if (ENHMOD==1)
1992  call calc_enhance_1()
1993 #elif (ENHMOD==2)
1994  call calc_enhance_2()
1995 #elif (ENHMOD==3)
1996  call calc_enhance_3(time_init)
1997 #elif (ENHMOD==4)
1998  call calc_enhance_4()
1999 #elif (ENHMOD==5)
2000  call calc_enhance_5()
2001 #else
2002  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
2003  call error(errormsg)
2004 #endif
2005 
2006 ! ------ Ice-free, relaxed bedrock
2007 
2008 #elif (ANF_DAT==2)
2009 
2010 call topography2(dxi, deta)
2011 
2012 #if (DISC>0) /* Ice discharge parameterization */
2013 call disc_param(dtime)
2014 call disc_fields()
2015 #endif
2016 
2017 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
2018 
2019 call boundary(time_init, dtime, dxi, deta, &
2020  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
2021 
2022 as_perp_apl = 0.0_dp
2023 
2024 smb_corr = 0.0_dp
2025 
2026 q_bm = 0.0_dp
2027 q_tld = 0.0_dp
2028 q_b_tot = 0.0_dp
2029 
2030 p_b_w = 0.0_dp
2031 h_w = 0.0_dp
2032 
2033 call init_temp_water_age_2()
2034 
2035 #if (ENHMOD==1)
2036  call calc_enhance_1()
2037 #elif (ENHMOD==2)
2038  call calc_enhance_2()
2039 #elif (ENHMOD==3)
2040  call calc_enhance_3(time_init)
2041 #elif (ENHMOD==4)
2042  call calc_enhance_4()
2043 #elif (ENHMOD==5)
2044  call calc_enhance_5()
2045 #else
2046  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
2047  call error(errormsg)
2048 #endif
2049 
2050 ! ------ Read initial state from output of previous simulation
2051 
2052 #elif (ANF_DAT==3)
2053 
2054 call topography3(dxi, deta, z_sl, anfdatname)
2055 
2056 #if (DISC>0) /* Ice discharge parameterization */
2057 call disc_param(dtime)
2058 call disc_fields()
2059 #endif
2060 
2061 call boundary(time_init, dtime, dxi, deta, &
2062  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
2063 
2064 where ((maske==0_i1b).or.(maske==3_i1b))
2065  ! grounded or floating ice
2067 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
2068  as_perp_apl = 0.0_dp
2069 end where
2070 
2071 smb_corr = 0.0_dp
2072 
2073 q_b_tot = q_bm + q_tld
2074 
2075 #endif
2076 
2077 !-------- Inner-point flag --------
2078 
2079 flag_inner_point = .true.
2080 
2081 flag_inner_point(0,:) = .false.
2082 flag_inner_point(jmax,:) = .false.
2083 
2084 flag_inner_point(:,0) = .false.
2085 flag_inner_point(:,imax) = .false.
2086 
2087 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
2088 
2089 #if (GRID==0 || GRID==1)
2090 
2091 do ir=-imax, imax
2092 do jr=-jmax, jmax
2093  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
2094  ! distortion due to stereographic projection not accounted for
2095 end do
2096 end do
2097 
2098 #endif
2099 
2100 !-------- Initial velocities --------
2101 
2102 call calc_temp_melt()
2103 call flag_update_gf_gl_cf()
2104 
2105 #if (DYNAMICS==1 || DYNAMICS==2)
2106 
2107 call calc_vxy_b_sia(time, z_sl)
2108 call calc_vxy_sia(dzeta_c, dzeta_t)
2109 
2110 #if (MARGIN==3 || DYNAMICS==2)
2111 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
2112 #endif
2113 
2114 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
2115 
2116 #if (MARGIN==3)
2117 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
2118 #endif
2119 
2120 #elif (DYNAMICS==0)
2121 
2122 call calc_vxy_static()
2123 call calc_vz_static()
2124 
2125 #else
2126 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
2127 call error(errormsg)
2128 #endif
2129 
2130 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
2131 
2132 !-------- Initial enthalpies --------
2133 
2134 #if (CALCMOD==0 || CALCMOD==-1)
2135 
2136 do i=0, imax
2137 do j=0, jmax
2138 
2139  do kc=0, kcmax
2140  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2141  end do
2142 
2143  do kt=0, ktmax
2144  enth_t(kt,j,i) = enth_c(0,j,i)
2145  end do
2146 
2147 end do
2148 end do
2149 
2150 #elif (CALCMOD==1)
2151 
2152 do i=0, imax
2153 do j=0, jmax
2154 
2155  do kc=0, kcmax
2156  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2157  end do
2158 
2159  if ( (maske(j,i) == 0_i1b).and.(n_cts(j,i) == 1_i1b) ) then
2160  do kt=0, ktmax
2161  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
2162  end do
2163  else
2164  do kt=0, ktmax
2165  enth_t(kt,j,i) = enth_c(0,j,i)
2166  end do
2167  end if
2168 
2169 end do
2170 end do
2171 
2172 #elif (CALCMOD==2 || CALCMOD==3)
2173 
2174 do i=0, imax
2175 do j=0, jmax
2176 
2177  do kc=0, kcmax
2178  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
2179  end do
2180 
2181  do kt=0, ktmax
2182  enth_t(kt,j,i) = enth_c(0,j,i)
2183  end do
2184 
2185 end do
2186 end do
2187 
2188 #else
2189 
2190 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
2191 call error(errormsg)
2192 
2193 #endif
2194 
2195 !-------- Initialize time-series files --------
2196 
2197 ! ------ Time-series file for the ice sheet on the whole
2198 
2199 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
2200 
2201 #if ( !defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD) ) /* Normal */
2202 open(12, iostat=ios, file=trim(filename_with_path), status='new')
2203 #else /* OpenAD */
2204 open(12, iostat=ios, file=trim(filename_with_path))
2205 #endif /* Normal vs. OpenAD */
2206 
2207 if (ios /= 0) then
2208  errormsg = ' >>> sico_init: Error when opening the ser file!'
2209  call error(errormsg)
2210 end if
2211 
2212 if (forcing_flag == 1) then
2213 
2214  write(12,1102)
2215  write(12,1103)
2216 
2217  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
2218  ' V(m^3) V_g(m^3) V_f(m^3)', &
2219  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2220  ' V_sle(m) V_t(m^3)', &
2221  ' A_t(m^2)',/, &
2222  ' H_max(m) H_t_max(m)', &
2223  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2224  1103 format('----------------------------------------------------', &
2225  '---------------------------------------')
2226 
2227 else if (forcing_flag == 2) then
2228 
2229  write(12,1112)
2230  write(12,1113)
2231 
2232  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
2233  ' V(m^3) V_g(m^3) V_f(m^3)', &
2234  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2235  ' V_sle(m) V_t(m^3)', &
2236  ' A_t(m^2)',/, &
2237  ' H_max(m) H_t_max(m)', &
2238  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2239  1113 format('----------------------------------------------------', &
2240  '---------------------------------------')
2241 
2242 else if (forcing_flag == 3) then
2243 
2244  write(12,1122)
2245  write(12,1123)
2246 
2247  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2248  ' V(m^3) V_g(m^3) V_f(m^3)', &
2249  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2250  ' V_sle(m) V_t(m^3)', &
2251  ' A_t(m^2)',/, &
2252  ' H_max(m) H_t_max(m)', &
2253  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2254  1123 format('----------------------------------------------------', &
2255  '---------------------------------------')
2256 
2257 end if
2258 
2259 ! ------ Time-series file for deep boreholes
2260 
2261 n_core = 7 ! GRIP, GISP2, Dye3, Camp Century (CC),
2262  ! NorthGRIP (NGRIP), NEEM, EastGRIP (EGRIP)
2263 
2264 allocate(lambda_core(n_core), phi_core(n_core), &
2266 
2267 ch_core(1) = 'GRIP'
2268 phi_core(1) = 72.58722_dp *pi_180 ! Geographical position of GRIP,
2269 lambda_core(1) = -37.64222_dp *pi_180 ! conversion deg -> rad
2270 
2271 ch_core(2) = 'GISP2'
2272 phi_core(2) = 72.58833_dp *pi_180 ! Geographical position of GISP2
2273 lambda_core(2) = -38.45750_dp *pi_180 ! conversion deg -> rad
2274 
2275 ch_core(3) = 'Dye3'
2276 phi_core(3) = 65.15139_dp *pi_180 ! Geographical position of Dye3,
2277 lambda_core(3) = -43.81722_dp *pi_180 ! conversion deg -> rad
2278 
2279 ch_core(4) = 'Camp Century'
2280 phi_core(4) = 77.17970_dp *pi_180 ! Geographical position of CC,
2281 lambda_core(4) = -61.10975_dp *pi_180 ! conversion deg -> rad
2282 
2283 ch_core(5) = 'NGRIP'
2284 phi_core(5) = 75.09694_dp *pi_180 ! Geographical position of NGRIP,
2285 lambda_core(5) = -42.31956_dp *pi_180 ! conversion deg -> rad
2286 
2287 ch_core(6) = 'NEEM'
2288 phi_core(6) = 77.5_dp *pi_180 ! Geographical position of NEEM,
2289 lambda_core(6) = -50.9_dp *pi_180 ! conversion deg -> rad
2290 
2291 ch_core(7) = 'EGRIP'
2292 phi_core(7) = 75.6299_dp *pi_180 ! Geographical position of EGRIP,
2293 lambda_core(7) = -35.9867_dp *pi_180 ! conversion deg -> rad
2294 
2295 #if (GRID==0 || GRID==1) /* Stereographic projection */
2296 
2297 do n=1, n_core
2299  a, b, lambda0, phi0, x_core(n), y_core(n))
2300 end do
2301 
2302 #elif (GRID==2) /* Geographical coordinates */
2303 
2305 y_core = phi_core
2306 
2307 #endif
2308 
2309 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
2310 
2311 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD) ) /* Normal */
2312 open(14, iostat=ios, file=trim(filename_with_path), status='new')
2313 #else /* OpenAD */
2314 open(14, iostat=ios, file=trim(filename_with_path))
2315 #endif /* Normal vs. OpenAD */
2316 
2317 if (ios /= 0) then
2318  errormsg = ' >>> sico_init: Error when opening the core file!'
2319  call error(errormsg)
2320 end if
2321 
2322 if (forcing_flag == 1) then
2323 
2324  write(14,1106)
2325  write(14,1107)
2326 
2327  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
2328  ' H_GR(m) H_G2(m) H_D3(m)', &
2329  ' H_CC(m) H_NG(m) H_NE(m) H_EG(m)',/, &
2330  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
2331  ' v_CC(m/a) v_NG(m/a) v_NE(m/a) v_EG(m/a)',/, &
2332  ' T_GR(C) T_G2(C) T_D3(C)', &
2333  ' T_CC(C) T_NG(C) T_NE(C) T_EG(C)')
2334  1107 format('----------------------------------------------------', &
2335  '----------------------------------------------------')
2336 
2337 else if (forcing_flag == 2) then
2338 
2339  write(14,1116)
2340  write(14,1117)
2341 
2342  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
2343  ' H_GR(m) H_G2(m) H_D3(m)', &
2344  ' H_CC(m) H_NG(m) H_NE(m) H_EG(m)',/, &
2345  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
2346  ' v_CC(m/a) v_NG(m/a) v_NE(m/a) v_EG(m/a)',/, &
2347  ' T_GR(C) T_G2(C) T_D3(C)', &
2348  ' T_CC(C) T_NG(C) T_NE(C) T_EG(C)')
2349  1117 format('----------------------------------------------------', &
2350  '----------------------------------------------------')
2351 
2352 else if (forcing_flag == 3) then
2353 
2354  write(14,1126)
2355  write(14,1127)
2356 
2357  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2358  ' H_GR(m) H_G2(m) H_D3(m)', &
2359  ' H_CC(m) H_NG(m) H_NE(m) H_EG(m)',/, &
2360  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
2361  ' v_CC(m/a) v_NG(m/a) v_NE(m/a) v_EG(m/a)',/, &
2362  ' T_GR(C) T_G2(C) T_D3(C)', &
2363  ' T_CC(C) T_NG(C) T_NE(C) T_EG(C)')
2364  1127 format('----------------------------------------------------', &
2365  '----------------------------------------------------')
2366 
2367 end if
2368 
2369 !-------- Output of the initial state --------
2370 
2371 #ifndef ALLOW_OPENAD /* Normal */
2372 
2373 #if (defined(OUTPUT_INIT))
2374 
2375 #if (OUTPUT_INIT==0)
2376  flag_init_output = .false.
2377 #elif (OUTPUT_INIT==1)
2378  flag_init_output = .true.
2379 #else
2380  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
2381  call error(errormsg)
2382 #endif
2383 
2384 #else
2385  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
2386 #endif
2387 
2388 #if (OUTPUT==1)
2389 
2390 #if (ERGDAT==0)
2391  flag_3d_output = .false.
2392 #elif (ERGDAT==1)
2393  flag_3d_output = .true.
2394 #else
2395  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
2396  call error(errormsg)
2397 #endif
2398 
2399  if (flag_init_output) &
2400  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2401  flag_3d_output, ndat2d, ndat3d)
2402 
2403 #elif (OUTPUT==2)
2404 
2405 if (time_output(1) <= time_init+eps) then
2406 
2407 #if (ERGDAT==0)
2408  flag_3d_output = .false.
2409 #elif (ERGDAT==1)
2410  flag_3d_output = .true.
2411 #else
2412  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
2413  call error(errormsg)
2414 #endif
2415 
2416  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2417  flag_3d_output, ndat2d, ndat3d)
2418 
2419 end if
2420 
2421 #elif (OUTPUT==3)
2422 
2423  flag_3d_output = .false.
2424 
2425  if (flag_init_output) &
2426  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2427  flag_3d_output, ndat2d, ndat3d)
2428 
2429 if (time_output(1) <= time_init+eps) then
2430 
2431  flag_3d_output = .true.
2432 
2433  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2434  flag_3d_output, ndat2d, ndat3d)
2435 
2436 end if
2437 
2438 #else
2439  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
2440  call error(errormsg)
2441 #endif
2442 
2443 if (flag_init_output) then
2444  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2445  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2446 end if
2447 
2448 #else /* OpenAD */
2449 
2450 print *, ' >>> sico_init: not producing initial, typical outputs'
2451 print *, ' in adjoint mode.'
2452 
2453 #endif /* Normal vs. OpenAD */
2454 
2455 #if (defined(EXEC_MAKE_C_DIS_0))
2456 
2457 #if (DISC>0)
2458 
2459 call calc_c_dis_0(dxi, deta)
2460 
2461 errormsg = ' >>> sico_init: Routine calc_c_dis_0 successfully completed,' &
2462  // end_of_line &
2463  //' c_dis_0 written on file out_runname.dat' &
2464  // end_of_line &
2465  //' (in directory specified by OUTPATH).' &
2466  // end_of_line &
2467  //' Execution of SICOPOLIS stopped.'
2468 call error(errormsg) ! actually not an error,
2469  ! just a regular stop with an info message
2470 
2471 #else
2472  errormsg = ' >>> sico_init: EXEC_MAKE_C_DIS_0 requires DISC>0!'
2473  call error(errormsg)
2474 #endif
2475 
2476 #endif
2477 
2478 end subroutine sico_init
2479 
2480 !-------------------------------------------------------------------------------
2481 !> Definition of the initial surface and bedrock topography
2482 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2483 !! For present-day initial topography.
2484 !<------------------------------------------------------------------------------
2485 subroutine topography1(dxi, deta)
2486 
2487 #if (GRID==0 || GRID==1)
2489 #endif
2490 
2491  use metric_m
2492  use topograd_m
2493 
2494 implicit none
2495 
2496 real(dp), intent(out) :: dxi, deta
2497 
2498 integer(i4b) :: i, j, n
2499 integer(i4b) :: ios, n_dummy
2500 real(dp) :: d_dummy
2501 real(dp) :: xi0, eta0
2502 real(dp) :: h_ice, freeboard_ratio
2503 character :: ch_dummy
2504 
2505 character(len= 8) :: ch_imax
2506 character(len=128) :: fmt4
2507 character(len=256) :: filename_with_path
2508 
2509 write(ch_imax, fmt='(i8)') imax
2510 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2511 
2512 !-------- Read topography --------
2513 
2514 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2515  trim(zs_present_file)
2516 
2517 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2518 
2519 if (ios /= 0) then
2520  errormsg = ' >>> topography1: Error when opening the zs file!'
2521  call error(errormsg)
2522 end if
2523 
2524 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
2525 
2526 do j=jmax, 0, -1
2527  read(21, fmt=*) (zs(j,i), i=0,imax)
2528 end do
2529 
2530 close(21, status='keep')
2531 
2532 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2533  trim(zl_present_file)
2534 
2535 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2536 
2537 if (ios /= 0) then
2538  errormsg = ' >>> topography1: Error when opening the zl file!'
2539  call error(errormsg)
2540 end if
2541 
2542 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
2543 
2544 do j=jmax, 0, -1
2545  read(22, fmt=*) (zl(j,i), i=0,imax)
2546 end do
2547 
2548 close(22, status='keep')
2549 
2550 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2551  trim(zl0_file)
2552 
2553 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2554 
2555 if (ios /= 0) then
2556  errormsg = ' >>> topography1: Error when opening the zl0 file!'
2557  call error(errormsg)
2558 end if
2559 
2560 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2561 
2562 do j=jmax, 0, -1
2563  read(23, fmt=*) (zl0(j,i), i=0,imax)
2564 end do
2565 
2566 close(23, status='keep')
2567 
2568 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2569  trim(mask_present_file)
2570 
2571 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2572 
2573 if (ios /= 0) then
2574  errormsg = ' >>> topography1: Error when opening the mask file!'
2575  call error(errormsg)
2576 end if
2577 
2578 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2579 
2580 do j=jmax, 0, -1
2581  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2582 end do
2583 
2584 close(24, status='keep')
2585 
2586 #if (defined(ZB_PRESENT_FILE))
2587 
2588 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2589  trim(zb_present_file)
2590 
2591 open(25, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2592 
2593 if (ios /= 0) then
2594  errormsg = ' >>> topography1: Error when opening the zb file!'
2595  call error(errormsg)
2596 end if
2597 
2598 do n=1, 6; read(25, fmt='(a)') ch_dummy; end do
2599 
2600 do j=jmax, 0, -1
2601  read(25, fmt=*) (zb(j,i), i=0,imax)
2602 end do
2603 
2604 close(25, status='keep')
2605 
2606 #else
2607 
2608 write(6, fmt='(a)') ' >>> topography1: ZB_PRESENT_FILE not defined,'
2609 write(6, fmt='(a)') ' thus zb = zl assumed.'
2610 
2611 zb = zl
2612 
2613 #endif
2614 
2615 !-------- Further stuff --------
2616 
2617 dxi = dx *1000.0_dp ! km -> m
2618 deta = dx *1000.0_dp ! km -> m
2619 
2620 xi0 = x0 *1000.0_dp ! km -> m
2621 eta0 = y0 *1000.0_dp ! km -> m
2622 
2623 freeboard_ratio = (rho_sw-rho)/rho_sw
2624 
2625 do i=0, imax
2626 do j=0, jmax
2627 
2628  if (maske(j,i) <= 1_i1b) then
2629 
2630  zb(j,i) = zl(j,i) ! ensure consistency
2631 
2632  else if (maske(j,i) == 2_i1b) then
2633 
2634 #if (MARGIN==1 || MARGIN==2)
2635  zs(j,i) = zl(j,i) ! ensure
2636  zb(j,i) = zl(j,i) ! consistency
2637 #elif (MARGIN==3)
2638  zs(j,i) = 0.0_dp ! present-day
2639  zb(j,i) = 0.0_dp ! sea level
2640 #endif
2641 
2642  else if (maske(j,i) == 3_i1b) then
2643 
2644 #if (MARGIN==1 || (MARGIN==2 && MARINE_ICE_FORMATION==1))
2645  maske(j,i) = 2_i1b ! floating ice cut off
2646  zs(j,i) = zl(j,i)
2647  zb(j,i) = zl(j,i)
2648 #elif (MARGIN==2 && MARINE_ICE_FORMATION==2)
2649  maske(j,i) = 0_i1b ! floating ice becomes "underwater ice"
2650  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2651  zs(j,i) = zl(j,i)+h_ice
2652  zb(j,i) = zl(j,i)
2653 #elif (MARGIN==3)
2654  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2655  zs(j,i) = freeboard_ratio*h_ice ! ensure properly
2656  zb(j,i) = zs(j,i)-h_ice ! floating ice
2657 #endif
2658 
2659  end if
2660 
2661  xi(i) = xi0 + real(i,dp)*dxi
2662  eta(j) = eta0 + real(j,dp)*deta
2663 
2664  zm(j,i) = zb(j,i)
2665  n_cts(j,i) = -1_i1b
2666  kc_cts(j,i) = 0
2667 
2668  h_c(j,i) = zs(j,i)-zm(j,i)
2669  h_t(j,i) = 0.0_dp
2670 
2671  dzs_dtau(j,i) = 0.0_dp
2672  dzm_dtau(j,i) = 0.0_dp
2673  dzb_dtau(j,i) = 0.0_dp
2674  dzl_dtau(j,i) = 0.0_dp
2675  dh_c_dtau(j,i) = 0.0_dp
2676  dh_t_dtau(j,i) = 0.0_dp
2677 
2678 end do
2679 end do
2680 
2681 maske_old = maske
2682 
2683 !-------- Geographic coordinates, metric tensor,
2684 ! gradients of the topography --------
2685 
2686 do i=0, imax
2687 do j=0, jmax
2688 
2689 #if (GRID==0 || GRID==1) /* Stereographic projection */
2690 
2691  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2692  lambda0, phi0, lambda(j,i), phi(j,i))
2693 
2694 #elif (GRID==2) /* Geographic coordinates */
2695 
2696  lambda(j,i) = xi(i)
2697  phi(j,i) = eta(j)
2698 
2699 #endif
2700 
2701 end do
2702 end do
2703 
2704 call metric()
2705 
2706 #if (TOPOGRAD==0)
2707 call topograd_1(dxi, deta, 1)
2708 #elif (TOPOGRAD==1)
2709 call topograd_2(dxi, deta, 1)
2710 #endif
2711 
2712 !-------- Corresponding area of grid points --------
2713 
2714 do i=0, imax
2715 do j=0, jmax
2716  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2717 end do
2718 end do
2719 
2720 end subroutine topography1
2721 
2722 !-------------------------------------------------------------------------------
2723 !> Definition of the initial surface and bedrock topography
2724 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2725 !! For ice-free initial topography with relaxed lithosphere.
2726 !<------------------------------------------------------------------------------
2727 subroutine topography2(dxi, deta)
2728 
2729 #if (GRID==0 || GRID==1)
2731 #endif
2732 
2733  use metric_m
2734  use topograd_m
2735 
2736 implicit none
2737 
2738 real(dp), intent(out) :: dxi, deta
2739 
2740 integer(i4b) :: i, j, n
2741 integer(i4b) :: ios
2742 real(dp) :: xi0, eta0
2743 character :: ch_dummy
2744 
2745 character(len= 8) :: ch_imax
2746 character(len=128) :: fmt4
2747 character(len=256) :: filename_with_path
2748 
2749 write(ch_imax, fmt='(i8)') imax
2750 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2751 
2752 !-------- Read topography --------
2753 
2754 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2755  trim(zl0_file)
2756 
2757 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2758 
2759 if (ios /= 0) then
2760  errormsg = ' >>> topography2: Error when opening the zl0 file!'
2761  call error(errormsg)
2762 end if
2763 
2764 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2765 
2766 do j=jmax, 0, -1
2767  read(23, fmt=*) (zl0(j,i), i=0,imax)
2768 end do
2769 
2770 close(23, status='keep')
2771 
2772 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2773  trim(mask_present_file)
2774 
2775 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2776 
2777 if (ios /= 0) then
2778  errormsg = ' >>> topography2: Error when opening the mask file!'
2779  call error(errormsg)
2780 end if
2781 
2782 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2783 
2784 do j=jmax, 0, -1
2785  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2786 end do
2787 
2788 close(24, status='keep')
2789 
2790 !-------- Further stuff --------
2791 
2792 dxi = dx *1000.0_dp ! km -> m
2793 deta = dx *1000.0_dp ! km -> m
2794 
2795 xi0 = x0 *1000.0_dp ! km -> m
2796 eta0 = y0 *1000.0_dp ! km -> m
2797 
2798 do i=0, imax
2799 do j=0, jmax
2800 
2801  if (maske(j,i) <= 1_i1b) then
2802  maske(j,i) = 1_i1b
2803  zs(j,i) = zl0(j,i)
2804  zb(j,i) = zl0(j,i)
2805  zl(j,i) = zl0(j,i)
2806  else ! (maske(j,i) >= 2_i1b)
2807  maske(j,i) = 2_i1b
2808 #if (MARGIN==1 || MARGIN==2)
2809  zs(j,i) = zl0(j,i)
2810  zb(j,i) = zl0(j,i)
2811 #elif (MARGIN==3)
2812  zs(j,i) = 0.0_dp ! present-day
2813  zb(j,i) = 0.0_dp ! sea level
2814 #endif
2815  zl(j,i) = zl0(j,i)
2816  end if
2817 
2818  xi(i) = xi0 + real(i,dp)*dxi
2819  eta(j) = eta0 + real(j,dp)*deta
2820 
2821  zm(j,i) = zb(j,i)
2822  n_cts(j,i) = -1_i1b
2823  kc_cts(j,i) = 0
2824 
2825  h_c(j,i) = 0.0_dp
2826  h_t(j,i) = 0.0_dp
2827 
2828  dzs_dtau(j,i) = 0.0_dp
2829  dzm_dtau(j,i) = 0.0_dp
2830  dzb_dtau(j,i) = 0.0_dp
2831  dzl_dtau(j,i) = 0.0_dp
2832  dh_c_dtau(j,i) = 0.0_dp
2833  dh_t_dtau(j,i) = 0.0_dp
2834 
2835 end do
2836 end do
2837 
2838 maske_old = maske
2839 
2840 !-------- Geographic coordinates, metric tensor,
2841 ! gradients of the topography --------
2842 
2843 do i=0, imax
2844 do j=0, jmax
2845 
2846 #if (GRID==0 || GRID==1) /* Stereographic projection */
2847 
2848  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2849  lambda0, phi0, lambda(j,i), phi(j,i))
2850 
2851 #elif (GRID==2) /* Geographic coordinates */
2852 
2853  lambda(j,i) = xi(i)
2854  phi(j,i) = eta(j)
2855 
2856 #endif
2857 
2858 end do
2859 end do
2860 
2861 call metric()
2862 
2863 #if (TOPOGRAD==0)
2864 call topograd_1(dxi, deta, 1)
2865 #elif (TOPOGRAD==1)
2866 call topograd_2(dxi, deta, 1)
2867 #endif
2868 
2869 !-------- Corresponding area of grid points --------
2870 
2871 do i=0, imax
2872 do j=0, jmax
2873  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2874 end do
2875 end do
2876 
2877 end subroutine topography2
2878 
2879 !-------------------------------------------------------------------------------
2880 !> Definition of the initial surface and bedrock topography
2881 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2882 !! For initial topography from previous simulation.
2883 !<------------------------------------------------------------------------------
2884 subroutine topography3(dxi, deta, z_sl, anfdatname)
2885 
2886  use read_m, only : read_erg_nc
2887 
2888 #if (GRID==0 || GRID==1)
2890 #endif
2891 
2892  use metric_m
2893  use topograd_m
2894 
2895 implicit none
2896 
2897 character(len=100), intent(in) :: anfdatname
2898 
2899 real(dp), intent(out) :: dxi, deta, z_sl
2900 
2901 integer(i4b) :: i, j, n
2902 integer(i4b) :: ios
2903 character(len=256) :: filename_with_path
2904 character :: ch_dummy
2905 
2906 !-------- Read data from time-slice file of previous simulation --------
2907 
2908 call read_erg_nc(z_sl, anfdatname)
2909 
2910 !-------- Read topography of the relaxed bedrock --------
2911 
2912 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2913  trim(zl0_file)
2914 
2915 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2916 
2917 if (ios /= 0) then
2918  errormsg = ' >>> topography3: Error when opening the zl0 file!'
2919  call error(errormsg)
2920 end if
2921 
2922 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2923 
2924 do j=jmax, 0, -1
2925  read(23, fmt=*) (zl0(j,i), i=0,imax)
2926 end do
2927 
2928 close(23, status='keep')
2929 
2930 !-------- Further stuff --------
2931 
2932 dxi = dx *1000.0_dp ! km -> m
2933 deta = dx *1000.0_dp ! km -> m
2934 
2935 !-------- Geographic coordinates, metric tensor,
2936 ! gradients of the topography --------
2937 
2938 do i=0, imax
2939 do j=0, jmax
2940 
2941 #if (GRID==0 || GRID==1) /* Stereographic projection */
2942 
2943  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2944  lambda0, phi0, lambda(j,i), phi(j,i))
2945 
2946 #elif (GRID==2) /* Geographic coordinates */
2947 
2948  lambda(j,i) = xi(i)
2949  phi(j,i) = eta(j)
2950 
2951 #endif
2952 
2953 end do
2954 end do
2955 
2956 call metric()
2957 
2958 #if (TOPOGRAD==0)
2959 call topograd_1(dxi, deta, 1)
2960 #elif (TOPOGRAD==1)
2961 call topograd_2(dxi, deta, 1)
2962 #endif
2963 
2964 !-------- Corresponding area of grid points --------
2965 
2966 do i=0, imax
2967 do j=0, jmax
2968  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2969 end do
2970 end do
2971 
2972 end subroutine topography3
2973 
2974 !-------------------------------------------------------------------------------
2975 
2976 end module sico_init_m
2977 !
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
integer(i1b), 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...
Update of the flags for the land-terminating grounded front, marine-terminating grounded front...
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:55
real(dp), dimension(:), allocatable phi_core
phi_core(n): Geographical latitude of the prescribed borehole positions
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:845
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(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
subroutine, public flag_update_gf_gl_cf()
Main subroutine of flag_update_gf_gl_cf_m: Update of the flags for the land-terminating grounded fron...
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:60
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 error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
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:60
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) temp_precip_time_stp
temp_precip_time_stp: Time step of the surface-temperature and precipitation data ...
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:903
real(dp) temp_precip_time_max
temp_precip_time_max: Maximum time of the surface-temperature and precipitation data ...
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)
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
integer(i1b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
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
subroutine, public calc_c_dis_0(dxi, deta)
Constant in ice discharge parameterization (Greenland). [Determine (amount of magnitude of) constant ...
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:57
real(dp), dimension(0:jmax, 0:imax) smb_corr_in
smb_corr_in(j,i): Prescribed SMB correction
Definition: sico_vars_m.F90:43
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 ...
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
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.
subroutine, public init_temp_water_age_1_2()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==2: present-day initial topogr...
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:5214
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:36
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.
subroutine, public disc_param(dtime)
Ice discharge parameters (Greenland). [Assign ice discharge parameters.].
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Initial temperature, water content and age.
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
integer(i4b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
integer(i1b), dimension(0:jmax, 0:imax) maske_sedi
maske_sedi(j,i): Sediment mask. 1: hard rock, 7: soft sediment, 2: ocean.
Definition: sico_vars_m.F90:46
real(dp), dimension(0:jmax, 0:imax) precip_ma_present
precip_ma_present(j,i): Present-day mean annual precipitation rate at the ice surface ...
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) temp_precip_time_min
temp_precip_time_min: Minimum time of the surface-temperature and precipitation data ...
real(dp), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
subroutine, public stereo_forw_ellipsoid(lambda_val, phi_val, A, B, lambda0, phi0, x_val, y_val)
Forward stereographic projection for an ellipsoidal planet.
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, public disc_fields()
Dependence of ice discharge coefficient on latitude (Greenland). [Determine dependence of ice dischar...
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:203
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
Writing of error messages and stopping execution.
Definition: error_m.F90:35
integer(i4b) ndata_temp_precip
ndata_temp_precip: Number of surface-temperature and precipitation data
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
real(dp), dimension(:), allocatable y_core
y_core(n): Coordinate eta (= y) of the prescribed borehole positions
real(dp), dimension(:), allocatable temp_precip_time
temp_precip_time(n): Times of the surface-temperature and precipitation data
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
real(dp), dimension(0:jmax, 0:imax) precip_ma_lgm_anom
precip_ma_lgm_anom(j,i): LGM anomaly (ratio LGM/present) of the mean annual precipitation rate at the...
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.
character, parameter end_of_line
end_of_line: End-of-line string
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.
integer(i1b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
NetCDF error capturing.
Definition: nc_check_m.F90:35
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
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(0:jmax, 0:imax) temp_mj_lgm_anom
temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean summer (northern hemisphere...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
integer(i1b), 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) 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:561
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, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
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:61
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
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) ...
real(dp), dimension(0:jmax, 0:imax) temp_ma_lgm_anom
temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean annual surface temperature ...
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:3869
integer(i1b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check_m.F90:48
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
integer(i4b) ncid_temp_precip
ncid_temp_precip: ID of the NetCDF file containing the surface-temperature and precipitation data as ...
character(len=16), dimension(:), allocatable ch_core
ch_core(n): Names of the prescribed borehole positions
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)
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(:), allocatable x_core
x_core(n): Coordinate xi (= x) of the prescribed borehole positions
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:57
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
real(dp) rho_sw
RHO_SW: Density of sea water.
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
character(len=256) errormsg
errormsg: Error-message string
Ice discharge parameterization for the Greenland ice sheet.
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
subroutine, public init_temp_water_age_1_4()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==4: present-day initial topogr...
subroutine, public init_temp_water_age_1_1()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==1: present-day initial topogr...
real(dp), dimension(:), allocatable lambda_core
lambda_core(n): Geographical longitude of the prescribed borehole positions
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
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
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_m.F90:734
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) 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:55
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...