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 character(len= 8) :: ch_imax
80 character(len=128) :: fmt4
81 
82 write(ch_imax, fmt='(i8)') imax
83 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
84 
85 !-------- Some initial values --------
86 
87 n_output = 0
88 
89 dtime = 0.0_dp
90 dtime_temp = 0.0_dp
91 dtime_wss = 0.0_dp
92 dtime_out = 0.0_dp
93 dtime_ser = 0.0_dp
94 
95 time = 0.0_dp
96 time_init = 0.0_dp
97 time_end = 0.0_dp
98 time_output = 0.0_dp
99 
100 !-------- Initialisation of the Library of Iterative Solvers Lis,
101 ! if required --------
102 
103 #if (CALCZS==4 || MARGIN==3)
104 
105 #include "lisf.h"
106 call lis_initialize(ierr)
107 
108 #endif
109 
110 !-------- Read physical parameters --------
111 
112 call phys_para()
113 
114 temp_min = temp_min
115 s_t = s_t *1.0e-09_dp ! K/km3 -> K/m3
116 x_hat = x_hat *1.0e+03_dp ! km -> m
117 y_hat = y_hat *1.0e+03_dp ! km -> m
118 rad = rad *1.0e+03_dp ! km -> m
119 b_min = b_min /year_sec ! m/a -> m/s
120 b_max = b_max /year_sec ! m/a -> m/s
121 
122 !-------- Compatibility check of the SICOPOLIS version with the header file
123 
124 if ( trim(version) /= trim(sico_version) ) &
125  stop ' sico_init: SICOPOLIS version not compatible with header file!'
126 
127 !-------- Compatibility check of the horizontal resolution with the
128 ! number of grid points --------
129 
130 #if GRID==0
131 
132 if (abs(dx-50.0_dp) < eps) then
133 
134  if ((imax /= 80).or.(jmax /= 80)) then
135  stop ' sico_init: IMAX and/or JMAX wrong!'
136  end if
137 
138 else if (abs(dx-25.0_dp) < eps) then
139 
140  if ((imax /= 160).or.(jmax /= 160)) then
141  stop ' sico_init: IMAX and/or JMAX wrong!'
142  end if
143 
144 else
145  stop ' sico_init: DX wrong!'
146 end if
147 
148 #elif GRID==1
149 
150 stop ' sico_init: GRID==1 not allowed for this application!'
151 
152 #elif GRID==2
153 
154 stop ' sico_init: GRID==2 not allowed for this application!'
155 
156 #endif
157 
158 !-------- Compatibility check of discretization schemes for the horizontal and
159 ! vertical advection terms in the temperature and age equations --------
160 
161 #if ADV_HOR==1
162 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
163 #endif
164 
165 !-------- Setting of forcing flag --------
166 
167 forcing_flag = 1 ! forcing by delta_ts
168 
169 !-------- Initialization of numerical time steps --------
170 
171 dtime0 = dtime0
172 dtime_temp0 = dtime_temp0
173 
174 !-------- Further initializations --------
175 
176 dzeta_c = 1.0_dp/real(kcmax,dp)
177 dzeta_t = 1.0_dp/real(ktmax,dp)
178 dzeta_r = 1.0_dp/real(krmax,dp)
179 
180 ndat2d = 1
181 ndat3d = 1
182 
183 !-------- Construction of a vector (with index n) from a 2-d array
184 ! (with indices i, j) by diagonal numbering --------
185 
186 n=1
187 
188 do m=0, imax+jmax
189  do i=m, 0, -1
190  j = m-i
191  if ((i <= imax).and.(j <= jmax)) then
192  ii(n) = i
193  jj(n) = j
194  nn(j,i) = n
195  n=n+1
196  end if
197  end do
198 end do
199 
200 !-------- Specification of current simulation --------
201 
202 runname = runname
203 anfdatname = anfdatname
204 
205 #if defined(YEAR_ZERO)
206 year_zero = year_zero
207 #else
208 year_zero = 2000.0_dp ! default value 2000 CE
209 #endif
210 
211 time_init0 = time_init0
212 time_end0 = time_end0
213 dtime_ser0 = dtime_ser0
214 
215 #if ( OUTPUT==1 || OUTPUT==3 )
216 dtime_out0 = dtime_out0
217 #endif
218 
219 #if ( OUTPUT==2 || OUTPUT==3 )
220 n_output = n_output
221 time_output0( 1) = time_out0_01
222 time_output0( 2) = time_out0_02
223 time_output0( 3) = time_out0_03
224 time_output0( 4) = time_out0_04
225 time_output0( 5) = time_out0_05
226 time_output0( 6) = time_out0_06
227 time_output0( 7) = time_out0_07
228 time_output0( 8) = time_out0_08
229 time_output0( 9) = time_out0_09
230 time_output0(10) = time_out0_10
231 time_output0(11) = time_out0_11
232 time_output0(12) = time_out0_12
233 time_output0(13) = time_out0_13
234 time_output0(14) = time_out0_14
235 time_output0(15) = time_out0_15
236 time_output0(16) = time_out0_16
237 time_output0(17) = time_out0_17
238 time_output0(18) = time_out0_18
239 time_output0(19) = time_out0_19
240 time_output0(20) = time_out0_20
241 #endif
242 
243 !-------- Write log file --------
244 
245 shell_command = 'if [ ! -d'
246 shell_command = trim(shell_command)//' '//outpath
247 shell_command = trim(shell_command)//' '//'] ; then mkdir'
248 shell_command = trim(shell_command)//' '//outpath
249 shell_command = trim(shell_command)//' '//'; fi'
250 call system(trim(shell_command))
251  ! Check whether directory OUTPATH exists. If not, it is created.
252 
253 open(10, iostat=ios, &
254  file=outpath//'/'//trim(runname)//'.log', &
255  status='new')
256 if (ios /= 0) stop ' sico_init: Error when opening the log file!'
257 
258 write(10, fmt=trim(fmt1)) 'Computational domain:'
259 #if defined(ANT)
260 write(10, fmt=trim(fmt1)) 'Antarctica'
261 #elif defined(GRL)
262 write(10, fmt=trim(fmt1)) 'Greenland'
263 #elif defined(SCAND)
264 write(10, fmt=trim(fmt1)) 'Scandinavia and Eurasia'
265 #elif defined(NHEM)
266 write(10, fmt=trim(fmt1)) 'Northern hemisphere'
267 #elif defined(EMTP2SGE)
268 write(10, fmt=trim(fmt1)) 'EISMINT Phase 2 Simplified Geometry Experiment'
269 #elif defined(HEINO)
270 write(10, fmt=trim(fmt1)) 'ISMIP HEINO'
271 #elif defined(NMARS)
272 write(10, fmt=trim(fmt1)) 'North polar cap of Mars'
273 #elif defined(SMARS)
274 write(10, fmt=trim(fmt1)) 'South polar cap of Mars'
275 #else
276 stop ' sico_init: Subroutine sico_init: No valid domain specified!'
277 #endif
278 write(10, fmt=trim(fmt1)) ' '
279 
280 write(10, fmt=trim(fmt2)) 'imax =', imax
281 write(10, fmt=trim(fmt2)) 'jmax =', jmax
282 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
283 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
284 write(10, fmt=trim(fmt2)) 'krmax =', krmax
285 write(10, fmt=trim(fmt1)) ' '
286 
287 write(10, fmt=trim(fmt3)) 'a =', deform
288 write(10, fmt=trim(fmt1)) ' '
289 
290 #if (GRID==0 || GRID==1)
291 write(10, fmt=trim(fmt3)) 'x0 =', x0
292 write(10, fmt=trim(fmt3)) 'y0 =', y0
293 write(10, fmt=trim(fmt3)) 'dx =', dx
294 #elif GRID==2
295 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
296 write(10, fmt=trim(fmt3)) 'dphi =', dphi
297 #endif
298 write(10, fmt=trim(fmt1)) ' '
299 
300 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
301 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
302 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
303 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
304 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
305 #if ( OUTPUT==1 || OUTPUT==3 )
306 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
307 #endif
308 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
309 write(10, fmt=trim(fmt1)) ' '
310 
311 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
312 write(10, fmt=trim(fmt1)) ' '
313 
314 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
315 write(10, fmt=trim(fmt1)) ' '
316 
317 #if (CALCZS==3 || CALCZS==4)
318 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
319 #if CALCZS==3
320 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
321 #endif
322 write(10, fmt=trim(fmt1)) ' '
323 #endif
324 
325 write(10, fmt=trim(fmt3)) 'temp_min =', temp_min
326 write(10, fmt=trim(fmt3)) 's_t =', s_t
327 write(10, fmt=trim(fmt3)) 'x_hat =', x_hat
328 write(10, fmt=trim(fmt3)) 'y_hat =', y_hat
329 write(10, fmt=trim(fmt3)) 'rad =', rad
330 write(10, fmt=trim(fmt3)) 'b_min =', b_min
331 write(10, fmt=trim(fmt3)) 'b_max =', b_max
332 write(10, fmt=trim(fmt1)) ' '
333 
334 #if TSURFACE==1
335 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
336 #elif TSURFACE==3
337 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
338 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
339 #elif TSURFACE==4
340 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
341 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
342 #endif
343 
344 #if SEA_LEVEL==1
345 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
346 #elif SEA_LEVEL==3
347 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
348 #endif
349 write(10, fmt=trim(fmt1)) ' '
350 
351 #if MARGIN==2
352 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
353 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
354 write(10, fmt=trim(fmt1)) ' '
355 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
356  || marine_ice_calving==6 || marine_ice_calving==7 )
357 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
358 write(10, fmt=trim(fmt1)) ' '
359 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
360 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
361 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
362 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
363 write(10, fmt=trim(fmt1)) ' '
364 #endif
365 #elif MARGIN==3
366 #if ICE_SHELF_CALVING==2
367 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
368 write(10, fmt=trim(fmt1)) ' '
369 #endif
370 #endif
371 
372 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
373 write(10, fmt=trim(fmt1)) ' '
374 
375 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
376 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
377 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
378 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
379 #if SLIDE_LAW==2
380 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
381 #endif
382 write(10, fmt=trim(fmt1)) ' '
383 
384 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
385 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
386 write(10, fmt=trim(fmt2)) 'p_weert_sedi =', p_weert_sedi
387 write(10, fmt=trim(fmt2)) 'q_weert_sedi =', q_weert_sedi
388 write(10, fmt=trim(fmt1)) ' '
389 
390 #if Q_GEO_MOD==1
391 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
392 write(10, fmt=trim(fmt1)) ' '
393 #endif
394 
395 #if REBOUND==1
396 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
397 write(10, fmt=trim(fmt1)) ' '
398 #endif
399 
400 #if FLOW_LAW==2
401 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
402 write(10, fmt=trim(fmt1)) ' '
403 #endif
404 #if FIN_VISC==2
405 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
406 write(10, fmt=trim(fmt1)) ' '
407 #endif
408 
409 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
410 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
411 #if ( ENHMOD==2 || ENHMOD==3 )
412 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
413 #endif
414 #if ENHMOD==2
415 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
416 #endif
417 #if ENHMOD==3
418 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
419 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
420 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
421 #endif
422 #if MARGIN==3
423 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
424 #endif
425 write(10, fmt=trim(fmt1)) ' '
426 
427 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
428 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
429 write(10, fmt=trim(fmt3)) 'H_min =', h_min
430 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
431 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
432 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
433 #if defined(QBM_MIN)
434 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
435 #elif defined(QB_MIN) /* obsolete */
436 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
437 #endif
438 #if defined(QBM_MAX)
439 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
440 #elif defined(QB_MAX) /* obsolete */
441 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
442 #endif
443 write(10, fmt=trim(fmt3)) 'age_min =', age_min
444 write(10, fmt=trim(fmt3)) 'age_max =', age_max
445 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
446 #if ADV_VERT==1
447 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
448 #endif
449 write(10, fmt=trim(fmt1)) ' '
450 
451 write(10, fmt=trim(fmt2)) 'GRID =', grid
452 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
453 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
454 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
455 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
456 #if MARGIN==2
457 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
458 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
459 #elif MARGIN==3
460 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
461 #endif
462 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
463 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
464 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
465 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
466 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
467 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
468 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
469 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
470 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
471 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
472 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
473 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
474 write(10, fmt=trim(fmt1)) ' '
475 
476 #if defined(OUT_TIMES)
477 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
478 #endif
479 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
480 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
481 #if defined(OUTSER_V_A)
482 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
483 #endif
484 #if ( OUTPUT==1 || OUTPUT==2 )
485 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
486 #endif
487 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
488 write(10, fmt=trim(fmt1)) ' '
489 #if ( OUTPUT==2 || OUTPUT==3 )
490 write(10, fmt=trim(fmt2)) 'n_output =', n_output
491 do n=1, n_output
492  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
493 end do
494 write(10, fmt=trim(fmt1)) ' '
495 #endif
496 
497 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
498 
499 close(10, status='keep')
500 
501 !-------- Conversion of time quantities --------
502 
503 year_zero = year_zero*year_sec ! a --> s
504 time_init = time_init0*year_sec ! a --> s
505 time_end = time_end0*year_sec ! a --> s
506 dtime = dtime0*year_sec ! a --> s
507 dtime_temp = dtime_temp0*year_sec ! a --> s
508 dtime_ser = dtime_ser0*year_sec ! a --> s
509 #if ( OUTPUT==1 || OUTPUT==3 )
510 dtime_out = dtime_out0*year_sec ! a --> s
511 #endif
512 #if ( OUTPUT==2 || OUTPUT==3 )
513 do n=1, n_output
514  time_output(n) = time_output0(n)*year_sec ! a --> s
515 end do
516 #endif
517 
518 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
519  stop ' sico_init: dtime_temp must be a multiple of dtime!'
520 
521 time = time_init
522 
523 !-------- Mean accumulation --------
524 
525 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
526 ! ! mm/a water equiv. --> m/s ice equiv.
527 
528 mean_accum_inv = 1.0_dp/mean_accum
529 
530 !-------- Read present topography mask --------
531 
532 open(24, iostat=ios, &
533  file=inpath//'/heino/'//mask_present_file, &
534  recl=8192, status='old')
535 
536 if (ios /= 0) stop ' sico_init: Error when opening the mask file!'
537 
538 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
539 
540 do j=jmax, 0, -1
541  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
542 end do
543 
544 close(24, status='keep')
545 
546 !-------- Read rock/sediment mask --------
547 
548 open(24, iostat=ios, &
549  file=inpath//'/heino/'//mask_sedi_file, &
550  recl=8192, status='old')
551 
552 if (ios /= 0) stop ' sico_init: Error when opening the rock/sediment mask file!'
553 
554 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
555 
556 do j=jmax, 0, -1
557  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
558 end do
559 
560 close(24, status='keep')
561 
562 !-------- Read data for delta_ts --------
563 
564 #if TSURFACE==4
565 
566 open(21, iostat=ios, &
567  file=inpath//'/general/'//grip_temp_file, &
568  status='old')
569 if (ios /= 0) stop ' sico_init: Error when opening the data file for delta_ts!'
570 
571 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
572 
573 if (ch_dummy /= '#') then
574  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
575  write(6, fmt=*) ' not defined in data file for delta_ts!'
576  stop
577 end if
578 
579 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
580 
581 allocate(griptemp(0:ndata_grip))
582 
583 do n=0, ndata_grip
584  read(21, fmt=*) d_dummy, griptemp(n)
585 end do
586 
587 close(21, status='keep')
588 
589 #endif
590 
591 !-------- Read data for z_sl --------
592 
593 #if SEA_LEVEL==3
594 
595 open(21, iostat=ios, &
596  file=inpath//'/general/'//sea_level_file, &
597  status='old')
598 if (ios /= 0) stop ' sico_init: Error when opening the data file for z_sl!'
599 
600 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
601 
602 if (ch_dummy /= '#') then
603  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
604  write(6, fmt=*) ' not defined in data file for z_sl!'
605  stop
606 end if
607 
608 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
609 
610 allocate(specmap_zsl(0:ndata_specmap))
611 
612 do n=0, ndata_specmap
613  read(21, fmt=*) d_dummy, specmap_zsl(n)
614 end do
615 
616 close(21, status='keep')
617 
618 #endif
619 
620 !-------- Determination of the geothermal heat flux --------
621 
622 #if Q_GEO_MOD==1
623 
624 ! ------ Constant value
625 
626 do i=0, imax
627 do j=0, jmax
628  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
629 end do
630 end do
631 
632 #elif Q_GEO_MOD==2
633 
634 stop ' sico_init: Option Q_GEO_MOD==2 not available for this application!'
635 
636 #endif
637 
638 !-------- Definition of initial values --------
639 
640 ! ------ Initial topography with a thin ice layer everywhere
641 
642 #if ANF_DAT==1
643 
644 call topography1(dxi, deta)
645 
646 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
647 
648 call boundary(time_init, dtime, dxi, deta, &
649  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
650 
651 do i=0, imax
652 do j=0, jmax
653 
654  q_bm(j,i) = 0.0_dp
655  q_tld(j,i) = 0.0_dp
656  q_b_tot(j,i) = 0.0_dp
657 
658  p_b_w(j,i) = 0.0_dp
659  h_w(j,i) = 0.0_dp
660 
661  do kc=0, kcmax
662  temp_c(kc,j,i) = temp_s(j,i)
663  age_c(kc,j,i) = 15000.0_dp*year_sec
664  end do
665 
666  do kt=0, ktmax
667  omega_t(kt,j,i) = 0.0_dp
668  age_t(kt,j,i) = 15000.0_dp*year_sec
669  end do
670 
671  do kr=0, krmax
672  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
673  *(1.0_dp-real(kr,dp)/real(krmax,dp))
674 ! (linear temperature distribution according to the
675 ! geothermal heat flux)
676  end do
677 
678 end do
679 end do
680 
681 call calc_enhance(time_init)
682 
683 ! ------ Ice-free, relaxed bedrock
684 
685 #elif ANF_DAT==2
686 
687 call topography2(dxi, deta)
688 
689 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
690 
691 call boundary(time_init, dtime, dxi, deta, &
692  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
693 
694 do i=0, imax
695 do j=0, jmax
696 
697  q_bm(j,i) = 0.0_dp
698  q_tld(j,i) = 0.0_dp
699  q_b_tot(j,i) = 0.0_dp
700 
701  p_b_w(j,i) = 0.0_dp
702  h_w(j,i) = 0.0_dp
703 
704  do kc=0, kcmax
705  temp_c(kc,j,i) = temp_s(j,i)
706  age_c(kc,j,i) = 15000.0_dp*year_sec
707  end do
708 
709  do kt=0, ktmax
710  omega_t(kt,j,i) = 0.0_dp
711  age_t(kt,j,i) = 15000.0_dp*year_sec
712  end do
713 
714  do kr=0, krmax
715  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
716  *(1.0_dp-real(kr,dp)/real(krmax,dp))
717 ! (linear temperature distribution according to the
718 ! geothermal heat flux)
719  end do
720 
721 end do
722 end do
723 
724 call calc_enhance(time_init)
725 
726 ! ------ Read initial state from output of previous simulation
727 
728 #elif ANF_DAT==3
729 
730 #if NETCDF==1
731 call topography3(dxi, deta, z_sl, anfdatname)
732 #elif NETCDF==2
733 call topography3_nc(dxi, deta, z_sl, anfdatname)
734 #else
735 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
736 #endif
737 
738 call boundary(time_init, dtime, dxi, deta, &
739  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
740 
741 q_b_tot = q_bm + q_tld
742 
743 call calc_enhance(time_init)
744 
745 #endif
746 
747 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
748 
749 #if (GRID==0 || GRID==1)
750 
751 do ir=-imax, imax
752 do jr=-jmax, jmax
753  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
754  ! distortion due to stereographic projection not accounted for
755 end do
756 end do
757 
758 #endif
759 
760 !-------- General abbreviations --------
761 
762 ea = exp(deform)
763 
764 do kc=0, kcmax
765  zeta_c(kc) = kc*dzeta_c
766  eaz_c(kc) = exp(deform*zeta_c(kc))
767  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
768 end do
769 
770 do kt=0, ktmax
771  zeta_t(kt) = kt*dzeta_t
772 end do
773 
774 !-------- Initial velocities --------
775 
776 call calc_temp_melt()
777 
778 call calc_vxy_b_sia(time, z_sl)
779 call calc_vxy_sia(dzeta_c, dzeta_t)
780 
781 #if MARGIN==3
782 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
783 #endif
784 
785 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
786 
787 #if MARGIN==3
788 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
789 #endif
790 
791 !-------- Initialize time-series files --------
792 
793 ! ------ Time-series file for the ice sheet on the whole
794 
795 open(12, iostat=ios, &
796  file=outpath//'/'//trim(runname)//'.ser', &
797  status='new')
798 if (ios /= 0) stop ' sico_init: Error when opening the ser file!'
799 
800 write(12,1102)
801 write(12,1103)
802 
803 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
804 
805  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
806  ' H_max(m) zs_max(m) V_g(m^3)', &
807  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
808  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
809  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
810  1103 format('----------------------------------------------------', &
811  '---------------------------------------')
812 
813 #elif OUTSER_V_A==2
814 
815  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
816  ' V(m^3) V_g(m^3) V_f(m^3)', &
817  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
818  ' H_max(m) zs_max(m)', &
819  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
820  ' A_t(m^2) V_bm(m^3/a)', &
821  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
822  1103 format('----------------------------------------------------', &
823  '---------------------------------------')
824 
825 #else
826  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
827 #endif
828 
829 ! ------ Time-series file for the sediment area ("Hudson Bay, Hudson Strait")
830 
831 open(15, iostat=ios, &
832  file=outpath//'/'//trim(runname)//'.sed', &
833  status='new')
834 if (ios /= 0) stop ' sico_init: Error when opening the sed file!'
835 
836 write(15,1108)
837 write(15,1109)
838 
839 1108 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
840  ' H_ave(m) Tbh_ave(C) Atb(m^2)')
841 1109 format('----------------------------------------------------')
842 
843 ! ------ Time-series file for the ice-thickness field
844 
845 #if OUTSER==2
846 
847 open(13, iostat=ios, &
848  file=outpath//'/'//trim(runname)//'.thk', &
849  status='new')
850 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
851 
852 write(13,1104)
853 write(13,1105)
854 
855 1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
856  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
857  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
858 1105 format('----------------------------------------------------')
859 
860 #endif
861 
862 ! ------ Time-series file for selected positions ("deep boreholes")
863 
864 #if OUTSER==3
865 
866 n_core = 7 ! Points P1 - P7
867 
868 allocate(lambda_core(n_core), phi_core(n_core), &
869  x_core(n_core), y_core(n_core))
870 
871 lambda_core(1) = 0.0_dp ! dummy
872 phi_core(1) = 0.0_dp ! dummy
873 x_core(1) = 3900.0_dp *1.0e+03_dp ! Point P1,
874 y_core(1) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
875 
876 lambda_core(2) = 0.0_dp ! dummy
877 phi_core(2) = 0.0_dp ! dummy
878 x_core(2) = 3800.0_dp *1.0e+03_dp ! Point P2,
879 y_core(2) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
880 
881 lambda_core(3) = 0.0_dp ! dummy
882 phi_core(3) = 0.0_dp ! dummy
883 x_core(3) = 3700.0_dp *1.0e+03_dp ! Point P3,
884 y_core(3) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
885 
886 lambda_core(4) = 0.0_dp ! dummy
887 phi_core(4) = 0.0_dp ! dummy
888 x_core(4) = 3500.0_dp *1.0e+03_dp ! Point P4,
889 y_core(4) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
890 
891 lambda_core(5) = 0.0_dp ! dummy
892 phi_core(5) = 0.0_dp ! dummy
893 x_core(5) = 3200.0_dp *1.0e+03_dp ! Point P5,
894 y_core(5) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
895 
896 lambda_core(6) = 0.0_dp ! dummy
897 phi_core(6) = 0.0_dp ! dummy
898 x_core(6) = 2900.0_dp *1.0e+03_dp ! Point P6,
899 y_core(6) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
900 
901 lambda_core(7) = 0.0_dp ! dummy
902 phi_core(7) = 0.0_dp ! dummy
903 x_core(7) = 2600.0_dp *1.0e+03_dp ! Point P7,
904 y_core(7) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
905 
906 open(14, iostat=ios, &
907  file=outpath//'/'//trim(runname)//'.core', &
908  status='new')
909 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
910 
911 if (forcing_flag == 1) then
912 
913  write(14,1106)
914  write(14,1107)
915 
916  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
917  ' H_P1(m) H_P2(m) H_P3(m) H_P4(m)', &
918  ' H_P5(m) H_P6(m) H_P7(m)',/, &
919  ' v_P1(m/a) v_P2(m/a) v_P3(m/a) v_P4(m/a)', &
920  ' v_P5(m/a) v_P6(m/a) v_P7(m/a)',/, &
921  ' T_P1(C) T_P2(C) T_P3(C) T_P4(C)', &
922  ' T_P5(C) T_P6(C) T_P7(C)',/, &
923  ' Rb_P1(W/m2) Rb_P2(W/m2) Rb_P3(W/m2) Rb_P4(W/m2)', &
924  ' Rb_P5(W/m2) Rb_P6(W/m2) Rb_P7(W/m2)')
925  1107 format('-----------------------------------------------------------------', &
926  '---------------------------------------')
927 
928 end if
929 
930 #endif
931 
932 !-------- Output of the initial state --------
933 
934 #if OUTPUT==1
935 
936 #if ERGDAT==0
937  flag_3d_output = .false.
938 #elif ERGDAT==1
939  flag_3d_output = .true.
940 #else
941  stop ' sico_init: ERGDAT must be either 0 or 1!'
942 #endif
943 
944 #if NETCDF==1
945  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
946  flag_3d_output, ndat2d, ndat3d)
947 #elif NETCDF==2
948  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
949  flag_3d_output, ndat2d, ndat3d)
950 #else
951  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
952 #endif
953 
954 #elif OUTPUT==2
955 
956 if (time_output(1) <= time_init+eps) then
957 
958 #if ERGDAT==0
959  flag_3d_output = .false.
960 #elif ERGDAT==1
961  flag_3d_output = .true.
962 #else
963  stop ' sico_init: ERGDAT must be either 0 or 1!'
964 #endif
965 
966 #if NETCDF==1
967  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
968  flag_3d_output, ndat2d, ndat3d)
969 #elif NETCDF==2
970  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
971  flag_3d_output, ndat2d, ndat3d)
972 #else
973  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
974 #endif
975 
976 end if
977 
978 #elif OUTPUT==3
979 
980  flag_3d_output = .false.
981 
982 #if NETCDF==1
983  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
984  flag_3d_output, ndat2d, ndat3d)
985 #elif NETCDF==2
986  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
987  flag_3d_output, ndat2d, ndat3d)
988 #else
989  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
990 #endif
991 
992 if (time_output(1) <= time_init+eps) then
993 
994  flag_3d_output = .true.
995 
996 #if NETCDF==1
997  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
998  flag_3d_output, ndat2d, ndat3d)
999 #elif NETCDF==2
1000  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1001  flag_3d_output, ndat2d, ndat3d)
1002 #else
1003  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1004 #endif
1005 
1006 end if
1007 
1008 #else
1009  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1010 #endif
1011 
1012 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1013 #if OUTSER==2
1014 call output3(time_init, delta_ts, glac_index, z_sl)
1015 #endif
1016 #if OUTSER==3
1017 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1018 #endif
1019 
1020 end subroutine sico_init
1021 !
Generated on Thu Sep 12 2013 09:44:59 for SICOPOLIS V3.1 by   doxygen 1.8.2