SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
sico_init.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : s i c o _ i n i t
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Initialisations for SICOPOLIS.
34 !<------------------------------------------------------------------------------
35 subroutine sico_init(delta_ts, glac_index, &
36  mean_accum, mean_accum_inv, &
37  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38  time, time_init, time_end, time_output, &
39  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40  z_sl, dzsl_dtau, z_mar, &
41  ii, jj, nn, &
42  ndat2d, ndat3d, n_output, &
43  runname)
44 
45 use sico_types
46 use sico_variables
47 #if (NETCDF > 1)
48 use netcdf
49 use nc_check
50 #endif
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, mean_accum_inv
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, ios1, ios2, ios3, ios4
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, target_topo_dat_name
75 character(len=256) :: shell_command
76 character :: ch_dummy
77 logical :: flag_3d_output
78 
79 #if (NETCDF > 1)
80 integer(i4b) :: dimid
81 ! dimid: Dimension ID
82 integer(i4b) :: ncv
83 ! ncv: Variable ID
84 #endif
85 
86 character(len=64), parameter :: thisroutine = 'sico_init'
87 
88 character(len=64), parameter :: fmt1 = '(a)', &
89  fmt2 = '(a,i4)', &
90  fmt3 = '(a,es12.4)'
91 
92 character(len= 8) :: ch_imax
93 character(len=128) :: fmt4
94 
95 write(ch_imax, fmt='(i8)') imax
96 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
97 
98 !-------- Some initial values --------
99 
100 n_output = 0
101 
102 dtime = 0.0_dp
103 dtime_temp = 0.0_dp
104 dtime_wss = 0.0_dp
105 dtime_out = 0.0_dp
106 dtime_ser = 0.0_dp
107 
108 time = 0.0_dp
109 time_init = 0.0_dp
110 time_end = 0.0_dp
111 time_output = 0.0_dp
112 
113 !-------- Initialisation of the Library of Iterative Solvers Lis,
114 ! if required --------
115 
116 #if (CALCZS==4 || MARGIN==3)
117 
118 #include "lisf.h"
119 call lis_initialize(ierr)
120 
121 #endif
122 
123 !-------- Read physical parameters --------
124 
125 call phys_para()
126 
127 !-------- Compatibility check of the SICOPOLIS version with the header file
128 
129 if ( trim(version) /= trim(sico_version) ) &
130  stop ' sico_init: SICOPOLIS version not compatible with header file!'
131 
132 !-------- Compatibility check of the horizontal resolution with the
133 ! number of grid points --------
134 
135 #if (GRID==0 || GRID==1)
136 
137 if (abs(dx-20.0_dp) < eps) then
138 
139  if ( (imax /= 75).or.(jmax /= 140) ) then
140  stop ' sico_init: IMAX and/or JMAX wrong!'
141  end if
142 
143 else if (abs(dx-10.0_dp) < eps) then
144 
145  if ( (imax /= 150).or.(jmax /= 280) ) then
146  stop ' sico_init: IMAX and/or JMAX wrong!'
147  end if
148 
149 else if (abs(dx-5.0_dp) < eps) then
150 
151  if ( (imax /= 300).or.(jmax /= 560) ) then
152  stop ' sico_init: IMAX and/or JMAX wrong!'
153  end if
154 
155 else
156  stop ' sico_init: DX wrong!'
157 end if
158 
159 #elif GRID==2
160 
161 stop ' sico_init: GRID==2 not allowed for this application!'
162 
163 #endif
164 
165 !-------- Compatibility check of surface-temperature
166 ! and precipitation switches --------
167 
168 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
169 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
170 #endif
171 
172 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
173 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
174 #endif
175 
176 #if (TSURFACE == 6)
177 #if (ACCSURFACE != 6)
178 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
179 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
180 stop
181 #elif (NETCDF != 2)
182 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
183 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
184 stop
185 #endif
186 #endif
187 
188 #if (ACCSURFACE == 6)
189 #if (TSURFACE != 6)
190 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
191 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
192 stop
193 #elif (NETCDF != 2)
194 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
195 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
196 stop
197 #endif
198 #endif
199 
200 !-------- Compatibility check of discretization schemes for the horizontal and
201 ! vertical advection terms in the temperature and age equations --------
202 
203 #if ADV_HOR==1
204 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
205 #endif
206 
207 !-------- Setting of forcing flag --------
208 
209 #if (TSURFACE <= 4)
210 
211 forcing_flag = 1 ! forcing by delta_ts
212 
213 #elif (TSURFACE == 5)
214 
215 forcing_flag = 2 ! forcing by glac_index
216 
217 #elif (TSURFACE == 6)
218 
219 forcing_flag = 3 ! forcing by time-dependent surface temperature
220  ! and precipitation data
221 
222 #endif
223 
224 !-------- Initialization of numerical time steps --------
225 
226 dtime0 = dtime0
227 dtime_temp0 = dtime_temp0
228 #if REBOUND==2
229 dtime_wss0 = dtime_wss0
230 #endif
231 
232 !-------- Further initializations --------
233 
234 dzeta_c = 1.0_dp/real(kcmax,dp)
235 dzeta_t = 1.0_dp/real(ktmax,dp)
236 dzeta_r = 1.0_dp/real(krmax,dp)
237 
238 ndat2d = 1
239 ndat3d = 1
240 
241 !-------- Construction of a vector (with index n) from a 2-d array
242 ! (with indices i, j) by diagonal numbering --------
243 
244 n=1
245 
246 do m=0, imax+jmax
247  do i=m, 0, -1
248  j = m-i
249  if ((i <= imax).and.(j <= jmax)) then
250  ii(n) = i
251  jj(n) = j
252  nn(j,i) = n
253  n=n+1
254  end if
255  end do
256 end do
257 
258 !-------- Specification of current simulation --------
259 
260 runname = runname
261 anfdatname = anfdatname
262 
263 #if defined(YEAR_ZERO)
264 year_zero = year_zero
265 #else
266 year_zero = 2000.0_dp ! default value 2000 CE
267 #endif
268 
269 time_init0 = time_init0
270 time_end0 = time_end0
271 dtime_ser0 = dtime_ser0
272 
273 #if ( OUTPUT==1 || OUTPUT==3 )
274 dtime_out0 = dtime_out0
275 #endif
276 
277 #if ( OUTPUT==2 || OUTPUT==3 )
278 n_output = n_output
279 time_output0( 1) = time_out0_01
280 time_output0( 2) = time_out0_02
281 time_output0( 3) = time_out0_03
282 time_output0( 4) = time_out0_04
283 time_output0( 5) = time_out0_05
284 time_output0( 6) = time_out0_06
285 time_output0( 7) = time_out0_07
286 time_output0( 8) = time_out0_08
287 time_output0( 9) = time_out0_09
288 time_output0(10) = time_out0_10
289 time_output0(11) = time_out0_11
290 time_output0(12) = time_out0_12
291 time_output0(13) = time_out0_13
292 time_output0(14) = time_out0_14
293 time_output0(15) = time_out0_15
294 time_output0(16) = time_out0_16
295 time_output0(17) = time_out0_17
296 time_output0(18) = time_out0_18
297 time_output0(19) = time_out0_19
298 time_output0(20) = time_out0_20
299 #endif
300 
301 !-------- Write log file --------
302 
303 shell_command = 'if [ ! -d'
304 shell_command = trim(shell_command)//' '//outpath
305 shell_command = trim(shell_command)//' '//'] ; then mkdir'
306 shell_command = trim(shell_command)//' '//outpath
307 shell_command = trim(shell_command)//' '//'; fi'
308 call system(trim(shell_command))
309  ! Check whether directory OUTPATH exists. If not, it is created.
310 
311 open(10, iostat=ios, &
312  file=outpath//'/'//trim(runname)//'.log', &
313  status='new')
314 if (ios /= 0) &
315  stop ' sico_init: Error when opening the log file!'
316 
317 write(10, fmt=trim(fmt1)) 'Computational domain:'
318 #if defined(ANT)
319 write(10, fmt=trim(fmt1)) 'Antarctica'
320 #elif defined(GRL)
321 write(10, fmt=trim(fmt1)) 'Greenland'
322 #elif defined(SCAND)
323 write(10, fmt=trim(fmt1)) 'Scandinavia and Eurasia'
324 #elif defined(NHEM)
325 write(10, fmt=trim(fmt1)) 'Northern hemisphere'
326 #elif defined(EMTP2SGE)
327 write(10, fmt=trim(fmt1)) 'EISMINT Phase 2 Simplified Geometry Experiment'
328 #elif defined(NMARS)
329 write(10, fmt=trim(fmt1)) 'North polar cap of Mars'
330 #elif defined(SMARS)
331 write(10, fmt=trim(fmt1)) 'South polar cap of Mars'
332 #else
333 stop ' sico_init: No valid domain specified!'
334 #endif
335 write(10, fmt=trim(fmt1)) ' '
336 
337 write(10, fmt=trim(fmt2)) 'imax =', imax
338 write(10, fmt=trim(fmt2)) 'jmax =', jmax
339 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
340 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
341 write(10, fmt=trim(fmt2)) 'krmax =', krmax
342 write(10, fmt=trim(fmt1)) ' '
343 
344 write(10, fmt=trim(fmt3)) 'a =', deform
345 write(10, fmt=trim(fmt1)) ' '
346 
347 #if (GRID==0 || GRID==1)
348 write(10, fmt=trim(fmt3)) 'x0 =', x0
349 write(10, fmt=trim(fmt3)) 'y0 =', y0
350 write(10, fmt=trim(fmt3)) 'dx =', dx
351 #elif GRID==2
352 stop ' sico_init: GRID==2 not allowed for this application!'
353 #endif
354 write(10, fmt=trim(fmt1)) ' '
355 
356 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
357 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
358 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
359 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
360 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
361 #if REBOUND==2
362 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
363 #endif
364 #if ( OUTPUT==1 || OUTPUT==3 )
365 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
366 #endif
367 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
368 write(10, fmt=trim(fmt1)) ' '
369 
370 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
371 #if ANF_DAT==1
372 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
373 #endif
374 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
375 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
376 #if ANF_DAT==3
377 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
378 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
379 #endif
380 write(10, fmt=trim(fmt1)) ' '
381 
382 #if ZS_EVOL==2
383 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
384 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
385 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
386 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
387 write(10, fmt=trim(fmt1)) ' '
388 #endif
389 
390 #if ZS_EVOL==3
391 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
392 write(10, fmt=trim(fmt1)) ' '
393 #endif
394 
395 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
396 write(10, fmt=trim(fmt1)) ' '
397 
398 #if (CALCZS==3 || CALCZS==4)
399 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
400 #if CALCZS==3
401 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
402 #endif
403 write(10, fmt=trim(fmt1)) ' '
404 #endif
405 
406 #if TSURFACE==1
407 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
408 #elif TSURFACE==3
409 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
410 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
411 #elif TSURFACE==4
412 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
413 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
414 #elif TSURFACE==5
415 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
416 write(10, fmt=trim(fmt1)) 'temp_ma_anom file = '//temp_ma_anom_file
417 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
418 write(10, fmt=trim(fmt1)) 'temp_mj_anom file = '//temp_mj_anom_file
419 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
420 #elif TSURFACE==6
421 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
422 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
423 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
424 #endif
425 
426 write(10, fmt=trim(fmt1)) 'precip_present file = '//precip_present_file
427 #if ACCSURFACE==1
428 write(10, fmt=trim(fmt3)) 'accfact =', accfact
429 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
430 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
431 #elif ( ACCSURFACE==5 )
432 write(10, fmt=trim(fmt1)) 'precip_anom file = '//precip_anom_file
433 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
434 #elif ( ACCSURFACE==6 )
435 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
436 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
437 #endif
438 #if (ACCSURFACE <= 3)
439 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
440 #if (ELEV_DESERT == 1)
441 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
442 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
443 #endif
444 #endif
445 
446 #if ABLSURFACE==3
447 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
448 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
449 #endif
450 
451 #if SEA_LEVEL==1
452 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
453 #elif SEA_LEVEL==3
454 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
455 #endif
456 write(10, fmt=trim(fmt1)) ' '
457 
458 #if MARGIN==2
459 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
460 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
461 write(10, fmt=trim(fmt1)) ' '
462 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
463  || marine_ice_calving==6 || marine_ice_calving==7 )
464 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
465 write(10, fmt=trim(fmt1)) ' '
466 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
467 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
468 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
469 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
470 write(10, fmt=trim(fmt1)) ' '
471 #endif
472 #elif MARGIN==3
473 #if ICE_SHELF_CALVING==2
474 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
475 write(10, fmt=trim(fmt1)) ' '
476 #endif
477 #endif
478 
479 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
480 write(10, fmt=trim(fmt3)) 'smw_coeff =', smw_coeff
481 write(10, fmt=trim(fmt2)) 'r_smw =', r_smw
482 write(10, fmt=trim(fmt2)) 's_smw =', s_smw
483 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
484 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
485 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
486 #if SLIDE_LAW==2
487 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
488 #endif
489 write(10, fmt=trim(fmt1)) ' '
490 
491 #if ICE_STREAM==2
492 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
493 write(10, fmt=trim(fmt1)) ' '
494 
495 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
496 write(10, fmt=trim(fmt3)) 'smw_coeff_sedi =', smw_coeff_sedi
497 write(10, fmt=trim(fmt2)) 'r_smw_sedi =', r_smw_sedi
498 write(10, fmt=trim(fmt2)) 's_smw_sedi =', s_smw_sedi
499 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
500 write(10, fmt=trim(fmt2)) 'p_weert_sedi =', p_weert_sedi
501 write(10, fmt=trim(fmt2)) 'q_weert_sedi =', q_weert_sedi
502 write(10, fmt=trim(fmt1)) ' '
503 #endif
504 
505 #if Q_GEO_MOD==1
506 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
507 #elif Q_GEO_MOD==2
508 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
509 #endif
510 write(10, fmt=trim(fmt1)) ' '
511 
512 #if defined(MARINE_ICE_BASAL_MELTING)
513 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
514 #if ( MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3 )
515 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
516 #endif
517 write(10, fmt=trim(fmt1)) ' '
518 #endif
519 
520 #if REBOUND==1
521 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
522 write(10, fmt=trim(fmt1)) ' '
523 #endif
524 
525 #if FLOW_LAW==2
526 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
527 write(10, fmt=trim(fmt1)) ' '
528 #endif
529 #if FIN_VISC==2
530 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
531 write(10, fmt=trim(fmt1)) ' '
532 #endif
533 
534 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
535 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
536 #if ( ENHMOD==2 || ENHMOD==3 )
537 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
538 #endif
539 #if ENHMOD==2
540 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
541 #endif
542 #if ENHMOD==3
543 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
544 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
545 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
546 #endif
547 #if MARGIN==3
548 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
549 #endif
550 write(10, fmt=trim(fmt1)) ' '
551 
552 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
553 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
554 write(10, fmt=trim(fmt3)) 'H_min =', h_min
555 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
556 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
557 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
558 #if defined(QBM_MIN)
559 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
560 #elif defined(QB_MIN) /* obsolete */
561 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
562 #endif
563 #if defined(QBM_MAX)
564 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
565 #elif defined(QB_MAX) /* obsolete */
566 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
567 #endif
568 write(10, fmt=trim(fmt3)) 'age_min =', age_min
569 write(10, fmt=trim(fmt3)) 'age_max =', age_max
570 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
571 #if ADV_VERT==1
572 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
573 #endif
574 write(10, fmt=trim(fmt1)) ' '
575 
576 write(10, fmt=trim(fmt2)) 'GRID =', grid
577 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
578 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
579 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
580 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
581 #if MARGIN==2
582 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
583 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
584 #elif MARGIN==3
585 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
586 #endif
587 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
588 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
589 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
590 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
591 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
592 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
593 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
594 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
595 write(10, fmt=trim(fmt1)) ' '
596 write(10, fmt=trim(fmt2)) 'TEMP_PRESENT_PARA =', temp_present_para
597 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
598 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
599 #if ACCSURFACE==5
600 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
601 #endif
602 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
603 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
604 #if ( ABLSURFACE==1 || ABLSURFACE==2 )
605 write(10, fmt=trim(fmt2)) 'PDD_MODIFIER =', pdd_modifier
606 #if PDD_MODIFIER==2
607 write(10, fmt=trim(fmt2)) 'N_PDD_MOD =', n_pdd_mod
608 write(10, fmt=trim(fmt3)) 'LON_W_E_SEP =', lon_w_e_sep
609 do n=1, n_pdd_mod
610  if (n== 1) then
611  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_01 =', pdd_mod_lat_01
612  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_01 =', pdd_mod_fac_w_01
613  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_01 =', pdd_mod_fac_e_01
614  end if
615  if (n== 2) then
616  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_02 =', pdd_mod_lat_02
617  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_02 =', pdd_mod_fac_w_02
618  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_02 =', pdd_mod_fac_e_02
619  end if
620  if (n== 3) then
621  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_03 =', pdd_mod_lat_03
622  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_03 =', pdd_mod_fac_w_03
623  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_03 =', pdd_mod_fac_e_03
624  end if
625  if (n== 4) then
626  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_04 =', pdd_mod_lat_04
627  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_04 =', pdd_mod_fac_w_04
628  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_04 =', pdd_mod_fac_e_04
629  end if
630  if (n== 5) then
631  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_05 =', pdd_mod_lat_05
632  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_05 =', pdd_mod_fac_w_05
633  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_05 =', pdd_mod_fac_e_05
634  end if
635  if (n== 6) then
636  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_06 =', pdd_mod_lat_06
637  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_06 =', pdd_mod_fac_w_06
638  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_06 =', pdd_mod_fac_e_06
639  end if
640  if (n== 7) then
641  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_07 =', pdd_mod_lat_07
642  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_07 =', pdd_mod_fac_w_07
643  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_07 =', pdd_mod_fac_e_07
644  end if
645  if (n== 8) then
646  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_08 =', pdd_mod_lat_08
647  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_08 =', pdd_mod_fac_w_08
648  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_08 =', pdd_mod_fac_e_08
649  end if
650  if (n== 9) then
651  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_09 =', pdd_mod_lat_09
652  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_09 =', pdd_mod_fac_w_09
653  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_09 =', pdd_mod_fac_e_09
654  end if
655  if (n==10) then
656  write(10, fmt=trim(fmt3)) 'PDD_MOD_LAT_10 =', pdd_mod_lat_10
657  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_W_10 =', pdd_mod_fac_w_10
658  write(10, fmt=trim(fmt3)) 'PDD_MOD_FAC_E_10 =', pdd_mod_fac_e_10
659  end if
660 end do
661 #endif
662 #endif
663 write(10, fmt=trim(fmt1)) ' '
664 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
665 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
666 write(10, fmt=trim(fmt2)) 'ICE_STREAM =', ice_stream
667 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
668 write(10, fmt=trim(fmt1)) ' '
669 
670 #if defined(OUT_TIMES)
671 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
672 #endif
673 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
674 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
675 #if defined(OUTSER_V_A)
676 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
677 #endif
678 #if ( OUTPUT==1 || OUTPUT==2 )
679 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
680 #endif
681 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
682 write(10, fmt=trim(fmt1)) ' '
683 #if ( OUTPUT==2 || OUTPUT==3 )
684 write(10, fmt=trim(fmt2)) 'n_output =', n_output
685 do n=1, n_output
686  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
687 end do
688 write(10, fmt=trim(fmt1)) ' '
689 #endif
690 
691 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
692 
693 close(10, status='keep')
694 
695 !-------- Conversion of time quantities --------
696 
697 year_zero = year_zero*year_sec ! a --> s
698 time_init = time_init0*year_sec ! a --> s
699 time_end = time_end0*year_sec ! a --> s
700 dtime = dtime0*year_sec ! a --> s
701 dtime_temp = dtime_temp0*year_sec ! a --> s
702 #if REBOUND==2
703 dtime_wss = dtime_wss0*year_sec ! a --> s
704 #endif
705 dtime_ser = dtime_ser0*year_sec ! a --> s
706 #if ( OUTPUT==1 || OUTPUT==3 )
707 dtime_out = dtime_out0*year_sec ! a --> s
708 #endif
709 #if ( OUTPUT==2 || OUTPUT==3 )
710 do n=1, n_output
711  time_output(n) = time_output0(n)*year_sec ! a --> s
712 end do
713 #endif
714 
715 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
716  stop ' sico_init: dtime_temp must be a multiple of dtime!'
717 
718 #if REBOUND==2
719 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
720  stop ' sico_init: dtime_wss must be a multiple of dtime!'
721 #endif
722 
723 #if ZS_EVOL==2
724 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
725 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
726 #endif
727 
728 time = time_init
729 
730 !-------- Reading of present mean-annual precipitation rate --------
731 
732 #if (GRID==0 || GRID==1)
733 
734 open(21, iostat=ios, &
735  file=inpath//'/grl/'//precip_present_file, &
736  recl=8192, status='old')
737 
738 #elif GRID==2
739 
740 stop ' sico_init: GRID==2 not allowed for this application!'
741 
742 #endif
743 
744 if (ios /= 0) stop ' sico_init: Error when opening the precip file!'
745 
746 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
747 
748 do j=jmax, 0, -1
749  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
750 end do
751 
752 close(21, status='keep')
753 
754 precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
755  ! mm/a water equiv. -> m/s ice equiv.
756 
757 !-------- Present monthly precipitation rates
758 ! (assumed to be equal to the mean annual precipitation rate) --------
759 
760 do n=1, 12 ! month counter
761 do i=0, imax
762 do j=0, jmax
763  precip_present(j,i,n) = precip_ma_present(j,i)
764 end do
765 end do
766 end do
767 
768 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
769 
770 #if ACCSURFACE==5
771 
772 #if (GRID==0 || GRID==1)
773 
774 open(21, iostat=ios, &
775  file=inpath//'/grl/'//precip_anom_file, &
776  recl=8192, status='old')
777 
778 #endif
779 
780 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
781 
782 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
783 
784 do j=jmax, 0, -1
785  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
786 end do
787 
788 close(21, status='keep')
789 
790 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
791 
792 !-------- LGM monthly precipitation-rate anomalies (assumed to be
793 ! equal to the mean annual precipitation-rate anomaly) --------
794 
795 do n=1, 12 ! month counter
796 do i=0, imax
797 do j=0, jmax
798 
799  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
800  ! positive values ensured
801 
802 #if (PRECIP_ANOM_INTERPOL==1)
803  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
804 #elif (PRECIP_ANOM_INTERPOL==2)
805  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
806 #else
807  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
808 #endif
809 
810 end do
811 end do
812 end do
813 
814 #endif
815 
816 !-------- Mean accumulation --------
817 
818 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
819  ! mm/a water equiv. -> m/s ice equiv.
820 
821 mean_accum_inv = 1.0_dp/mean_accum
822 
823 !-------- Read present topography mask --------
824 
825 #if (GRID==0 || GRID==1)
826 
827 open(24, iostat=ios, &
828  file=inpath//'/grl/'//mask_present_file, &
829  recl=8192, status='old')
830 
831 #elif GRID==2
832 
833 stop ' sico_init: GRID==2 not allowed for this application!'
834 
835 #endif
836 
837 if (ios /= 0) stop ' sico_init: Error when opening the mask file!'
838 
839 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
840 
841 do j=jmax, 0, -1
842  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
843 end do
844 
845 close(24, status='keep')
846 
847 !-------- Read rock/sediment mask --------
848 
849 #if ICE_STREAM==2
850 
851 open(24, iostat=ios, &
852  file=inpath//'/grl/'//mask_sedi_file, &
853  recl=8192, status='old')
854 
855 if (ios /= 0) stop ' sico_init: Error when opening the rock/sediment mask file!'
856 
857 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
858 
859 do j=jmax, 0, -1
860  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
861 end do
862 
863 close(24, status='keep')
864 
865 #endif
866 
867 !-------- Present reference elevation (for precipitation data) --------
868 
869 #if (GRID==0 || GRID==1)
870 
871 open(21, iostat=ios, &
872  file=inpath//'/grl/'//zs_present_file, &
873  recl=8192, status='old')
874 
875 #elif GRID==2
876 
877 stop ' sico_init: GRID==2 not allowed for this application!'
878 
879 #endif
880 
881 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
882 
883 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
884 
885 do j=jmax, 0, -1
886  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
887 end do
888 
889 close(21, status='keep')
890 
891 ! ------ Reset bathymetry data (sea floor elevation) to the
892 ! sea surface
893 
894 do i=0, imax
895 do j=0, jmax
896  if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
897 end do
898 end do
899 
900 !-------- Reading of the prescribed target topography --------
901 
902 #if ZS_EVOL==2
903 
904 target_topo_dat_name = trim(target_topo_dat_name)
905 
906 #if NETCDF==1
907 write(6, fmt='(a)') ' sico_init: Reading of target topography'
908 write(6, fmt='(a)') ' only implemented for NetCDF files (NETCDF==2)!'
909 stop
910 #elif NETCDF==2
911 call read_target_topo_nc(target_topo_dat_name)
912 #else
913 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
914 #endif
915 
916 #endif
917 
918 !-------- Reading of the maximum ice extent mask --------
919 
920 #if ZS_EVOL==3
921 
922 #if (GRID==0 || GRID==1)
923 
924 open(24, iostat=ios, &
925  file=inpath//'/grl/'//mask_maxextent_file, &
926  recl=8192, status='old')
927 
928 #elif GRID==2
929 
930 stop ' sico_init: GRID==2 not allowed for this application!'
931 
932 #endif
933 
934 if (ios /= 0) &
935  stop ' sico_init: Error when opening the maximum ice extent mask file!'
936 
937 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
938 
939 do j=jmax, 0, -1
940  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
941 end do
942 
943 close(24, status='keep')
944 
945 #endif
946 
947 !-------- Reading of LGM mean-annual and mean-July surface-temperature
948 ! anomalies --------
949 
950 #if TSURFACE==5
951 
952 #if (GRID==0 || GRID==1)
953 
954 open(21, iostat=ios, &
955  file=inpath//'/grl/'//temp_ma_anom_file, &
956  recl=8192, status='old')
957 
958 #endif
959 
960 if (ios /= 0) stop ' sico_init: Error when opening the temp_ma anomaly file!'
961 
962 #if (GRID==0 || GRID==1)
963 
964 open(22, iostat=ios, &
965  file=inpath//'/grl/'//temp_mj_anom_file, &
966  recl=8192, status='old')
967 
968 #endif
969 
970 if (ios /= 0) stop ' sico_init: Error when opening the temp_mj anomaly file!'
971 
972 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
973 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
974 
975 do j=jmax, 0, -1
976  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
977  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
978 end do
979 
980 close(21, status='keep')
981 close(22, status='keep')
982 
983 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
984 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
985 
986 #endif
987 
988 !-------- Read data for delta_ts --------
989 
990 #if TSURFACE==4
991 
992 open(21, iostat=ios, &
993  file=inpath//'/general/'//grip_temp_file, &
994  status='old')
995 if (ios /= 0) &
996  stop ' sico_init: Error when opening the data file for delta_ts!'
997 
998 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
999 
1000 if (ch_dummy /= '#') then
1001  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
1002  write(6, fmt=*) ' not defined in data file for delta_ts!'
1003  stop
1004 end if
1005 
1006 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
1007 
1008 allocate(griptemp(0:ndata_grip))
1009 
1010 do n=0, ndata_grip
1011  read(21, fmt=*) d_dummy, griptemp(n)
1012 end do
1013 
1014 close(21, status='keep')
1015 
1016 #endif
1017 
1018 !-------- Read data for the glacial index --------
1019 
1020 #if TSURFACE==5
1021 
1022 open(21, iostat=ios, &
1023  file=inpath//'/general/'//glac_ind_file, status='old')
1024 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
1025 
1026 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1027 
1028 if (ch_dummy /= '#') then
1029  write(6, fmt=*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
1030  write(6, fmt=*) ' not defined in glacial-index file!'
1031  stop
1032 end if
1033 
1034 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
1035 
1036 allocate(glacial_index(0:ndata_gi))
1037 
1038 do n=0, ndata_gi
1039  read(21, fmt=*) d_dummy, glacial_index(n)
1040 end do
1041 
1042 close(21, status='keep')
1043 
1044 #endif
1045 
1046 !-------- Reading of time-dependent surface-temperature
1047 ! and precipitation anomaly data --------
1048 
1049 #if ( (TSURFACE==6) && (ACCSURFACE==6) )
1050 
1051 call check( nf90_open(inpath//'/grl/'//temp_precip_anom_file, &
1052  nf90_nowrite, ncid_temp_precip), thisroutine )
1053 
1054 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1055  thisroutine )
1056 
1057 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1058  len = ndata_temp_precip), thisroutine )
1059 
1060 ndata_temp_precip = ndata_temp_precip-1
1061 
1062 allocate(temp_precip_time(0:ndata_temp_precip))
1063 
1064 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1065 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1066  ! in a
1067 
1068 temp_precip_time_min = temp_precip_time(0) ! in a
1069 temp_precip_time_stp = temp_precip_time(1)-temp_precip_time(0)
1070  ! in a; constant time step assumed
1071 temp_precip_time_max = temp_precip_time(ndata_temp_precip) ! in a
1072 
1073 #endif
1074 
1075 !-------- Read data for z_sl --------
1076 
1077 #if SEA_LEVEL==3
1078 
1079 open(21, iostat=ios, &
1080  file=inpath//'/general/'//sea_level_file, &
1081  status='old')
1082 if (ios /= 0) &
1083  stop ' sico_init: Error when opening the data file for z_sl!'
1084 
1085 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1086 
1087 if (ch_dummy /= '#') then
1088  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1089  write(6, fmt=*) ' not defined in data file for z_sl!'
1090  stop
1091 end if
1092 
1093 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
1094 
1095 allocate(specmap_zsl(0:ndata_specmap))
1096 
1097 do n=0, ndata_specmap
1098  read(21, fmt=*) d_dummy, specmap_zsl(n)
1099 end do
1100 
1101 close(21, status='keep')
1102 
1103 #endif
1104 
1105 !-------- Determination of the geothermal heat flux --------
1106 
1107 #if Q_GEO_MOD==1
1108 
1109 ! ------ Constant value
1110 
1111 do i=0, imax
1112 do j=0, jmax
1113  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1114 end do
1115 end do
1116 
1117 #elif Q_GEO_MOD==2
1118 
1119 ! ------ Read data from file
1120 
1121 open(21, iostat=ios, &
1122  file=inpath//'/grl/'//q_geo_file, &
1123  recl=8192, status='old')
1124 
1125 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
1126 
1127 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1128 
1129 do j=jmax, 0, -1
1130  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1131 end do
1132 
1133 close(21, status='keep')
1134 
1135 do i=0, imax
1136 do j=0, jmax
1137  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1138 end do
1139 end do
1140 
1141 #endif
1142 
1143 !-------- Reading of tabulated kei function--------
1144 
1145 #if (REBOUND==0 || REBOUND==1)
1146 
1147 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1148  ! dummy values
1149 #elif REBOUND==2
1150 
1151 call read_kei()
1152 
1153 #endif
1154 
1155 !-------- Definition of initial values --------
1156 
1157 ! ------ Present topography
1158 
1159 #if ANF_DAT==1
1160 
1161 call topography1(dxi, deta)
1162 
1163 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1164 
1165 call boundary(time_init, dtime, dxi, deta, &
1166  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1167 
1168 do i=0, imax
1169 do j=0, jmax
1170 
1171  q_bm(j,i) = 0.0_dp
1172  q_tld(j,i) = 0.0_dp
1173  q_b_tot(j,i) = 0.0_dp
1174 
1175  p_b_w(j,i) = 0.0_dp
1176  h_w(j,i) = 0.0_dp
1177 
1178  do kc=0, kcmax
1179  temp_c(kc,j,i) = -10.0_dp
1180  age_c(kc,j,i) = 15000.0_dp*year_sec
1181  end do
1182 
1183  do kt=0, ktmax
1184  omega_t(kt,j,i) = 0.0_dp
1185  age_t(kt,j,i) = 15000.0_dp*year_sec
1186  end do
1187 
1188  do kr=0, krmax
1189  temp_r(kr,j,i) = -10.0_dp+q_geo(j,i)*h_r/kappa_r &
1190  *(1.0_dp-real(kr,dp)/real(krmax,dp))
1191 ! (linear temperature distribution according to the
1192 ! geothermal heat flux)
1193  end do
1194 
1195 end do
1196 end do
1197 
1198 call calc_enhance(time_init)
1199 
1200 ! ------ Ice-free, relaxed bedrock
1201 
1202 #elif ANF_DAT==2
1203 
1204 call topography2(dxi, deta)
1205 
1206 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1207 
1208 call boundary(time_init, dtime, dxi, deta, &
1209  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1210 
1211 do i=0, imax
1212 do j=0, jmax
1213 
1214  q_bm(j,i) = 0.0_dp
1215  q_tld(j,i) = 0.0_dp
1216  q_b_tot(j,i) = 0.0_dp
1217 
1218  p_b_w(j,i) = 0.0_dp
1219  h_w(j,i) = 0.0_dp
1220 
1221  do kc=0, kcmax
1222  temp_c(kc,j,i) = temp_s(j,i)
1223  age_c(kc,j,i) = 15000.0_dp*year_sec
1224  end do
1225 
1226  do kt=0, ktmax
1227  omega_t(kt,j,i) = 0.0_dp
1228  age_t(kt,j,i) = 15000.0_dp*year_sec
1229  end do
1230 
1231  do kr=0, krmax
1232  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
1233  *(1.0_dp-real(kr,dp)/real(krmax,dp))
1234 ! (linear temperature distribution according to the
1235 ! geothermal heat flux)
1236  end do
1237 
1238 end do
1239 end do
1240 
1241 call calc_enhance(time_init)
1242 
1243 ! ------ Read initial state from output of previous simulation
1244 
1245 #elif ANF_DAT==3
1246 
1247 #if NETCDF==1
1248 call topography3(dxi, deta, z_sl, anfdatname)
1249 #elif NETCDF==2
1250 call topography3_nc(dxi, deta, z_sl, anfdatname)
1251 #else
1252 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1253 #endif
1254 
1255 call boundary(time_init, dtime, dxi, deta, &
1256  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1257 
1258 q_b_tot = q_bm + q_tld
1259 
1260 call calc_enhance(time_init)
1261 
1262 #endif
1263 
1264 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1265 
1266 #if (GRID==0 || GRID==1)
1267 
1268 do ir=-imax, imax
1269 do jr=-jmax, jmax
1270  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1271  ! distortion due to stereographic projection not accounted for
1272 end do
1273 end do
1274 
1275 #endif
1276 
1277 !-------- General abbreviations --------
1278 
1279 ea = exp(deform)
1280 
1281 do kc=0, kcmax
1282  zeta_c(kc) = kc*dzeta_c
1283  eaz_c(kc) = exp(deform*zeta_c(kc))
1284  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
1285 end do
1286 
1287 do kt=0, ktmax
1288  zeta_t(kt) = kt*dzeta_t
1289 end do
1290 
1291 !-------- Initial velocities --------
1292 
1293 call calc_temp_melt()
1294 
1295 call calc_vxy_b_sia(time, z_sl)
1296 call calc_vxy_sia(dzeta_c, dzeta_t)
1297 
1298 #if MARGIN==3
1299 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
1300 #endif
1301 
1302 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1303 
1304 #if MARGIN==3
1305 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1306 #endif
1307 
1308 !-------- Initialize time-series files --------
1309 
1310 ! ------ Time-series file for the ice sheet on the whole
1311 
1312 open(12, iostat=ios, &
1313  file=outpath//'/'//trim(runname)//'.ser', &
1314  status='new')
1315 if (ios /= 0) &
1316  stop ' sico_init: Error when opening the ser file!'
1317 
1318 if (forcing_flag == 1) then
1319 
1320  write(12,1102)
1321  write(12,1103)
1322 
1323 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1324 
1325  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1326  ' H_max(m) zs_max(m) V_g(m^3)', &
1327  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1328  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1329  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1330  1103 format('----------------------------------------------------', &
1331  '---------------------------------------')
1332 
1333 #elif OUTSER_V_A==2
1334 
1335  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1336  ' V(m^3) V_g(m^3) V_f(m^3)', &
1337  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1338  ' H_max(m) zs_max(m)', &
1339  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1340  ' A_t(m^2) V_bm(m^3/a)', &
1341  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1342  1103 format('----------------------------------------------------', &
1343  '---------------------------------------')
1344 
1345 #else
1346  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1347 #endif
1348 
1349 else if (forcing_flag == 2) then
1350 
1351  write(12,1112)
1352  write(12,1113)
1353 
1354 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1355 
1356  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1357  ' H_max(m) zs_max(m) V_g(m^3)', &
1358  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1359  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1360  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1361  1113 format('----------------------------------------------------', &
1362  '---------------------------------------')
1363 
1364 #elif OUTSER_V_A==2
1365 
1366  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1367  ' V(m^3) V_g(m^3) V_f(m^3)', &
1368  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1369  ' H_max(m) zs_max(m)', &
1370  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1371  ' A_t(m^2) V_bm(m^3/a)', &
1372  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1373  1113 format('----------------------------------------------------', &
1374  '---------------------------------------')
1375 
1376 #else
1377  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1378 #endif
1379 
1380 else if (forcing_flag == 3) then
1381 
1382  write(12,1122)
1383  write(12,1123)
1384 
1385 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1386 
1387  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1388  ' H_max(m) zs_max(m) V_g(m^3)', &
1389  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1390  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1391  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1392  1123 format('----------------------------------------------------', &
1393  '---------------------------------------')
1394 
1395 #elif OUTSER_V_A==2
1396 
1397  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1398  ' V(m^3) V_g(m^3) V_f(m^3)', &
1399  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1400  ' H_max(m) zs_max(m)', &
1401  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1402  ' A_t(m^2) V_bm(m^3/a)', &
1403  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1404  1123 format('----------------------------------------------------', &
1405  '---------------------------------------')
1406 
1407 #else
1408  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1409 #endif
1410 
1411 end if
1412 
1413 ! ------ Time-series file for the ice-thickness field
1414 
1415 #if OUTSER==2
1416 
1417 open(13, iostat=ios, &
1418  file=outpath//'/'//trim(runname)//'.thk', &
1419  status='new')
1420 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1421 
1422 if (forcing_flag == 1) then
1423 
1424  write(13,1104)
1425  write(13,1105)
1426 
1427  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1428  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1429  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1430  1105 format('----------------------------------------------------')
1431 
1432 else if (forcing_flag == 2) then
1433 
1434  write(13,1114)
1435  write(13,1115)
1436 
1437  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1438  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1439  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1440  1115 format('----------------------------------------------------')
1441 
1442 else if (forcing_flag == 3) then
1443 
1444  write(13,1124)
1445  write(13,1125)
1446 
1447  1124 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1448  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1449  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1450  1125 format('----------------------------------------------------')
1451 
1452 end if
1453 
1454 #endif
1455 
1456 ! ------ Time-series file for deep boreholes
1457 
1458 #if OUTSER==3
1459 
1460 n_core = 6 ! GRIP, GISP2, Dye3, Camp Century (CC), NorthGRIP (NGRIP), NEEM
1461 
1462 allocate(lambda_core(n_core), phi_core(n_core), &
1463  x_core(n_core), y_core(n_core))
1464 
1465 phi_core(1) = 72.58722_dp *pi_180 ! Geographical position of GRIP,
1466 lambda_core(1) = -37.64222_dp *pi_180 ! conversion deg -> rad
1467 call num_coord(lambda_core(1), phi_core(1), x_core(1), y_core(1))
1468 
1469 phi_core(2) = 72.58833_dp *pi_180 ! Geographical position of GISP2
1470 lambda_core(2) = -38.45750_dp *pi_180 ! conversion deg -> rad
1471 call num_coord(lambda_core(2), phi_core(2), x_core(2), y_core(2))
1472 
1473 phi_core(3) = 65.15139_dp *pi_180 ! Geographical position of Dye3,
1474 lambda_core(3) = -43.81722_dp *pi_180 ! conversion deg -> rad
1475 call num_coord(lambda_core(3), phi_core(3), x_core(3), y_core(3))
1476 
1477 phi_core(4) = 77.17970_dp *pi_180 ! Geographical position of CC,
1478 lambda_core(4) = -61.10975_dp *pi_180 ! conversion deg -> rad
1479 call num_coord(lambda_core(4), phi_core(4), x_core(4), y_core(4))
1480 
1481 phi_core(5) = 75.09694_dp *pi_180 ! Geographical position of NGRIP,
1482 lambda_core(5) = -42.31956_dp *pi_180 ! conversion deg -> rad
1483 call num_coord(lambda_core(5), phi_core(5), x_core(5), y_core(5))
1484 
1485 phi_core(6) = 77.5_dp *pi_180 ! Geographical position of NEEM,
1486 lambda_core(6) = -50.9_dp *pi_180 ! conversion deg -> rad
1487 call num_coord(lambda_core(6), phi_core(6), x_core(6), y_core(6))
1488 
1489 open(14, iostat=ios, &
1490  file=outpath//'/'//trim(runname)//'.core', &
1491  status='new')
1492 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
1493 
1494 if (forcing_flag == 1) then
1495 
1496  write(14,1106)
1497  write(14,1107)
1498 
1499  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1500  ' H_GR(m) H_G2(m) H_D3(m)', &
1501  ' H_CC(m) H_NG(m) H_NE(m)',/, &
1502  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
1503  ' v_CC(m/a) v_NG(m/a) v_NE(m/a)',/, &
1504  ' T_GR(C) T_G2(C) T_D3(C)', &
1505  ' T_CC(C) T_NG(C) T_NE(C)',/, &
1506  ' Rb_GR(W/m2) Rb_G2(W/m2) Rb_D3(W/m2)', &
1507  ' Rb_CC(W/m2) Rb_NG(W/m2) Rb_NE(W/m2)')
1508  1107 format('----------------------------------------------------', &
1509  '---------------------------------------')
1510 
1511 else if (forcing_flag == 2) then
1512 
1513  write(14,1116)
1514  write(14,1117)
1515 
1516  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
1517  ' H_GR(m) H_G2(m) H_D3(m)', &
1518  ' H_CC(m) H_NG(m) H_NE(m)',/, &
1519  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
1520  ' v_CC(m/a) v_NG(m/a) v_NE(m/a)',/, &
1521  ' T_GR(C) T_G2(C) T_D3(C)', &
1522  ' T_CC(C) T_NG(C) T_NE(C)',/, &
1523  ' Rb_GR(W/m2) Rb_G2(W/m2) Rb_D3(W/m2)', &
1524  ' Rb_CC(W/m2) Rb_NG(W/m2) Rb_NE(W/m2)')
1525  1117 format('----------------------------------------------------', &
1526  '---------------------------------------')
1527 
1528 else if (forcing_flag == 3) then
1529 
1530  write(14,1126)
1531  write(14,1127)
1532 
1533  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1534  ' H_GR(m) H_G2(m) H_D3(m)', &
1535  ' H_CC(m) H_NG(m) H_NE(m)',/, &
1536  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
1537  ' v_CC(m/a) v_NG(m/a) v_NE(m/a)',/, &
1538  ' T_GR(C) T_G2(C) T_D3(C)', &
1539  ' T_CC(C) T_NG(C) T_NE(C)',/, &
1540  ' Rb_GR(W/m2) Rb_G2(W/m2) Rb_D3(W/m2)', &
1541  ' Rb_CC(W/m2) Rb_NG(W/m2) Rb_NE(W/m2)')
1542  1127 format('----------------------------------------------------', &
1543  '---------------------------------------')
1544 
1545 end if
1546 
1547 #endif
1548 
1549 !-------- Output of the initial state --------
1550 
1551 #if OUTPUT==1
1552 
1553 #if ERGDAT==0
1554  flag_3d_output = .false.
1555 #elif ERGDAT==1
1556  flag_3d_output = .true.
1557 #else
1558  stop ' sico_init: ERGDAT must be either 0 or 1!'
1559 #endif
1560 
1561 #if NETCDF==1
1562  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1563  flag_3d_output, ndat2d, ndat3d)
1564 #elif NETCDF==2
1565  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1566  flag_3d_output, ndat2d, ndat3d)
1567 #else
1568  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1569 #endif
1570 
1571 #elif OUTPUT==2
1572 
1573 if (time_output(1) <= time_init+eps) then
1574 
1575 #if ERGDAT==0
1576  flag_3d_output = .false.
1577 #elif ERGDAT==1
1578  flag_3d_output = .true.
1579 #else
1580  stop ' sico_init: ERGDAT must be either 0 or 1!'
1581 #endif
1582 
1583 #if NETCDF==1
1584  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1585  flag_3d_output, ndat2d, ndat3d)
1586 #elif NETCDF==2
1587  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1588  flag_3d_output, ndat2d, ndat3d)
1589 #else
1590  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1591 #endif
1592 
1593 end if
1594 
1595 #elif OUTPUT==3
1596 
1597  flag_3d_output = .false.
1598 
1599 #if NETCDF==1
1600  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1601  flag_3d_output, ndat2d, ndat3d)
1602 #elif NETCDF==2
1603  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1604  flag_3d_output, ndat2d, ndat3d)
1605 #else
1606  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1607 #endif
1608 
1609 if (time_output(1) <= time_init+eps) then
1610 
1611  flag_3d_output = .true.
1612 
1613 #if NETCDF==1
1614  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1615  flag_3d_output, ndat2d, ndat3d)
1616 #elif NETCDF==2
1617  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1618  flag_3d_output, ndat2d, ndat3d)
1619 #else
1620  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1621 #endif
1622 
1623 end if
1624 
1625 #else
1626  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1627 #endif
1628 
1629 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1630 #if OUTSER==2
1631 call output3(time_init, delta_ts, glac_index, z_sl)
1632 #endif
1633 #if OUTSER==3
1634 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1635 #endif
1636 
1637 end subroutine sico_init
1638 !
Generated on Thu Sep 12 2013 09:44:59 for SICOPOLIS V3.1 by   doxygen 1.8.2