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