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