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