SICOPOLIS V3.2
 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-2016 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, &
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 use sico_vars
48 use ice_material_quantities, only : ice_mat_eqs_pars
49 use enth_temp_omega, only : calc_c_int_table, calc_c_int_inv_table, &
50  enth_fct_temp_omega
51 
52 implicit none
53 
54 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
55  jj((imax+1)*(jmax+1)), &
56  nn(0:jmax,0:imax)
57 integer(i4b), intent(out) :: ndat2d, ndat3d
58 integer(i4b), intent(out) :: n_output
59 real(dp), intent(out) :: delta_ts, glac_index
60 real(dp), intent(out) :: mean_accum
61 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
62  dtime_out, dtime_ser
63 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
64 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
65 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
66 character(len=100), intent(out) :: runname
67 
68 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
69 integer(i4b) :: ios
70 integer(i4b) :: ierr
71 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
72 real(dp) :: time_init0, time_end0, time_output0(100)
73 real(dp) :: d_dummy
74 character (len=100) :: anfdatname
75 character (len=256) :: shell_command
76 character :: ch_dummy
77 logical :: flag_3d_output
78 
79 character(len=64), parameter :: fmt1 = '(a)', &
80  fmt2 = '(a,i4)', &
81  fmt3 = '(a,es12.4)'
82 
83 character(len= 8) :: ch_imax
84 character(len=128) :: fmt4
85 
86 write(ch_imax, fmt='(i8)') imax
87 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
88 
89 !-------- Name of the computational domain --------
90 
91 #if defined(ANT)
92 ch_domain_long = 'Antarctica'
93 ch_domain_short = 'ant'
94 
95 #elif defined(ASF)
96 ch_domain_long = 'Austfonna'
97 ch_domain_short = 'asf'
98 
99 #elif defined(EMTP2SGE)
100 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
101 ch_domain_short = 'emtp2sge'
102 
103 #elif defined(GRL)
104 ch_domain_long = 'Greenland'
105 ch_domain_short = 'grl'
106 
107 #elif defined(NHEM)
108 ch_domain_long = 'Northern hemisphere'
109 ch_domain_short = 'nhem'
110 
111 #elif defined(SCAND)
112 ch_domain_long = 'Scandinavia and Eurasia'
113 ch_domain_short = 'scand'
114 
115 #elif defined(TIBET)
116 ch_domain_long = 'Tibet'
117 ch_domain_short = 'tibet'
118 
119 #elif defined(NMARS)
120 ch_domain_long = 'North polar cap of Mars'
121 ch_domain_short = 'nmars'
122 
123 #elif defined(SMARS)
124 ch_domain_long = 'South polar cap of Mars'
125 ch_domain_short = 'smars'
126 
127 #elif (defined(XYZ))
128 ch_domain_long = 'XYZ'
129 ch_domain_short = 'xyz'
130 #if (defined(HEINO))
131 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
132 #endif
133 
134 #else
135 stop ' sico_init: No valid domain specified!'
136 #endif
137 
138 !-------- Some initial values --------
139 
140 n_output = 0
141 
142 dtime = 0.0_dp
143 dtime_temp = 0.0_dp
144 dtime_wss = 0.0_dp
145 dtime_out = 0.0_dp
146 dtime_ser = 0.0_dp
147 
148 time = 0.0_dp
149 time_init = 0.0_dp
150 time_end = 0.0_dp
151 time_output = 0.0_dp
152 
153 !-------- Initialisation of the Library of Iterative Solvers Lis,
154 ! if required --------
155 
156 #if (CALCZS==4 || MARGIN==3 || DYNAMICS==2)
157 
158 #include "lisf.h"
159 call lis_initialize(ierr)
160 
161 #endif
162 
163 !-------- Read physical parameters --------
164 
165 call phys_para()
166 
167 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
168 
169 ! ------ Some auxiliary quantities required for the enthalpy method
170 
171 call calc_c_int_table(c, -190, 10, l)
172 call calc_c_int_inv_table()
173 
174 !-------- Compatibility check of the SICOPOLIS version with the header file
175 
176 if ( trim(version) /= trim(sico_version) ) &
177  stop ' sico_init: SICOPOLIS version not compatible with header file!'
178 
179 !-------- Check whether the dynamics and thermodynamics modes are defined
180 
181 #if ( !defined(DYNAMICS) )
182 stop ' sico_init: DYNAMICS not defined in the header file!'
183 #endif
184 
185 #if ( !defined(CALCMOD) )
186 stop ' sico_init: CALCMOD not defined in the header file!'
187 #endif
188 
189 #if ( defined(ENTHMOD) )
190 write(6, fmt='(a)') ' sico_init: ENTHMOD must not be defined any more.'
191 write(6, fmt='(a)') ' Please update your header file!'
192 stop
193 #endif
194 
195 !-------- Compatibility check of the horizontal resolution with the
196 ! number of grid points --------
197 
198 #if (GRID==0 || GRID==1)
199 
200  stop ' sico_init: GRID==0 or 1 not allowed for the Tibet application!'
201 
202 #elif GRID==2
203 
204 if ( (abs(dlambda-1.0_dp/3.0_dp) < eps) &
205  .and.(abs(dphi- 1.0_dp/3.0_dp) < eps) ) then
206 
207  if ((imax /= 135).or.(jmax /= 51)) then
208  stop ' sico_init: IMAX and/or JMAX wrong!'
209  end if
210 
211 else if ( (abs(dlambda-1.0_dp/6.0_dp) < eps) &
212  .and.(abs(dphi -1.0_dp/6.0_dp) < eps) ) then
213 
214  if ((imax /= 270).or.(jmax /= 102)) then
215  stop ' sico_init: IMAX and/or JMAX wrong!'
216  end if
217 
218 else if ( (abs(dlambda-1.0_dp/12.0_dp) < eps) &
219  .and.(abs(dphi -1.0_dp/12.0_dp) < eps) ) then
220 
221  if ((imax /= 540).or.(jmax /= 204)) then
222  stop ' sico_init: IMAX and/or JMAX wrong!'
223  end if
224 
225 else
226  stop ' sico_init: DLAMBDA / DPHI wrong!'
227 end if
228 
229 #endif
230 
231 !-------- Compatibility check of the thermodynamics mode
232 ! (cold vs. polythermal vs. enthalpy method)
233 ! and the number of grid points in the lower (kt) ice domain --------
234 
235 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
236 
237 if (ktmax /= 1) then
238  write(6, fmt='(a)') ' sico_init: For options CALCMOD==0, 2, 3 or -1,'
239  write(6, fmt='(a)') ' the separate kt domain is redundant.'
240  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 1.'
241  write(6, fmt='(a)') ' '
242 end if
243 
244 #endif
245 
246 !-------- Compatibility check of surface-temperature and precipitation
247 ! determination by interpolation between present and LGM values
248 ! with a glacial index --------
249 
250 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
251 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
252 #endif
253 
254 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
255 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
256 #endif
257 
258 !-------- Compatibility check of discretization schemes for the horizontal and
259 ! vertical advection terms in the temperature and age equations --------
260 
261 #if (ADV_HOR==1)
262 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
263 #endif
264 
265 !-------- Check whether for the shallow shelf
266 ! or shelfy stream approximation
267 ! the chosen grid is Cartesian coordinates
268 ! without distortion correction (GRID==0) --------
269 
270 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
271 #if (GRID != 0)
272 write(6, fmt='(a)') ' sico_init: Distortion correction not implemented'
273 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
274 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
275 write(6, fmt='(a)') ' -> GRID==0 required!'
276 stop
277 #endif
278 #endif
279 
280 !-------- Setting of forcing flag --------
281 
282 #if (TSURFACE <= 4)
283 
284 forcing_flag = 1 ! forcing by delta_ts
285 
286 #elif (TSURFACE == 5)
287 
288 forcing_flag = 2 ! forcing by glac_index
289 
290 #endif
291 
292 !-------- Initialization of numerical time steps --------
293 
294 dtime0 = dtime0
295 dtime_temp0 = dtime_temp0
296 #if (REBOUND==2)
297 dtime_wss0 = dtime_wss0
298 #endif
299 
300 !-------- Further initializations --------
301 
302 dzeta_c = 1.0_dp/real(kcmax,dp)
303 dzeta_t = 1.0_dp/real(ktmax,dp)
304 dzeta_r = 1.0_dp/real(krmax,dp)
305 
306 ndat2d = 1
307 ndat3d = 1
308 
309 !-------- General abbreviations --------
310 
311 ! ------ kc domain
312 
313 if (deform >= eps) then
314 
315  flag_aa_nonzero = .true. ! non-equidistant grid
316 
317  aa = deform
318  ea = exp(aa)
319 
320  kc=0
321  zeta_c(kc) = 0.0_dp
322  eaz_c(kc) = 1.0_dp
323  eaz_c_quotient(kc) = 0.0_dp
324 
325  do kc=1, kcmax-1
326  zeta_c(kc) = kc*dzeta_c
327  eaz_c(kc) = exp(aa*zeta_c(kc))
328  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
329  end do
330 
331  kc=kcmax
332  zeta_c(kc) = 1.0_dp
333  eaz_c(kc) = exp(aa)
334  eaz_c_quotient(kc) = 1.0_dp
335 
336 else
337 
338  flag_aa_nonzero = .false. ! equidistant grid
339 
340  aa = 0.0_dp
341  ea = 1.0_dp
342 
343  kc=0
344  zeta_c(kc) = 0.0_dp
345  eaz_c(kc) = 1.0_dp
346  eaz_c_quotient(kc) = 0.0_dp
347 
348  do kc=1, kcmax-1
349  zeta_c(kc) = kc*dzeta_c
350  eaz_c(kc) = 1.0_dp
351  eaz_c_quotient(kc) = zeta_c(kc)
352  end do
353 
354  kc=kcmax
355  zeta_c(kc) = 1.0_dp
356  eaz_c(kc) = 1.0_dp
357  eaz_c_quotient(kc) = 1.0_dp
358 
359 end if
360 
361 ! ------ kt domain
362 
363 kt=0
364 zeta_t(kt) = 0.0_dp
365 
366 do kt=1, ktmax-1
367  zeta_t(kt) = kt*dzeta_t
368 end do
369 
370 kt=ktmax
371 zeta_t(kt) = 1.0_dp
372 
373 ! ------ kr domain
374 
375 kr=0
376 zeta_r(kr) = 0.0_dp
377 
378 do kr=1, krmax-1
379  zeta_r(kr) = kr*dzeta_r
380 end do
381 
382 kr=krmax
383 zeta_r(kr) = 1.0_dp
384 
385 !-------- Construction of a vector (with index n) from a 2-d array
386 ! (with indices i, j) by diagonal numbering --------
387 
388 n=1
389 
390 do m=0, imax+jmax
391  do i=m, 0, -1
392  j = m-i
393  if ((i <= imax).and.(j <= jmax)) then
394  ii(n) = i
395  jj(n) = j
396  nn(j,i) = n
397  n=n+1
398  end if
399  end do
400 end do
401 
402 !-------- Specification of current simulation --------
403 
404 runname = runname
405 anfdatname = anfdatname
406 
407 #if defined(YEAR_ZERO)
408 year_zero = year_zero
409 #else
410 year_zero = 2000.0_dp ! default value 2000 CE
411 #endif
412 
413 time_init0 = time_init0
414 time_end0 = time_end0
415 dtime_ser0 = dtime_ser0
416 
417 #if ( OUTPUT==1 || OUTPUT==3 )
418 dtime_out0 = dtime_out0
419 #endif
420 
421 #if ( OUTPUT==2 || OUTPUT==3 )
422 n_output = n_output
423 time_output0( 1) = time_out0_01
424 time_output0( 2) = time_out0_02
425 time_output0( 3) = time_out0_03
426 time_output0( 4) = time_out0_04
427 time_output0( 5) = time_out0_05
428 time_output0( 6) = time_out0_06
429 time_output0( 7) = time_out0_07
430 time_output0( 8) = time_out0_08
431 time_output0( 9) = time_out0_09
432 time_output0(10) = time_out0_10
433 time_output0(11) = time_out0_11
434 time_output0(12) = time_out0_12
435 time_output0(13) = time_out0_13
436 time_output0(14) = time_out0_14
437 time_output0(15) = time_out0_15
438 time_output0(16) = time_out0_16
439 time_output0(17) = time_out0_17
440 time_output0(18) = time_out0_18
441 time_output0(19) = time_out0_19
442 time_output0(20) = time_out0_20
443 #endif
444 
445 !-------- Write log file --------
446 
447 shell_command = 'if [ ! -d'
448 shell_command = trim(shell_command)//' '//outpath
449 shell_command = trim(shell_command)//' '//'] ; then mkdir'
450 shell_command = trim(shell_command)//' '//outpath
451 shell_command = trim(shell_command)//' '//'; fi'
452 call system(trim(shell_command))
453  ! Check whether directory OUTPATH exists. If not, it is created.
454 
455 open(10, iostat=ios, &
456  file=outpath//'/'//trim(runname)//'.log', &
457  status='new')
458 if (ios /= 0) &
459  stop ' sico_init: Error when opening the log file!'
460 
461 write(10, fmt=trim(fmt1)) 'Computational domain:'
462 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
463 write(10, fmt=trim(fmt1)) ' '
464 
465 write(10, fmt=trim(fmt2)) 'imax =', imax
466 write(10, fmt=trim(fmt2)) 'jmax =', jmax
467 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
468 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
469 write(10, fmt=trim(fmt2)) 'krmax =', krmax
470 write(10, fmt=trim(fmt1)) ' '
471 
472 write(10, fmt=trim(fmt3)) 'a =', aa
473 write(10, fmt=trim(fmt1)) ' '
474 
475 #if (GRID==0 || GRID==1)
476 stop ' sico_init: GRID==0 or 1 not allowed for this application!'
477 #elif GRID==2
478 write(10, fmt=trim(fmt3)) 'lambda0 =', lambda_0
479 write(10, fmt=trim(fmt3)) 'phi0 =', phi_0
480 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
481 write(10, fmt=trim(fmt3)) 'dphi =', dphi
482 #endif
483 write(10, fmt=trim(fmt1)) ' '
484 
485 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
486 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
487 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
488 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
489 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
490 #if (REBOUND==2)
491 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
492 #endif
493 #if ( OUTPUT==1 || OUTPUT==3 )
494 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
495 #endif
496 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
497 write(10, fmt=trim(fmt1)) ' '
498 
499 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
500 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
501 #if (ANF_DAT==1)
502 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
503 #endif
504 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
505 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
506 #if (ANF_DAT==1 && defined(TEMP_INIT))
507 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
508 #endif
509 #if (ANF_DAT==3)
510 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
511 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
512 #endif
513 write(10, fmt=trim(fmt1)) ' '
514 
515 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
516 write(10, fmt=trim(fmt1)) ' '
517 
518 #if (CALCZS==3 || CALCZS==4)
519 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
520 #if CALCZS==3
521 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
522 #endif
523 write(10, fmt=trim(fmt1)) ' '
524 #endif
525 
526 write(10, fmt=trim(fmt1)) 'temp_mm_present file = '//temp_mm_present_file
527 #if TSURFACE==1
528 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
529 #elif TSURFACE==3
530 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
531 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
532 #elif TSURFACE==4
533 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
534 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
535 #elif TSURFACE==5
536 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
537 write(10, fmt=trim(fmt1)) 'temp_mm_anom file = '//temp_mm_anom_file
538 write(10, fmt=trim(fmt3)) 'temp_mm_anom fact = ', temp_mm_anom_fact
539 #endif
540 
541 write(10, fmt=trim(fmt1)) 'precip_mm_present file = '//precip_mm_present_file
542 #if ACCSURFACE==1
543 write(10, fmt=trim(fmt3)) 'accfact =', accfact
544 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
545 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
546 #elif ( ACCSURFACE==5 )
547 write(10, fmt=trim(fmt1)) 'precip_mm_anom file = '//precip_mm_anom_file
548 write(10, fmt=trim(fmt3)) 'precip_mm_anom fact = ', precip_mm_anom_fact
549 #endif
550 #if (ACCSURFACE <= 3)
551 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
552 #if (ELEV_DESERT == 1)
553 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
554 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
555 #endif
556 #endif
557 
558 #if ABLSURFACE==3
559 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
560 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
561 #endif
562 
563 #if (SEA_LEVEL==1)
564 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
565 #elif (SEA_LEVEL==3)
566 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
567 #endif
568 write(10, fmt=trim(fmt1)) ' '
569 
570 #if (MARGIN==2)
571 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
572 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
573 write(10, fmt=trim(fmt1)) ' '
574 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
575  || marine_ice_calving==6 || marine_ice_calving==7 )
576 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
577 write(10, fmt=trim(fmt1)) ' '
578 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
579 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
580 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
581 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
582 write(10, fmt=trim(fmt1)) ' '
583 #endif
584 #elif (MARGIN==3)
585 #if ICE_SHELF_CALVING==2
586 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
587 write(10, fmt=trim(fmt1)) ' '
588 #endif
589 #endif
590 
591 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
592 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
593 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
594 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
595 #if SLIDE_LAW==2
596 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
597 #endif
598 write(10, fmt=trim(fmt1)) ' '
599 
600 #if Q_GEO_MOD==1
601 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
602 #elif Q_GEO_MOD==2
603 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
604 #endif
605 write(10, fmt=trim(fmt1)) ' '
606 
607 #if (defined(MARINE_ICE_BASAL_MELTING))
608 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
609 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
610 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
611 #endif
612 write(10, fmt=trim(fmt1)) ' '
613 #endif
614 
615 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
616 #if (REBOUND==1)
617 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
618 #endif
619 #if (REBOUND==1 || REBOUND==2)
620 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
621 #if (TIME_LAG_MOD==1)
622 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
623 #elif (TIME_LAG_MOD==2)
624 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
625 #else
626 stop ' sico_init: TIME_LAG_MOD must be either 1 or 2!'
627 #endif
628 #endif
629 #if (REBOUND==2)
630 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
631 #if (FLEX_RIG_MOD==1)
632 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
633 #elif (FLEX_RIG_MOD==2)
634 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
635 #else
636 stop ' sico_init: FLEX_RIG_MOD must be either 1 or 2!'
637 #endif
638 #endif
639 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
640 write(10, fmt=trim(fmt1)) ' '
641 
642 #if FLOW_LAW==2
643 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
644 write(10, fmt=trim(fmt1)) ' '
645 #endif
646 #if FIN_VISC==2
647 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
648 write(10, fmt=trim(fmt1)) ' '
649 #endif
650 
651 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
652 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
653 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
654 #endif
655 #if (ENHMOD==2 || ENHMOD==3)
656 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
657 #endif
658 #if (ENHMOD==2)
659 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
660 #endif
661 #if (ENHMOD==3)
662 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
663 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
664 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
665 #endif
666 #if (ENHMOD==4 || ENHMOD==5)
667 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
668 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
669 #endif
670 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
671 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
672 #endif
673 write(10, fmt=trim(fmt1)) ' '
674 
675 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
676 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
677 write(10, fmt=trim(fmt3)) 'H_min =', h_min
678 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
679 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
680 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
681 #if defined(QBM_MIN)
682 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
683 #elif defined(QB_MIN) /* obsolete */
684 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
685 #endif
686 #if defined(QBM_MAX)
687 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
688 #elif defined(QB_MAX) /* obsolete */
689 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
690 #endif
691 write(10, fmt=trim(fmt3)) 'age_min =', age_min
692 write(10, fmt=trim(fmt3)) 'age_max =', age_max
693 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
694 #if ADV_VERT==1
695 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
696 #endif
697 write(10, fmt=trim(fmt1)) ' '
698 
699 write(10, fmt=trim(fmt2)) 'GRID =', grid
700 write(10, fmt=trim(fmt2)) 'DYNAMICS =', dynamics
701 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
702 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
703 #endif
704 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
705 #if ( CALCMOD==-1 && defined(TEMP_CONST) )
706 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
707 #endif
708 #if ( CALCMOD==-1 && defined(AGE_CONST) )
709 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
710 #endif
711 #if ( CALCMOD==1 && defined(CTS_MELTING_FREEZING) )
712 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
713 #endif
714 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
715 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
716 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
717 #if (MARGIN==2)
718 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
719 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
720 #elif (MARGIN==3)
721 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
722 #endif
723 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
724 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
725 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
726 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
727 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
728 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
729 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
730 #if ACCSURFACE==5
731 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
732 #endif
733 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
734 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
735 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
736 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
737 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
738 write(10, fmt=trim(fmt1)) ' '
739 
740 #if defined(OUT_TIMES)
741 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
742 #endif
743 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
744 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
745 #if defined(OUTSER_V_A)
746 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
747 #endif
748 #if ( OUTPUT==1 || OUTPUT==2 )
749 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
750 #endif
751 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
752 write(10, fmt=trim(fmt1)) ' '
753 #if ( OUTPUT==2 || OUTPUT==3 )
754 write(10, fmt=trim(fmt2)) 'n_output =', n_output
755 do n=1, n_output
756  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
757 end do
758 write(10, fmt=trim(fmt1)) ' '
759 #endif
760 
761 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
762 
763 close(10, status='keep')
764 
765 !-------- Conversion of time quantities --------
766 
767 year_zero = year_zero*year_sec ! a --> s
768 time_init = time_init0*year_sec ! a --> s
769 time_end = time_end0*year_sec ! a --> s
770 dtime = dtime0*year_sec ! a --> s
771 dtime_temp = dtime_temp0*year_sec ! a --> s
772 #if (REBOUND==2)
773 dtime_wss = dtime_wss0*year_sec ! a --> s
774 #endif
775 dtime_ser = dtime_ser0*year_sec ! a --> s
776 #if ( OUTPUT==1 || OUTPUT==3 )
777 dtime_out = dtime_out0*year_sec ! a --> s
778 #endif
779 #if ( OUTPUT==2 || OUTPUT==3 )
780 do n=1, n_output
781  time_output(n) = time_output0(n)*year_sec ! a --> s
782 end do
783 #endif
784 
785 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
786  stop ' sico_init: dtime_temp must be a multiple of dtime!'
787 
788 #if (REBOUND==2)
789 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
790  stop ' sico_init: dtime_wss must be a multiple of dtime!'
791 #endif
792 
793 time = time_init
794 
795 !-------- Reading of measurements for present monthly-mean precipitation --------
796 
797 #if (GRID==0 || GRID==1)
798 
799 stop ' sico_init: GRID==0 or 1 not allowed for the Tibet application!'
800 
801 #elif GRID==2
802 
803 open(21, iostat=ios, &
804  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_mm_present_file, &
805  recl=16384, status='old')
806 
807 #endif
808 
809 if (ios /= 0) stop ' sico_init: Error when opening the precipitation file!'
810 
811 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
812 
813 do n=1, 12 ! month counter
814  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
815  do j=jmax, 0, -1
816  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
817  end do
818 end do
819 
820 close(21, status='keep')
821 
822 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
823 
824 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
825  ! mm/a water equiv. -> m/s ice equiv.
826 
827 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
828 
829 #if ACCSURFACE==5
830 
831 #if (GRID==2)
832 
833 open(21, iostat=ios, &
834  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_mm_anom_file, &
835  recl=16384, status='old')
836 
837 #endif
838 
839 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
840 
841 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
842 
843 do n=1, 12 ! month counter
844  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
845  do j=jmax, 0, -1
846  read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
847  end do
848 end do
849 
850 close(21, status='keep')
851 
852 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
853 
854 do i=0, imax
855 do j=0, jmax
856 
857 #if (PRECIP_ANOM_INTERPOL==1)
858  do n=1, 12 ! month counter
859  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
860  end do
861 #elif (PRECIP_ANOM_INTERPOL==2)
862  do n=1, 12 ! month counter
863  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
864  end do
865 #else
866  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
867 #endif
868 
869 end do
870 end do
871 
872 #endif
873 
874 !-------- Mean accumulation --------
875 
876 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
877  ! mm/a water equiv. -> m/s ice equiv.
878 
879 !-------- Reading of present topography mask --------
880 
881 #if (GRID==0 || GRID==1)
882 
883 stop ' sico_init: GRID==0 or 1 not allowed for the Tibet application!'
884 
885 #elif GRID==2
886 
887 open(24, iostat=ios, &
888  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_present_file, &
889  recl=2048, status='old')
890 
891 #endif
892 
893 if (ios /= 0) &
894  stop ' sico_init: Error when opening the mask file!'
895 
896 do m=1, 6; read(24, fmt='(a)') ch_dummy; end do
897 
898 do j=jmax, 0, -1
899  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
900 end do
901 
902 close(24, status='keep')
903 
904 !-------- Reading of data for present monthly-mean surface temperature --------
905 
906 #if (GRID==2)
907 
908 open(21, iostat=ios, &
909  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_mm_present_file, &
910  recl=16384, status='old')
911 
912 #endif
913 
914 if (ios /= 0) stop ' sico_init: Error when opening the temperature file!'
915 
916 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
917 
918 do n=1, 12 ! month counter
919  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
920  do j=jmax, 0, -1
921  read(21, fmt=*) (temp_mm_present(j,i,n), i=0,imax)
922  end do
923 end do
924 
925 close(21, status='keep')
926 
927 !-------- Reading of LGM monthly-mean surface-temperature anomalies --------
928 
929 #if TSURFACE==5
930 
931 #if (GRID==2)
932 
933 open(21, iostat=ios, &
934  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_mm_anom_file, &
935  recl=16384, status='old')
936 
937 #endif
938 
939 if (ios /= 0) stop ' sico_init: Error when opening the temperature anomaly file!'
940 
941 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
942 
943 do n=1, 12 ! month counter
944  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
945  do j=jmax, 0, -1
946  read(21, fmt=*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
947  end do
948 end do
949 
950 close(21, status='keep')
951 
952 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
953 
954 #endif
955 
956 !-------- Present reference elevation
957 ! (for precipitation and surface-temperature data) --------
958 
959 #if (GRID==0 || GRID==1)
960 
961 stop ' sico_init: GRID==0 or 1 not allowed for the Tibet application!'
962 
963 #elif GRID==2
964 
965 open(21, iostat=ios, &
966  file=inpath//'/'//trim(ch_domain_short)//'/'//zs_present_file, &
967  recl=8192, status='old')
968 
969 #endif
970 
971 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
972 
973 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
974 
975 do j=jmax, 0, -1
976  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
977 end do
978 
979 close(21, status='keep')
980 
981 ! ------ Reset bathymetry data (sea floor elevation) to the
982 ! sea surface
983 
984 do i=0, imax
985 do j=0, jmax
986  if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
987 end do
988 end do
989 
990 !-------- Read data for delta_ts --------
991 
992 #if TSURFACE==4
993 
994 open(21, iostat=ios, &
995  file=inpath//'/general/'//grip_temp_file, &
996  status='old')
997 if (ios /= 0) &
998  stop ' sico_init: Error when opening the data file for delta_ts!'
999 
1000 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1001 
1002 if (ch_dummy /= '#') then
1003  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
1004  write(6, fmt=*) ' not defined in data file for delta_ts!'
1005  stop
1006 end if
1007 
1008 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
1009 
1010 allocate(griptemp(0:ndata_grip))
1011 
1012 do n=0, ndata_grip
1013  read(21, fmt=*) d_dummy, griptemp(n)
1014 end do
1015 
1016 close(21, status='keep')
1017 
1018 #endif
1019 
1020 !-------- Read data for the glacial index --------
1021 
1022 #if TSURFACE==5
1023 
1024 open(21, iostat=ios, &
1025  file=inpath//'/general/'//glac_ind_file, status='old')
1026 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
1027 
1028 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1029 
1030 if (ch_dummy /= '#') then
1031  write(6, fmt=*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
1032  write(6, fmt=*) ' not defined in glacial-index file!'
1033  stop
1034 end if
1035 
1036 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
1037 
1038 allocate(glacial_index(0:ndata_gi))
1039 
1040 do n=0, ndata_gi
1041  read(21, fmt=*) d_dummy, glacial_index(n)
1042 end do
1043 
1044 close(21, status='keep')
1045 
1046 #endif
1047 
1048 !-------- Read data for z_sl --------
1049 
1050 #if (SEA_LEVEL==3)
1051 
1052 open(21, iostat=ios, &
1053  file=inpath//'/general/'//sea_level_file, &
1054  status='old')
1055 if (ios /= 0) &
1056  stop ' sico_init: Error when opening the data file for z_sl!'
1057 
1058 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1059 
1060 if (ch_dummy /= '#') then
1061  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1062  write(6, fmt=*) ' not defined in data file for z_sl!'
1063  stop
1064 end if
1065 
1066 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
1067 
1068 allocate(specmap_zsl(0:ndata_specmap))
1069 
1070 do n=0, ndata_specmap
1071  read(21, fmt=*) d_dummy, specmap_zsl(n)
1072 end do
1073 
1074 close(21, status='keep')
1075 
1076 #endif
1077 
1078 !-------- Determination of the geothermal heat flux --------
1079 
1080 #if Q_GEO_MOD==1
1081 
1082 ! ------ Constant value
1083 
1084 do i=0, imax
1085 do j=0, jmax
1086  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1087 end do
1088 end do
1089 
1090 #elif Q_GEO_MOD==2
1091 
1092 ! ------ Read data from file
1093 
1094 open(21, iostat=ios, &
1095  file=inpath//'/'//trim(ch_domain_short)//'/'//q_geo_file, &
1096  recl=8192, status='old')
1097 
1098 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
1099 
1100 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1101 
1102 do j=jmax, 0, -1
1103  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1104 end do
1105 
1106 close(21, status='keep')
1107 
1108 do i=0, imax
1109 do j=0, jmax
1110  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1111 end do
1112 end do
1113 
1114 #endif
1115 
1116 !-------- Reading of tabulated kei function--------
1117 
1118 #if (REBOUND==0 || REBOUND==1)
1119 
1120 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1121  ! dummy values
1122 #elif (REBOUND==2)
1123 
1124 call read_kei()
1125 
1126 #endif
1127 
1128 !-------- Determination of the time lag
1129 ! of the relaxing asthenosphere --------
1130 
1131 #if (REBOUND==1 || REBOUND==2)
1132 
1133 #if (TIME_LAG_MOD==1)
1134 
1135 time_lag_asth = time_lag*year_sec ! a -> s
1136 
1137 #elif (TIME_LAG_MOD==2)
1138 
1139 open(29, iostat=ios, &
1140  file=inpath//'/'//trim(ch_domain_short)//'/'//time_lag_file, &
1141  recl=8192, status='old')
1142 
1143 if (ios /= 0) stop ' sico_init: Error when opening the time-lag file!'
1144 
1145 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1146 
1147 do j=jmax, 0, -1
1148  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1149 end do
1150 
1151 close(29, status='keep')
1152 
1153 time_lag_asth = time_lag_asth*year_sec ! a -> s
1154 
1155 #endif
1156 
1157 #elif (REBOUND==0)
1158 
1159 time_lag_asth = 0.0_dp ! dummy values
1160 
1161 #endif
1162 
1163 !-------- Determination of the flexural rigidity of the lithosphere --------
1164 
1165 #if (REBOUND==2)
1166 
1167 #if (FLEX_RIG_MOD==1)
1168 
1169 flex_rig_lith = flex_rig
1170 
1171 #elif (FLEX_RIG_MOD==2)
1172 
1173 open(29, iostat=ios, &
1174  file=inpath//'/'//trim(ch_domain_short)//'/'//flex_rig_file, &
1175  recl=8192, status='old')
1176 
1177 if (ios /= 0) stop ' sico_init: Error when opening the flex-rig file!'
1178 
1179 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1180 
1181 do j=jmax, 0, -1
1182  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1183 end do
1184 
1185 close(29, status='keep')
1186 
1187 #endif
1188 
1189 #elif (REBOUND==0 || REBOUND==1)
1190 
1191 flex_rig_lith = 0.0_dp ! dummy values
1192 
1193 #endif
1194 
1195 !-------- Definition of initial values --------
1196 
1197 ! ------ Present topography
1198 
1199 #if (ANF_DAT==1)
1200 
1201 stop ' sico_init: ANF_DAT==1 not allowed for Tibet application!'
1202 
1203 ! ------ Ice-free, relaxed bedrock
1204 
1205 #elif (ANF_DAT==2)
1206 
1207 call topography2(dxi, deta)
1208 
1209 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1210 
1211 call boundary(time_init, dtime, dxi, deta, &
1212  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1213 
1214 q_bm = 0.0_dp
1215 q_tld = 0.0_dp
1216 q_b_tot = 0.0_dp
1217 
1218 p_b_w = 0.0_dp
1219 h_w = 0.0_dp
1220 
1221 call init_temp_water_age()
1222 
1223 call calc_enhance(time_init)
1224 
1225 ! ------ Read initial state from output of previous simulation
1226 
1227 #elif (ANF_DAT==3)
1228 
1229 #if (NETCDF==1)
1230 call topography3(dxi, deta, z_sl, anfdatname)
1231 #elif (NETCDF==2)
1232 call topography3_nc(dxi, deta, z_sl, anfdatname)
1233 #else
1234 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1235 #endif
1236 
1237 call boundary(time_init, dtime, dxi, deta, &
1238  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1239 
1240 q_b_tot = q_bm + q_tld
1241 
1242 #endif
1243 
1244 !-------- Grounding line and calving front flags --------
1245 
1246 flag_grounding_line_1 = .false.
1247 flag_grounding_line_2 = .false.
1248 
1249 flag_calving_front_1 = .false.
1250 flag_calving_front_2 = .false.
1251 
1252 #if (MARGIN==1 || MARGIN==2)
1253 
1254 continue
1255 
1256 #elif (MARGIN==3)
1257 
1258 do i=1, imax-1
1259 do j=1, jmax-1
1260 
1261  if ( (maske(j,i)==0_i2b) & ! grounded ice
1262  .and. &
1263  ( (maske(j,i+1)==3_i2b) & ! with
1264  .or.(maske(j,i-1)==3_i2b) & ! one
1265  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1266  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1267  ) &
1268  flag_grounding_line_1(j,i) = .true.
1269 
1270  if ( (maske(j,i)==3_i2b) & ! floating ice
1271  .and. &
1272  ( (maske(j,i+1)==0_i2b) & ! with
1273  .or.(maske(j,i-1)==0_i2b) & ! one
1274  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1275  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1276  ) &
1277  flag_grounding_line_2(j,i) = .true.
1278 
1279  if ( (maske(j,i)==3_i2b) & ! floating ice
1280  .and. &
1281  ( (maske(j,i+1)==2_i2b) & ! with
1282  .or.(maske(j,i-1)==2_i2b) & ! one
1283  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1284  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1285  ) &
1286  flag_calving_front_1(j,i) = .true.
1287 
1288  if ( (maske(j,i)==2_i2b) & ! ocean
1289  .and. &
1290  ( (maske(j,i+1)==3_i2b) & ! with
1291  .or.(maske(j,i-1)==3_i2b) & ! one
1292  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1293  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1294  ) &
1295  flag_calving_front_2(j,i) = .true.
1296 
1297 end do
1298 end do
1299 
1300 do i=1, imax-1
1301 
1302  j=0
1303  if ( (maske(j,i)==2_i2b) & ! ocean
1304  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1305  ) & ! floating ice point
1306  flag_calving_front_2(j,i) = .true.
1307 
1308  j=jmax
1309  if ( (maske(j,i)==2_i2b) & ! ocean
1310  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1311  ) & ! floating ice point
1312  flag_calving_front_2(j,i) = .true.
1313 
1314 end do
1315 
1316 do j=1, jmax-1
1317 
1318  i=0
1319  if ( (maske(j,i)==2_i2b) & ! ocean
1320  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1321  ) & ! floating ice point
1322  flag_calving_front_2(j,i) = .true.
1323 
1324  i=imax
1325  if ( (maske(j,i)==2_i2b) & ! ocean
1326  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1327  ) & ! floating ice point
1328  flag_calving_front_2(j,i) = .true.
1329 
1330 end do
1331 
1332 #else
1333 stop ' sico_init: MARGIN must be either 1, 2 or 3!'
1334 #endif
1335 
1336 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1337 
1338 #if (GRID==2)
1339 
1340 do ir=-imax, imax
1341 do jr=-jmax, jmax
1342 
1343  dist_dxdy(jr,ir) = sqrt( (sq_g11_g(jmax/2,imax/2)*real(ir,dp)*dxi)**2 &
1344  + (sq_g22_g(jmax/2,imax/2)*real(jr,dp)*deta)**2 )
1345 
1346  ! This uses the metric tensor in the center of the Tibet
1347  ! domain for the entire domain; DIRTY TRICK !!!
1348 
1349 end do
1350 end do
1351 
1352 #endif
1353 
1354 !-------- Initial velocities --------
1355 
1356 call calc_temp_melt()
1357 
1358 #if (DYNAMICS==1 || DYNAMICS==2)
1359 
1360 call calc_vxy_b_sia(time, z_sl)
1361 call calc_vxy_sia(dzeta_c, dzeta_t)
1362 
1363 #if (MARGIN==3 || DYNAMICS==2)
1364 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1365 #endif
1366 
1367 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1368 
1369 #if (MARGIN==3)
1370 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1371 #endif
1372 
1373 #elif (DYNAMICS==0)
1374 
1375 call calc_vxy_static()
1376 call calc_vz_static()
1377 
1378 #else
1379  stop ' sico_init: DYNAMICS must be either 0, 1 or 2!'
1380 #endif
1381 
1382 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1383 
1384 !-------- Initial enthalpies --------
1385 
1386 #if (CALCMOD==0 || CALCMOD==-1)
1387 
1388 do i=0, imax
1389 do j=0, jmax
1390 
1391  do kc=0, kcmax
1392  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1393  end do
1394 
1395  do kt=0, ktmax
1396  enth_t(kt,j,i) = enth_c(0,j,i)
1397  end do
1398 
1399 end do
1400 end do
1401 
1402 #elif (CALCMOD==1)
1403 
1404 do i=0, imax
1405 do j=0, jmax
1406 
1407  do kc=0, kcmax
1408  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1409  end do
1410 
1411  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1412  do kt=0, ktmax
1413  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1414  end do
1415  else
1416  do kt=0, ktmax
1417  enth_t(kt,j,i) = enth_c(0,j,i)
1418  end do
1419  end if
1420 
1421 end do
1422 end do
1423 
1424 #elif (CALCMOD==2 || CALCMOD==3)
1425 
1426 do i=0, imax
1427 do j=0, jmax
1428 
1429  do kc=0, kcmax
1430  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1431  end do
1432 
1433  do kt=0, ktmax
1434  enth_t(kt,j,i) = enth_c(0,j,i)
1435  end do
1436 
1437 end do
1438 end do
1439 
1440 #else
1441 
1442 stop ' sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1443 
1444 #endif
1445 
1446 !-------- Initialize time-series files --------
1447 
1448 ! ------ Time-series file for the ice sheet on the whole
1449 
1450 open(12, iostat=ios, &
1451  file=outpath//'/'//trim(runname)//'.ser', &
1452  status='new')
1453 if (ios /= 0) &
1454  stop ' sico_init: Error when opening the ser file!'
1455 
1456 if (forcing_flag == 1) then
1457 
1458  write(12,1102)
1459  write(12,1103)
1460 
1461 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1462 
1463  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1464  ' H_max(m) zs_max(m) V_g(m^3)', &
1465  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1466  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1467  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1468  1103 format('----------------------------------------------------', &
1469  '---------------------------------------')
1470 
1471 #elif OUTSER_V_A==2
1472 
1473  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1474  ' V(m^3) V_g(m^3) V_f(m^3)', &
1475  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1476  ' H_max(m) zs_max(m)', &
1477  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1478  ' A_t(m^2) V_bm(m^3/a)', &
1479  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1480  1103 format('----------------------------------------------------', &
1481  '---------------------------------------')
1482 
1483 #else
1484  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1485 #endif
1486 
1487 else if (forcing_flag == 2) then
1488 
1489  write(12,1112)
1490  write(12,1113)
1491 
1492 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1493 
1494  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1495  ' H_max(m) zs_max(m) V_g(m^3)', &
1496  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1497  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1498  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1499  1113 format('----------------------------------------------------', &
1500  '---------------------------------------')
1501 
1502 #elif OUTSER_V_A==2
1503 
1504  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1505  ' V(m^3) V_g(m^3) V_f(m^3)', &
1506  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1507  ' H_max(m) zs_max(m)', &
1508  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1509  ' A_t(m^2) V_bm(m^3/a)', &
1510  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1511  1113 format('----------------------------------------------------', &
1512  '---------------------------------------')
1513 
1514 #else
1515  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1516 #endif
1517 
1518 end if
1519 
1520 ! ------ Time-series file for the ice-thickness field
1521 
1522 #if OUTSER==2
1523 
1524 open(13, iostat=ios, &
1525  file=outpath//'/'//trim(runname)//'.thk', &
1526  status='new')
1527 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1528 
1529 if (forcing_flag == 1) then
1530 
1531  write(13,1104)
1532  write(13,1105)
1533 
1534  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1535  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1536  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1537  1105 format('----------------------------------------------------')
1538 
1539 else if (forcing_flag == 2) then
1540 
1541  write(13,1114)
1542  write(13,1115)
1543 
1544  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1545  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1546  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1547  1115 format('----------------------------------------------------')
1548 
1549 end if
1550 
1551 #endif
1552 
1553 !-------- Output of the initial state --------
1554 
1555 #if OUTPUT==1
1556 
1557 #if ERGDAT==0
1558  flag_3d_output = .false.
1559 #elif ERGDAT==1
1560  flag_3d_output = .true.
1561 #else
1562  stop ' sico_init: ERGDAT must be either 0 or 1!'
1563 #endif
1564 
1565 #if NETCDF==1
1566  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1567  flag_3d_output, ndat2d, ndat3d)
1568 #elif NETCDF==2
1569  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1570  flag_3d_output, ndat2d, ndat3d)
1571 #else
1572  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1573 #endif
1574 
1575 #elif OUTPUT==2
1576 
1577 if (time_output(1) <= time_init+eps) then
1578 
1579 #if ERGDAT==0
1580  flag_3d_output = .false.
1581 #elif ERGDAT==1
1582  flag_3d_output = .true.
1583 #else
1584  stop ' sico_init: ERGDAT must be either 0 or 1!'
1585 #endif
1586 
1587 #if NETCDF==1
1588  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1589  flag_3d_output, ndat2d, ndat3d)
1590 #elif NETCDF==2
1591  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1592  flag_3d_output, ndat2d, ndat3d)
1593 #else
1594  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1595 #endif
1596 
1597 end if
1598 
1599 #elif OUTPUT==3
1600 
1601  flag_3d_output = .false.
1602 
1603 #if NETCDF==1
1604  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1605  flag_3d_output, ndat2d, ndat3d)
1606 #elif NETCDF==2
1607  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1608  flag_3d_output, ndat2d, ndat3d)
1609 #else
1610  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1611 #endif
1612 
1613 if (time_output(1) <= time_init+eps) then
1614 
1615  flag_3d_output = .true.
1616 
1617 #if NETCDF==1
1618  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1619  flag_3d_output, ndat2d, ndat3d)
1620 #elif NETCDF==2
1621  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1622  flag_3d_output, ndat2d, ndat3d)
1623 #else
1624  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1625 #endif
1626 
1627 end if
1628 
1629 #else
1630  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1631 #endif
1632 
1633 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1634 #if OUTSER==2
1635 call output3(time_init, delta_ts, glac_index, z_sl)
1636 #endif
1637 
1638 end subroutine sico_init
1639 !
calc_vz_ssa
subroutine calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz in the shallow shelf approximation.
Definition: calc_vz_ssa.F90:35
calc_vxy_static
subroutine calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
Definition: calc_vxy_static.F90:37
calc_vxy_sia
subroutine calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_sia.F90:41
ice_material_quantities::ice_mat_eqs_pars
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_temp_min, n_temp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
Definition: ice_material_quantities.F90:74
output2
subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format.
Definition: output2.F90:35
phys_para
subroutine phys_para()
Reading of physical parameters.
Definition: phys_para.F90:35
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
topography3
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography3.F90:39
output3
subroutine output3(time, delta_ts, glac_index, z_sl)
Writing of time-series data of the ice thickness field on file in ASCII format.
Definition: output3.F90:38
topography3_nc
subroutine topography3_nc(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography3_nc.F90:39
calc_temp_melt
subroutine calc_temp_melt()
Computation of the melting temperatures.
Definition: calc_temp_melt.F90:35
output_nc
subroutine output_nc(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in NetCDF format.
Definition: output_nc.F90:35
calc_vxy_ssa
subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Definition: calc_vxy_ssa.F90:39
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
topography2
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography2.F90:39
calc_vz_static
subroutine calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_static.F90:35
calc_dxyz
subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz.F90:37
calc_enhance
subroutine calc_enhance(time)
Computation of the flow enhancement factor.
Definition: calc_enhance.F90:35
enth_temp_omega::calc_c_int_table
subroutine, public calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
Computation of the temperature integral of the specific heat of ice as a table (c_int_table). Further, definition of the latent heat of ice.
Definition: enth_temp_omega.F90:84
calc_vxy_b_sia
subroutine calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
Definition: calc_vxy_b_sia.F90:37
sico_init
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Initialisations for SICOPOLIS.
Definition: sico_init.F90:35
boundary
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary.F90:37
read_kei
subroutine read_kei()
Reading of the tabulated kei function.
Definition: read_kei.F90:35
init_temp_water_age
subroutine init_temp_water_age()
Initial temperature, water content and age.
Definition: init_temp_water_age.F90:35
ice_material_quantities
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Definition: ice_material_quantities.F90:37
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
enth_temp_omega::calc_c_int_inv_table
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
Definition: enth_temp_omega.F90:141
enth_temp_omega
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
Definition: enth_temp_omega.F90:37
calc_vz_sia
subroutine calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz in the shallow ice approximation.
Definition: calc_vz_sia.F90:35
output1
subroutine output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in binary format.
Definition: output1.F90:35
enth_temp_omega::enth_fct_temp_omega
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
Definition: enth_temp_omega.F90:259
Generated on Wed Mar 2 2016 11:26:50 for SICOPOLIS V3.2 by   doxygen 1.8.5