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