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