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