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