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
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, ios1, ios2, ios3, ios4
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 !-------- Some initial values --------
80 
81 n_output = 0
82 
83 dtime = 0.0_dp
84 dtime_temp = 0.0_dp
85 dtime_wss = 0.0_dp
86 dtime_out = 0.0_dp
87 dtime_ser = 0.0_dp
88 
89 time = 0.0_dp
90 time_init = 0.0_dp
91 time_end = 0.0_dp
92 time_output = 0.0_dp
93 
94 !-------- Initialisation of the Library of Iterative Solvers Lis,
95 ! if required --------
96 
97 #if (CALCZS==4 || MARGIN==3)
98 
99 #include "lisf.h"
100 call lis_initialize(ierr)
101 
102 #endif
103 
104 !-------- Read physical parameters --------
105 
106 call phys_para()
107 
108 temp_min = temp_min
109 s_t = s_t *1.0e-03_dp ! K/km -> K/m
110 x_hat = x_hat *1.0e+03_dp ! km -> m
111 y_hat = y_hat *1.0e+03_dp ! km -> m
112 b_max = b_max /year_sec ! m/a -> m/s
113 s_b = s_b *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
114 eld = eld *1.0e+03_dp ! km -> m
115 
116 !-------- Compatibility check of the SICOPOLIS version with the header file
117 
118 if ( trim(version) /= trim(sico_version) ) &
119  stop ' sico_init: SICOPOLIS version not compatible with header file!'
120 
121 !-------- Compatibility check of the horizontal resolution with the
122 ! number of grid points --------
123 
124 #if GRID==0
125 
126 if (abs(dx-25.0_dp) < eps) then
127 
128  if ((imax /= 60).or.(jmax /= 60)) then
129  stop ' sico_init: IMAX and/or JMAX wrong!'
130  end if
131 
132 else if (abs(dx-75.0_dp) < eps) then
133 
134  if ((imax /= 20).or.(jmax /= 20)) then
135  stop ' sico_init: IMAX and/or JMAX wrong!'
136  end if
137 
138 else
139  stop ' sico_init: DX wrong!'
140 end if
141 
142 #elif GRID==1
143 
144 stop ' sico_init: GRID==1 not allowed for this application!'
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 discretization schemes for the horizontal and
153 ! vertical advection terms in the temperature and age equations --------
154 
155 #if ADV_HOR==1
156 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
157 #endif
158 
159 !-------- Setting of forcing flag --------
160 
161 forcing_flag = 1 ! forcing by delta_ts
162 
163 !-------- Initialization of numerical time steps --------
164 
165 dtime0 = dtime0
166 dtime_temp0 = dtime_temp0
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 <= imax).and.(j <= 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, fmt=trim(fmt1)) 'Computational domain:'
254 #if defined(ANT)
255 write(10, fmt=trim(fmt1)) 'Antarctica'
256 #elif defined(GRL)
257 write(10, fmt=trim(fmt1)) 'Greenland'
258 #elif defined(SCAND)
259 write(10, fmt=trim(fmt1)) 'Scandinavia and Eurasia'
260 #elif defined(NHEM)
261 write(10, fmt=trim(fmt1)) 'Northern hemisphere'
262 #elif defined(EMTP2SGE)
263 write(10, fmt=trim(fmt1)) 'EISMINT Phase 2 Simplified Geometry Experiment'
264 #elif defined(NMARS)
265 write(10, fmt=trim(fmt1)) 'North polar cap of Mars'
266 #elif defined(SMARS)
267 write(10, fmt=trim(fmt1)) 'South polar cap of Mars'
268 #else
269 stop ' sico_init: No valid domain specified!'
270 #endif
271 write(10, fmt=trim(fmt1)) ' '
272 
273 write(10, fmt=trim(fmt2)) 'imax =', imax
274 write(10, fmt=trim(fmt2)) 'jmax =', jmax
275 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
276 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
277 write(10, fmt=trim(fmt2)) 'krmax =', krmax
278 write(10, fmt=trim(fmt1)) ' '
279 
280 write(10, fmt=trim(fmt3)) 'a =', deform
281 write(10, fmt=trim(fmt1)) ' '
282 
283 #if (GRID==0 || GRID==1)
284 write(10, fmt=trim(fmt3)) 'x0 =', x0
285 write(10, fmt=trim(fmt3)) 'y0 =', y0
286 write(10, fmt=trim(fmt3)) 'dx =', dx
287 #elif GRID==2
288 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
289 write(10, fmt=trim(fmt3)) 'dphi =', dphi
290 #endif
291 write(10, fmt=trim(fmt1)) ' '
292 
293 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
294 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
295 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
296 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
297 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
298 #if ( OUTPUT==1 || OUTPUT==3 )
299 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
300 #endif
301 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
302 write(10, fmt=trim(fmt1)) ' '
303 
304 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
305 write(10, fmt=trim(fmt1)) ' '
306 
307 #if (CALCZS==3 || CALCZS==4)
308 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
309 #if CALCZS==3
310 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
311 #endif
312 write(10, fmt=trim(fmt1)) ' '
313 #endif
314 
315 write(10, fmt=trim(fmt3)) 'temp_min =', temp_min
316 write(10, fmt=trim(fmt3)) 's_t =', s_t
317 write(10, fmt=trim(fmt3)) 'x_hat =', x_hat
318 write(10, fmt=trim(fmt3)) 'y_hat =', y_hat
319 write(10, fmt=trim(fmt3)) 'b_max =', b_max
320 write(10, fmt=trim(fmt3)) 's_b =', s_b
321 write(10, fmt=trim(fmt3)) 'eld =', eld
322 write(10, fmt=trim(fmt1)) ' '
323 
324 #if TSURFACE==1
325 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
326 #elif TSURFACE==3
327 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
328 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
329 #elif TSURFACE==4
330 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
331 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
332 #endif
333 
334 #if SEA_LEVEL==1
335 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
336 #elif SEA_LEVEL==3
337 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
338 #endif
339 write(10, fmt=trim(fmt1)) ' '
340 
341 #if MARGIN==2
342 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
343 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
344 write(10, fmt=trim(fmt1)) ' '
345 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
346  || marine_ice_calving==6 || marine_ice_calving==7 )
347 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
348 write(10, fmt=trim(fmt1)) ' '
349 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
350 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
351 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
352 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
353 write(10, fmt=trim(fmt1)) ' '
354 #endif
355 #elif MARGIN==3
356 #if ICE_SHELF_CALVING==2
357 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
358 write(10, fmt=trim(fmt1)) ' '
359 #endif
360 #endif
361 
362 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
363 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
364 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
365 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
366 #if SLIDE_LAW==2
367 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
368 #endif
369 write(10, fmt=trim(fmt1)) ' '
370 
371 #if Q_GEO_MOD==1
372 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
373 write(10, fmt=trim(fmt1)) ' '
374 #endif
375 
376 #if REBOUND==1
377 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
378 write(10, fmt=trim(fmt1)) ' '
379 #endif
380 
381 #if FLOW_LAW==2
382 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
383 write(10, fmt=trim(fmt1)) ' '
384 #endif
385 #if FIN_VISC==2
386 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
387 write(10, fmt=trim(fmt1)) ' '
388 #endif
389 
390 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
391 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
392 #if ( ENHMOD==2 || ENHMOD==3 )
393 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
394 #endif
395 #if ENHMOD==2
396 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
397 #endif
398 #if ENHMOD==3
399 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
400 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
401 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
402 #endif
403 #if MARGIN==3
404 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
405 #endif
406 write(10, fmt=trim(fmt1)) ' '
407 
408 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
409 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
410 write(10, fmt=trim(fmt3)) 'H_min =', h_min
411 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
412 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
413 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
414 #if defined(QBM_MIN)
415 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
416 #elif defined(QB_MIN) /* obsolete */
417 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
418 #endif
419 #if defined(QBM_MAX)
420 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
421 #elif defined(QB_MAX) /* obsolete */
422 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
423 #endif
424 write(10, fmt=trim(fmt3)) 'age_min =', age_min
425 write(10, fmt=trim(fmt3)) 'age_max =', age_max
426 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
427 #if ADV_VERT==1
428 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
429 #endif
430 write(10, fmt=trim(fmt1)) ' '
431 
432 write(10, fmt=trim(fmt2)) 'GRID =', grid
433 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
434 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
435 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
436 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
437 #if MARGIN==2
438 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
439 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
440 #elif MARGIN==3
441 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
442 #endif
443 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
444 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
445 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
446 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
447 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
448 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
449 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
450 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
451 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
452 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
453 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
454 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
455 write(10, fmt=trim(fmt1)) ' '
456 
457 #if defined(OUT_TIMES)
458 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
459 #endif
460 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
461 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
462 #if defined(OUTSER_V_A)
463 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
464 #endif
465 #if ( OUTPUT==1 || OUTPUT==2 )
466 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
467 #endif
468 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
469 write(10, fmt=trim(fmt1)) ' '
470 #if ( OUTPUT==2 || OUTPUT==3 )
471 write(10, fmt=trim(fmt2)) 'n_output =', n_output
472 do n=1, n_output
473  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
474 end do
475 write(10, fmt=trim(fmt1)) ' '
476 #endif
477 
478 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
479 
480 close(10, status='keep')
481 
482 !-------- Conversion of time quantities --------
483 
484 year_zero = year_zero*year_sec ! a --> s
485 time_init = time_init0*year_sec ! a --> s
486 time_end = time_end0*year_sec ! a --> s
487 dtime = dtime0*year_sec ! a --> s
488 dtime_temp = dtime_temp0*year_sec ! a --> s
489 dtime_ser = dtime_ser0*year_sec ! a --> s
490 #if ( OUTPUT==1 || OUTPUT==3 )
491 dtime_out = dtime_out0*year_sec ! a --> s
492 #endif
493 #if ( OUTPUT==2 || OUTPUT==3 )
494 do n=1, n_output
495  time_output(n) = time_output0(n)*year_sec ! a --> s
496 end do
497 #endif
498 
499 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
500  stop ' sico_init: dtime_temp must be a multiple of dtime!'
501 
502 time = time_init
503 
504 !-------- Mean accumulation --------
505 
506 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
507 ! ! mm/a water equiv. --> m/s ice equiv.
508 
509 mean_accum_inv = 1.0_dp/mean_accum
510 
511 !-------- Set present topography mask --------
512 
513 maske_help = 1
514 
515 maske_help(0,:) = 2
516 maske_help(jmax,:) = 2
517 
518 maske_help(:,0) = 2
519 maske_help(:,imax) = 2
520 
521 !-------- Read data for delta_ts --------
522 
523 #if TSURFACE==4
524 
525 open(21, iostat=ios, &
526  file=inpath//'/general/'//grip_temp_file, &
527  status='old')
528 if (ios /= 0) &
529  stop ' sico_init: Error when opening the data file for delta_ts!'
530 
531 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
532 
533 if (ch_dummy /= '#') then
534  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
535  write(6, fmt=*) ' not defined in data file for delta_ts!'
536  stop
537 end if
538 
539 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
540 
541 allocate(griptemp(0:ndata_grip))
542 
543 do n=0, ndata_grip
544  read(21, fmt=*) d_dummy, griptemp(n)
545 end do
546 
547 close(21, status='keep')
548 
549 #endif
550 
551 !-------- Read data for z_sl --------
552 
553 #if SEA_LEVEL==3
554 
555 open(21, iostat=ios, &
556  file=inpath//'/general/'//sea_level_file, &
557  status='old')
558 if (ios /= 0) &
559  stop ' sico_init: Error when opening the data file for z_sl!'
560 
561 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
562 
563 if (ch_dummy /= '#') then
564  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
565  write(6, fmt=*) ' not defined in data file for z_sl!'
566  stop
567 end if
568 
569 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
570 
571 allocate(specmap_zsl(0:ndata_specmap))
572 
573 do n=0, ndata_specmap
574  read(21, fmt=*) d_dummy, specmap_zsl(n)
575 end do
576 
577 close(21, status='keep')
578 
579 #endif
580 
581 !-------- Determination of the geothermal heat flux --------
582 
583 #if Q_GEO_MOD==1
584 
585 ! ------ Constant value
586 
587 do i=0, imax
588 do j=0, jmax
589  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
590 end do
591 end do
592 
593 #elif Q_GEO_MOD==2
594 
595 stop ' sico_init: Option Q_GEO_MOD==2 not available for this application!'
596 
597 #endif
598 
599 !-------- Definition of initial values --------
600 
601 ! ------ Present topography
602 
603 #if ANF_DAT==1
604 
605 stop ' sico_init: ANF_DAT==1 not allowed for this application!'
606 
607 ! ------ Ice-free, relaxed bedrock
608 
609 #elif ANF_DAT==2
610 
611 call topography2(dxi, deta)
612 
613 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
614 
615 call boundary(time_init, dtime, dxi, deta, &
616  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
617 
618 do i=0, imax
619 do j=0, jmax
620 
621  q_bm(j,i) = 0.0_dp
622  q_tld(j,i) = 0.0_dp
623  q_b_tot(j,i) = 0.0_dp
624 
625  p_b_w(j,i) = 0.0_dp
626  h_w(j,i) = 0.0_dp
627 
628  do kc=0, kcmax
629  temp_c(kc,j,i) = temp_s(j,i)
630  age_c(kc,j,i) = 15000.0_dp*year_sec
631  end do
632 
633  do kt=0, ktmax
634  omega_t(kt,j,i) = 0.0_dp
635  age_t(kt,j,i) = 15000.0_dp*year_sec
636  end do
637 
638  do kr=0, krmax
639  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
640  *(1.0_dp-real(kr,dp)/real(krmax,dp))
641 ! (linear temperature distribution according to the
642 ! geothermal heat flux)
643  end do
644 
645 end do
646 end do
647 
648 call calc_enhance(time_init)
649 
650 ! ------ Read initial state from output of previous simulation
651 
652 #elif ANF_DAT==3
653 
654 #if NETCDF==1
655 call topography3(dxi, deta, z_sl, anfdatname)
656 #elif NETCDF==2
657 call topography3_nc(dxi, deta, z_sl, anfdatname)
658 #else
659 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
660 #endif
661 
662 call boundary(time_init, dtime, dxi, deta, &
663  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
664 
665 q_b_tot = q_bm + q_tld
666 
667 call calc_enhance(time_init)
668 
669 #endif
670 
671 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
672 
673 #if (GRID==0 || GRID==1)
674 
675 do ir=-imax, imax
676 do jr=-jmax, jmax
677  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
678  ! distortion due to stereographic projection not accounted for
679 end do
680 end do
681 
682 #endif
683 
684 !-------- General abbreviations --------
685 
686 ea = exp(deform)
687 
688 do kc=0, kcmax
689  zeta_c(kc) = kc*dzeta_c
690  eaz_c(kc) = exp(deform*zeta_c(kc))
691  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
692 end do
693 
694 do kt=0, ktmax
695  zeta_t(kt) = kt*dzeta_t
696 end do
697 
698 !-------- Initial velocities --------
699 
700 call calc_temp_melt()
701 
702 call calc_vxy_b_sia(time, z_sl)
703 call calc_vxy_sia(dzeta_c, dzeta_t)
704 
705 #if MARGIN==3
706 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
707 #endif
708 
709 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
710 
711 #if MARGIN==3
712 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
713 #endif
714 
715 !-------- Initialize time-series files --------
716 
717 ! ------ Time-series file for the ice sheet on the whole
718 
719 open(12, iostat=ios, &
720  file=outpath//'/'//trim(runname)//'.ser', &
721  status='new')
722 if (ios /= 0) &
723  stop ' sico_init: Error when opening the ser file!'
724 
725 write(12,1102)
726 write(12,1103)
727 
728 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
729 
730  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
731  ' H_max(m) zs_max(m) V_g(m^3)', &
732  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
733  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
734  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
735  1103 format('----------------------------------------------------', &
736  '---------------------------------------')
737 
738 #elif OUTSER_V_A==2
739 
740  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
741  ' V(m^3) V_g(m^3) V_f(m^3)', &
742  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
743  ' H_max(m) zs_max(m)', &
744  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
745  ' A_t(m^2) V_bm(m^3/a)', &
746  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
747  1103 format('----------------------------------------------------', &
748  '---------------------------------------')
749 
750 #else
751  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
752 #endif
753 
754 ! ------ Time-series file for the ice-thickness field
755 
756 #if OUTSER==2
757 
758 open(13, iostat=ios, &
759  file=outpath//'/'//trim(runname)//'.thk', &
760  status='new')
761 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
762 
763 write(13,1104)
764 write(13,1105)
765 
766 1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
767  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
768  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
769 1105 format('----------------------------------------------------')
770 
771 #endif
772 
773 ! ------ Time-series file for selected positions ("deep boreholes")
774 
775 #if OUTSER==3
776 
777 n_core = 2 ! central dome, position halfway to coast
778 
779 allocate(lambda_core(n_core), phi_core(n_core), &
780  x_core(n_core), y_core(n_core))
781 
782 lambda_core(1) = 0.0_dp ! dummy
783 phi_core(1) = 0.0_dp ! dummy
784 x_core(1) = 750.0_dp *1.0e+03_dp ! Position of the central dome
785 y_core(1) = 750.0_dp *1.0e+03_dp ! (750 km, 750 km),
786  ! conversion km -> m
787 
788 lambda_core(2) = 0.0_dp ! dummy
789 phi_core(2) = 0.0_dp ! dummy
790 x_core(2) = 750.0_dp *1.0e+03_dp ! Position halfway to the coast
791 y_core(2) = 1125.0_dp *1.0e+03_dp ! (750 km, 1125 km),
792  ! conversion km -> m
793 
794 open(14, iostat=ios, &
795  file=outpath//'/'//trim(runname)//'.core', &
796  status='new')
797 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
798 
799 if (forcing_flag == 1) then
800 
801  write(14,1106)
802  write(14,1107)
803 
804  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
805  ' H_P1(m) H_P2(m)',/, &
806  ' v_P1(m/a) v_P2(m/a)',/, &
807  ' T_P1(C) T_P2(C)',/, &
808  ' Rb_P1(W/m2) Rb_P2(W/m2)')
809  1107 format('---------------------------------------')
810 
811 end if
812 
813 #endif
814 
815 !-------- Output of the initial state --------
816 
817 #if OUTPUT==1
818 
819 #if ERGDAT==0
820  flag_3d_output = .false.
821 #elif ERGDAT==1
822  flag_3d_output = .true.
823 #else
824  stop ' sico_init: ERGDAT must be either 0 or 1!'
825 #endif
826 
827 #if NETCDF==1
828  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
829  flag_3d_output, ndat2d, ndat3d)
830 #elif NETCDF==2
831  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
832  flag_3d_output, ndat2d, ndat3d)
833 #else
834  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
835 #endif
836 
837 #elif OUTPUT==2
838 
839 if (time_output(1) <= time_init+eps) then
840 
841 #if ERGDAT==0
842  flag_3d_output = .false.
843 #elif ERGDAT==1
844  flag_3d_output = .true.
845 #else
846  stop ' sico_init: ERGDAT must be either 0 or 1!'
847 #endif
848 
849 #if NETCDF==1
850  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
851  flag_3d_output, ndat2d, ndat3d)
852 #elif NETCDF==2
853  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
854  flag_3d_output, ndat2d, ndat3d)
855 #else
856  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
857 #endif
858 
859 end if
860 
861 #elif OUTPUT==3
862 
863  flag_3d_output = .false.
864 
865 #if NETCDF==1
866  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
867  flag_3d_output, ndat2d, ndat3d)
868 #elif NETCDF==2
869  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
870  flag_3d_output, ndat2d, ndat3d)
871 #else
872  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
873 #endif
874 
875 if (time_output(1) <= time_init+eps) then
876 
877  flag_3d_output = .true.
878 
879 #if NETCDF==1
880  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
881  flag_3d_output, ndat2d, ndat3d)
882 #elif NETCDF==2
883  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
884  flag_3d_output, ndat2d, ndat3d)
885 #else
886  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
887 #endif
888 
889 end if
890 
891 #else
892  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
893 #endif
894 
895 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
896 #if OUTSER==2
897 call output3(time_init, delta_ts, glac_index, z_sl)
898 #endif
899 #if OUTSER==3
900 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
901 #endif
902 
903 end subroutine sico_init
904 !
Generated on Thu Sep 12 2013 09:44:59 for SICOPOLIS V3.1 by   doxygen 1.8.2