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