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