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 if (abs(dx-40.0_dp) < eps) then
201 
202  if ((imax /= 150).or.(jmax /= 70)) then
203  stop ' sico_init: IMAX and/or JMAX wrong!'
204  end if
205 
206 else if (abs(dx-20.0_dp) < eps) then
207 
208  if ((imax /= 300).or.(jmax /= 140)) then
209  stop ' sico_init: IMAX and/or JMAX wrong!'
210  end if
211 
212 else if (abs(dx-10.0_dp) < eps) then
213 
214  if ((imax /= 600).or.(jmax /= 280)) then
215  stop ' sico_init: IMAX and/or JMAX wrong!'
216  end if
217 
218 else
219  stop ' sico_init: DX wrong!'
220 end if
221 
222 #elif GRID==2
223 
224  stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
225 
226 #endif
227 
228 !-------- Compatibility check of the thermodynamics mode
229 ! (cold vs. polythermal vs. enthalpy method)
230 ! and the number of grid points in the lower (kt) ice domain --------
231 
232 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
233 
234 if (ktmax /= 1) then
235  write(6, fmt='(a)') ' sico_init: For options CALCMOD==0, 2, 3 or -1,'
236  write(6, fmt='(a)') ' the separate kt domain is redundant.'
237  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 1.'
238  write(6, fmt='(a)') ' '
239 end if
240 
241 #endif
242 
243 !-------- Compatibility check of surface-temperature and precipitation
244 ! determination by interpolation between present and LGM values
245 ! with a glacial index --------
246 
247 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
248 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
249 #endif
250 
251 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
252 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
253 #endif
254 
255 !-------- Compatibility check of discretization schemes for the horizontal and
256 ! vertical advection terms in the temperature and age equations --------
257 
258 #if (ADV_HOR==1)
259 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
260 #endif
261 
262 !-------- Check whether for the shallow shelf
263 ! or shelfy stream approximation
264 ! the chosen grid is Cartesian coordinates
265 ! without distortion correction (GRID==0) --------
266 
267 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
268 #if (GRID != 0)
269 write(6, fmt='(a)') ' sico_init: Distortion correction not implemented'
270 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
271 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
272 write(6, fmt='(a)') ' -> GRID==0 required!'
273 stop
274 #endif
275 #endif
276 
277 !-------- Setting of forcing flag --------
278 
279 #if (TSURFACE <= 4)
280 
281 forcing_flag = 1 ! forcing by delta_ts
282 
283 #elif (TSURFACE == 5)
284 
285 forcing_flag = 2 ! forcing by glac_index
286 
287 #endif
288 
289 !-------- Initialization of numerical time steps --------
290 
291 dtime0 = dtime0
292 dtime_temp0 = dtime_temp0
293 #if (REBOUND==2)
294 dtime_wss0 = dtime_wss0
295 #endif
296 
297 !-------- Further initializations --------
298 
299 dzeta_c = 1.0_dp/real(kcmax,dp)
300 dzeta_t = 1.0_dp/real(ktmax,dp)
301 dzeta_r = 1.0_dp/real(krmax,dp)
302 
303 ndat2d = 1
304 ndat3d = 1
305 
306 !-------- General abbreviations --------
307 
308 ! ------ kc domain
309 
310 if (deform >= eps) then
311 
312  flag_aa_nonzero = .true. ! non-equidistant grid
313 
314  aa = deform
315  ea = exp(aa)
316 
317  kc=0
318  zeta_c(kc) = 0.0_dp
319  eaz_c(kc) = 1.0_dp
320  eaz_c_quotient(kc) = 0.0_dp
321 
322  do kc=1, kcmax-1
323  zeta_c(kc) = kc*dzeta_c
324  eaz_c(kc) = exp(aa*zeta_c(kc))
325  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
326  end do
327 
328  kc=kcmax
329  zeta_c(kc) = 1.0_dp
330  eaz_c(kc) = exp(aa)
331  eaz_c_quotient(kc) = 1.0_dp
332 
333 else
334 
335  flag_aa_nonzero = .false. ! equidistant grid
336 
337  aa = 0.0_dp
338  ea = 1.0_dp
339 
340  kc=0
341  zeta_c(kc) = 0.0_dp
342  eaz_c(kc) = 1.0_dp
343  eaz_c_quotient(kc) = 0.0_dp
344 
345  do kc=1, kcmax-1
346  zeta_c(kc) = kc*dzeta_c
347  eaz_c(kc) = 1.0_dp
348  eaz_c_quotient(kc) = zeta_c(kc)
349  end do
350 
351  kc=kcmax
352  zeta_c(kc) = 1.0_dp
353  eaz_c(kc) = 1.0_dp
354  eaz_c_quotient(kc) = 1.0_dp
355 
356 end if
357 
358 ! ------ kt domain
359 
360 kt=0
361 zeta_t(kt) = 0.0_dp
362 
363 do kt=1, ktmax-1
364  zeta_t(kt) = kt*dzeta_t
365 end do
366 
367 kt=ktmax
368 zeta_t(kt) = 1.0_dp
369 
370 ! ------ kr domain
371 
372 kr=0
373 zeta_r(kr) = 0.0_dp
374 
375 do kr=1, krmax-1
376  zeta_r(kr) = kr*dzeta_r
377 end do
378 
379 kr=krmax
380 zeta_r(kr) = 1.0_dp
381 
382 !-------- Construction of a vector (with index n) from a 2-d array
383 ! (with indices i, j) by diagonal numbering --------
384 
385 n=1
386 
387 do m=0, imax+jmax
388  do i=m, 0, -1
389  j = m-i
390  if ((i <= imax).and.(j <= jmax)) then
391  ii(n) = i
392  jj(n) = j
393  nn(j,i) = n
394  n=n+1
395  end if
396  end do
397 end do
398 
399 !-------- Specification of current simulation --------
400 
401 runname = runname
402 anfdatname = anfdatname
403 
404 #if defined(YEAR_ZERO)
405 year_zero = year_zero
406 #else
407 year_zero = 2000.0_dp ! default value 2000 CE
408 #endif
409 
410 time_init0 = time_init0
411 time_end0 = time_end0
412 dtime_ser0 = dtime_ser0
413 
414 #if ( OUTPUT==1 || OUTPUT==3 )
415 dtime_out0 = dtime_out0
416 #endif
417 
418 #if ( OUTPUT==2 || OUTPUT==3 )
419 n_output = n_output
420 time_output0( 1) = time_out0_01
421 time_output0( 2) = time_out0_02
422 time_output0( 3) = time_out0_03
423 time_output0( 4) = time_out0_04
424 time_output0( 5) = time_out0_05
425 time_output0( 6) = time_out0_06
426 time_output0( 7) = time_out0_07
427 time_output0( 8) = time_out0_08
428 time_output0( 9) = time_out0_09
429 time_output0(10) = time_out0_10
430 time_output0(11) = time_out0_11
431 time_output0(12) = time_out0_12
432 time_output0(13) = time_out0_13
433 time_output0(14) = time_out0_14
434 time_output0(15) = time_out0_15
435 time_output0(16) = time_out0_16
436 time_output0(17) = time_out0_17
437 time_output0(18) = time_out0_18
438 time_output0(19) = time_out0_19
439 time_output0(20) = time_out0_20
440 #endif
441 
442 !-------- Write log file --------
443 
444 shell_command = 'if [ ! -d'
445 shell_command = trim(shell_command)//' '//outpath
446 shell_command = trim(shell_command)//' '//'] ; then mkdir'
447 shell_command = trim(shell_command)//' '//outpath
448 shell_command = trim(shell_command)//' '//'; fi'
449 call system(trim(shell_command))
450  ! Check whether directory OUTPATH exists. If not, it is created.
451 
452 open(10, iostat=ios, &
453  file=outpath//'/'//trim(runname)//'.log', &
454  status='new')
455 if (ios /= 0) &
456  stop 'Error when opening the log file!'
457 
458 write(10, fmt=trim(fmt1)) 'Computational domain:'
459 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
460 write(10, fmt=trim(fmt1)) ' '
461 
462 write(10, fmt=trim(fmt2)) 'imax =', imax
463 write(10, fmt=trim(fmt2)) 'jmax =', jmax
464 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
465 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
466 write(10, fmt=trim(fmt2)) 'krmax =', krmax
467 write(10, fmt=trim(fmt1)) ' '
468 
469 write(10, fmt=trim(fmt3)) 'a =', aa
470 write(10, fmt=trim(fmt1)) ' '
471 
472 #if (GRID==0 || GRID==1)
473 write(10, fmt=trim(fmt3)) 'x0 =', x0
474 write(10, fmt=trim(fmt3)) 'y0 =', y0
475 write(10, fmt=trim(fmt3)) 'dx =', dx
476 #elif GRID==2
477 stop ' sico_init: GRID==2 not allowed for this application!'
478 #endif
479 write(10, fmt=trim(fmt1)) ' '
480 
481 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
482 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
483 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
484 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
485 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
486 #if (REBOUND==2)
487 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
488 #endif
489 #if ( OUTPUT==1 || OUTPUT==3 )
490 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
491 #endif
492 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
493 write(10, fmt=trim(fmt1)) ' '
494 
495 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
496 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
497 #if (ANF_DAT==1)
498 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
499 #endif
500 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
501 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
502 #if (ANF_DAT==1 && defined(TEMP_INIT))
503 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
504 #endif
505 #if (ANF_DAT==3)
506 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
507 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
508 #endif
509 write(10, fmt=trim(fmt1)) ' '
510 
511 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
512 write(10, fmt=trim(fmt1)) ' '
513 
514 #if (CALCZS==3 || CALCZS==4)
515 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
516 #if CALCZS==3
517 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
518 #endif
519 write(10, fmt=trim(fmt1)) ' '
520 #endif
521 
522 write(10, fmt=trim(fmt1)) 'temp_mm_present file = '//temp_mm_present_file
523 #if TSURFACE==1
524 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
525 #elif TSURFACE==3
526 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
527 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
528 #elif TSURFACE==4
529 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
530 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
531 #elif TSURFACE==5
532 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
533 write(10, fmt=trim(fmt1)) 'temp_mm_anom file = '//temp_mm_anom_file
534 write(10, fmt=trim(fmt3)) 'temp_mm_anom fact = ', temp_mm_anom_fact
535 #endif
536 
537 write(10, fmt=trim(fmt1)) 'precip_mm_present file = '//precip_mm_present_file
538 #if ACCSURFACE==1
539 write(10, fmt=trim(fmt3)) 'accfact =', accfact
540 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
541 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
542 #elif ( ACCSURFACE==5 )
543 write(10, fmt=trim(fmt1)) 'precip_mm_anom file = '//precip_mm_anom_file
544 write(10, fmt=trim(fmt3)) 'precip_mm_anom fact = ', precip_mm_anom_fact
545 #endif
546 #if (ACCSURFACE <= 3)
547 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
548 #if (ELEV_DESERT == 1)
549 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
550 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
551 #endif
552 #endif
553 
554 #if ABLSURFACE==3
555 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
556 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
557 #endif
558 
559 #if (SEA_LEVEL==1)
560 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
561 #elif (SEA_LEVEL==3)
562 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
563 #endif
564 write(10, fmt=trim(fmt1)) ' '
565 
566 #if (MARGIN==2)
567 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
568 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
569 write(10, fmt=trim(fmt1)) ' '
570 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
571  || marine_ice_calving==6 || marine_ice_calving==7 )
572 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
573 write(10, fmt=trim(fmt1)) ' '
574 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
575 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
576 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
577 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
578 write(10, fmt=trim(fmt1)) ' '
579 #endif
580 #elif (MARGIN==3)
581 #if ICE_SHELF_CALVING==2
582 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
583 write(10, fmt=trim(fmt1)) ' '
584 #endif
585 #endif
586 
587 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
588 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
589 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
590 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
591 #if SLIDE_LAW==2
592 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
593 #endif
594 write(10, fmt=trim(fmt1)) ' '
595 
596 #if Q_GEO_MOD==1
597 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
598 #elif Q_GEO_MOD==2
599 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
600 #endif
601 write(10, fmt=trim(fmt1)) ' '
602 
603 #if (defined(MARINE_ICE_BASAL_MELTING))
604 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
605 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
606 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
607 #endif
608 write(10, fmt=trim(fmt1)) ' '
609 #endif
610 
611 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
612 #if (REBOUND==1)
613 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
614 #endif
615 #if (REBOUND==1 || REBOUND==2)
616 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
617 #if (TIME_LAG_MOD==1)
618 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
619 #elif (TIME_LAG_MOD==2)
620 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
621 #else
622 stop ' sico_init: TIME_LAG_MOD must be either 1 or 2!'
623 #endif
624 #endif
625 #if (REBOUND==2)
626 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
627 #if (FLEX_RIG_MOD==1)
628 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
629 #elif (FLEX_RIG_MOD==2)
630 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
631 #else
632 stop ' sico_init: FLEX_RIG_MOD must be either 1 or 2!'
633 #endif
634 #endif
635 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
636 write(10, fmt=trim(fmt1)) ' '
637 
638 #if FLOW_LAW==2
639 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
640 write(10, fmt=trim(fmt1)) ' '
641 #endif
642 #if FIN_VISC==2
643 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
644 write(10, fmt=trim(fmt1)) ' '
645 #endif
646 
647 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
648 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
649 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
650 #endif
651 #if (ENHMOD==2 || ENHMOD==3)
652 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
653 #endif
654 #if (ENHMOD==2)
655 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
656 #endif
657 #if (ENHMOD==3)
658 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
659 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
660 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
661 #endif
662 #if (ENHMOD==4 || ENHMOD==5)
663 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
664 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
665 #endif
666 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
667 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
668 #endif
669 write(10, fmt=trim(fmt1)) ' '
670 
671 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
672 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
673 write(10, fmt=trim(fmt3)) 'H_min =', h_min
674 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
675 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
676 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
677 #if defined(QBM_MIN)
678 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
679 #elif defined(QB_MIN) /* obsolete */
680 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
681 #endif
682 #if defined(QBM_MAX)
683 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
684 #elif defined(QB_MAX) /* obsolete */
685 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
686 #endif
687 write(10, fmt=trim(fmt3)) 'age_min =', age_min
688 write(10, fmt=trim(fmt3)) 'age_max =', age_max
689 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
690 #if ADV_VERT==1
691 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
692 #endif
693 write(10, fmt=trim(fmt1)) ' '
694 
695 write(10, fmt=trim(fmt2)) 'GRID =', grid
696 write(10, fmt=trim(fmt2)) 'DYNAMICS =', dynamics
697 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
698 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
699 #endif
700 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
701 #if ( CALCMOD==-1 && defined(TEMP_CONST) )
702 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
703 #endif
704 #if ( CALCMOD==-1 && defined(AGE_CONST) )
705 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
706 #endif
707 #if ( CALCMOD==1 && defined(CTS_MELTING_FREEZING) )
708 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
709 #endif
710 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
711 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
712 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
713 #if (MARGIN==2)
714 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
715 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
716 #elif (MARGIN==3)
717 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
718 #endif
719 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
720 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
721 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
722 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
723 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
724 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
725 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
726 #if ACCSURFACE==5
727 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
728 #endif
729 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
730 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
731 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
732 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
733 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
734 write(10, fmt=trim(fmt1)) ' '
735 
736 #if defined(OUT_TIMES)
737 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
738 #endif
739 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
740 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
741 #if defined(OUTSER_V_A)
742 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
743 #endif
744 #if ( OUTPUT==1 || OUTPUT==2 )
745 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
746 #endif
747 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
748 write(10, fmt=trim(fmt1)) ' '
749 #if ( OUTPUT==2 || OUTPUT==3 )
750 write(10, fmt=trim(fmt2)) 'n_output =', n_output
751 do n=1, n_output
752  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
753 end do
754 write(10, fmt=trim(fmt1)) ' '
755 #endif
756 
757 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
758 
759 close(10, status='keep')
760 
761 !-------- Conversion of time quantities --------
762 
763 year_zero = year_zero*year_sec ! a --> s
764 time_init = time_init0*year_sec ! a --> s
765 time_end = time_end0*year_sec ! a --> s
766 dtime = dtime0*year_sec ! a --> s
767 dtime_temp = dtime_temp0*year_sec ! a --> s
768 #if (REBOUND==2)
769 dtime_wss = dtime_wss0*year_sec ! a --> s
770 #endif
771 dtime_ser = dtime_ser0*year_sec ! a --> s
772 #if ( OUTPUT==1 || OUTPUT==3 )
773 dtime_out = dtime_out0*year_sec ! a --> s
774 #endif
775 #if ( OUTPUT==2 || OUTPUT==3 )
776 do n=1, n_output
777  time_output(n) = time_output0(n)*year_sec ! a --> s
778 end do
779 #endif
780 
781 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
782  stop ' sico_init: dtime_temp must be a multiple of dtime!'
783 
784 #if (REBOUND==2)
785 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
786  stop ' sico_init: dtime_wss must be a multiple of dtime!'
787 #endif
788 
789 time = time_init
790 
791 !-------- Reading of measurements for present monthly-mean precipitation --------
792 
793 #if (GRID==0 || GRID==1)
794 
795 open(21, iostat=ios, &
796  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_mm_present_file, &
797  recl=16384, status='old')
798 
799 #elif GRID==2
800 
801 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
802 
803 #endif
804 
805 if (ios /= 0) stop ' sico_init: Error when opening the precipitation file!'
806 
807 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
808 
809 do n=1, 12 ! month counter
810  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
811  do j=jmax, 0, -1
812  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
813  end do
814 end do
815 
816 close(21, status='keep')
817 
818 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
819 
820 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
821  ! mm/a water equiv. -> m/s ice equiv.
822 
823 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
824 
825 #if ACCSURFACE==5
826 
827 #if (GRID==0 || GRID==1)
828 
829 open(21, iostat=ios, &
830  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_mm_anom_file, &
831  recl=16384, status='old')
832 
833 #endif
834 
835 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
836 
837 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
838 
839 do n=1, 12 ! month counter
840  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
841  do j=jmax, 0, -1
842  read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
843  end do
844 end do
845 
846 close(21, status='keep')
847 
848 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
849 
850 do i=0, imax
851 do j=0, jmax
852 
853 #if (PRECIP_ANOM_INTERPOL==1)
854  do n=1, 12 ! month counter
855  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
856  end do
857 #elif (PRECIP_ANOM_INTERPOL==2)
858  do n=1, 12 ! month counter
859  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
860  end do
861 #else
862  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
863 #endif
864 
865 end do
866 end do
867 
868 #endif
869 
870 !-------- Mean accumulation --------
871 
872 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
873  ! mm/a water equiv. -> m/s ice equiv.
874 
875 !-------- Reading of present topography mask --------
876 
877 #if (GRID==0 || GRID==1)
878 
879 open(24, iostat=ios, &
880  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_present_file, &
881  recl=2048, status='old')
882 
883 #elif GRID==2
884 
885 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
886 
887 #endif
888 
889 if (ios /= 0) &
890  stop ' sico_init: Error when opening the mask file!'
891 
892 do m=1, 6; read(24, fmt='(a)') ch_dummy; end do
893 
894 do j=jmax, 0, -1
895  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
896 end do
897 
898 close(24, status='keep')
899 
900 !-------- Reading of data for present monthly-mean surface temperature --------
901 
902 #if (GRID==0 || GRID==1)
903 
904 open(21, iostat=ios, &
905  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_mm_present_file, &
906  recl=16384, status='old')
907 
908 #elif GRID==2
909 
910 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
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==0 || GRID==1)
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 open(21, iostat=ios, &
962  file=inpath//'/'//trim(ch_domain_short)//'/'//zs_present_file, &
963  recl=8192, status='old')
964 
965 #elif GRID==2
966 
967 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
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 scand 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==0 || GRID==1)
1339 
1340 do ir=-imax, imax
1341 do jr=-jmax, jmax
1342  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1343  ! distortion due to stereographic projection not accounted for
1344 end do
1345 end do
1346 
1347 #endif
1348 
1349 !-------- Initial velocities --------
1350 
1351 call calc_temp_melt()
1352 
1353 #if (DYNAMICS==1 || DYNAMICS==2)
1354 
1355 call calc_vxy_b_sia(time, z_sl)
1356 call calc_vxy_sia(dzeta_c, dzeta_t)
1357 
1358 #if (MARGIN==3 || DYNAMICS==2)
1359 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1360 #endif
1361 
1362 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1363 
1364 #if (MARGIN==3)
1365 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1366 #endif
1367 
1368 #elif (DYNAMICS==0)
1369 
1370 call calc_vxy_static()
1371 call calc_vz_static()
1372 
1373 #else
1374  stop ' sico_init: DYNAMICS must be either 0, 1 or 2!'
1375 #endif
1376 
1377 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1378 
1379 !-------- Initial enthalpies --------
1380 
1381 #if (CALCMOD==0 || CALCMOD==-1)
1382 
1383 do i=0, imax
1384 do j=0, jmax
1385 
1386  do kc=0, kcmax
1387  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1388  end do
1389 
1390  do kt=0, ktmax
1391  enth_t(kt,j,i) = enth_c(0,j,i)
1392  end do
1393 
1394 end do
1395 end do
1396 
1397 #elif (CALCMOD==1)
1398 
1399 do i=0, imax
1400 do j=0, jmax
1401 
1402  do kc=0, kcmax
1403  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1404  end do
1405 
1406  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1407  do kt=0, ktmax
1408  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1409  end do
1410  else
1411  do kt=0, ktmax
1412  enth_t(kt,j,i) = enth_c(0,j,i)
1413  end do
1414  end if
1415 
1416 end do
1417 end do
1418 
1419 #elif (CALCMOD==2 || CALCMOD==3)
1420 
1421 do i=0, imax
1422 do j=0, jmax
1423 
1424  do kc=0, kcmax
1425  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1426  end do
1427 
1428  do kt=0, ktmax
1429  enth_t(kt,j,i) = enth_c(0,j,i)
1430  end do
1431 
1432 end do
1433 end do
1434 
1435 #else
1436 
1437 stop ' sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1438 
1439 #endif
1440 
1441 !-------- Initialize time-series files --------
1442 
1443 ! ------ Time-series file for the ice sheet on the whole
1444 
1445 open(12, iostat=ios, &
1446  file=outpath//'/'//trim(runname)//'.ser', &
1447  status='new')
1448 if (ios /= 0) &
1449  stop ' sico_init: Error when opening the ser file!'
1450 
1451 if (forcing_flag == 1) then
1452 
1453  write(12,1102)
1454  write(12,1103)
1455 
1456 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1457 
1458  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1459  ' H_max(m) zs_max(m) V_g(m^3)', &
1460  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1461  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1462  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1463  1103 format('----------------------------------------------------', &
1464  '---------------------------------------')
1465 
1466 #elif OUTSER_V_A==2
1467 
1468  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1469  ' V(m^3) V_g(m^3) V_f(m^3)', &
1470  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1471  ' H_max(m) zs_max(m)', &
1472  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1473  ' A_t(m^2) V_bm(m^3/a)', &
1474  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1475  1103 format('----------------------------------------------------', &
1476  '---------------------------------------')
1477 
1478 #else
1479  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1480 #endif
1481 
1482 else if (forcing_flag == 2) then
1483 
1484  write(12,1112)
1485  write(12,1113)
1486 
1487 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1488 
1489  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1490  ' H_max(m) zs_max(m) V_g(m^3)', &
1491  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1492  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1493  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1494  1113 format('----------------------------------------------------', &
1495  '---------------------------------------')
1496 
1497 #elif OUTSER_V_A==2
1498 
1499  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1500  ' V(m^3) V_g(m^3) V_f(m^3)', &
1501  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1502  ' H_max(m) zs_max(m)', &
1503  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1504  ' A_t(m^2) V_bm(m^3/a)', &
1505  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1506  1113 format('----------------------------------------------------', &
1507  '---------------------------------------')
1508 
1509 #else
1510  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1511 #endif
1512 
1513 end if
1514 
1515 ! ------ Time-series file for the ice-thickness field
1516 
1517 #if OUTSER==2
1518 
1519 open(13, iostat=ios, &
1520  file=outpath//'/'//trim(runname)//'.thk', &
1521  status='new')
1522 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1523 
1524 if (forcing_flag == 1) then
1525 
1526  write(13,1104)
1527  write(13,1105)
1528 
1529  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1530  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1531  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1532  1105 format('----------------------------------------------------')
1533 
1534 else if (forcing_flag == 2) then
1535 
1536  write(13,1114)
1537  write(13,1115)
1538 
1539  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1540  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1541  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1542  1115 format('----------------------------------------------------')
1543 
1544 end if
1545 
1546 #endif
1547 
1548 !-------- Output of the initial state --------
1549 
1550 #if OUTPUT==1
1551 
1552 #if ERGDAT==0
1553  flag_3d_output = .false.
1554 #elif ERGDAT==1
1555  flag_3d_output = .true.
1556 #else
1557  stop ' sico_init: ERGDAT must be either 0 or 1!'
1558 #endif
1559 
1560 #if NETCDF==1
1561  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1562  flag_3d_output, ndat2d, ndat3d)
1563 #elif NETCDF==2
1564  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1565  flag_3d_output, ndat2d, ndat3d)
1566 #else
1567  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1568 #endif
1569 
1570 #elif OUTPUT==2
1571 
1572 if (time_output(1) <= time_init+eps) then
1573 
1574 #if ERGDAT==0
1575  flag_3d_output = .false.
1576 #elif ERGDAT==1
1577  flag_3d_output = .true.
1578 #else
1579  stop ' sico_init: ERGDAT must be either 0 or 1!'
1580 #endif
1581 
1582 #if NETCDF==1
1583  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1584  flag_3d_output, ndat2d, ndat3d)
1585 #elif NETCDF==2
1586  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1587  flag_3d_output, ndat2d, ndat3d)
1588 #else
1589  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1590 #endif
1591 
1592 end if
1593 
1594 #elif OUTPUT==3
1595 
1596  flag_3d_output = .false.
1597 
1598 #if NETCDF==1
1599  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1600  flag_3d_output, ndat2d, ndat3d)
1601 #elif NETCDF==2
1602  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1603  flag_3d_output, ndat2d, ndat3d)
1604 #else
1605  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1606 #endif
1607 
1608 if (time_output(1) <= time_init+eps) then
1609 
1610  flag_3d_output = .true.
1611 
1612 #if NETCDF==1
1613  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1614  flag_3d_output, ndat2d, ndat3d)
1615 #elif NETCDF==2
1616  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1617  flag_3d_output, ndat2d, ndat3d)
1618 #else
1619  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1620 #endif
1621 
1622 end if
1623 
1624 #else
1625  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1626 #endif
1627 
1628 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1629 #if OUTSER==2
1630 call output3(time_init, delta_ts, glac_index, z_sl)
1631 #endif
1632 
1633 end subroutine sico_init
1634 !
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