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