SICOPOLIS V5-dev  Revision 1279
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 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
862 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
863 #endif
864 #if (defined(SEDI_SLIDE))
865 write(10, fmt=trim(fmt2a)) 'SEDI_SLIDE = ', sedi_slide
866 #endif
867 #if (defined(SEDI_SLIDE) && SEDI_SLIDE==2)
868 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
869 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
870 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
871 write(10, fmt=trim(fmt2a)) 'p_weert_sedi = ', p_weert_sedi
872 write(10, fmt=trim(fmt2a)) 'q_weert_sedi = ', q_weert_sedi
873 #if (defined(TRANS_LENGTH_SL))
874 write(10, fmt=trim(fmt3)) 'trans_length_sl =', trans_length_sl
875 #endif
876 #endif
877 write(10, fmt=trim(fmt1)) ' '
878 
879 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
880 #if (Q_GEO_MOD==1)
881 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
882 #elif (Q_GEO_MOD==2)
883 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
884 #endif
885 write(10, fmt=trim(fmt1)) ' '
886 
887 #if (defined(MARINE_ICE_BASAL_MELTING))
888 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
889 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
890 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
891 #endif
892 write(10, fmt=trim(fmt1)) ' '
893 #endif
894 
895 #if (MARGIN==3)
896 write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
897 #if (FLOATING_ICE_BASAL_MELTING<=3)
898 write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
899 write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
900 #endif
901 write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
902 write(10, fmt=trim(fmt3)) 'z_abyss =', z_abyss
903 #if (FLOATING_ICE_BASAL_MELTING==4)
904 write(10, fmt=trim(fmt3)) 'temp_ocean =', temp_ocean
905 write(10, fmt=trim(fmt3)) 'Omega_qbm =', omega_qbm
906 write(10, fmt=trim(fmt3)) 'alpha_qbm =', alpha_qbm
907 #endif
908 #if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
909 write(10, fmt=trim(fmt3)) 'H_w_0 =', h_w_0
910 #endif
911 #if (defined(INITMIP_BMB_ANOM_FILE))
912 if ( trim(adjustl(initmip_bmb_anom_file)) /= 'none' ) then
913  write(10, fmt=trim(fmt1)) 'initmip_bmb_anom file = '//initmip_bmb_anom_file
914 end if
915 #endif
916 #if (defined(LARMIP_REGIONS_FILE))
917 write(10, fmt=trim(fmt1)) 'larmip_regions_file = '//larmip_regions_file
918 if ( trim(adjustl(larmip_regions_file)) /= 'none' ) then
919  write(10, fmt=trim(fmt3)) 'larmip_qbm_anom_1 =', larmip_qbm_anom_1
920  write(10, fmt=trim(fmt3)) 'larmip_qbm_anom_2 =', larmip_qbm_anom_2
921  write(10, fmt=trim(fmt3)) 'larmip_qbm_anom_3 =', larmip_qbm_anom_3
922  write(10, fmt=trim(fmt3)) 'larmip_qbm_anom_4 =', larmip_qbm_anom_4
923  write(10, fmt=trim(fmt3)) 'larmip_qbm_anom_5 =', larmip_qbm_anom_5
924 end if
925 #endif
926 write(10, fmt=trim(fmt1)) ' '
927 #endif
928 
929 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
930 #if (REBOUND==1)
931 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
932 #endif
933 #if (REBOUND==1 || REBOUND==2)
934 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
935 #if (TIME_LAG_MOD==1)
936 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
937 #elif (TIME_LAG_MOD==2)
938 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
939 #else
940 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
941 call error(errormsg)
942 #endif
943 #endif
944 #if (REBOUND==2)
945 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
946 #if (FLEX_RIG_MOD==1)
947 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
948 #elif (FLEX_RIG_MOD==2)
949 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
950 #else
951 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
952 call error(errormsg)
953 #endif
954 #endif
955 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
956 write(10, fmt=trim(fmt1)) ' '
957 
958 #if (FLOW_LAW==2)
959 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
960 write(10, fmt=trim(fmt1)) ' '
961 #endif
962 #if (FIN_VISC==2)
963 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
964 write(10, fmt=trim(fmt1)) ' '
965 #endif
966 
967 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
968 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
969 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
970 #endif
971 #if (ENHMOD==2 || ENHMOD==3)
972 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
973 #endif
974 #if (ENHMOD==2)
975 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
976 #endif
977 #if (ENHMOD==3)
978 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
979 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
980 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
981 #endif
982 #if (ENHMOD==4 || ENHMOD==5)
983 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
984 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
985 #endif
986 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
987 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
988 #endif
989 write(10, fmt=trim(fmt1)) ' '
990 
991 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
992 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
993 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
994 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
995 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
996 #if (defined(QBM_MIN))
997 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
998 #elif (defined(QB_MIN)) /* obsolete */
999 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
1000 #endif
1001 #if (defined(QBM_MAX))
1002 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
1003 #elif (defined(QB_MAX)) /* obsolete */
1004 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
1005 #endif
1006 write(10, fmt=trim(fmt3)) 'age_min =', age_min
1007 write(10, fmt=trim(fmt3)) 'age_max =', age_max
1008 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
1009 #if (ADV_VERT==1)
1010 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
1011 #endif
1012 write(10, fmt=trim(fmt1)) ' '
1013 
1014 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
1015 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
1016 #if (defined(LIS_OPTS))
1017 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
1018 #endif
1019 #if (defined(N_ITER_SSA))
1020 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
1021 #endif
1022 #if (defined(RELAX_FACT_SSA))
1023 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
1024 #endif
1025 #endif
1026 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
1027 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
1028 #endif
1029 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
1030 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
1031 #endif
1032 write(10, fmt=trim(fmt1)) ' '
1033 
1034 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
1035 #if (CALCMOD==-1 && defined(TEMP_CONST))
1036 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
1037 #endif
1038 #if (CALCMOD==-1 && defined(AGE_CONST))
1039 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
1040 #endif
1041 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
1042 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
1043 #endif
1044 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
1045 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
1046 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
1047 #if (MARGIN==2)
1048 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
1049 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
1050 #elif (MARGIN==3)
1051 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
1052 #endif
1053 #if (defined(THK_EVOL))
1054 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
1055 #else
1056 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
1057 call error(errormsg)
1058 #endif
1059 #if (defined(CALCTHK))
1060 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
1061 #else
1062 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
1063 call error(errormsg)
1064 #endif
1065 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
1066 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
1067 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
1068 write(10, fmt=trim(fmt1)) ' '
1069 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
1070 #if (ACCSURFACE==5)
1071 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
1072 #endif
1073 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
1074 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
1075 #if (defined(MB_ACCOUNT))
1076 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
1077 #endif
1078 write(10, fmt=trim(fmt1)) ' '
1079 
1080 #if (defined(OUT_TIMES))
1081 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
1082 #endif
1083 #if (defined(OUTPUT_INIT))
1084 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
1085 #endif
1086 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
1087 #if (OUTPUT==1 || OUTPUT==3)
1088 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
1089 #endif
1090 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
1091 #if (OUTPUT==1 || OUTPUT==2)
1092 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
1093 #endif
1094 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
1095 #if (OUTPUT==2 || OUTPUT==3)
1096 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
1097 do n=1, n_output
1098  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
1099 end do
1100 #endif
1101 write(10, fmt=trim(fmt1)) ' '
1102 
1103 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
1104 
1105 close(10, status='keep')
1106 
1107 #else /* OpenAD */
1108 
1109 print *, ' >>> sico_init: not writing to the log file in adjoint!'
1110 
1111 #endif /* Normal vs. OpenAD */
1112 
1113 !-------- Conversion of time quantities --------
1114 
1115 year_zero = year_zero*year_sec ! a --> s
1116 time_init = time_init0*year_sec ! a --> s
1117 time_end = time_end0*year_sec ! a --> s
1118 dtime = dtime0*year_sec ! a --> s
1119 dtime_temp = dtime_temp0*year_sec ! a --> s
1120 #if (REBOUND==2)
1121 dtime_wss = dtime_wss0*year_sec ! a --> s
1122 #endif
1123 dtime_ser = dtime_ser0*year_sec ! a --> s
1124 #if (OUTPUT==1 || OUTPUT==3)
1125 dtime_out = dtime_out0*year_sec ! a --> s
1126 #endif
1127 #if (OUTPUT==2 || OUTPUT==3)
1128 do n=1, n_output
1129  time_output(n) = time_output0(n)*year_sec ! a --> s
1130 end do
1131 #endif
1132 
1133 #ifndef ALLOW_OPENAD /* Normal */
1134 
1135 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
1136  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
1137  call error(errormsg)
1138 end if
1139 
1140 #if (REBOUND==2)
1141 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) then
1142  errormsg = ' >>> sico_init: dtime_wss must be a multiple of dtime!'
1143  call error(errormsg)
1144 end if
1145 #endif
1146 
1147 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
1148  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
1149  call error(errormsg)
1150 end if
1151 
1152 #if (OUTPUT==1 || OUTPUT==3)
1153 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
1154  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
1155  call error(errormsg)
1156 end if
1157 #endif
1158 
1159 #else /* OpenAD */
1160 
1161 print *, ' >>> sico_init: compare_float not used in adjoint applications!'
1162 print *, ' double check your time step sizes are multples'
1163 print *, ' of each other.'
1164 
1165 #endif /* Normal vs. OpenAD */
1166 
1167 #if (THK_EVOL==2)
1168 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
1169 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
1170 #endif
1171 
1172 #if (THK_EVOL==3)
1173 target_topo_tau = target_topo_tau0 *year_sec ! a --> s
1174 #endif
1175 
1176 
1177 time = time_init
1178 
1179 !-------- Reading of present mean-annual precipitation rate --------
1180 
1181 #if (GRID==0 || GRID==1)
1182 
1183 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1184  trim(precip_present_file)
1185 
1186 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1187 
1188 #elif (GRID==2)
1189 
1190 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1191 call error(errormsg)
1192 
1193 #endif
1194 
1195 if (ios /= 0) then
1196  errormsg = ' >>> sico_init: Error when opening the precip file!'
1197  call error(errormsg)
1198 end if
1199 
1200 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1201 
1202 do j=jmax, 0, -1
1203  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
1204 end do
1205 
1206 close(21, status='keep')
1207 
1208 precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
1209  ! mm/a water equiv. -> m/s ice equiv.
1210 
1211 !-------- Present monthly precipitation rates
1212 ! (assumed to be equal to the mean annual precipitation rate) --------
1213 
1214 do n=1, 12 ! month counter
1215 do i=0, imax
1216 do j=0, jmax
1217  precip_present(j,i,n) = precip_ma_present(j,i)
1218 end do
1219 end do
1220 end do
1221 
1222 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
1223 
1224 #if (ACCSURFACE==5)
1225 
1226 #if (GRID==0 || GRID==1)
1227 
1228 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1229  trim(precip_anom_file)
1230 
1231 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1232 
1233 #endif
1234 
1235 if (ios /= 0) then
1236  errormsg = ' >>> sico_init: Error when opening the precip anomaly file!'
1237  call error(errormsg)
1238 end if
1239 
1240 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1241 
1242 do j=jmax, 0, -1
1243  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
1244 end do
1245 
1246 close(21, status='keep')
1247 
1248 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
1249 
1250 !-------- LGM monthly precipitation-rate anomalies (assumed to be
1251 ! equal to the mean annual precipitation-rate anomaly) --------
1252 
1253 do n=1, 12 ! month counter
1254 do i=0, imax
1255 do j=0, jmax
1256 
1257  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
1258  ! positive values ensured
1259 
1260 #if (PRECIP_ANOM_INTERPOL==1)
1261  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
1262 #elif (PRECIP_ANOM_INTERPOL==2)
1263  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1264 #else
1265  errormsg = ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1266  call error(errormsg)
1267 #endif
1268 
1269 end do
1270 end do
1271 end do
1272 
1273 #endif
1274 
1275 !-------- Mean accumulation --------
1276 
1277 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1278  ! mm/a water equiv. -> m/s ice equiv.
1279 
1280 !-------- Read present topography mask --------
1281 
1282 #if (GRID==0 || GRID==1)
1283 
1284 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1285  trim(mask_present_file)
1286 
1287 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1288 
1289 #elif (GRID==2)
1290 
1291 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1292 call error(errormsg)
1293 
1294 #endif
1295 
1296 if (ios /= 0) then
1297  errormsg = ' >>> sico_init: Error when opening the mask file!'
1298  call error(errormsg)
1299 end if
1300 
1301 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1302 
1303 do j=jmax, 0, -1
1304  read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
1305 end do
1306 
1307 close(24, status='keep')
1308 
1309 !-------- Read rock/sediment mask --------
1310 
1311 #if (defined(SEDI_SLIDE) && SEDI_SLIDE==2)
1312 
1313 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1314  trim(mask_sedi_file)
1315 
1316 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1317 
1318 if (ios /= 0) then
1319  errormsg = ' >>> sico_init: Error when opening the rock/sediment mask file!'
1320  call error(errormsg)
1321 end if
1322 
1323 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1324 
1325 do j=jmax, 0, -1
1326  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
1327 end do
1328 
1329 close(24, status='keep')
1330 
1331 do i=0, imax
1332 do j=0, jmax
1333  if (maske_sedi(j,i) == 2_i2b) maske_sedi(j,i) = 7_i2b
1334  ! ocean to be treated like soft sediment
1335 end do
1336 end do
1337 
1338 #endif
1339 
1340 !-------- Initialisation of the marker for the different sectors for
1341 ! ice shelf basal melting --------
1342 
1343 n_sector = 0_i2b
1344 
1345 !-------- Present reference elevation (for precipitation data) --------
1346 
1347 #if (GRID==0 || GRID==1)
1348 
1349 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1350  trim(zs_present_file)
1351 
1352 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1353 
1354 #elif (GRID==2)
1355 
1356 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1357 call error(errormsg)
1358 
1359 #endif
1360 
1361 if (ios /= 0) then
1362  errormsg = ' >>> sico_init: Error when opening the zs file!'
1363  call error(errormsg)
1364 end if
1365 
1366 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1367 
1368 do j=jmax, 0, -1
1369  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1370 end do
1371 
1372 close(21, status='keep')
1373 
1374 ! ------ Reset bathymetry data (sea floor elevation) to the
1375 ! sea surface
1376 
1377 do i=0, imax
1378 do j=0, jmax
1379  if (maske_ref(j,i) >= 2_i2b) zs_ref(j,i) = 0.0_dp
1380 end do
1381 end do
1382 
1383 !-------- Reading of the prescribed target topography --------
1384 
1385 #if ((THK_EVOL==2) || (THK_EVOL==3))
1386 
1387 target_topo_dat_name = trim(target_topo_dat_name)
1388 
1389 #if (NETCDF==1)
1390 errormsg = ' >>> sico_init: Reading of target topography' &
1391  // end_of_line &
1392  //' only implemented for NetCDF files (NETCDF==2)!'
1393 call error(errormsg)
1394 #elif (NETCDF==2)
1395 call read_target_topo_nc(target_topo_dat_name)
1396 #else
1397 errormsg = ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1398 call error(errormsg)
1399 #endif
1400 
1401 #endif
1402 
1403 !-------- Reading of the maximum ice extent mask --------
1404 
1405 #if (THK_EVOL==4)
1406 
1407 #if (GRID==0 || GRID==1)
1408 
1409 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1410  trim(mask_maxextent_file)
1411 
1412 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1413 
1414 #elif (GRID==2)
1415 
1416 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1417 call error(errormsg)
1418 
1419 #endif
1420 
1421 if (ios /= 0) then
1422  errormsg = ' >>> sico_init: ' &
1423  //'Error when opening the maximum ice extent mask file!'
1424  call error(errormsg)
1425 end if
1426 
1427 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1428 
1429 do j=jmax, 0, -1
1430  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1431 end do
1432 
1433 close(24, status='keep')
1434 
1435 #endif
1436 
1437 !-------- Reading of LGM mean-annual and mean-January (summer)
1438 ! surface-temperature anomalies --------
1439 
1440 #if (TSURFACE==5)
1441 
1442 #if (GRID==0 || GRID==1)
1443 
1444 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1445  trim(temp_ma_anom_file)
1446 
1447 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1448 
1449 #endif
1450 
1451 if (ios /= 0) then
1452  errormsg = ' >>> sico_init: Error when opening the temp_ma anomaly file!'
1453  call error(errormsg)
1454 end if
1455 
1456 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1457 
1458 do j=jmax, 0, -1
1459  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
1460 end do
1461 
1462 close(21, status='keep')
1463 
1464 #if (GRID==0 || GRID==1)
1465 
1466 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1467  trim(temp_mj_anom_file)
1468 
1469 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1470 
1471 #endif
1472 
1473 if (ios /= 0) then
1474  errormsg = ' >>> sico_init: Error when opening the temp_mj anomaly file!'
1475  call error(errormsg)
1476 end if
1477 
1478 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1479 
1480 do j=jmax, 0, -1
1481  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
1482 end do
1483 
1484 close(22, status='keep')
1485 
1486 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
1487 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
1488 
1489 #endif
1490 
1491 !-------- Read data for delta_ts --------
1492 
1493 #if (TSURFACE==4)
1494 
1495 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1496 
1497 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1498 
1499 if (ios /= 0) then
1500  errormsg = ' >>> sico_init: Error when opening the data file for delta_ts!'
1501  call error(errormsg)
1502 end if
1503 
1504 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1505 
1506 if (ch_dummy /= '#') then
1507  errormsg = ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' &
1508  // end_of_line &
1509  //' not defined in data file for delta_ts!'
1510  call error(errormsg)
1511 end if
1512 
1514 
1515 allocate(griptemp(0:ndata_grip))
1516 
1517 do n=0, ndata_grip
1518  read(21, fmt=*) d_dummy, griptemp(n)
1519 end do
1520 
1521 close(21, status='keep')
1522 
1523 #endif
1524 
1525 !-------- Read data for the glacial index --------
1526 
1527 #if (TSURFACE==5)
1528 
1529 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1530 
1531 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1532 
1533 if (ios /= 0) then
1534  errormsg = ' >>> sico_init: Error when opening the glacial-index file!'
1535  call error(errormsg)
1536 end if
1537 
1538 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1539 
1540 if (ch_dummy /= '#') then
1541  errormsg = ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max' &
1542  // end_of_line &
1543  //' not defined in glacial-index file!'
1544  call error(errormsg)
1545 end if
1546 
1548 
1549 allocate(glacial_index(0:ndata_gi))
1550 
1551 do n=0, ndata_gi
1552  read(21, fmt=*) d_dummy, glacial_index(n)
1553 end do
1554 
1555 close(21, status='keep')
1556 
1557 #endif
1558 
1559 !-------- Reading of time-dependent surface-temperature
1560 ! and precipitation anomaly data --------
1561 
1562 #if (TSURFACE==6 && ACCSURFACE==6)
1563 
1564 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1565  '/'//trim(temp_precip_anom_file)
1566 
1567 call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
1568  ncid_temp_precip), thisroutine )
1569 
1570 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1571  thisroutine )
1572 
1573 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1574  len = ndata_temp_precip), thisroutine )
1575 
1577 
1579 
1580 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1581 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1582  ! in a
1583 
1586  ! in a; constant time step assumed
1588 
1589 #endif
1590 
1591 !-------- Prescribed surface mass balance correction --------
1592 
1593 smb_corr_in = 0.0_dp
1594 
1595 #if (defined(SMB_CORR_FILE))
1596 
1597 if (trim(adjustl(smb_corr_file)) /= 'none') then
1598 
1599  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1600  trim(smb_corr_file)
1601 
1602  filename_with_path = adjustr(filename_with_path)
1603  n = len(filename_with_path)
1604  ch_nc_test = filename_with_path(n-2:n)
1605  filename_with_path = adjustl(filename_with_path)
1606 
1607  if (ch_nc_test == '.nc') then ! NetCDF file
1608 
1609 #if (NETCDF==1)
1610 
1611  errormsg = ' >>> sico_init: Reading of NetCDF SMB_CORR_FILE' &
1612  // end_of_line &
1613  //' requires NETCDF==2!'
1614  call error(errormsg)
1615 
1616 #elif (NETCDF==2)
1617 
1618  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1619  thisroutine )
1620 
1621  call check( nf90_inq_varid(ncid, 'DSMB', ncv) )
1622  call check( nf90_get_var(ncid, ncv, smb_corr_in_conv) )
1623 
1624  call check( nf90_close(ncid) )
1625 
1626  do i=0, imax
1627  do j=0, jmax
1628  smb_corr_in(j,i) = smb_corr_in_conv(i,j) /year_sec
1629  ! m/a ice equiv. -> m/s ice equiv.
1630  end do
1631  end do
1632 
1633 #else
1634  errormsg = ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1635  call error(errormsg)
1636 #endif
1637 
1638  else ! ASCII file
1639 
1640  open(21, iostat=ios, &
1641  file=trim(filename_with_path), recl=rcl1, status='old')
1642 
1643  if (ios /= 0) then
1644  errormsg = ' >>> sico_init: ' &
1645  //'Error when opening the SMB correction file!'
1646  call error(errormsg)
1647  end if
1648 
1649  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1650 
1651  do j=jmax, 0, -1
1652  read(21, fmt=*) (smb_corr_in(j,i), i=0,imax)
1653  end do
1654 
1655  close(21, status='keep')
1656 
1657  smb_corr_in = smb_corr_in /year_sec
1658  ! m/a ice equiv. -> m/s ice equiv.
1659 
1660  end if
1661 
1662 end if
1663 
1664 #endif
1665 
1666 !-------- Reading of ISMIP6 SMB and BMB anomaly data --------
1667 
1668 ! ------ SMB (InitMIP)
1669 
1670 #if (defined(INITMIP_SMB_ANOM_FILE))
1671 
1672 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
1673 
1674  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1675  '/'//trim(initmip_smb_anom_file)
1676 
1677  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1678  thisroutine )
1679 
1680  call check( nf90_inq_varid(ncid, 'asmb', ncv) )
1681  call check( nf90_get_var(ncid, ncv, smb_anom_initmip_conv) )
1682 
1683  call check( nf90_close(ncid) )
1684 
1685  do i=0, imax
1686  do j=0, jmax
1687  smb_anom_initmip(j,i) = real(smb_anom_initmip_conv(i,j),dp) /YEAR_SEC
1688  ! m/a ice equiv. -> m/s ice equiv.
1689  end do
1690  end do
1691 
1692 else
1693  smb_anom_initmip = 0.0_dp
1694 end if
1695 
1696 #endif
1697 
1698 ! ------ BMB (InitMIP)
1699 
1700 #if (defined(INITMIP_BMB_ANOM_FILE))
1701 
1702 if ( trim(adjustl(initmip_bmb_anom_file)) /= 'none' ) then
1703 
1704  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1705  '/'//trim(initmip_bmb_anom_file)
1706 
1707  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1708  thisroutine )
1709 
1710  call check( nf90_inq_varid(ncid, 'abmb', ncv) )
1711  call check( nf90_get_var(ncid, ncv, ab_anom_initmip_conv) )
1712 
1713  call check( nf90_close(ncid) )
1714 
1715  do i=0, imax
1716  do j=0, jmax
1717  ab_anom_initmip(j,i) = real(ab_anom_initmip_conv(i,j),dp) /YEAR_SEC
1718  ! m/a ice equiv. -> m/s ice equiv.
1719  end do
1720  end do
1721 
1722 else
1723  ab_anom_initmip = 0.0_dp
1724 end if
1725 
1726 #endif
1727 
1728 ! ------ BMB (LARMIP)
1729 
1730 #if (defined(LARMIP_REGIONS_FILE))
1731 
1732 if ( trim(adjustl(larmip_regions_file)) /= 'none' ) then
1733 
1734  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1735  '/'//trim(larmip_regions_file)
1736 
1737  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1738  thisroutine )
1739 
1740  call check( nf90_inq_varid(ncid, 'regions', ncv) )
1741  call check( nf90_get_var(ncid, ncv, regions_larmip_conv) )
1742 
1743  call check( nf90_close(ncid) )
1744 
1745  do i=0, imax
1746  do j=0, jmax
1747  n_regions_larmip(j,i) = nint(regions_larmip_conv(i,j),i1b)
1748  end do
1749  end do
1750 
1751  ab_anom_larmip = 0.0_dp
1752 
1753  ab_anom_larmip(1) = larmip_qbm_anom_1 /year_sec
1754  ! m/a ice equiv. -> m/s ice equiv.
1755  ab_anom_larmip(2) = larmip_qbm_anom_2 /year_sec
1756  ! m/a ice equiv. -> m/s ice equiv.
1757  ab_anom_larmip(3) = larmip_qbm_anom_3 /year_sec
1758  ! m/a ice equiv. -> m/s ice equiv.
1759  ab_anom_larmip(4) = larmip_qbm_anom_4 /year_sec
1760  ! m/a ice equiv. -> m/s ice equiv.
1761  ab_anom_larmip(5) = larmip_qbm_anom_5 /year_sec
1762  ! m/a ice equiv. -> m/s ice equiv.
1763 else
1764  n_regions_larmip = 0_i1b
1765  ab_anom_larmip = 0.0_dp
1766 end if
1767 
1768 #endif
1769 
1770 !-------- Read data for z_sl --------
1771 
1772 #if (SEA_LEVEL==3)
1773 
1774 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1775 
1776 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1777 
1778 if (ios /= 0) then
1779  errormsg = ' >>> sico_init: Error when opening the data file for z_sl!'
1780  call error(errormsg)
1781 end if
1782 
1783 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1784 
1785 if (ch_dummy /= '#') then
1786  errormsg = ' >>> sico_init:' &
1787  // end_of_line &
1788  //' specmap_time_min, specmap_time_stp, specmap_time_max' &
1789  // end_of_line &
1790  //' not defined in data file for z_sl!'
1791  call error(errormsg)
1792 end if
1793 
1795 
1796 allocate(specmap_zsl(0:ndata_specmap))
1797 
1798 do n=0, ndata_specmap
1799  read(21, fmt=*) d_dummy, specmap_zsl(n)
1800 end do
1801 
1802 close(21, status='keep')
1803 
1804 #endif
1805 
1806 !-------- Determination of the geothermal heat flux --------
1807 
1808 #if (Q_GEO_MOD==1)
1809 
1810 ! ------ Constant value
1811 
1812 do i=0, imax
1813 do j=0, jmax
1814  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1815 end do
1816 end do
1817 
1818 #elif (Q_GEO_MOD==2)
1819 
1820 ! ------ Read data from file
1821 
1822 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1823  trim(q_geo_file)
1824 
1825 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1826 
1827 if (ios /= 0) then
1828  errormsg = ' >>> sico_init: Error when opening the qgeo file!'
1829  call error(errormsg)
1830 end if
1831 
1832 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1833 
1834 do j=jmax, 0, -1
1835  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1836 end do
1837 
1838 close(21, status='keep')
1839 
1840 do i=0, imax
1841 do j=0, jmax
1842  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1843 end do
1844 end do
1845 
1846 #endif
1847 
1848 !-------- Reading of tabulated kei function--------
1849 
1850 #if (REBOUND==0 || REBOUND==1)
1851 
1852 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1853  ! dummy values
1854 #elif (REBOUND==2)
1855 
1856 call read_kei()
1857 
1858 #endif
1859 
1860 !-------- Determination of the time lag
1861 ! of the relaxing asthenosphere --------
1862 
1863 #if (REBOUND==1 || REBOUND==2)
1864 
1865 #if (TIME_LAG_MOD==1)
1866 
1867 time_lag_asth = time_lag*year_sec ! a -> s
1868 
1869 #elif (TIME_LAG_MOD==2)
1870 
1871 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1872  trim(time_lag_file)
1873 
1874 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1875 
1876 if (ios /= 0) then
1877  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
1878  call error(errormsg)
1879 end if
1880 
1881 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1882 
1883 do j=jmax, 0, -1
1884  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1885 end do
1886 
1887 close(29, status='keep')
1888 
1889 time_lag_asth = time_lag_asth*year_sec ! a -> s
1890 
1891 #endif
1892 
1893 #elif (REBOUND==0)
1894 
1895 time_lag_asth = 0.0_dp ! dummy values
1896 
1897 #endif
1898 
1899 !-------- Determination of the flexural rigidity of the lithosphere --------
1900 
1901 #if (REBOUND==2)
1902 
1903 #if (FLEX_RIG_MOD==1)
1904 
1905 flex_rig_lith = flex_rig
1906 
1907 #elif (FLEX_RIG_MOD==2)
1908 
1909 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1910  trim(flex_rig_file)
1911 
1912 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1913 
1914 if (ios /= 0) then
1915  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1916  call error(errormsg)
1917 end if
1918 
1919 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1920 
1921 do j=jmax, 0, -1
1922  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1923 end do
1924 
1925 close(29, status='keep')
1926 
1927 #endif
1928 
1929 #elif (REBOUND==0 || REBOUND==1)
1930 
1931 flex_rig_lith = 0.0_dp ! dummy values
1932 
1933 #endif
1934 
1935 !-------- Definition of initial values --------
1936 
1937 ! ------ Present topography
1938 
1939 #if (ANF_DAT==1)
1940 
1941 call topography1(dxi, deta)
1942 
1943 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1944 
1945 call boundary(time_init, dtime, dxi, deta, &
1946  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1947 
1948 #if ( !defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD) ) /* Normal */
1949 
1950 where ((maske==0_i2b).or.(maske==3_i2b))
1951  ! grounded or floating ice
1953 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1954  as_perp_apl = 0.0_dp
1955 end where
1956 
1957 #else /* OpenAD */
1958 
1959 do j=0,jmax
1960 do i=0,imax
1961  if ((maske(j,i)==0_i2b).or.(maske(j,i)==3_i2b)) then
1962  as_perp_apl(j,i) = as_perp(j,i)
1963  else
1964  as_perp_apl(j,i) = 0.0_dp
1965  end if
1966 end do
1967 end do
1968 
1969 #endif /* Normal vs. OpenAD */
1970 
1971 smb_corr = 0.0_dp
1972 
1973 q_bm = 0.0_dp
1974 q_tld = 0.0_dp
1975 q_b_tot = 0.0_dp
1976 
1977 p_b_w = 0.0_dp
1978 h_w = 0.0_dp
1979 
1980 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1982 #elif (TEMP_INIT==2)
1984 #elif (TEMP_INIT==3)
1986 #elif (TEMP_INIT==4)
1988 #else
1989  errormsg = ' >>> sico_init: TEMP_INIT must be set to either 1, 2, 3 or 4!'
1990  call error(errormsg)
1991 #endif
1992 
1993 #if (ENHMOD==1)
1994  call calc_enhance_1()
1995 #elif (ENHMOD==2)
1996  call calc_enhance_2()
1997 #elif (ENHMOD==3)
1998  call calc_enhance_3(time_init)
1999 #elif (ENHMOD==4)
2000  call calc_enhance_4()
2001 #elif (ENHMOD==5)
2002  call calc_enhance_5()
2003 #else
2004  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
2005  call error(errormsg)
2006 #endif
2007 
2008 ! ------ Ice-free, relaxed bedrock
2009 
2010 #elif (ANF_DAT==2)
2011 
2012 call topography2(dxi, deta)
2013 
2014 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
2015 
2016 call boundary(time_init, dtime, dxi, deta, &
2017  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
2018 
2019 as_perp_apl = 0.0_dp
2020 
2021 smb_corr = 0.0_dp
2022 
2023 q_bm = 0.0_dp
2024 q_tld = 0.0_dp
2025 q_b_tot = 0.0_dp
2026 
2027 p_b_w = 0.0_dp
2028 h_w = 0.0_dp
2029 
2030 call init_temp_water_age_2()
2031 
2032 #if (ENHMOD==1)
2033  call calc_enhance_1()
2034 #elif (ENHMOD==2)
2035  call calc_enhance_2()
2036 #elif (ENHMOD==3)
2037  call calc_enhance_3(time_init)
2038 #elif (ENHMOD==4)
2039  call calc_enhance_4()
2040 #elif (ENHMOD==5)
2041  call calc_enhance_5()
2042 #else
2043  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
2044  call error(errormsg)
2045 #endif
2046 
2047 ! ------ Read initial state from output of previous simulation
2048 
2049 #elif (ANF_DAT==3)
2050 
2051 call topography3(dxi, deta, z_sl, anfdatname)
2052 
2053 call boundary(time_init, dtime, dxi, deta, &
2054  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
2055 
2056 #if ( !defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD) ) /* Normal */
2057 
2058 where ((maske==0_i2b).or.(maske==3_i2b))
2059  ! grounded or floating ice
2061 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
2062  as_perp_apl = 0.0_dp
2063 end where
2064 
2065 #else /* OpenAD */
2066 
2067 do j=0,jmax
2068 do i=0,imax
2069  if ((maske(j,i)==0_i2b).or.(maske(j,i)==3_i2b)) then
2070  as_perp_apl(j,i) = as_perp(j,i)
2071  else
2072  as_perp_apl(j,i) = 0.0_dp
2073  end if
2074 end do
2075 end do
2076 
2077 #endif /* Normal vs. OpenAD */
2078 
2079 smb_corr = 0.0_dp
2080 
2081 q_b_tot = q_bm + q_tld
2082 
2083 #endif
2084 
2085 !-------- Inner-point flag --------
2086 
2087 flag_inner_point = .true.
2088 
2089 flag_inner_point(0,:) = .false.
2090 flag_inner_point(jmax,:) = .false.
2091 
2092 flag_inner_point(:,0) = .false.
2093 flag_inner_point(:,imax) = .false.
2094 
2095 !-------- Grounding line and calving front flags --------
2096 
2097 flag_grounding_line_1 = .false.
2098 flag_grounding_line_2 = .false.
2099 
2100 flag_calving_front_1 = .false.
2101 flag_calving_front_2 = .false.
2102 
2103 #if (MARGIN==1 || MARGIN==2)
2104 
2105 continue
2106 
2107 #elif (MARGIN==3)
2108 
2109 do i=1, imax-1
2110 do j=1, jmax-1
2111 
2112  if ( (maske(j,i)==0_i2b) & ! grounded ice
2113  .and. &
2114  ( (maske(j,i+1)==3_i2b) & ! with
2115  .or.(maske(j,i-1)==3_i2b) & ! one
2116  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
2117  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
2118  ) &
2119  flag_grounding_line_1(j,i) = .true.
2120 
2121  if ( (maske(j,i)==3_i2b) & ! floating ice
2122  .and. &
2123  ( (maske(j,i+1)==0_i2b) & ! with
2124  .or.(maske(j,i-1)==0_i2b) & ! one
2125  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
2126  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
2127  ) &
2128  flag_grounding_line_2(j,i) = .true.
2129 
2130  if ( (maske(j,i)==3_i2b) & ! floating ice
2131  .and. &
2132  ( (maske(j,i+1)==2_i2b) & ! with
2133  .or.(maske(j,i-1)==2_i2b) & ! one
2134  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
2135  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
2136  ) &
2137  flag_calving_front_1(j,i) = .true.
2138 
2139  if ( (maske(j,i)==2_i2b) & ! ocean
2140  .and. &
2141  ( (maske(j,i+1)==3_i2b) & ! with
2142  .or.(maske(j,i-1)==3_i2b) & ! one
2143  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
2144  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
2145  ) &
2146  flag_calving_front_2(j,i) = .true.
2147 
2148 end do
2149 end do
2150 
2151 do i=1, imax-1
2152 
2153  j=0
2154  if ( (maske(j,i)==2_i2b) & ! ocean
2155  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
2156  ) & ! floating ice point
2157  flag_calving_front_2(j,i) = .true.
2158 
2159  j=jmax
2160  if ( (maske(j,i)==2_i2b) & ! ocean
2161  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
2162  ) & ! floating ice point
2163  flag_calving_front_2(j,i) = .true.
2164 
2165 end do
2166 
2167 do j=1, jmax-1
2168 
2169  i=0
2170  if ( (maske(j,i)==2_i2b) & ! ocean
2171  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
2172  ) & ! floating ice point
2173  flag_calving_front_2(j,i) = .true.
2174 
2175  i=imax
2176  if ( (maske(j,i)==2_i2b) & ! ocean
2177  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
2178  ) & ! floating ice point
2179  flag_calving_front_2(j,i) = .true.
2180 
2181 end do
2182 
2183 #else
2184 errormsg = ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
2185 call error(errormsg)
2186 #endif
2187 
2188 !-------- Rock/sediment mask: smoothing over the transition length --------
2189 
2190 #if (defined(SEDI_SLIDE) && SEDI_SLIDE==2 && defined(TRANS_LENGTH_SL))
2191 
2192 do i=0, imax
2193 do j=0, jmax
2194  if (maske_sedi(j,i) /= 7_i2b) then ! no soft sediment
2195  r_mask_sedi(j,i) = 0.0_dp
2196  else ! maske_sedi(j,i) == 7_i2b, soft sediment
2197  r_mask_sedi(j,i) = 1.0_dp
2198  end if
2199 end do
2200 end do
2201 
2202 if (trans_length_sl > eps) then
2203 
2204  if ((p_weert /= p_weert_sedi).or.(q_weert /= q_weert_sedi)) then
2205  errormsg = ' >>> sico_init: Specifying a transition length' &
2206  // end_of_line &
2207  //' TRANS_LENGTH_SL > 0 for the basal sliding' &
2208  // end_of_line &
2209  //' regimes requires P_WEERT==P_WEERT_SEDI' &
2210  // end_of_line &
2211  //' and Q_WEERT==Q_WEERT_SEDI!'
2212  call error(errormsg)
2213  end if
2214 
2215  transition_length_sliding = trans_length_sl *1.0e+03_dp ! km -> m
2216 
2217  quasi_time = 0.25_dp*transition_length_sliding**2
2218  ! Time scale resulting from the fundamental solution
2219  ! of the diffusion equation for diffusitivity = 1
2220  d_quasi_time = 0.1_dp*quasi_time
2221  ! First guess for the time step
2222 
2223 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
2224 
2225  do
2226  if ( d_quasi_time < 0.1_dp*min(dxi,deta)**2 ) exit
2227  ! CFL condition fulfilled (with some leeway)
2228  d_quasi_time = 0.5_dp*d_quasi_time
2229  end do
2230 
2231 #else /* OpenAD */
2232 
2233  do while ( d_quasi_time >= 0.1_dp*min(dxi,deta)**2 )
2234  ! CFL condition UN!!fulfilled (with some leeway)
2235  d_quasi_time = 0.5_dp*d_quasi_time
2236  end do
2237 
2238 #endif /* Normal vs. OpenAD */
2239 
2240  m = nint(quasi_time/d_quasi_time)
2241 
2242  do n=1, m
2243 
2244  r_mask_sedi_new = r_mask_sedi
2245 
2246  do i=1, imax-1
2247  do j=1, jmax-1
2248  r_mask_sedi_new(j,i) = r_mask_sedi(j,i) &
2249  + d_quasi_time/dxi**2 &
2250  * (r_mask_sedi(j,i+1)-2.0_dp*r_mask_sedi(j,i)+r_mask_sedi(j,i-1)) &
2251  + d_quasi_time/deta**2 &
2252  * (r_mask_sedi(j+1,i)-2.0_dp*r_mask_sedi(j,i)+r_mask_sedi(j-1,i))
2253  ! Solving the diffusion equation for diffusitivity = 1
2254  end do
2255  end do
2256 
2257  r_mask_sedi = r_mask_sedi_new
2258 
2259  end do
2260 
2261 end if
2262 
2263 #endif
2264 
2265 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
2266 
2267 #if (GRID==0 || GRID==1)
2268 
2269 do ir=-imax, imax
2270 do jr=-jmax, jmax
2271  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
2272  ! distortion due to stereographic projection not accounted for
2273 end do
2274 end do
2275 
2276 #endif
2277 
2278 !-------- Initial velocities --------
2279 
2280 call calc_temp_melt()
2281 
2282 #if (DYNAMICS==1 || DYNAMICS==2)
2283 
2284 call calc_vxy_b_sia(time, z_sl)
2285 call calc_vxy_sia(dzeta_c, dzeta_t)
2286 
2287 #if (MARGIN==3 || DYNAMICS==2)
2288 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
2289 #endif
2290 
2291 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
2292 
2293 #if (MARGIN==3)
2294 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
2295 #endif
2296 
2297 #elif (DYNAMICS==0)
2298 
2299 call calc_vxy_static()
2300 call calc_vz_static()
2301 
2302 #else
2303  errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
2304  call error(errormsg)
2305 #endif
2306 
2307 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
2308 
2309 !-------- Initial enthalpies --------
2310 
2311 #if (CALCMOD==0 || CALCMOD==-1)
2312 
2313 do i=0, imax
2314 do j=0, jmax
2315 
2316  do kc=0, kcmax
2317  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2318  end do
2319 
2320  do kt=0, ktmax
2321  enth_t(kt,j,i) = enth_c(0,j,i)
2322  end do
2323 
2324 end do
2325 end do
2326 
2327 #elif (CALCMOD==1)
2328 
2329 do i=0, imax
2330 do j=0, jmax
2331 
2332  do kc=0, kcmax
2333  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2334  end do
2335 
2336  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
2337  do kt=0, ktmax
2338  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
2339  end do
2340  else
2341  do kt=0, ktmax
2342  enth_t(kt,j,i) = enth_c(0,j,i)
2343  end do
2344  end if
2345 
2346 end do
2347 end do
2348 
2349 #elif (CALCMOD==2 || CALCMOD==3)
2350 
2351 do i=0, imax
2352 do j=0, jmax
2353 
2354  do kc=0, kcmax
2355  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
2356  end do
2357 
2358  do kt=0, ktmax
2359  enth_t(kt,j,i) = enth_c(0,j,i)
2360  end do
2361 
2362 end do
2363 end do
2364 
2365 #else
2366 
2367 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
2368 call error(errormsg)
2369 
2370 #endif
2371 
2372 !-------- Initialize time-series files --------
2373 
2374 ! ------ Time-series file for the ice sheet on the whole
2375 
2376 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
2377 
2378 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
2379 
2380 open(12, iostat=ios, file=trim(filename_with_path), status='new')
2381 
2382 if (ios /= 0) then
2383  errormsg = ' >>> sico_init: Error when opening the ser file!'
2384  call error(errormsg)
2385 end if
2386 
2387 if (forcing_flag == 1) then
2388 
2389  write(12,1102)
2390  write(12,1103)
2391 
2392  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
2393  ' V(m^3) V_g(m^3) V_f(m^3)', &
2394  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2395  ' V_sle(m) V_t(m^3)', &
2396  ' A_t(m^2)',/, &
2397  ' H_max(m) H_t_max(m)', &
2398  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2399  1103 format('----------------------------------------------------', &
2400  '---------------------------------------')
2401 
2402 else if (forcing_flag == 2) then
2403 
2404  write(12,1112)
2405  write(12,1113)
2406 
2407  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
2408  ' V(m^3) V_g(m^3) V_f(m^3)', &
2409  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2410  ' V_sle(m) V_t(m^3)', &
2411  ' A_t(m^2)',/, &
2412  ' H_max(m) H_t_max(m)', &
2413  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2414  1113 format('----------------------------------------------------', &
2415  '---------------------------------------')
2416 
2417 else if (forcing_flag == 3) then
2418 
2419  write(12,1122)
2420  write(12,1123)
2421 
2422  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2423  ' V(m^3) V_g(m^3) V_f(m^3)', &
2424  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2425  ' V_sle(m) V_t(m^3)', &
2426  ' A_t(m^2)',/, &
2427  ' H_max(m) H_t_max(m)', &
2428  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2429  1123 format('----------------------------------------------------', &
2430  '---------------------------------------')
2431 
2432 end if
2433 
2434 #else /* OpenAD */
2435 
2436 print *, ' >>> sico_init: not writing to the ser file in '
2437 print *, ' adjoint applications!'
2438 
2439 #endif /* Normal vs. OpenAD */
2440 
2441 ! ------ Time-series file for deep boreholes
2442 
2443 #ifndef ALLOW_OPENAD /* Normal */
2444 
2445 n_core = 6 ! Vostok, Dome A, Dome C, Dome F, Kohnen, Byrd
2446 
2447 allocate(lambda_core(n_core), phi_core(n_core), &
2449 
2450 ch_core(1) = 'Vostok'
2451 phi_core(1) = -78.467_dp *pi_180 ! Geographical position of Vostok,
2452 lambda_core(1) = 106.8_dp *pi_180 ! conversion deg -> rad
2453 
2454 ch_core(2) = 'Dome A'
2455 phi_core(2) = -80.367_dp *pi_180 ! Geographical position of Dome A,
2456 lambda_core(2) = 77.35_dp *pi_180 ! conversion deg -> rad
2457 
2458 ch_core(3) = 'Dome C'
2459 phi_core(3) = -75.1_dp *pi_180 ! Geographical position of Dome C,
2460 lambda_core(3) = 123.4_dp *pi_180 ! conversion deg -> rad
2461 
2462 ch_core(4) = 'Dome F'
2463 phi_core(4) = -77.317_dp *pi_180 ! Geographical position of Dome F,
2464 lambda_core(4) = 39.7_dp *pi_180 ! conversion deg -> rad
2465 
2466 ch_core(5) = 'Kohnen'
2467 phi_core(5) = -75.0_dp *pi_180 ! Geographical position of Kohnen,
2468 lambda_core(5) = 0.067_dp *pi_180 ! conversion deg -> rad
2469 
2470 ch_core(6) = 'Byrd'
2471 phi_core(6) = -80.017_dp *pi_180 ! Geographical position of Byrd,
2472 lambda_core(6) = -119.517_dp *pi_180 ! conversion deg -> rad
2473 
2474 #if (GRID==0 || GRID==1) /* Stereographic projection */
2475 
2476 do n=1, n_core
2478  a, b, lambda0, phi0, x_core(n), y_core(n))
2479 end do
2480 
2481 #elif (GRID==2) /* Geographical coordinates */
2482 
2484 y_core = phi_core
2485 
2486 #endif
2487 
2488 #endif /* Normal (OpenAD: No core data for adjoint) */
2489 
2490 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
2491 
2492 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD) ) /* Normal */
2493 open(14, iostat=ios, file=trim(filename_with_path), status='new')
2494 #else /* OpenAD */
2495 open(14, iostat=ios, file=trim(filename_with_path))
2496 #endif /* Normal vs. OpenAD */
2497 
2498 if (ios /= 0) then
2499  errormsg = ' >>> sico_init: Error when opening the core file!'
2500  call error(errormsg)
2501 end if
2502 
2503 if (forcing_flag == 1) then
2504 
2505  write(14,1106)
2506  write(14,1107)
2507 
2508  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
2509  ' H_Vo(m) H_DA(m) H_DC(m)', &
2510  ' H_DF(m) H_Ko(m) H_By(m)',/, &
2511  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
2512  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
2513  ' T_Vo(C) T_DA(C) T_DC(C)', &
2514  ' T_DF(C) T_Ko(C) T_By(C)')
2515  1107 format('----------------------------------------------------', &
2516  '---------------------------------------')
2517 
2518 else if (forcing_flag == 2) then
2519 
2520  write(14,1116)
2521  write(14,1117)
2522 
2523  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
2524  ' H_Vo(m) H_DA(m) H_DC(m)', &
2525  ' H_DF(m) H_Ko(m) H_By(m)',/, &
2526  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
2527  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
2528  ' T_Vo(C) T_DA(C) T_DC(C)', &
2529  ' T_DF(C) T_Ko(C) T_By(C)')
2530  1117 format('----------------------------------------------------', &
2531  '---------------------------------------')
2532 
2533 else if (forcing_flag == 3) then
2534 
2535  write(14,1126)
2536  write(14,1127)
2537 
2538  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2539  ' H_Vo(m) H_DA(m) H_DC(m)', &
2540  ' H_DF(m) H_Ko(m) H_By(m)',/, &
2541  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
2542  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
2543  ' T_Vo(C) T_DA(C) T_DC(C)', &
2544  ' T_DF(C) T_Ko(C) T_By(C)')
2545  1127 format('----------------------------------------------------', &
2546  '---------------------------------------')
2547 
2548 end if
2549 
2550 !-------- Output of the initial state --------
2551 
2552 #ifndef ALLOW_OPENAD /* Normal */
2553 
2554 #if (defined(OUTPUT_INIT))
2555 
2556 #if (OUTPUT_INIT==0)
2557  flag_init_output = .false.
2558 #elif (OUTPUT_INIT==1)
2559  flag_init_output = .true.
2560 #else
2561  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
2562  call error(errormsg)
2563 #endif
2564 
2565 #else
2566  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
2567 #endif
2568 
2569 #if (OUTPUT==1)
2570 
2571 #if (ERGDAT==0)
2572  flag_3d_output = .false.
2573 #elif (ERGDAT==1)
2574  flag_3d_output = .true.
2575 #else
2576  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
2577  call error(errormsg)
2578 #endif
2579 
2580  if (flag_init_output) &
2581  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2582  flag_3d_output, ndat2d, ndat3d)
2583 
2584 #elif (OUTPUT==2)
2585 
2586 if (time_output(1) <= time_init+eps) then
2587 
2588 #if (ERGDAT==0)
2589  flag_3d_output = .false.
2590 #elif (ERGDAT==1)
2591  flag_3d_output = .true.
2592 #else
2593  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
2594  call error(errormsg)
2595 #endif
2596 
2597  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2598  flag_3d_output, ndat2d, ndat3d)
2599 
2600 end if
2601 
2602 #elif (OUTPUT==3)
2603 
2604  flag_3d_output = .false.
2605 
2606  if (flag_init_output) &
2607  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2608  flag_3d_output, ndat2d, ndat3d)
2609 
2610 if (time_output(1) <= time_init+eps) then
2611 
2612  flag_3d_output = .true.
2613 
2614  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2615  flag_3d_output, ndat2d, ndat3d)
2616 
2617 end if
2618 
2619 #else
2620  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
2621  call error(errormsg)
2622 #endif
2623 
2624 if (flag_init_output) then
2625  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2626  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2627 end if
2628 
2629 #else /* OpenAD */
2630 
2631 print *, ' >>> sico_init: not producing initial, typical outputs'
2632 print *, ' in adjoint mode.'
2633 
2634 #endif /* Normal vs. OpenAD */
2635 
2636 end subroutine sico_init
2637 
2638 !-------------------------------------------------------------------------------
2639 !> Definition of the initial surface and bedrock topography
2640 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2641 !! For present-day initial topography.
2642 !<------------------------------------------------------------------------------
2643 subroutine topography1(dxi, deta)
2645 #if (GRID==0 || GRID==1)
2647 #endif
2648 
2649  use metric_m
2650  use topograd_m
2651 
2652 implicit none
2653 
2654 real(dp), intent(out) :: dxi, deta
2655 
2656 integer(i4b) :: i, j, n
2657 integer(i4b) :: ios, n_dummy
2658 real(dp) :: d_dummy
2659 real(dp) :: xi0, eta0
2660 real(dp) :: H_ice, freeboard_ratio
2661 character :: ch_dummy
2662 
2663 character(len= 8) :: ch_imax
2664 character(len=128) :: fmt4
2665 character(len=256) :: filename_with_path
2666 
2667 write(ch_imax, fmt='(i8)') imax
2668 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2669 
2670 !-------- Read topography --------
2671 
2672 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2673  trim(zs_present_file)
2674 
2675 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2676 
2677 if (ios /= 0) then
2678  errormsg = ' >>> topography1: Error when opening the zs file!'
2679  call error(errormsg)
2680 end if
2681 
2682 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
2683 
2684 do j=jmax, 0, -1
2685  read(21, fmt=*) (zs(j,i), i=0,imax)
2686 end do
2687 
2688 close(21, status='keep')
2689 
2690 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2691  trim(zl_present_file)
2692 
2693 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2694 
2695 if (ios /= 0) then
2696  errormsg = ' >>> topography1: Error when opening the zl file!'
2697  call error(errormsg)
2698 end if
2699 
2700 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
2701 
2702 do j=jmax, 0, -1
2703  read(22, fmt=*) (zl(j,i), i=0,imax)
2704 end do
2705 
2706 close(22, status='keep')
2707 
2708 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2709  trim(zl0_file)
2710 
2711 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2712 
2713 if (ios /= 0) then
2714  errormsg = ' >>> topography1: Error when opening the zl0 file!'
2715  call error(errormsg)
2716 end if
2717 
2718 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2719 
2720 do j=jmax, 0, -1
2721  read(23, fmt=*) (zl0(j,i), i=0,imax)
2722 end do
2723 
2724 close(23, status='keep')
2725 
2726 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2727  trim(mask_present_file)
2728 
2729 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2730 
2731 if (ios /= 0) then
2732  errormsg = ' >>> topography1: Error when opening the mask file!'
2733  call error(errormsg)
2734 end if
2735 
2736 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2737 
2738 do j=jmax, 0, -1
2739  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2740 end do
2741 
2742 close(24, status='keep')
2743 
2744 #if (defined(ZB_PRESENT_FILE))
2745 
2746 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2747  trim(zb_present_file)
2748 
2749 open(25, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2750 
2751 if (ios /= 0) then
2752  errormsg = ' >>> topography1: Error when opening the zb file!'
2753  call error(errormsg)
2754 end if
2755 
2756 do n=1, 6; read(25, fmt='(a)') ch_dummy; end do
2757 
2758 do j=jmax, 0, -1
2759  read(25, fmt=*) (zb(j,i), i=0,imax)
2760 end do
2761 
2762 close(25, status='keep')
2763 
2764 #else
2765 
2766 write(6, fmt='(a)') ' >>> topography1: ZB_PRESENT_FILE not defined,'
2767 write(6, fmt='(a)') ' thus zb = zl assumed.'
2768 
2769 zb = zl
2770 
2771 #endif
2772 
2773 !-------- Further stuff --------
2774 
2775 dxi = dx *1000.0_dp ! km -> m
2776 deta = dx *1000.0_dp ! km -> m
2777 
2778 xi0 = x0 *1000.0_dp ! km -> m
2779 eta0 = y0 *1000.0_dp ! km -> m
2780 
2781 freeboard_ratio = (rho_sw-rho)/rho_sw
2782 
2783 do i=0, imax
2784 do j=0, jmax
2785 
2786  if (maske(j,i) <= 1_i2b) then
2787 
2788  zb(j,i) = zl(j,i) ! ensure consistency
2789 
2790  if (zs(j,i) < zb(j,i)) then
2791  zs(j,i) = zb(j,i)
2792  maske(j,i) = 1_i2b
2793  end if
2794 
2795  else if (maske(j,i) == 2_i2b) then
2796 
2797 #if (MARGIN==1 || MARGIN==2)
2798  zs(j,i) = zl(j,i) ! ensure
2799  zb(j,i) = zl(j,i) ! consistency
2800 #elif (MARGIN==3)
2801  zs(j,i) = 0.0_dp ! present-day
2802  zb(j,i) = 0.0_dp ! sea level
2803 #endif
2804 
2805  else if (maske(j,i) == 3_i2b) then
2806 
2807 #if (MARGIN==1 || (MARGIN==2 && MARINE_ICE_FORMATION==1))
2808  maske(j,i) = 2_i2b ! floating ice cut off
2809  zs(j,i) = zl(j,i)
2810  zb(j,i) = zl(j,i)
2811 #elif (MARGIN==2 && MARINE_ICE_FORMATION==2)
2812  maske(j,i) = 0_i2b ! floating ice becomes "underwater ice"
2813  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2814  if (h_ice < 0.0_dp) then
2815  h_ice = 0.0_dp
2816  maske(j,i) = 2_i2b
2817  end if
2818  zs(j,i) = zl(j,i)+h_ice
2819  zb(j,i) = zl(j,i)
2820 #elif (MARGIN==3)
2821  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2822  if (h_ice < 0.0_dp) then
2823  h_ice = 0.0_dp
2824  maske(j,i) = 2_i2b
2825  end if
2826  zs(j,i) = freeboard_ratio*h_ice ! ensure properly
2827  zb(j,i) = zs(j,i)-h_ice ! floating ice
2828 #endif
2829 
2830  end if
2831 
2832  xi(i) = xi0 + real(i,dp)*dxi
2833  eta(j) = eta0 + real(j,dp)*deta
2834 
2835  zm(j,i) = zb(j,i)
2836  n_cts(j,i) = -1_i2b
2837  kc_cts(j,i) = 0_i2b
2838 
2839  h_c(j,i) = zs(j,i)-zm(j,i)
2840  h_t(j,i) = 0.0_dp
2841 
2842  dzs_dtau(j,i) = 0.0_dp
2843  dzm_dtau(j,i) = 0.0_dp
2844  dzb_dtau(j,i) = 0.0_dp
2845  dzl_dtau(j,i) = 0.0_dp
2846  dh_c_dtau(j,i) = 0.0_dp
2847  dh_t_dtau(j,i) = 0.0_dp
2848 
2849 end do
2850 end do
2851 
2852 maske_old = maske
2853 
2854 !-------- Geographic coordinates, metric tensor,
2855 ! gradients of the topography --------
2856 
2857 do i=0, imax
2858 do j=0, jmax
2859 
2860 #if (GRID==0 || GRID==1) /* Stereographic projection */
2861 
2862  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2863  lambda0, phi0, lambda(j,i), phi(j,i))
2864 
2865 #elif (GRID==2) /* Geographic coordinates */
2866 
2867  lambda(j,i) = xi(i)
2868  phi(j,i) = eta(j)
2869 
2870 #endif
2871 
2872 end do
2873 end do
2874 
2875 call metric()
2876 
2877 #if (TOPOGRAD==0)
2878 call topograd_1(dxi, deta, 1)
2879 #elif (TOPOGRAD==1)
2880 call topograd_2(dxi, deta, 1)
2881 #endif
2882 
2883 !-------- Corresponding area of grid points --------
2884 
2885 do i=0, imax
2886 do j=0, jmax
2887  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2888 end do
2889 end do
2890 
2891 end subroutine topography1
2892 
2893 !-------------------------------------------------------------------------------
2894 !> Definition of the initial surface and bedrock topography
2895 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2896 !! For ice-free initial topography with relaxed lithosphere.
2897 !<------------------------------------------------------------------------------
2898 subroutine topography2(dxi, deta)
2900 #if (GRID==0 || GRID==1)
2902 #endif
2903 
2904  use metric_m
2905  use topograd_m
2906 
2907 implicit none
2908 
2909 real(dp), intent(out) :: dxi, deta
2910 
2911 integer(i4b) :: i, j, n
2912 integer(i4b) :: ios
2913 real(dp) :: xi0, eta0
2914 character :: ch_dummy
2915 
2916 character(len= 8) :: ch_imax
2917 character(len=128) :: fmt4
2918 character(len=256) :: filename_with_path
2919 
2920 write(ch_imax, fmt='(i8)') imax
2921 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2922 
2923 !-------- Read topography --------
2924 
2925 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2926  trim(zl0_file)
2927 
2928 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2929 
2930 if (ios /= 0) then
2931  errormsg = ' >>> topography2: Error when opening the zl0 file!'
2932  call error(errormsg)
2933 end if
2934 
2935 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2936 
2937 do j=jmax, 0, -1
2938  read(23, fmt=*) (zl0(j,i), i=0,imax)
2939 end do
2940 
2941 close(23, status='keep')
2942 
2943 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2944  trim(mask_present_file)
2945 
2946 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2947 
2948 if (ios /= 0) then
2949  errormsg = ' >>> topography2: Error when opening the mask file!'
2950  call error(errormsg)
2951 end if
2952 
2953 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2954 
2955 do j=jmax, 0, -1
2956  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2957 end do
2958 
2959 close(24, status='keep')
2960 
2961 !-------- Further stuff --------
2962 
2963 dxi = dx *1000.0_dp ! km -> m
2964 deta = dx *1000.0_dp ! km -> m
2965 
2966 xi0 = x0 *1000.0_dp ! km -> m
2967 eta0 = y0 *1000.0_dp ! km -> m
2968 
2969 do i=0, imax
2970 do j=0, jmax
2971 
2972  if (maske(j,i) <= 1_i2b) then
2973  maske(j,i) = 1_i2b
2974  zs(j,i) = zl0(j,i)
2975  zb(j,i) = zl0(j,i)
2976  zl(j,i) = zl0(j,i)
2977  else ! (maske(j,i) >= 2_i2b)
2978  maske(j,i) = 2_i2b
2979 #if (MARGIN==1 || MARGIN==2)
2980  zs(j,i) = zl0(j,i)
2981  zb(j,i) = zl0(j,i)
2982 #elif (MARGIN==3)
2983  zs(j,i) = 0.0_dp ! present-day
2984  zb(j,i) = 0.0_dp ! sea level
2985 #endif
2986  zl(j,i) = zl0(j,i)
2987  end if
2988 
2989  xi(i) = xi0 + real(i,dp)*dxi
2990  eta(j) = eta0 + real(j,dp)*deta
2991 
2992  zm(j,i) = zb(j,i)
2993  n_cts(j,i) = -1_i2b
2994  kc_cts(j,i) = 0_i2b
2995 
2996  h_c(j,i) = 0.0_dp
2997  h_t(j,i) = 0.0_dp
2998 
2999  dzs_dtau(j,i) = 0.0_dp
3000  dzm_dtau(j,i) = 0.0_dp
3001  dzb_dtau(j,i) = 0.0_dp
3002  dzl_dtau(j,i) = 0.0_dp
3003  dh_c_dtau(j,i) = 0.0_dp
3004  dh_t_dtau(j,i) = 0.0_dp
3005 
3006 end do
3007 end do
3008 
3009 maske_old = maske
3010 
3011 !-------- Geographic coordinates, metric tensor,
3012 ! gradients of the topography --------
3013 
3014 do i=0, imax
3015 do j=0, jmax
3016 
3017 #if (GRID==0 || GRID==1) /* Stereographic projection */
3018 
3019  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
3020  lambda0, phi0, lambda(j,i), phi(j,i))
3021 
3022 #elif (GRID==2) /* Geographic coordinates */
3023 
3024  lambda(j,i) = xi(i)
3025  phi(j,i) = eta(j)
3026 
3027 #endif
3028 
3029 end do
3030 end do
3031 
3032 call metric()
3033 
3034 #if (TOPOGRAD==0)
3035 call topograd_1(dxi, deta, 1)
3036 #elif (TOPOGRAD==1)
3037 call topograd_2(dxi, deta, 1)
3038 #endif
3039 
3040 !-------- Corresponding area of grid points --------
3041 
3042 do i=0, imax
3043 do j=0, jmax
3044  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
3045 end do
3046 end do
3047 
3048 end subroutine topography2
3049 
3050 !-------------------------------------------------------------------------------
3051 !> Definition of the initial surface and bedrock topography
3052 !! (including gradients) and of the horizontal grid spacings dxi, deta.
3053 !! For initial topography from previous simulation.
3054 !<------------------------------------------------------------------------------
3055 subroutine topography3(dxi, deta, z_sl, anfdatname)
3057  use read_m, only : read_erg_nc
3058 
3059 #if (GRID==0 || GRID==1)
3061 #endif
3062 
3063  use metric_m
3064  use topograd_m
3065 
3066 implicit none
3067 
3068 character(len=100), intent(in) :: anfdatname
3069 
3070 real(dp), intent(out) :: dxi, deta, z_sl
3071 
3072 integer(i4b) :: i, j, n
3073 integer(i4b) :: ios
3074 character(len=256) :: filename_with_path
3075 character :: ch_dummy
3076 
3077 !-------- Read data from time-slice file of previous simulation --------
3078 
3079 call read_erg_nc(z_sl, anfdatname)
3080 
3081 !-------- Read topography of the relaxed bedrock --------
3082 
3083 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
3084  trim(zl0_file)
3085 
3086 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
3087 
3088 if (ios /= 0) then
3089  errormsg = ' >>> topography3: Error when opening the zl0 file!'
3090  call error(errormsg)
3091 end if
3092 
3093 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
3094 
3095 do j=jmax, 0, -1
3096  read(23, fmt=*) (zl0(j,i), i=0,imax)
3097 end do
3098 
3099 close(23, status='keep')
3100 
3101 !-------- Further stuff --------
3102 
3103 dxi = dx *1000.0_dp ! km -> m
3104 deta = dx *1000.0_dp ! km -> m
3105 
3106 !-------- Geographic coordinates, metric tensor,
3107 ! gradients of the topography --------
3108 
3109 do i=0, imax
3110 do j=0, jmax
3111 
3112 #if (GRID==0 || GRID==1) /* Stereographic projection */
3113 
3114  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
3115  lambda0, phi0, lambda(j,i), phi(j,i))
3116 
3117 #elif (GRID==2) /* Geographic coordinates */
3118 
3119  lambda(j,i) = xi(i)
3120  phi(j,i) = eta(j)
3121 
3122 #endif
3123 
3124 end do
3125 end do
3126 
3127 call metric()
3128 
3129 #if (TOPOGRAD==0)
3130 call topograd_1(dxi, deta, 1)
3131 #elif (TOPOGRAD==1)
3132 call topograd_2(dxi, deta, 1)
3133 #endif
3134 
3135 !-------- Corresponding area of grid points --------
3136 
3137 do i=0, imax
3138 do j=0, jmax
3139  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
3140 end do
3141 end do
3142 
3143 end subroutine topography3
3144 
3145 !-------------------------------------------------------------------------------
3146 
3147 end module sico_init_m
3148 !
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:497
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...