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