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