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