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