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 #if (NETCDF > 1)
53 use netcdf
54 use nc_check
55 #endif
56 
57 implicit none
58 
59 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
60  jj((imax+1)*(jmax+1)), &
61  nn(0:jmax,0:imax)
62 integer(i4b), intent(out) :: ndat2d, ndat3d
63 integer(i4b), intent(out) :: n_output
64 real(dp), intent(out) :: delta_ts, glac_index
65 real(dp), intent(out) :: mean_accum
66 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
67  dtime_out, dtime_ser
68 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
69 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
70 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
71 character(len=100), intent(out) :: runname
72 
73 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
74 integer(i4b) :: ios, ios1, ios2, ios3, ios4
75 integer(i4b) :: ierr
76 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
77 real(dp) :: time_init0, time_end0, time_output0(100)
78 real(dp), dimension(0:JMAX,0:IMAX) :: precip_factor
79 real(dp) :: d_dummy
80 character(len=100) :: anfdatname, target_topo_dat_name
81 character(len=256) :: filename_with_path
82 character(len=256) :: shell_command
83 character :: ch_dummy
84 logical :: flag_3d_output
85 
86 #if (NETCDF > 1)
87 integer(i4b) :: dimid
88 ! dimid: Dimension ID
89 integer(i4b) :: ncv
90 ! ncv: Variable ID
91 #endif
92 
93 character(len=64), parameter :: thisroutine = 'sico_init'
94 
95 character(len=64), parameter :: fmt1 = '(a)', &
96  fmt2 = '(a,i4)', &
97  fmt3 = '(a,es12.4)'
98 
99 character(len= 8) :: ch_imax
100 character(len=128) :: fmt4
101 
102 write(ch_imax, fmt='(i8)') imax
103 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
104 
105 !-------- Name of the computational domain --------
106 
107 #if defined(ANT)
108 ch_domain_long = 'Antarctica'
109 ch_domain_short = 'ant'
110 
111 #elif defined(ASF)
112 ch_domain_long = 'Austfonna'
113 ch_domain_short = 'asf'
114 
115 #elif defined(EMTP2SGE)
116 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
117 ch_domain_short = 'emtp2sge'
118 
119 #elif defined(GRL)
120 ch_domain_long = 'Greenland'
121 ch_domain_short = 'grl'
122 
123 #elif defined(NHEM)
124 ch_domain_long = 'Northern hemisphere'
125 ch_domain_short = 'nhem'
126 
127 #elif defined(SCAND)
128 ch_domain_long = 'Scandinavia and Eurasia'
129 ch_domain_short = 'scand'
130 
131 #elif defined(TIBET)
132 ch_domain_long = 'Tibet'
133 ch_domain_short = 'tibet'
134 
135 #elif defined(NMARS)
136 ch_domain_long = 'North polar cap of Mars'
137 ch_domain_short = 'nmars'
138 
139 #elif defined(SMARS)
140 ch_domain_long = 'South polar cap of Mars'
141 ch_domain_short = 'smars'
142 
143 #elif (defined(XYZ))
144 ch_domain_long = 'XYZ'
145 ch_domain_short = 'xyz'
146 #if (defined(HEINO))
147 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
148 #endif
149 
150 #else
151 stop ' sico_init: No valid domain specified!'
152 #endif
153 
154 !-------- Some initial values --------
155 
156 n_output = 0
157 
158 dtime = 0.0_dp
159 dtime_temp = 0.0_dp
160 dtime_wss = 0.0_dp
161 dtime_out = 0.0_dp
162 dtime_ser = 0.0_dp
163 
164 time = 0.0_dp
165 time_init = 0.0_dp
166 time_end = 0.0_dp
167 time_output = 0.0_dp
168 
169 !-------- Initialisation of the Library of Iterative Solvers Lis,
170 ! if required --------
171 
172 #if (CALCZS==4 || MARGIN==3 || DYNAMICS==2)
173 
174 #include "lisf.h"
175 call lis_initialize(ierr)
176 
177 #endif
178 
179 !-------- Read physical parameters --------
180 
181 call phys_para()
182 
183 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
184 
185 ! ------ Some auxiliary quantities required for the enthalpy method
186 
187 call calc_c_int_table(c, -190, 10, l)
188 call calc_c_int_inv_table()
189 
190 !-------- Compatibility check of the SICOPOLIS version with the header file
191 
192 if ( trim(version) /= trim(sico_version) ) &
193  stop ' sico_init: SICOPOLIS version not compatible with header file!'
194 
195 !-------- Check whether the dynamics and thermodynamics modes are defined
196 
197 #if ( !defined(DYNAMICS) )
198 stop ' sico_init: DYNAMICS not defined in the header file!'
199 #endif
200 
201 #if ( !defined(CALCMOD) )
202 stop ' sico_init: CALCMOD not defined in the header file!'
203 #endif
204 
205 #if ( defined(ENTHMOD) )
206 write(6, fmt='(a)') ' sico_init: ENTHMOD must not be defined any more.'
207 write(6, fmt='(a)') ' Please update your header file!'
208 stop
209 #endif
210 
211 !-------- Compatibility check of the horizontal resolution with the
212 ! number of grid points --------
213 
214 #if (GRID==0 || GRID==1)
215 
216 if (abs(dx-20.0_dp) < eps) then
217 
218  if ( (imax /= 75).or.(jmax /= 140) ) then
219  stop ' sico_init: IMAX and/or JMAX wrong!'
220  end if
221 
222 else if (abs(dx-10.0_dp) < eps) then
223 
224  if ( (imax /= 150).or.(jmax /= 280) ) then
225  stop ' sico_init: IMAX and/or JMAX wrong!'
226  end if
227 
228 else if (abs(dx-5.0_dp) < eps) then
229 
230  if ( (imax /= 300).or.(jmax /= 560) ) then
231  stop ' sico_init: IMAX and/or JMAX wrong!'
232  end if
233 
234 else
235  stop ' sico_init: DX wrong!'
236 end if
237 
238 #elif GRID==2
239 
240 stop ' sico_init: GRID==2 not allowed for this application!'
241 
242 #endif
243 
244 !-------- Compatibility check of the thermodynamics mode
245 ! (cold vs. polythermal vs. enthalpy method)
246 ! and the number of grid points in the lower (kt) ice domain --------
247 
248 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
249 
250 if (ktmax /= 1) then
251  write(6, fmt='(a)') ' sico_init: For options CALCMOD==0, 2, 3 or -1,'
252  write(6, fmt='(a)') ' the separate kt domain is redundant.'
253  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 1.'
254  write(6, fmt='(a)') ' '
255 end if
256 
257 #endif
258 
259 !-------- Compatibility check of surface-temperature
260 ! and precipitation switches --------
261 
262 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
263 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
264 #endif
265 
266 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
267 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
268 #endif
269 
270 #if (TSURFACE == 6)
271 #if (ACCSURFACE != 6)
272 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
273 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
274 stop
275 #elif (NETCDF != 2)
276 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
277 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
278 stop
279 #endif
280 #endif
281 
282 #if (ACCSURFACE == 6)
283 #if (TSURFACE != 6)
284 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
285 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
286 stop
287 #elif (NETCDF != 2)
288 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
289 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
290 stop
291 #endif
292 #endif
293 
294 !-------- Compatibility check of discretization schemes for the horizontal and
295 ! vertical advection terms in the temperature and age equations --------
296 
297 #if (ADV_HOR==1)
298 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
299 #endif
300 
301 !-------- Check whether for the shallow shelf
302 ! or shelfy stream approximation
303 ! the chosen grid is Cartesian coordinates
304 ! without distortion correction (GRID==0) --------
305 
306 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
307 #if (GRID != 0)
308 write(6, fmt='(a)') ' sico_init: Distortion correction not implemented'
309 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
310 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
311 write(6, fmt='(a)') ' -> GRID==0 required!'
312 stop
313 #endif
314 #endif
315 
316 !-------- Setting of forcing flag --------
317 
318 #if (TSURFACE <= 4)
319 
320 forcing_flag = 1 ! forcing by delta_ts
321 
322 #elif (TSURFACE == 5)
323 
324 forcing_flag = 2 ! forcing by glac_index
325 
326 #elif (TSURFACE == 6)
327 
328 forcing_flag = 3 ! forcing by time-dependent surface temperature
329  ! and precipitation data
330 
331 #endif
332 
333 !-------- Initialization of numerical time steps --------
334 
335 dtime0 = dtime0
336 dtime_temp0 = dtime_temp0
337 #if (REBOUND==2)
338 dtime_wss0 = dtime_wss0
339 #endif
340 
341 !-------- Further initializations --------
342 
343 dzeta_c = 1.0_dp/real(kcmax,dp)
344 dzeta_t = 1.0_dp/real(ktmax,dp)
345 dzeta_r = 1.0_dp/real(krmax,dp)
346 
347 ndat2d = 1
348 ndat3d = 1
349 
350 !-------- General abbreviations --------
351 
352 ! ------ kc domain
353 
354 if (deform >= eps) then
355 
356  flag_aa_nonzero = .true. ! non-equidistant grid
357 
358  aa = deform
359  ea = exp(aa)
360 
361  kc=0
362  zeta_c(kc) = 0.0_dp
363  eaz_c(kc) = 1.0_dp
364  eaz_c_quotient(kc) = 0.0_dp
365 
366  do kc=1, kcmax-1
367  zeta_c(kc) = kc*dzeta_c
368  eaz_c(kc) = exp(aa*zeta_c(kc))
369  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
370  end do
371 
372  kc=kcmax
373  zeta_c(kc) = 1.0_dp
374  eaz_c(kc) = exp(aa)
375  eaz_c_quotient(kc) = 1.0_dp
376 
377 else
378 
379  flag_aa_nonzero = .false. ! equidistant grid
380 
381  aa = 0.0_dp
382  ea = 1.0_dp
383 
384  kc=0
385  zeta_c(kc) = 0.0_dp
386  eaz_c(kc) = 1.0_dp
387  eaz_c_quotient(kc) = 0.0_dp
388 
389  do kc=1, kcmax-1
390  zeta_c(kc) = kc*dzeta_c
391  eaz_c(kc) = 1.0_dp
392  eaz_c_quotient(kc) = zeta_c(kc)
393  end do
394 
395  kc=kcmax
396  zeta_c(kc) = 1.0_dp
397  eaz_c(kc) = 1.0_dp
398  eaz_c_quotient(kc) = 1.0_dp
399 
400 end if
401 
402 ! ------ kt domain
403 
404 kt=0
405 zeta_t(kt) = 0.0_dp
406 
407 do kt=1, ktmax-1
408  zeta_t(kt) = kt*dzeta_t
409 end do
410 
411 kt=ktmax
412 zeta_t(kt) = 1.0_dp
413 
414 ! ------ kr domain
415 
416 kr=0
417 zeta_r(kr) = 0.0_dp
418 
419 do kr=1, krmax-1
420  zeta_r(kr) = kr*dzeta_r
421 end do
422 
423 kr=krmax
424 zeta_r(kr) = 1.0_dp
425 
426 !-------- Construction of a vector (with index n) from a 2-d array
427 ! (with indices i, j) by diagonal numbering --------
428 
429 n=1
430 
431 do m=0, imax+jmax
432  do i=m, 0, -1
433  j = m-i
434  if ((i <= imax).and.(j <= jmax)) then
435  ii(n) = i
436  jj(n) = j
437  nn(j,i) = n
438  n=n+1
439  end if
440  end do
441 end do
442 
443 !-------- Specification of current simulation --------
444 
445 runname = runname
446 anfdatname = anfdatname
447 
448 #if defined(YEAR_ZERO)
449 year_zero = year_zero
450 #else
451 year_zero = 2000.0_dp ! default value 2000 CE
452 #endif
453 
454 time_init0 = time_init0
455 time_end0 = time_end0
456 dtime_ser0 = dtime_ser0
457 
458 #if ( OUTPUT==1 || OUTPUT==3 )
459 dtime_out0 = dtime_out0
460 #endif
461 
462 #if ( OUTPUT==2 || OUTPUT==3 )
463 n_output = n_output
464 time_output0( 1) = time_out0_01
465 time_output0( 2) = time_out0_02
466 time_output0( 3) = time_out0_03
467 time_output0( 4) = time_out0_04
468 time_output0( 5) = time_out0_05
469 time_output0( 6) = time_out0_06
470 time_output0( 7) = time_out0_07
471 time_output0( 8) = time_out0_08
472 time_output0( 9) = time_out0_09
473 time_output0(10) = time_out0_10
474 time_output0(11) = time_out0_11
475 time_output0(12) = time_out0_12
476 time_output0(13) = time_out0_13
477 time_output0(14) = time_out0_14
478 time_output0(15) = time_out0_15
479 time_output0(16) = time_out0_16
480 time_output0(17) = time_out0_17
481 time_output0(18) = time_out0_18
482 time_output0(19) = time_out0_19
483 time_output0(20) = time_out0_20
484 #endif
485 
486 !-------- Write log file --------
487 
488 shell_command = 'if [ ! -d'
489 shell_command = trim(shell_command)//' '//outpath
490 shell_command = trim(shell_command)//' '//'] ; then mkdir'
491 shell_command = trim(shell_command)//' '//outpath
492 shell_command = trim(shell_command)//' '//'; fi'
493 call system(trim(shell_command))
494  ! Check whether directory OUTPATH exists. If not, it is created.
495 
496 open(10, iostat=ios, &
497  file=outpath//'/'//trim(runname)//'.log', &
498  status='new')
499 if (ios /= 0) &
500  stop ' sico_init: Error when opening the log file!'
501 
502 write(10, fmt=trim(fmt1)) 'Computational domain:'
503 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
504 write(10, fmt=trim(fmt1)) ' '
505 
506 write(10, fmt=trim(fmt2)) 'imax =', imax
507 write(10, fmt=trim(fmt2)) 'jmax =', jmax
508 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
509 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
510 write(10, fmt=trim(fmt2)) 'krmax =', krmax
511 write(10, fmt=trim(fmt1)) ' '
512 
513 write(10, fmt=trim(fmt3)) 'a =', aa
514 write(10, fmt=trim(fmt1)) ' '
515 
516 #if (GRID==0 || GRID==1)
517 write(10, fmt=trim(fmt3)) 'x0 =', x0
518 write(10, fmt=trim(fmt3)) 'y0 =', y0
519 write(10, fmt=trim(fmt3)) 'dx =', dx
520 #elif GRID==2
521 stop ' sico_init: GRID==2 not allowed for this application!'
522 #endif
523 write(10, fmt=trim(fmt1)) ' '
524 
525 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
526 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
527 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
528 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
529 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
530 #if (REBOUND==2)
531 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
532 #endif
533 #if ( OUTPUT==1 || OUTPUT==3 )
534 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
535 #endif
536 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
537 write(10, fmt=trim(fmt1)) ' '
538 
539 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
540 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
541 #if (ANF_DAT==1)
542 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
543 #endif
544 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
545 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
546 #if (ANF_DAT==1 && defined(TEMP_INIT))
547 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
548 #endif
549 #if (ANF_DAT==3)
550 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
551 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
552 #endif
553 write(10, fmt=trim(fmt1)) ' '
554 
555 #if ZS_EVOL==2
556 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
557 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
558 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
559 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
560 write(10, fmt=trim(fmt1)) ' '
561 #endif
562 
563 #if ZS_EVOL==3
564 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
565 write(10, fmt=trim(fmt1)) ' '
566 #endif
567 
568 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
569 write(10, fmt=trim(fmt1)) ' '
570 
571 #if (CALCZS==3 || CALCZS==4)
572 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
573 #if (CALCZS==3)
574 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
575 #endif
576 write(10, fmt=trim(fmt1)) ' '
577 #endif
578 
579 #if (TSURFACE==1)
580 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
581 #elif (TSURFACE==3)
582 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
583 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
584 #elif (TSURFACE==4)
585 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
586 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
587 #elif (TSURFACE==5)
588 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
589 write(10, fmt=trim(fmt1)) 'temp_ma_anom file = '//temp_ma_anom_file
590 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
591 write(10, fmt=trim(fmt1)) 'temp_mj_anom file = '//temp_mj_anom_file
592 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
593 #elif (TSURFACE==6)
594 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
595 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
596 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
597 #endif
598 
599 write(10, fmt=trim(fmt1)) 'precip_present_file = '//precip_present_file
600 #if (defined(PRECIP_FACTOR_FILE))
601 if ( trim(adjustl(precip_factor_file)) /= 'none' ) &
602  write(10, fmt=trim(fmt1)) 'precip_factor_file = '//precip_factor_file
603 #endif
604 #if (ACCSURFACE==1)
605 write(10, fmt=trim(fmt3)) 'accfact =', accfact
606 #elif (ACCSURFACE==2 || ACCSURFACE==3)
607 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
608 #elif (ACCSURFACE==5)
609 write(10, fmt=trim(fmt1)) 'precip_anom file = '//precip_anom_file
610 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
611 #elif (ACCSURFACE==6)
612 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
613 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
614 #endif
615 #if (ACCSURFACE <= 3)
616 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
617 #if (ELEV_DESERT == 1)
618 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
619 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
620 #endif
621 #endif
622 
623 #if (ABLSURFACE==3)
624 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
625 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
626 #endif
627 
628 #if (SEA_LEVEL==1)
629 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
630 #elif (SEA_LEVEL==3)
631 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
632 #endif
633 write(10, fmt=trim(fmt1)) ' '
634 
635 #if (MARGIN==2)
636 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
637 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
638 write(10, fmt=trim(fmt1)) ' '
639 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
640  || marine_ice_calving==6 || marine_ice_calving==7 )
641 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
642 write(10, fmt=trim(fmt1)) ' '
643 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
644 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
645 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
646 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
647 write(10, fmt=trim(fmt1)) ' '
648 #endif
649 #elif (MARGIN==3)
650 #if ICE_SHELF_CALVING==2
651 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
652 write(10, fmt=trim(fmt1)) ' '
653 #endif
654 #endif
655 
656 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
657 write(10, fmt=trim(fmt3)) 'smw_coeff =', smw_coeff
658 write(10, fmt=trim(fmt2)) 'r_smw =', r_smw
659 write(10, fmt=trim(fmt2)) 's_smw =', s_smw
660 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
661 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
662 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
663 #if SLIDE_LAW==2
664 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
665 #endif
666 write(10, fmt=trim(fmt1)) ' '
667 
668 #if ICE_STREAM==2
669 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
670 write(10, fmt=trim(fmt1)) ' '
671 
672 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
673 write(10, fmt=trim(fmt3)) 'smw_coeff_sedi =', smw_coeff_sedi
674 write(10, fmt=trim(fmt2)) 'r_smw_sedi =', r_smw_sedi
675 write(10, fmt=trim(fmt2)) 's_smw_sedi =', s_smw_sedi
676 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
677 write(10, fmt=trim(fmt2)) 'p_weert_sedi =', p_weert_sedi
678 write(10, fmt=trim(fmt2)) 'q_weert_sedi =', q_weert_sedi
679 write(10, fmt=trim(fmt1)) ' '
680 #endif
681 
682 #if Q_GEO_MOD==1
683 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
684 #elif Q_GEO_MOD==2
685 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
686 #endif
687 write(10, fmt=trim(fmt1)) ' '
688 
689 #if (defined(MARINE_ICE_BASAL_MELTING))
690 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
691 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
692 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
693 #endif
694 write(10, fmt=trim(fmt1)) ' '
695 #endif
696 
697 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
698 #if (REBOUND==1)
699 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
700 #endif
701 #if (REBOUND==1 || REBOUND==2)
702 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
703 #if (TIME_LAG_MOD==1)
704 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
705 #elif (TIME_LAG_MOD==2)
706 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
707 #else
708 stop ' sico_init: TIME_LAG_MOD must be either 1 or 2!'
709 #endif
710 #endif
711 #if (REBOUND==2)
712 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
713 #if (FLEX_RIG_MOD==1)
714 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
715 #elif (FLEX_RIG_MOD==2)
716 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
717 #else
718 stop ' sico_init: FLEX_RIG_MOD must be either 1 or 2!'
719 #endif
720 #endif
721 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
722 write(10, fmt=trim(fmt1)) ' '
723 
724 #if FLOW_LAW==2
725 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
726 write(10, fmt=trim(fmt1)) ' '
727 #endif
728 #if FIN_VISC==2
729 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
730 write(10, fmt=trim(fmt1)) ' '
731 #endif
732 
733 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
734 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
735 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
736 #endif
737 #if (ENHMOD==2 || ENHMOD==3)
738 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
739 #endif
740 #if (ENHMOD==2)
741 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
742 #endif
743 #if (ENHMOD==3)
744 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
745 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
746 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
747 #endif
748 #if (ENHMOD==4 || ENHMOD==5)
749 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
750 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
751 #endif
752 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
753 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
754 #endif
755 write(10, fmt=trim(fmt1)) ' '
756 
757 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
758 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
759 write(10, fmt=trim(fmt3)) 'H_min =', h_min
760 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
761 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
762 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
763 #if defined(QBM_MIN)
764 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
765 #elif defined(QB_MIN) /* obsolete */
766 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
767 #endif
768 #if defined(QBM_MAX)
769 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
770 #elif defined(QB_MAX) /* obsolete */
771 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
772 #endif
773 write(10, fmt=trim(fmt3)) 'age_min =', age_min
774 write(10, fmt=trim(fmt3)) 'age_max =', age_max
775 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
776 #if ADV_VERT==1
777 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
778 #endif
779 write(10, fmt=trim(fmt1)) ' '
780 
781 write(10, fmt=trim(fmt2)) 'GRID =', grid
782 write(10, fmt=trim(fmt2)) 'DYNAMICS =', dynamics
783 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
784 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
785 #endif
786 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
787 #if ( CALCMOD==-1 && defined(TEMP_CONST) )
788 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
789 #endif
790 #if ( CALCMOD==-1 && defined(AGE_CONST) )
791 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
792 #endif
793 #if ( CALCMOD==1 && defined(CTS_MELTING_FREEZING) )
794 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
795 #endif
796 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
797 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
798 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
799 #if (MARGIN==2)
800 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
801 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
802 #elif (MARGIN==3)
803 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
804 #endif
805 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
806 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
807 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
808 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
809 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
810 write(10, fmt=trim(fmt1)) ' '
811 write(10, fmt=trim(fmt2)) 'TEMP_PRESENT_PARA =', temp_present_para
812 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
813 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
814 #if ACCSURFACE==5
815 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
816 #endif
817 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
818 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
819 #if ( ABLSURFACE==1 || ABLSURFACE==2 )
820 write(10, fmt=trim(fmt2)) 'PDD_MODIFIER =', pdd_modifier
821 #if PDD_MODIFIER==2
822 write(10, fmt=trim(fmt2)) 'N_PDD_MOD =', n_pdd_mod
823 write(10, fmt=trim(fmt3)) 'LON_W_E_SEP =', lon_w_e_sep
824 do n=1, n_pdd_mod
825  if (n== 1) then
826  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_01 =', pdd_mod_lat_01
827  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_01 =', pdd_mod_fac_w_01
828  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_01 =', pdd_mod_fac_e_01
829  end if
830  if (n== 2) then
831  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_02 =', pdd_mod_lat_02
832  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_02 =', pdd_mod_fac_w_02
833  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_02 =', pdd_mod_fac_e_02
834  end if
835  if (n== 3) then
836  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_03 =', pdd_mod_lat_03
837  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_03 =', pdd_mod_fac_w_03
838  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_03 =', pdd_mod_fac_e_03
839  end if
840  if (n== 4) then
841  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_04 =', pdd_mod_lat_04
842  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_04 =', pdd_mod_fac_w_04
843  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_04 =', pdd_mod_fac_e_04
844  end if
845  if (n== 5) then
846  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_05 =', pdd_mod_lat_05
847  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_05 =', pdd_mod_fac_w_05
848  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_05 =', pdd_mod_fac_e_05
849  end if
850  if (n== 6) then
851  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_06 =', pdd_mod_lat_06
852  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_06 =', pdd_mod_fac_w_06
853  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_06 =', pdd_mod_fac_e_06
854  end if
855  if (n== 7) then
856  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_07 =', pdd_mod_lat_07
857  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_07 =', pdd_mod_fac_w_07
858  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_07 =', pdd_mod_fac_e_07
859  end if
860  if (n== 8) then
861  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_08 =', pdd_mod_lat_08
862  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_08 =', pdd_mod_fac_w_08
863  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_08 =', pdd_mod_fac_e_08
864  end if
865  if (n== 9) then
866  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_09 =', pdd_mod_lat_09
867  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_09 =', pdd_mod_fac_w_09
868  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_09 =', pdd_mod_fac_e_09
869  end if
870  if (n==10) then
871  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_10 =', pdd_mod_lat_10
872  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_10 =', pdd_mod_fac_w_10
873  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_10 =', pdd_mod_fac_e_10
874  end if
875 end do
876 #endif
877 #endif
878 write(10, fmt=trim(fmt1)) ' '
879 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
880 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
881 write(10, fmt=trim(fmt2)) 'ICE_STREAM =', ice_stream
882 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
883 write(10, fmt=trim(fmt1)) ' '
884 
885 #if defined(OUT_TIMES)
886 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
887 #endif
888 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
889 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
890 #if defined(OUTSER_V_A)
891 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
892 #endif
893 #if ( OUTPUT==1 || OUTPUT==2 )
894 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
895 #endif
896 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
897 write(10, fmt=trim(fmt1)) ' '
898 #if ( OUTPUT==2 || OUTPUT==3 )
899 write(10, fmt=trim(fmt2)) 'n_output =', n_output
900 do n=1, n_output
901  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
902 end do
903 write(10, fmt=trim(fmt1)) ' '
904 #endif
905 
906 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
907 
908 close(10, status='keep')
909 
910 !-------- Conversion of time quantities --------
911 
912 year_zero = year_zero*year_sec ! a --> s
913 time_init = time_init0*year_sec ! a --> s
914 time_end = time_end0*year_sec ! a --> s
915 dtime = dtime0*year_sec ! a --> s
916 dtime_temp = dtime_temp0*year_sec ! a --> s
917 #if (REBOUND==2)
918 dtime_wss = dtime_wss0*year_sec ! a --> s
919 #endif
920 dtime_ser = dtime_ser0*year_sec ! a --> s
921 #if ( OUTPUT==1 || OUTPUT==3 )
922 dtime_out = dtime_out0*year_sec ! a --> s
923 #endif
924 #if ( OUTPUT==2 || OUTPUT==3 )
925 do n=1, n_output
926  time_output(n) = time_output0(n)*year_sec ! a --> s
927 end do
928 #endif
929 
930 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
931  stop ' sico_init: dtime_temp must be a multiple of dtime!'
932 
933 #if (REBOUND==2)
934 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
935  stop ' sico_init: dtime_wss must be a multiple of dtime!'
936 #endif
937 
938 #if ZS_EVOL==2
939 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
940 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
941 #endif
942 
943 time = time_init
944 
945 !-------- Reading of present mean-annual precipitation rate --------
946 
947 #if (GRID==0 || GRID==1)
948 
949 open(21, iostat=ios, &
950  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_present_file, &
951  recl=8192, status='old')
952 
953 #elif (GRID==2)
954 
955 stop ' sico_init: GRID==2 not allowed for this application!'
956 
957 #endif
958 
959 if (ios /= 0) stop ' sico_init: Error when opening the precip file!'
960 
961 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
962 
963 do j=jmax, 0, -1
964  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
965 end do
966 
967 close(21, status='keep')
968 
969 precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
970  ! mm/a water equiv. -> m/s ice equiv.
971 
972 precip_factor = 1.0_dp
973 
974 #if (defined(PRECIP_FACTOR_FILE))
975 
976 if ( trim(adjustl(precip_factor_file)) /= 'none' ) then
977 
978  open(21, iostat=ios, &
979  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_factor_file, &
980  recl=8192, status='old')
981 
982  if (ios /= 0) stop ' sico_init: Error when opening the precip factor file!'
983 
984  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
985 
986  do j=jmax, 0, -1
987  read(21, fmt=*) (precip_factor(j,i), i=0,imax)
988  end do
989 
990  close(21, status='keep')
991 
992 end if
993 
994 #endif
995 
996 precip_ma_present = precip_ma_present * precip_factor
997 
998 !-------- Present monthly precipitation rates
999 ! (assumed to be equal to the mean annual precipitation rate) --------
1000 
1001 do n=1, 12 ! month counter
1002 do i=0, imax
1003 do j=0, jmax
1004  precip_present(j,i,n) = precip_ma_present(j,i)
1005 end do
1006 end do
1007 end do
1008 
1009 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
1010 
1011 #if ACCSURFACE==5
1012 
1013 #if (GRID==0 || GRID==1)
1014 
1015 open(21, iostat=ios, &
1016  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_anom_file, &
1017  recl=8192, status='old')
1018 
1019 #endif
1020 
1021 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
1022 
1023 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1024 
1025 do j=jmax, 0, -1
1026  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
1027 end do
1028 
1029 close(21, status='keep')
1030 
1031 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
1032 
1033 !-------- LGM monthly precipitation-rate anomalies (assumed to be
1034 ! equal to the mean annual precipitation-rate anomaly) --------
1035 
1036 do n=1, 12 ! month counter
1037 do i=0, imax
1038 do j=0, jmax
1039 
1040  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
1041  ! positive values ensured
1042 
1043 #if (PRECIP_ANOM_INTERPOL==1)
1044  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
1045 #elif (PRECIP_ANOM_INTERPOL==2)
1046  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1047 #else
1048  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1049 #endif
1050 
1051 end do
1052 end do
1053 end do
1054 
1055 #endif
1056 
1057 !-------- Mean accumulation --------
1058 
1059 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1060  ! mm/a water equiv. -> m/s ice equiv.
1061 
1062 !-------- Read present topography mask --------
1063 
1064 #if (GRID==0 || GRID==1)
1065 
1066 open(24, iostat=ios, &
1067  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_present_file, &
1068  recl=2048, status='old')
1069 
1070 #elif GRID==2
1071 
1072 stop ' sico_init: GRID==2 not allowed for this application!'
1073 
1074 #endif
1075 
1076 if (ios /= 0) stop ' sico_init: Error when opening the mask file!'
1077 
1078 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1079 
1080 do j=jmax, 0, -1
1081  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
1082 end do
1083 
1084 close(24, status='keep')
1085 
1086 !-------- Read rock/sediment mask --------
1087 
1088 #if ICE_STREAM==2
1089 
1090 open(24, iostat=ios, &
1091  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_sedi_file, &
1092  recl=2048, status='old')
1093 
1094 if (ios /= 0) stop ' sico_init: Error when opening the rock/sediment mask file!'
1095 
1096 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1097 
1098 do j=jmax, 0, -1
1099  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
1100 end do
1101 
1102 close(24, status='keep')
1103 
1104 #endif
1105 
1106 !-------- Present reference elevation (for precipitation data) --------
1107 
1108 #if (GRID==0 || GRID==1)
1109 
1110 open(21, iostat=ios, &
1111  file=inpath//'/'//trim(ch_domain_short)//'/'//zs_present_file, &
1112  recl=8192, status='old')
1113 
1114 #elif GRID==2
1115 
1116 stop ' sico_init: GRID==2 not allowed for this application!'
1117 
1118 #endif
1119 
1120 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
1121 
1122 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1123 
1124 do j=jmax, 0, -1
1125  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1126 end do
1127 
1128 close(21, status='keep')
1129 
1130 ! ------ Reset bathymetry data (sea floor elevation) to the
1131 ! sea surface
1132 
1133 do i=0, imax
1134 do j=0, jmax
1135  if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
1136 end do
1137 end do
1138 
1139 !-------- Reading of the prescribed target topography --------
1140 
1141 #if ZS_EVOL==2
1142 
1143 target_topo_dat_name = trim(target_topo_dat_name)
1144 
1145 #if NETCDF==1
1146 write(6, fmt='(a)') ' sico_init: Reading of target topography'
1147 write(6, fmt='(a)') ' only implemented for NetCDF files (NETCDF==2)!'
1148 stop
1149 #elif NETCDF==2
1150 call read_target_topo_nc(target_topo_dat_name)
1151 #else
1152 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1153 #endif
1154 
1155 #endif
1156 
1157 !-------- Reading of the maximum ice extent mask --------
1158 
1159 #if ZS_EVOL==3
1160 
1161 #if (GRID==0 || GRID==1)
1162 
1163 open(24, iostat=ios, &
1164  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_maxextent_file, &
1165  recl=2048, status='old')
1166 
1167 #elif GRID==2
1168 
1169 stop ' sico_init: GRID==2 not allowed for this application!'
1170 
1171 #endif
1172 
1173 if (ios /= 0) &
1174  stop ' sico_init: Error when opening the maximum ice extent mask file!'
1175 
1176 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1177 
1178 do j=jmax, 0, -1
1179  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1180 end do
1181 
1182 close(24, status='keep')
1183 
1184 #endif
1185 
1186 !-------- Reading of LGM mean-annual and mean-July surface-temperature
1187 ! anomalies --------
1188 
1189 #if TSURFACE==5
1190 
1191 #if (GRID==0 || GRID==1)
1192 
1193 open(21, iostat=ios, &
1194  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_ma_anom_file, &
1195  recl=8192, status='old')
1196 
1197 #endif
1198 
1199 if (ios /= 0) stop ' sico_init: Error when opening the temp_ma anomaly file!'
1200 
1201 #if (GRID==0 || GRID==1)
1202 
1203 open(22, iostat=ios, &
1204  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_mj_anom_file, &
1205  recl=8192, status='old')
1206 
1207 #endif
1208 
1209 if (ios /= 0) stop ' sico_init: Error when opening the temp_mj anomaly file!'
1210 
1211 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1212 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1213 
1214 do j=jmax, 0, -1
1215  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
1216  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
1217 end do
1218 
1219 close(21, status='keep')
1220 close(22, status='keep')
1221 
1222 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
1223 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
1224 
1225 #endif
1226 
1227 !-------- Read data for delta_ts --------
1228 
1229 #if TSURFACE==4
1230 
1231 open(21, iostat=ios, &
1232  file=inpath//'/general/'//grip_temp_file, &
1233  status='old')
1234 if (ios /= 0) &
1235  stop ' sico_init: Error when opening the data file for delta_ts!'
1236 
1237 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1238 
1239 if (ch_dummy /= '#') then
1240  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
1241  write(6, fmt=*) ' not defined in data file for delta_ts!'
1242  stop
1243 end if
1244 
1245 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
1246 
1247 allocate(griptemp(0:ndata_grip))
1248 
1249 do n=0, ndata_grip
1250  read(21, fmt=*) d_dummy, griptemp(n)
1251 end do
1252 
1253 close(21, status='keep')
1254 
1255 #endif
1256 
1257 !-------- Read data for the glacial index --------
1258 
1259 #if TSURFACE==5
1260 
1261 open(21, iostat=ios, &
1262  file=inpath//'/general/'//glac_ind_file, status='old')
1263 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
1264 
1265 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1266 
1267 if (ch_dummy /= '#') then
1268  write(6, fmt=*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
1269  write(6, fmt=*) ' not defined in glacial-index file!'
1270  stop
1271 end if
1272 
1273 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
1274 
1275 allocate(glacial_index(0:ndata_gi))
1276 
1277 do n=0, ndata_gi
1278  read(21, fmt=*) d_dummy, glacial_index(n)
1279 end do
1280 
1281 close(21, status='keep')
1282 
1283 #endif
1284 
1285 !-------- Reading of time-dependent surface-temperature
1286 ! and precipitation anomaly data --------
1287 
1288 #if ( (TSURFACE==6) && (ACCSURFACE==6) )
1289 
1290 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1291  '/'//trim(temp_precip_anom_file)
1292 
1293 call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
1294  ncid_temp_precip), thisroutine )
1295 
1296 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1297  thisroutine )
1298 
1299 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1300  len = ndata_temp_precip), thisroutine )
1301 
1302 ndata_temp_precip = ndata_temp_precip-1
1303 
1304 allocate(temp_precip_time(0:ndata_temp_precip))
1305 
1306 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1307 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1308  ! in a
1309 
1310 temp_precip_time_min = temp_precip_time(0) ! in a
1311 temp_precip_time_stp = temp_precip_time(1)-temp_precip_time(0)
1312  ! in a; constant time step assumed
1313 temp_precip_time_max = temp_precip_time(ndata_temp_precip) ! in a
1314 
1315 #endif
1316 
1317 !-------- Read data for z_sl --------
1318 
1319 #if (SEA_LEVEL==3)
1320 
1321 open(21, iostat=ios, &
1322  file=inpath//'/general/'//sea_level_file, &
1323  status='old')
1324 if (ios /= 0) &
1325  stop ' sico_init: Error when opening the data file for z_sl!'
1326 
1327 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1328 
1329 if (ch_dummy /= '#') then
1330  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1331  write(6, fmt=*) ' not defined in data file for z_sl!'
1332  stop
1333 end if
1334 
1335 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
1336 
1337 allocate(specmap_zsl(0:ndata_specmap))
1338 
1339 do n=0, ndata_specmap
1340  read(21, fmt=*) d_dummy, specmap_zsl(n)
1341 end do
1342 
1343 close(21, status='keep')
1344 
1345 #endif
1346 
1347 !-------- Determination of the geothermal heat flux --------
1348 
1349 #if Q_GEO_MOD==1
1350 
1351 ! ------ Constant value
1352 
1353 do i=0, imax
1354 do j=0, jmax
1355  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1356 end do
1357 end do
1358 
1359 #elif Q_GEO_MOD==2
1360 
1361 ! ------ Read data from file
1362 
1363 open(21, iostat=ios, &
1364  file=inpath//'/'//trim(ch_domain_short)//'/'//q_geo_file, &
1365  recl=8192, status='old')
1366 
1367 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
1368 
1369 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1370 
1371 do j=jmax, 0, -1
1372  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1373 end do
1374 
1375 close(21, status='keep')
1376 
1377 do i=0, imax
1378 do j=0, jmax
1379  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1380 end do
1381 end do
1382 
1383 #endif
1384 
1385 !-------- Reading of tabulated kei function--------
1386 
1387 #if (REBOUND==0 || REBOUND==1)
1388 
1389 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1390  ! dummy values
1391 #elif (REBOUND==2)
1392 
1393 call read_kei()
1394 
1395 #endif
1396 
1397 !-------- Determination of the time lag
1398 ! of the relaxing asthenosphere --------
1399 
1400 #if (REBOUND==1 || REBOUND==2)
1401 
1402 #if (TIME_LAG_MOD==1)
1403 
1404 time_lag_asth = time_lag*year_sec ! a -> s
1405 
1406 #elif (TIME_LAG_MOD==2)
1407 
1408 open(29, iostat=ios, &
1409  file=inpath//'/'//trim(ch_domain_short)//'/'//time_lag_file, &
1410  recl=8192, status='old')
1411 
1412 if (ios /= 0) stop ' sico_init: Error when opening the time-lag file!'
1413 
1414 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1415 
1416 do j=jmax, 0, -1
1417  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1418 end do
1419 
1420 close(29, status='keep')
1421 
1422 time_lag_asth = time_lag_asth*year_sec ! a -> s
1423 
1424 #endif
1425 
1426 #elif (REBOUND==0)
1427 
1428 time_lag_asth = 0.0_dp ! dummy values
1429 
1430 #endif
1431 
1432 !-------- Determination of the flexural rigidity of the lithosphere --------
1433 
1434 #if (REBOUND==2)
1435 
1436 #if (FLEX_RIG_MOD==1)
1437 
1438 flex_rig_lith = flex_rig
1439 
1440 #elif (FLEX_RIG_MOD==2)
1441 
1442 open(29, iostat=ios, &
1443  file=inpath//'/'//trim(ch_domain_short)//'/'//flex_rig_file, &
1444  recl=8192, status='old')
1445 
1446 if (ios /= 0) stop ' sico_init: Error when opening the flex-rig file!'
1447 
1448 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1449 
1450 do j=jmax, 0, -1
1451  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1452 end do
1453 
1454 close(29, status='keep')
1455 
1456 #endif
1457 
1458 #elif (REBOUND==0 || REBOUND==1)
1459 
1460 flex_rig_lith = 0.0_dp ! dummy values
1461 
1462 #endif
1463 
1464 !-------- Definition of initial values --------
1465 
1466 ! ------ Present topography
1467 
1468 #if (ANF_DAT==1)
1469 
1470 call topography1(dxi, deta)
1471 
1472 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1473 
1474 call boundary(time_init, dtime, dxi, deta, &
1475  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1476 
1477 q_bm = 0.0_dp
1478 q_tld = 0.0_dp
1479 q_b_tot = 0.0_dp
1480 
1481 p_b_w = 0.0_dp
1482 h_w = 0.0_dp
1483 
1484 call init_temp_water_age()
1485 
1486 call calc_enhance(time_init)
1487 
1488 ! ------ Ice-free, relaxed bedrock
1489 
1490 #elif (ANF_DAT==2)
1491 
1492 call topography2(dxi, deta)
1493 
1494 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1495 
1496 call boundary(time_init, dtime, dxi, deta, &
1497  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1498 
1499 q_bm = 0.0_dp
1500 q_tld = 0.0_dp
1501 q_b_tot = 0.0_dp
1502 
1503 p_b_w = 0.0_dp
1504 h_w = 0.0_dp
1505 
1506 call init_temp_water_age()
1507 
1508 call calc_enhance(time_init)
1509 
1510 ! ------ Read initial state from output of previous simulation
1511 
1512 #elif (ANF_DAT==3)
1513 
1514 #if (NETCDF==1)
1515 call topography3(dxi, deta, z_sl, anfdatname)
1516 #elif (NETCDF==2)
1517 call topography3_nc(dxi, deta, z_sl, anfdatname)
1518 #else
1519 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1520 #endif
1521 
1522 call boundary(time_init, dtime, dxi, deta, &
1523  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1524 
1525 q_b_tot = q_bm + q_tld
1526 
1527 #endif
1528 
1529 !-------- Grounding line and calving front flags --------
1530 
1531 flag_grounding_line_1 = .false.
1532 flag_grounding_line_2 = .false.
1533 
1534 flag_calving_front_1 = .false.
1535 flag_calving_front_2 = .false.
1536 
1537 #if (MARGIN==1 || MARGIN==2)
1538 
1539 continue
1540 
1541 #elif (MARGIN==3)
1542 
1543 do i=1, imax-1
1544 do j=1, jmax-1
1545 
1546  if ( (maske(j,i)==0_i2b) & ! grounded ice
1547  .and. &
1548  ( (maske(j,i+1)==3_i2b) & ! with
1549  .or.(maske(j,i-1)==3_i2b) & ! one
1550  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1551  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1552  ) &
1553  flag_grounding_line_1(j,i) = .true.
1554 
1555  if ( (maske(j,i)==3_i2b) & ! floating ice
1556  .and. &
1557  ( (maske(j,i+1)==0_i2b) & ! with
1558  .or.(maske(j,i-1)==0_i2b) & ! one
1559  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1560  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1561  ) &
1562  flag_grounding_line_2(j,i) = .true.
1563 
1564  if ( (maske(j,i)==3_i2b) & ! floating ice
1565  .and. &
1566  ( (maske(j,i+1)==2_i2b) & ! with
1567  .or.(maske(j,i-1)==2_i2b) & ! one
1568  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1569  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1570  ) &
1571  flag_calving_front_1(j,i) = .true.
1572 
1573  if ( (maske(j,i)==2_i2b) & ! ocean
1574  .and. &
1575  ( (maske(j,i+1)==3_i2b) & ! with
1576  .or.(maske(j,i-1)==3_i2b) & ! one
1577  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1578  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1579  ) &
1580  flag_calving_front_2(j,i) = .true.
1581 
1582 end do
1583 end do
1584 
1585 do i=1, imax-1
1586 
1587  j=0
1588  if ( (maske(j,i)==2_i2b) & ! ocean
1589  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1590  ) & ! floating ice point
1591  flag_calving_front_2(j,i) = .true.
1592 
1593  j=jmax
1594  if ( (maske(j,i)==2_i2b) & ! ocean
1595  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1596  ) & ! floating ice point
1597  flag_calving_front_2(j,i) = .true.
1598 
1599 end do
1600 
1601 do j=1, jmax-1
1602 
1603  i=0
1604  if ( (maske(j,i)==2_i2b) & ! ocean
1605  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1606  ) & ! floating ice point
1607  flag_calving_front_2(j,i) = .true.
1608 
1609  i=imax
1610  if ( (maske(j,i)==2_i2b) & ! ocean
1611  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1612  ) & ! floating ice point
1613  flag_calving_front_2(j,i) = .true.
1614 
1615 end do
1616 
1617 #else
1618 stop ' sico_init: MARGIN must be either 1, 2 or 3!'
1619 #endif
1620 
1621 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1622 
1623 #if (GRID==0 || GRID==1)
1624 
1625 do ir=-imax, imax
1626 do jr=-jmax, jmax
1627  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1628  ! distortion due to stereographic projection not accounted for
1629 end do
1630 end do
1631 
1632 #endif
1633 
1634 !-------- Initial velocities --------
1635 
1636 call calc_temp_melt()
1637 
1638 #if (DYNAMICS==1 || DYNAMICS==2)
1639 
1640 call calc_vxy_b_sia(time, z_sl)
1641 call calc_vxy_sia(dzeta_c, dzeta_t)
1642 
1643 #if (MARGIN==3 || DYNAMICS==2)
1644 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1645 #endif
1646 
1647 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1648 
1649 #if (MARGIN==3)
1650 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1651 #endif
1652 
1653 #elif (DYNAMICS==0)
1654 
1655 call calc_vxy_static()
1656 call calc_vz_static()
1657 
1658 #else
1659  stop ' sico_init: DYNAMICS must be either 0, 1 or 2!'
1660 #endif
1661 
1662 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1663 
1664 !-------- Initial enthalpies --------
1665 
1666 #if (CALCMOD==0 || CALCMOD==-1)
1667 
1668 do i=0, imax
1669 do j=0, jmax
1670 
1671  do kc=0, kcmax
1672  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1673  end do
1674 
1675  do kt=0, ktmax
1676  enth_t(kt,j,i) = enth_c(0,j,i)
1677  end do
1678 
1679 end do
1680 end do
1681 
1682 #elif (CALCMOD==1)
1683 
1684 do i=0, imax
1685 do j=0, jmax
1686 
1687  do kc=0, kcmax
1688  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1689  end do
1690 
1691  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1692  do kt=0, ktmax
1693  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1694  end do
1695  else
1696  do kt=0, ktmax
1697  enth_t(kt,j,i) = enth_c(0,j,i)
1698  end do
1699  end if
1700 
1701 end do
1702 end do
1703 
1704 #elif (CALCMOD==2 || CALCMOD==3)
1705 
1706 do i=0, imax
1707 do j=0, jmax
1708 
1709  do kc=0, kcmax
1710  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1711  end do
1712 
1713  do kt=0, ktmax
1714  enth_t(kt,j,i) = enth_c(0,j,i)
1715  end do
1716 
1717 end do
1718 end do
1719 
1720 #else
1721 
1722 stop ' sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1723 
1724 #endif
1725 
1726 !-------- Initialize time-series files --------
1727 
1728 ! ------ Time-series file for the ice sheet on the whole
1729 
1730 open(12, iostat=ios, &
1731  file=outpath//'/'//trim(runname)//'.ser', &
1732  status='new')
1733 if (ios /= 0) &
1734  stop ' sico_init: Error when opening the ser file!'
1735 
1736 if (forcing_flag == 1) then
1737 
1738  write(12,1102)
1739  write(12,1103)
1740 
1741 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1742 
1743  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1744  ' H_max(m) zs_max(m) V_g(m^3)', &
1745  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1746  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1747  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1748  1103 format('----------------------------------------------------', &
1749  '---------------------------------------')
1750 
1751 #elif OUTSER_V_A==2
1752 
1753  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1754  ' V(m^3) V_g(m^3) V_f(m^3)', &
1755  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1756  ' H_max(m) zs_max(m)', &
1757  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1758  ' A_t(m^2) V_bm(m^3/a)', &
1759  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1760  1103 format('----------------------------------------------------', &
1761  '---------------------------------------')
1762 
1763 #else
1764  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1765 #endif
1766 
1767 else if (forcing_flag == 2) then
1768 
1769  write(12,1112)
1770  write(12,1113)
1771 
1772 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1773 
1774  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1775  ' H_max(m) zs_max(m) V_g(m^3)', &
1776  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1777  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1778  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1779  1113 format('----------------------------------------------------', &
1780  '---------------------------------------')
1781 
1782 #elif OUTSER_V_A==2
1783 
1784  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1785  ' V(m^3) V_g(m^3) V_f(m^3)', &
1786  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1787  ' H_max(m) zs_max(m)', &
1788  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1789  ' A_t(m^2) V_bm(m^3/a)', &
1790  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1791  1113 format('----------------------------------------------------', &
1792  '---------------------------------------')
1793 
1794 #else
1795  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1796 #endif
1797 
1798 else if (forcing_flag == 3) then
1799 
1800  write(12,1122)
1801  write(12,1123)
1802 
1803 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1804 
1805  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1806  ' H_max(m) zs_max(m) V_g(m^3)', &
1807  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1808  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1809  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1810  1123 format('----------------------------------------------------', &
1811  '---------------------------------------')
1812 
1813 #elif OUTSER_V_A==2
1814 
1815  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1816  ' V(m^3) V_g(m^3) V_f(m^3)', &
1817  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1818  ' H_max(m) zs_max(m)', &
1819  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1820  ' A_t(m^2) V_bm(m^3/a)', &
1821  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1822  1123 format('----------------------------------------------------', &
1823  '---------------------------------------')
1824 
1825 #else
1826  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1827 #endif
1828 
1829 end if
1830 
1831 ! ------ Time-series file for the ice-thickness field
1832 
1833 #if OUTSER==2
1834 
1835 open(13, iostat=ios, &
1836  file=outpath//'/'//trim(runname)//'.thk', &
1837  status='new')
1838 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1839 
1840 if (forcing_flag == 1) then
1841 
1842  write(13,1104)
1843  write(13,1105)
1844 
1845  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1846  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1847  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1848  1105 format('----------------------------------------------------')
1849 
1850 else if (forcing_flag == 2) then
1851 
1852  write(13,1114)
1853  write(13,1115)
1854 
1855  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1856  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1857  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1858  1115 format('----------------------------------------------------')
1859 
1860 else if (forcing_flag == 3) then
1861 
1862  write(13,1124)
1863  write(13,1125)
1864 
1865  1124 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1866  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1867  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1868  1125 format('----------------------------------------------------')
1869 
1870 end if
1871 
1872 #endif
1873 
1874 ! ------ Time-series file for deep boreholes
1875 
1876 #if OUTSER==3
1877 
1878 n_core = 6 ! GRIP, GISP2, Dye3, Camp Century (CC), NorthGRIP (NGRIP), NEEM
1879 
1880 allocate(lambda_core(n_core), phi_core(n_core), &
1881  x_core(n_core), y_core(n_core))
1882 
1883 phi_core(1) = 72.58722_dp *pi_180 ! Geographical position of GRIP,
1884 lambda_core(1) = -37.64222_dp *pi_180 ! conversion deg -> rad
1885 call num_coord(lambda_core(1), phi_core(1), x_core(1), y_core(1))
1886 
1887 phi_core(2) = 72.58833_dp *pi_180 ! Geographical position of GISP2
1888 lambda_core(2) = -38.45750_dp *pi_180 ! conversion deg -> rad
1889 call num_coord(lambda_core(2), phi_core(2), x_core(2), y_core(2))
1890 
1891 phi_core(3) = 65.15139_dp *pi_180 ! Geographical position of Dye3,
1892 lambda_core(3) = -43.81722_dp *pi_180 ! conversion deg -> rad
1893 call num_coord(lambda_core(3), phi_core(3), x_core(3), y_core(3))
1894 
1895 phi_core(4) = 77.17970_dp *pi_180 ! Geographical position of CC,
1896 lambda_core(4) = -61.10975_dp *pi_180 ! conversion deg -> rad
1897 call num_coord(lambda_core(4), phi_core(4), x_core(4), y_core(4))
1898 
1899 phi_core(5) = 75.09694_dp *pi_180 ! Geographical position of NGRIP,
1900 lambda_core(5) = -42.31956_dp *pi_180 ! conversion deg -> rad
1901 call num_coord(lambda_core(5), phi_core(5), x_core(5), y_core(5))
1902 
1903 phi_core(6) = 77.5_dp *pi_180 ! Geographical position of NEEM,
1904 lambda_core(6) = -50.9_dp *pi_180 ! conversion deg -> rad
1905 call num_coord(lambda_core(6), phi_core(6), x_core(6), y_core(6))
1906 
1907 open(14, iostat=ios, &
1908  file=outpath//'/'//trim(runname)//'.core', &
1909  status='new')
1910 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
1911 
1912 if (forcing_flag == 1) then
1913 
1914  write(14,1106)
1915  write(14,1107)
1916 
1917  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1918  ' H_GR(m) H_G2(m) H_D3(m)', &
1919  ' H_CC(m) H_NG(m) H_NE(m)',/, &
1920  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
1921  ' v_CC(m/a) v_NG(m/a) v_NE(m/a)',/, &
1922  ' T_GR(C) T_G2(C) T_D3(C)', &
1923  ' T_CC(C) T_NG(C) T_NE(C)',/, &
1924  ' Rb_GR(W/m2) Rb_G2(W/m2) Rb_D3(W/m2)', &
1925  ' Rb_CC(W/m2) Rb_NG(W/m2) Rb_NE(W/m2)')
1926  1107 format('----------------------------------------------------', &
1927  '---------------------------------------')
1928 
1929 else if (forcing_flag == 2) then
1930 
1931  write(14,1116)
1932  write(14,1117)
1933 
1934  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
1935  ' H_GR(m) H_G2(m) H_D3(m)', &
1936  ' H_CC(m) H_NG(m) H_NE(m)',/, &
1937  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
1938  ' v_CC(m/a) v_NG(m/a) v_NE(m/a)',/, &
1939  ' T_GR(C) T_G2(C) T_D3(C)', &
1940  ' T_CC(C) T_NG(C) T_NE(C)',/, &
1941  ' Rb_GR(W/m2) Rb_G2(W/m2) Rb_D3(W/m2)', &
1942  ' Rb_CC(W/m2) Rb_NG(W/m2) Rb_NE(W/m2)')
1943  1117 format('----------------------------------------------------', &
1944  '---------------------------------------')
1945 
1946 else if (forcing_flag == 3) then
1947 
1948  write(14,1126)
1949  write(14,1127)
1950 
1951  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1952  ' H_GR(m) H_G2(m) H_D3(m)', &
1953  ' H_CC(m) H_NG(m) H_NE(m)',/, &
1954  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
1955  ' v_CC(m/a) v_NG(m/a) v_NE(m/a)',/, &
1956  ' T_GR(C) T_G2(C) T_D3(C)', &
1957  ' T_CC(C) T_NG(C) T_NE(C)',/, &
1958  ' Rb_GR(W/m2) Rb_G2(W/m2) Rb_D3(W/m2)', &
1959  ' Rb_CC(W/m2) Rb_NG(W/m2) Rb_NE(W/m2)')
1960  1127 format('----------------------------------------------------', &
1961  '---------------------------------------')
1962 
1963 end if
1964 
1965 #endif
1966 
1967 !-------- Output of the initial state --------
1968 
1969 #if OUTPUT==1
1970 
1971 #if ERGDAT==0
1972  flag_3d_output = .false.
1973 #elif ERGDAT==1
1974  flag_3d_output = .true.
1975 #else
1976  stop ' sico_init: ERGDAT must be either 0 or 1!'
1977 #endif
1978 
1979 #if NETCDF==1
1980  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1981  flag_3d_output, ndat2d, ndat3d)
1982 #elif NETCDF==2
1983  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1984  flag_3d_output, ndat2d, ndat3d)
1985 #else
1986  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1987 #endif
1988 
1989 #elif OUTPUT==2
1990 
1991 if (time_output(1) <= time_init+eps) then
1992 
1993 #if ERGDAT==0
1994  flag_3d_output = .false.
1995 #elif ERGDAT==1
1996  flag_3d_output = .true.
1997 #else
1998  stop ' sico_init: ERGDAT must be either 0 or 1!'
1999 #endif
2000 
2001 #if NETCDF==1
2002  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2003  flag_3d_output, ndat2d, ndat3d)
2004 #elif NETCDF==2
2005  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2006  flag_3d_output, ndat2d, ndat3d)
2007 #else
2008  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2009 #endif
2010 
2011 end if
2012 
2013 #elif OUTPUT==3
2014 
2015  flag_3d_output = .false.
2016 
2017 #if NETCDF==1
2018  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2019  flag_3d_output, ndat2d, ndat3d)
2020 #elif NETCDF==2
2021  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2022  flag_3d_output, ndat2d, ndat3d)
2023 #else
2024  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2025 #endif
2026 
2027 if (time_output(1) <= time_init+eps) then
2028 
2029  flag_3d_output = .true.
2030 
2031 #if NETCDF==1
2032  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2033  flag_3d_output, ndat2d, ndat3d)
2034 #elif NETCDF==2
2035  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2036  flag_3d_output, ndat2d, ndat3d)
2037 #else
2038  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2039 #endif
2040 
2041 end if
2042 
2043 #else
2044  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
2045 #endif
2046 
2047 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2048 #if OUTSER==2
2049 call output3(time_init, delta_ts, glac_index, z_sl)
2050 #endif
2051 #if OUTSER==3
2052 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2053 #endif
2054 
2055 end subroutine sico_init
2056 !
nc_check
NetCDF error capturing.
Definition: nc_check.F90:35
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
topography1
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography1.F90:39
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
num_coord
subroutine num_coord(lambda_val, phi_val, x_val, y_val)
Computation of position (x,y) in the numerical domain for longitude lambda and latitude phi (in rad)...
Definition: num_coord.F90:37
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
nc_check::check
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check.F90:45
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
output4
subroutine output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format.
Definition: output4.F90:35
read_target_topo_nc
subroutine read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_target_topo_nc.F90:35
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