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