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