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