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 /= 90).or.(jmax /= 90)) 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 /= 180).or.(jmax /= 180)) 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 /= 360).or.(jmax /= 360)) 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 nmars 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_90N =', temp0_ma_90n
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_90n_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//'/nmars/'//mask_present_file, &
623  recl=1024, status='old')
624 
625 if (ios /= 0) &
626  stop ' sico_init: Error when opening the mask file!'
627 
628 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
629 
630 do j=jmax, 0, -1
631  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
632 end do
633 
634 close(24, status='keep')
635 
636 !-------- Read chasm mask --------
637 
638 #if CHASM==2
639 
640 open(24, iostat=ios, &
641  file=inpath//'/nmars/'//mask_chasm_file, &
642  recl=1024, status='old')
643 
644 if (ios /= 0) stop ' sico_init: Error when opening the chasm mask file!'
645 
646 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
647 
648 do j=jmax, 0, -1
649  read(24, fmt=trim(fmt4)) (maske_chasm(j,i), i=0,imax)
650 end do
651 
652 close(24, status='keep')
653 
654 #endif
655 
656 !-------- Read data for z_sl --------
657 
658 #if SEA_LEVEL==3
659 
660 stop ' sico_init: SEA_LEVEL==3 not allowed for nmars application!'
661 
662 #endif
663 
664 !-------- Reading of insolation data --------
665 
666 insol_ma_90 = 0.0_dp ! Assignment of dummy values
667 obl_data = 0.0_dp
668 ecc_data = 0.0_dp
669 ave_data = 0.0_dp
670 cp_data = 0.0_dp
671 
672 #if ( TSURFACE==5 || TSURFACE==6 )
673 
674 open(21, iostat=ios, &
675  file=inpath//'/nmars/'//insol_ma_90n_file, &
676  status='old')
677 if (ios /= 0) stop ' sico_init: Error when opening the insolation-data file!'
678 
679 read(21, fmt=*) ch_dummy, insol_time_min, insol_time_stp, insol_time_max
680 
681 if (ch_dummy /= '#') then
682  write(6, fmt=*) ' sico_init: insol_time_min, insol_time_stp, insol_time_max'
683  write(6, fmt=*) ' not defined in insolation-data file!'
684  stop
685 end if
686 
687 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
688 
689 if (ndata_insol > 100000) &
690  stop ' sico_init: Too many data in insolation-data file!'
691 
692 do n=0, ndata_insol
693  read(21, fmt=*) d_dummy, ecc_data(n), obl_data(n), cp_data(n), ave_data(n), insol_ma_90(n)
694  obl_data(n) = obl_data(n) *pi_180
695  ave_data(n) = ave_data(n) *pi_180
696 end do
697 
698 close(21, status='keep')
699 
700 #endif
701 
702 !-------- Determination of the normal geothermal heat flux
703 ! (without the possible contribution from an active chasm area) --------
704 
705 #if Q_GEO_MOD==1
706 
707 ! ------ Constant value
708 
709 do i=0, imax
710 do j=0, jmax
711  q_geo_normal(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
712 end do
713 end do
714 
715 #elif Q_GEO_MOD==2
716 
717 ! ------ Read data from file
718 
719 open(21, iostat=ios, &
720  file=inpath//'/nmars/'//q_geo_file, &
721  recl=8192, status='old')
722 
723 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
724 
725 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
726 
727 do j=jmax, 0, -1
728  read(21, fmt=*) (q_geo_normal(j,i), i=0,imax)
729 end do
730 
731 close(21, status='keep')
732 
733 do i=0, imax
734 do j=0, jmax
735  q_geo_normal(j,i) = q_geo_normal(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
736 end do
737 end do
738 
739 #endif
740 
741 !-------- Reading of tabulated kei function--------
742 
743 #if (REBOUND==0 || REBOUND==1)
744 
745 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
746  ! dummy values
747 #elif REBOUND==2
748 
749 call read_kei()
750 
751 #endif
752 
753 !-------- Definition of initial values --------
754 
755 ! ------ Present topography
756 
757 #if ANF_DAT==1
758 
759 call topography1(dxi, deta)
760 
761 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
762 
763 call boundary(time_init, dtime, dxi, deta, &
764  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
765 
766 do i=0, imax
767 do j=0, jmax
768 
769  q_bm(j,i) = 0.0_dp
770  q_tld(j,i) = 0.0_dp
771  q_b_tot(j,i) = 0.0_dp
772 
773  p_b_w(j,i) = 0.0_dp
774  h_w(j,i) = 0.0_dp
775 
776  do kc=0, kcmax
777  temp_c(kc,j,i) = temp_s(j,i)
778  age_c(kc,j,i) = 15000.0_dp*year_sec
779  end do
780 
781  do kt=0, ktmax
782  omega_t(kt,j,i) = 0.0_dp
783  age_t(kt,j,i) = 15000.0_dp*year_sec
784  end do
785 
786  do kr=0, krmax
787  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
788  *(1.0_dp-real(kr,dp)/real(krmax,dp))
789 ! (linear temperature distribution according to the
790 ! geothermal heat flux)
791  end do
792 
793 end do
794 end do
795 
796 call calc_enhance(time_init)
797 
798 ! ------ Ice-free, relaxed bedrock
799 
800 #elif ANF_DAT==2
801 
802 call topography2(dxi, deta)
803 
804 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
805 
806 call boundary(time_init, dtime, dxi, deta, &
807  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
808 
809 do i=0, imax
810 do j=0, jmax
811 
812  q_bm(j,i) = 0.0_dp
813  q_tld(j,i) = 0.0_dp
814  q_b_tot(j,i) = 0.0_dp
815 
816  p_b_w(j,i) = 0.0_dp
817  h_w(j,i) = 0.0_dp
818 
819  do kc=0, kcmax
820  temp_c(kc,j,i) = temp_s(j,i)
821  age_c(kc,j,i) = 15000.0_dp*year_sec
822  end do
823 
824  do kt=0, ktmax
825  omega_t(kt,j,i) = 0.0_dp
826  age_t(kt,j,i) = 15000.0_dp*year_sec
827  end do
828 
829  do kr=0, krmax
830  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
831  *(1.0_dp-real(kr,dp)/real(krmax,dp))
832 ! (linear temperature distribution according to the
833 ! geothermal heat flux)
834  end do
835 
836 end do
837 end do
838 
839 call calc_enhance(time_init)
840 
841 ! ------ Read initial state from output of previous simulation
842 
843 #elif ANF_DAT==3
844 
845 #if NETCDF==1
846 call topography3(dxi, deta, z_sl, anfdatname)
847 #elif NETCDF==2
848 call topography3_nc(dxi, deta, z_sl, anfdatname)
849 #else
850 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
851 #endif
852 
853 call boundary(time_init, dtime, dxi, deta, &
854  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
855 
856 q_b_tot = q_bm + q_tld
857 
858 call calc_enhance(time_init)
859 
860 #endif
861 
862 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
863 
864 #if (GRID==0 || GRID==1)
865 
866 do ir=-imax, imax
867 do jr=-jmax, jmax
868  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
869  ! distortion due to stereographic projection not accounted for
870 end do
871 end do
872 
873 #endif
874 
875 !-------- General abbreviations --------
876 
877 ea = exp(deform)
878 
879 do kc=0, kcmax
880  zeta_c(kc) = kc*dzeta_c
881  eaz_c(kc) = exp(deform*zeta_c(kc))
882  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
883 end do
884 
885 do kt=0, ktmax
886  zeta_t(kt) = kt*dzeta_t
887 end do
888 
889 !-------- Initial velocities --------
890 
891 call calc_temp_melt()
892 
893 call calc_vxy_b_sia(time, z_sl)
894 call calc_vxy_sia(dzeta_c, dzeta_t)
895 
896 #if MARGIN==3
897 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
898 #endif
899 
900 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
901 
902 #if MARGIN==3
903 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
904 #endif
905 
906 !-------- Initialize time-series files --------
907 
908 ! ------ Time-series file for the ice sheet on the whole
909 
910 open(12, iostat=ios, &
911  file=outpath//'/'//trim(runname)//'.ser', &
912  status='new')
913 if (ios /= 0) &
914  stop ' sico_init: Error when opening the ser file!'
915 
916 write(12,1102)
917 write(12,1103)
918 
919 1102 format(' t(a) D_Ts(C) z_sl(m)',/, &
920  ' H_max(m) zs_max(m) V(m^3)', &
921  ' V_t(m^3) V_fw(m^3/a) Tbh_max(C)',/, &
922  ' A(m^2) A_t(m^2) V_bm(m^3/a)', &
923  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
924 1103 format('----------------------------------------------------', &
925  '---------------------------------------')
926 
927 ! ------ Time-series file for the ice-thickness field
928 
929 #if OUTSER==2
930 
931 open(13, iostat=ios, &
932  file=outpath//'/'//trim(runname)//'.thk', &
933  status='new')
934 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
935 
936 write(13,1104)
937 write(13,1105)
938 
939 1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
940  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
941  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
942 1105 format('----------------------------------------------------')
943 
944 #endif
945 
946 ! ------ Time-series file for selected positions ("deep boreholes")
947 
948 #if OUTSER==3
949 
950 n_core = 3 ! Points NP (north pole), C1, C2 (in Chasma Borealis)
951 
952 allocate(lambda_core(n_core), phi_core(n_core), &
953  x_core(n_core), y_core(n_core))
954 
955 lambda_core(1) = 0.0_dp ! dummy
956 phi_core(1) = 0.0_dp ! dummy
957 x_core(1) = 0.0_dp *1.0e+03_dp ! North pole,
958 y_core(1) = 0.0_dp *1.0e+03_dp ! conversion km -> m
959 
960 lambda_core(2) = 0.0_dp ! dummy
961 phi_core(2) = 0.0_dp ! dummy
962 x_core(2) = -150.0_dp *1.0e+03_dp ! Point C1,
963 y_core(2) = -290.0_dp *1.0e+03_dp ! conversion km -> m
964 
965 lambda_core(3) = 0.0_dp ! dummy
966 phi_core(3) = 0.0_dp ! dummy
967 x_core(3) = -300.0_dp *1.0e+03_dp ! Point C2,
968 y_core(3) = -280.0_dp *1.0e+03_dp ! conversion km -> m
969 
970 open(14, iostat=ios, &
971  file=outpath//'/'//trim(runname)//'.core', &
972  status='new')
973 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
974 
975 if (forcing_flag == 1) then
976 
977  write(14,1106)
978  write(14,1107)
979 
980  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
981  ' H_NP(m) H_C1(m) H_C2(m)',/, &
982  ' v_NP(m/a) v_C1(m/a) v_C2(m/a)',/, &
983  ' T_NP(C) T_C1(C) T_C2(C)',/, &
984  ' Rb_NP(W/m2) Rb_C1(W/m2) Rb_C2(W/m2)')
985  1107 format('----------------------------------------------------')
986 
987 end if
988 
989 #endif
990 
991 !-------- Output of the initial state --------
992 
993 #if OUTPUT==1
994 
995 #if ERGDAT==0
996  flag_3d_output = .false.
997 #elif ERGDAT==1
998  flag_3d_output = .true.
999 #else
1000  stop ' sico_init: ERGDAT must be either 0 or 1!'
1001 #endif
1002 
1003 #if NETCDF==1
1004  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1005  flag_3d_output, ndat2d, ndat3d)
1006 #elif NETCDF==2
1007  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1008  flag_3d_output, ndat2d, ndat3d)
1009 #else
1010  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1011 #endif
1012 
1013 #elif OUTPUT==2
1014 
1015 if (time_output(1) <= time_init+eps) then
1016 
1017 #if ERGDAT==0
1018  flag_3d_output = .false.
1019 #elif ERGDAT==1
1020  flag_3d_output = .true.
1021 #else
1022  stop ' sico_init: ERGDAT must be either 0 or 1!'
1023 #endif
1024 
1025 #if NETCDF==1
1026  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1027  flag_3d_output, ndat2d, ndat3d)
1028 #elif NETCDF==2
1029  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1030  flag_3d_output, ndat2d, ndat3d)
1031 #else
1032  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1033 #endif
1034 
1035 end if
1036 
1037 #elif OUTPUT==3
1038 
1039  flag_3d_output = .false.
1040 
1041 #if NETCDF==1
1042  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1043  flag_3d_output, ndat2d, ndat3d)
1044 #elif NETCDF==2
1045  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1046  flag_3d_output, ndat2d, ndat3d)
1047 #else
1048  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1049 #endif
1050 
1051 if (time_output(1) <= time_init+eps) then
1052 
1053  flag_3d_output = .true.
1054 
1055 #if NETCDF==1
1056  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1057  flag_3d_output, ndat2d, ndat3d)
1058 #elif NETCDF==2
1059  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1060  flag_3d_output, ndat2d, ndat3d)
1061 #else
1062  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1063 #endif
1064 
1065 end if
1066 
1067 #else
1068  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1069 #endif
1070 
1071 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1072 #if OUTSER==2
1073 call output3(time_init, delta_ts, glac_index, z_sl)
1074 #endif
1075 #if OUTSER==3
1076 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1077 #endif
1078 
1079 end subroutine sico_init
1080 !
Generated on Thu Sep 12 2013 09:44:59 for SICOPOLIS V3.1 by   doxygen 1.8.2