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