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