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 
48 implicit none
49 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
50 integer(i4b) :: ios
51 integer(i4b) :: ndat2d, ndat3d
52 integer(i4b) :: n_output
53 integer(i4b) :: ierr
54 integer(i4b) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1)), nn(0:jmax,0:imax)
55 real(dp) :: delta_ts, glac_index
56 real(dp) :: sum_accum, mean_accum, mean_accum_inv
57 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
58 real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
59 real(dp) :: time_init0, time_end0, time_output0(100)
60 real(dp) :: time, time_init, time_end, time_output(100)
61 real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
62 real(dp) :: z_sl, dzsl_dtau, z_mar
63 real(dp) :: d_dummy
64 character (len=100) :: runname, anfdatname
65 character (len=256) :: shell_command
66 character :: ch_dummy
67 logical :: flag_3d_output
68 
69 !-------- Initialisation of the Library of Iterative Solvers Lis,
70 ! if required --------
71 
72 #if (CALCZS==4 || MARGIN==3)
73 
74 #include "lisf.h"
75 call lis_initialize(ierr)
76 
77 #endif
78 
79 !-------- Read physical parameters --------
80 
81 call phys_para
82 
83 !-------- Compatibility check of the SICOPOLIS version with the header file
84 
85 if ( trim(version) /= trim(sico_version) ) &
86  stop ' sico_init: SICOPOLIS version not compatible with header file!'
87 
88 !-------- Compatibility check of the horizontal resolution with the
89 ! number of grid points --------
90 
91 #if (GRID==0 || GRID==1)
92 
93 if (abs(dx-40.0_dp) < eps) then
94 
95  if ((imax /= 150).or.(jmax /= 70)) then
96  stop ' sico_init: IMAX and/or JMAX wrong!'
97  end if
98 
99 else if (abs(dx-20.0_dp) < eps) then
100 
101  if ((imax /= 300).or.(jmax /= 140)) then
102  stop ' sico_init: IMAX and/or JMAX wrong!'
103  end if
104 
105 else if (abs(dx-10.0_dp) < eps) then
106 
107  if ((imax /= 600).or.(jmax /= 280)) then
108  stop ' sico_init: IMAX and/or JMAX wrong!'
109  end if
110 
111 else
112  stop ' sico_init: DX wrong!'
113 end if
114 
115 #elif GRID==2
116 
117  stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
118 
119 #endif
120 
121 !-------- Compatibility check of surface-temperature and precipitation
122 ! determination by interpolation between present and LGM values
123 ! with a glacial index --------
124 
125 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
126 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
127 #endif
128 
129 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
130 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
131 #endif
132 
133 !-------- Compatibility check of discretization schemes for the horizontal and
134 ! vertical advection terms in the temperature and age equations --------
135 
136 #if ADV_HOR==1
137 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
138 #endif
139 
140 #if ((ADV_HOR == 3) && (ADV_VERT != 3))
141 stop ' sico_init: Options ADV_HOR==3 and ADV_VERT==3 must be used together!'
142 #endif
143 
144 #if ((ADV_HOR != 3) && (ADV_VERT == 3))
145 stop ' sico_init: Options ADV_HOR==3 and ADV_VERT==3 must be used together!'
146 #endif
147 
148 !-------- Setting of forcing flag --------
149 
150 #if (TSURFACE <= 4)
151 
152 forcing_flag = 1 ! forcing by delta_ts
153 
154 #elif (TSURFACE == 5)
155 
156 forcing_flag = 2 ! forcing by glac_index
157 
158 #endif
159 
160 !-------- Initialization of numerical time steps --------
161 
162 dtime0 = dtime0
163 dtime_temp0 = dtime_temp0
164 #if REBOUND==2
165 dtime_wss0 = dtime_wss0
166 #endif
167 
168 !-------- Further initializations --------
169 
170 dzeta_c = 1.0_dp/real(kcmax,dp)
171 dzeta_t = 1.0_dp/real(ktmax,dp)
172 dzeta_r = 1.0_dp/real(krmax,dp)
173 
174 ndat2d = 1
175 ndat3d = 1
176 
177 !-------- Construction of a vector (with index n) from a 2-d array
178 ! (with indices i, j) by diagonal numbering --------
179 
180 n=1
181 
182 do m=0, imax+jmax
183  do i=m, 0, -1
184  j = m-i
185  if ((i.le.imax).and.(j.le.jmax)) then
186  ii(n) = i
187  jj(n) = j
188  nn(j,i) = n
189  n=n+1
190  end if
191  end do
192 end do
193 
194 !-------- Specification of current simulation --------
195 
196 runname = runname
197 anfdatname = anfdatname
198 
199 #if defined(YEAR_ZERO)
200 year_zero = year_zero
201 #else
202 year_zero = 2000.0_dp ! default value 2000 CE
203 #endif
204 
205 time_init0 = time_init0
206 time_end0 = time_end0
207 dtime_ser0 = dtime_ser0
208 
209 #if ( OUTPUT==1 || OUTPUT==3 )
210 dtime_out0 = dtime_out0
211 #endif
212 
213 #if ( OUTPUT==2 || OUTPUT==3 )
214 n_output = n_output
215 time_output0( 1) = time_out0_01
216 time_output0( 2) = time_out0_02
217 time_output0( 3) = time_out0_03
218 time_output0( 4) = time_out0_04
219 time_output0( 5) = time_out0_05
220 time_output0( 6) = time_out0_06
221 time_output0( 7) = time_out0_07
222 time_output0( 8) = time_out0_08
223 time_output0( 9) = time_out0_09
224 time_output0(10) = time_out0_10
225 time_output0(11) = time_out0_11
226 time_output0(12) = time_out0_12
227 time_output0(13) = time_out0_13
228 time_output0(14) = time_out0_14
229 time_output0(15) = time_out0_15
230 time_output0(16) = time_out0_16
231 time_output0(17) = time_out0_17
232 time_output0(18) = time_out0_18
233 time_output0(19) = time_out0_19
234 time_output0(20) = time_out0_20
235 #endif
236 
237 !-------- Write log file --------
238 
239 shell_command = 'if [ ! -d'
240 shell_command = trim(shell_command)//' '//outpath
241 shell_command = trim(shell_command)//' '//'] ; then mkdir'
242 shell_command = trim(shell_command)//' '//outpath
243 shell_command = trim(shell_command)//' '//'; fi'
244 call system(trim(shell_command))
245  ! Check whether directory OUTPATH exists. If not, it is created.
246 
247 open(10, iostat=ios, &
248  file=outpath//'/'//trim(runname)//'.log', &
249  status='new')
250 if (ios /= 0) &
251  stop 'Error when opening the log file!'
252 
253 write(10,1000) 'Computational domain:'
254 #if defined(ANT)
255 write(10,1000) 'Antarctica'
256 #elif defined(GRL)
257 write(10,1000) 'Greenland'
258 #elif defined(SCAND)
259 write(10,1000) 'Scandinavia and Eurasia'
260 #elif defined(NHEM)
261 write(10,1000) 'Northern hemisphere'
262 #elif defined(EMTP2SGE)
263 write(10,1000) 'EISMINT Phase 2 Simplified Geometry Experiment'
264 #elif defined(NMARS)
265 write(10,1000) 'North polar cap of Mars'
266 #elif defined(SMARS)
267 write(10,1000) 'South polar cap of Mars'
268 #else
269 stop ' sico_init: Subroutine sico_init: No valid domain specified!'
270 #endif
271 write(10,1000) ' '
272 
273 write(10,1001) 'imax =', imax
274 write(10,1001) 'jmax =', jmax
275 write(10,1001) 'kcmax =', kcmax
276 write(10,1001) 'ktmax =', ktmax
277 write(10,1001) 'krmax =', krmax
278 write(10,1000) ' '
279 
280 write(10,1002) 'a =', deform
281 write(10,1000) ' '
282 
283 #if (GRID==0 || GRID==1)
284 write(10,1002) 'x0 =', x0
285 write(10,1002) 'y0 =', y0
286 write(10,1002) 'dx =', dx
287 #elif GRID==2
288 stop ' sico_init: GRID==2 not allowed for this application!'
289 #endif
290 write(10,1000) ' '
291 
292 write(10,1002) 'year_zero =', year_zero
293 write(10,1002) 'time_init =', time_init0
294 write(10,1002) 'time_end =', time_end0
295 write(10,1002) 'dtime =', dtime0
296 write(10,1002) 'dtime_temp =', dtime_temp0
297 #if REBOUND==2
298 write(10,1002) 'dtime_wss =', dtime_wss0
299 #endif
300 #if ( OUTPUT==1 || OUTPUT==3 )
301 write(10,1002) 'dtime_out =', dtime_out0
302 #endif
303 write(10,1002) 'dtime_ser =', dtime_ser0
304 write(10,1000) ' '
305 
306 write(10,1000) 'zs_present file = '//zs_present_file
307 #if ANF_DAT==1
308 write(10,1000) 'zl_present file = '//zl_present_file
309 #endif
310 write(10,1000) 'zl0 file = '//zl0_file
311 write(10,1000) 'mask_present file = '//mask_present_file
312 write(10,1000) ' '
313 
314 write(10,1000) 'Physical-parameter file = '//phys_para_file
315 write(10,1000) ' '
316 
317 #if (CALCZS==3 || CALCZS==4)
318 write(10,1002) 'ovi_weight =', ovi_weight
319 #if CALCZS==3
320 write(10,1002) 'omega_sor =', omega_sor
321 #endif
322 write(10,1000) ' '
323 #endif
324 
325 write(10,1000) 'temp_mm_present file = '//temp_mm_present_file
326 #if TSURFACE==1
327 write(10,1002) 'delta_ts0 =', delta_ts0
328 #elif TSURFACE==3
329 write(10,1002) 'sine_amplit =', sine_amplit
330 write(10,1002) 'sine_period =', sine_period
331 #elif TSURFACE==4
332 write(10,1000) 'GRIP file = '//grip_temp_file
333 write(10,1002) 'grip_temp_fact =', grip_temp_fact
334 #elif TSURFACE==5
335 write(10,1000) 'Glacial-index file = '//glac_ind_file
336 write(10,1000) 'temp_mm_anom file = '//temp_mm_anom_file
337 write(10,1002) 'temp_mm_anom fact = ', temp_mm_anom_fact
338 #endif
339 
340 write(10,1000) 'precip_mm_present file = '//precip_mm_present_file
341 #if ACCSURFACE==1
342 write(10,1002) 'accfact =', accfact
343 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
344 write(10,1002) 'gamma_s =', gamma_s
345 #elif ( ACCSURFACE==5 )
346 write(10,1000) 'precip_mm_anom file = '//precip_mm_anom_file
347 write(10,1002) 'precip_mm_anom fact = ', precip_mm_anom_fact
348 #endif
349 #if (ACCSURFACE <= 3)
350 write(10,1001) 'ELEV_DESERT =', elev_desert
351 #if (ELEV_DESERT == 1)
352 write(10,1002) 'gamma_p =', gamma_p
353 write(10,1002) 'zs_thresh =', zs_thresh
354 #endif
355 #endif
356 
357 #if ABLSURFACE==3
358 write(10,1002) 'lambda_lti =', lambda_lti
359 write(10,1002) 'temp_lti =', temp_lti
360 #endif
361 
362 #if SEA_LEVEL==1
363 write(10,1002) 'z_sl0 =', z_sl0
364 #elif SEA_LEVEL==3
365 write(10,1000) 'sea-level file = '//sea_level_file
366 #endif
367 write(10,1000) ' '
368 
369 #if MARGIN==2
370 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
371 write(10,1002) 'z_mar =', z_mar
372 write(10,1000) ' '
373 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
374  || marine_ice_calving==6 || marine_ice_calving==7 )
375 write(10,1002) 'fact_z_mar =', fact_z_mar
376 write(10,1000) ' '
377 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
378 write(10,1002) 'calv_uw_coeff =', calv_uw_coeff
379 write(10,1002) 'r1_calv_uw =', r1_calv_uw
380 write(10,1002) 'r2_calv_uw =', r2_calv_uw
381 write(10,1000) ' '
382 #endif
383 #elif MARGIN==3
384 #if ICE_SHELF_CALVING==2
385 write(10,1002) 'H_calv =', h_calv
386 write(10,1000) ' '
387 #endif
388 #endif
389 
390 write(10,1002) 'c_slide =', c_slide
391 write(10,1002) 'gamma_slide =', gamma_slide
392 write(10,1001) 'p_weert =', p_weert
393 write(10,1001) 'q_weert =', q_weert
394 #if SLIDE_LAW==2
395 write(10,1002) 'red_pres_limit_fact =', red_pres_limit_fact
396 #endif
397 write(10,1000) ' '
398 
399 #if Q_GEO_MOD==1
400 write(10,1002) 'q_geo =', q_geo
401 #elif Q_GEO_MOD==2
402 write(10,1000) 'q_geo file = '//q_geo_file
403 #endif
404 write(10,1000) ' '
405 
406 #if defined(MARINE_ICE_BASAL_MELTING)
407 write(10,1001) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
408 #if ( MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3 )
409 write(10,1002) 'qbm_marine =', qbm_marine
410 #endif
411 write(10,1000) ' '
412 #endif
413 
414 #if REBOUND==1
415 write(10,1002) 'frac_llra =', frac_llra
416 write(10,1000) ' '
417 #endif
418 
419 #if FLOW_LAW==2
420 write(10,1002) 'gr_size =', gr_size
421 write(10,1000) ' '
422 #endif
423 #if FIN_VISC==2
424 write(10,1002) 'sigma_res =', sigma_res
425 write(10,1000) ' '
426 #endif
427 
428 write(10,1001) 'ENHMOD =', enhmod
429 write(10,1002) 'enh_fact =', enh_fact
430 #if ( ENHMOD==2 || ENHMOD==3 )
431 write(10,1002) 'enh_intg =', enh_intg
432 #endif
433 #if ENHMOD==2
434 write(10,1002) 'age_trans =', age_trans_0
435 #endif
436 #if ENHMOD==3
437 write(10,1002) 'date_trans1 =', date_trans1_0
438 write(10,1002) 'date_trans2 =', date_trans2_0
439 write(10,1002) 'date_trans3 =', date_trans3_0
440 #endif
441 #if MARGIN==3
442 write(10,1002) 'enh_shelf =', enh_shelf
443 #endif
444 write(10,1000) ' '
445 
446 write(10,1002) 'numdiff_H_t =', numdiff_h_t
447 write(10,1002) 'tau_cts =', tau_cts
448 write(10,1002) 'H_min =', h_min
449 write(10,1002) 'vh_max =', vh_max
450 write(10,1002) 'hd_min =', hd_min
451 write(10,1002) 'hd_max =', hd_max
452 #if defined(QBM_MIN)
453 write(10,1002) 'qbm_min =', qbm_min
454 #elif defined(QB_MIN) /* obsolete */
455 write(10,1002) 'qb_min =', qb_min
456 #endif
457 #if defined(QBM_MAX)
458 write(10,1002) 'qbm_max =', qbm_max
459 #elif defined(QB_MAX) /* obsolete */
460 write(10,1002) 'qb_max =', qb_max
461 #endif
462 write(10,1002) 'age_min =', age_min
463 write(10,1002) 'age_max =', age_max
464 write(10,1002) 'mean_accum =', mean_accum
465 #if ADV_VERT==1
466 write(10,1002) 'age_diff =', agediff
467 write(10,1002) 'm_age =', m_age
468 #endif
469 write(10,1000) ' '
470 
471 write(10,1001) 'GRID =', grid
472 write(10,1001) 'CALCMOD =', calcmod
473 write(10,1001) 'FLOW_LAW =', flow_law
474 write(10,1001) 'FIN_VISC =', fin_visc
475 write(10,1001) 'MARGIN =', margin
476 #if MARGIN==2
477 write(10,1001) 'MARINE_ICE_FORMATION =', marine_ice_formation
478 write(10,1001) 'MARINE_ICE_CALVING =', marine_ice_calving
479 #elif MARGIN==3
480 write(10,1001) 'ICE_SHELF_CALVING =', ice_shelf_calving
481 #endif
482 write(10,1001) 'ANF_DAT =', anf_dat
483 write(10,1001) 'REBOUND =', rebound
484 write(10,1001) 'Q_LITHO =', q_litho
485 write(10,1001) 'ZS_EVOL =', zs_evol
486 write(10,1001) 'CALCZS =', calczs
487 write(10,1001) 'ADV_HOR =', adv_hor
488 write(10,1001) 'ADV_VERT =', adv_vert
489 write(10,1001) 'TOPOGRAD =', topograd
490 write(10,1001) 'TSURFACE =', tsurface
491 write(10,1001) 'ACCSURFACE =', accsurface
492 #if ACCSURFACE==5
493 write(10,1001) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
494 #endif
495 write(10,1001) 'SOLID_PRECIP =', solid_precip
496 write(10,1001) 'ABLSURFACE =', ablsurface
497 write(10,1001) 'SEA_LEVEL =', sea_level
498 write(10,1001) 'SLIDE_LAW =', slide_law
499 write(10,1001) 'Q_GEO_MOD =', q_geo_mod
500 write(10,1000) ' '
501 
502 #if defined(OUT_TIMES)
503 write(10,1001) 'OUT_TIMES =', out_times
504 #endif
505 write(10,1001) 'OUTPUT =', output
506 write(10,1001) 'OUTSER =', outser
507 #if defined(OUTSER_V_A)
508 write(10,1001) 'OUTSER_V_A =', outser_v_a
509 #endif
510 #if ( OUTPUT==1 || OUTPUT==2 )
511 write(10,1001) 'ERGDAT =', ergdat
512 #endif
513 write(10,1001) 'NETCDF =', netcdf
514 write(10,1000) ' '
515 #if ( OUTPUT==2 || OUTPUT==3 )
516 write(10,1001) 'n_output =', n_output
517 do n=1, n_output
518  write(10,1002) 'time_output =', time_output0(n)
519 end do
520 write(10,1000) ' '
521 #endif
522 
523 write(10,1000) 'Program version and date: '//version//' / '//date
524 
525  1000 format(a)
526  1001 format(a,i4)
527  1002 format(a,es12.4)
528 
529 close(10, status='keep')
530 
531 !-------- Conversion of time quantities --------
532 
533 year_zero = year_zero*year_sec ! a --> s
534 time_init = time_init0*year_sec ! a --> s
535 time_end = time_end0*year_sec ! a --> s
536 dtime = dtime0*year_sec ! a --> s
537 dtime_temp = dtime_temp0*year_sec ! a --> s
538 #if REBOUND==2
539 dtime_wss = dtime_wss0*year_sec ! a --> s
540 #endif
541 dtime_ser = dtime_ser0*year_sec ! a --> s
542 #if ( OUTPUT==1 || OUTPUT==3 )
543 dtime_out = dtime_out0*year_sec ! a --> s
544 #endif
545 #if ( OUTPUT==2 || OUTPUT==3 )
546 do n=1, n_output
547  time_output(n) = time_output0(n)*year_sec ! a --> s
548 end do
549 #endif
550 
551 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
552  stop ' sico_init: dtime_temp must be a multiple of dtime!'
553 
554 #if REBOUND==2
555 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
556  stop ' sico_init: dtime_wss must be a multiple of dtime!'
557 #endif
558 
559 !-------- Reading of measurements for present monthly-mean precipitation --------
560 
561 #if (GRID==0 || GRID==1)
562 
563 open(21, iostat=ios, &
564  file=inpath//'/scand/'//precip_mm_present_file, &
565  recl=8192, status='old')
566 
567 #elif GRID==2
568 
569 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
570 
571 #endif
572 
573 if (ios /= 0) stop ' sico_init: Error when opening the precipitation file!'
574 
575 do m=1, 6; read(21,'(a)') ch_dummy; end do
576 
577 do n=1, 12 ! month counter
578  do m=1, 3; read(21,'(a)') ch_dummy; end do
579  do j=jmax, 0, -1
580  read(21,*) (precip_present(j,i,n), i=0,imax)
581  end do
582 end do
583 
584 close(21, status='keep')
585 
586 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
587 
588 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
589  ! mm/a water equiv. -> m/s ice equiv.
590 
591 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
592 
593 #if ACCSURFACE==5
594 
595 #if (GRID==0 || GRID==1)
596 
597 open(21, iostat=ios, &
598  file=inpath//'/scand/'//precip_mm_anom_file, &
599  recl=8192, status='old')
600 
601 #endif
602 
603 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
604 
605 do m=1, 6; read(21,'(a)') ch_dummy; end do
606 
607 do n=1, 12 ! month counter
608  do m=1, 3; read(21,'(a)') ch_dummy; end do
609  do j=jmax, 0, -1
610  read(21,*) (precip_lgm_anom(j,i,n), i=0,imax)
611  end do
612 end do
613 
614 close(21, status='keep')
615 
616 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
617 
618 do i=0, imax
619 do j=0, jmax
620 
621 #if (PRECIP_ANOM_INTERPOL==1)
622  do n=1, 12 ! month counter
623  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
624  end do
625 #elif (PRECIP_ANOM_INTERPOL==2)
626  do n=1, 12 ! month counter
627  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
628  end do
629 #else
630  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
631 #endif
632 
633 end do
634 end do
635 
636 #endif
637 
638 !-------- Mean accumulation --------
639 
640 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
641  ! mm/a water equiv. -> m/s ice equiv.
642 
643 mean_accum_inv = 1.0_dp/mean_accum
644 
645 !-------- Reading of present topography mask --------
646 
647 #if (GRID==0 || GRID==1)
648 
649 open(24, iostat=ios, &
650  file=inpath//'/scand/'//mask_present_file, &
651  status='old')
652 
653 #elif GRID==2
654 
655 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
656 
657 #endif
658 
659 if (ios /= 0) &
660  stop ' sico_init: Error when opening the mask file!'
661 
662 do m=1, 6; read(24,'(a)') ch_dummy; end do
663 
664 do j=jmax, 0, -1
665  read(24,2300) (maske_help(j,i), i=0,imax)
666 end do
667 
668 close(24, status='keep')
669 
670 2300 format(imax(i1),i1)
671 
672 !-------- Reading of data for present monthly-mean surface temperature --------
673 
674 #if (GRID==0 || GRID==1)
675 
676 open(21, iostat=ios, &
677  file=inpath//'/scand/'//temp_mm_present_file, &
678  recl=8192, status='old')
679 
680 #elif GRID==2
681 
682 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
683 
684 #endif
685 
686 if (ios /= 0) stop ' sico_init: Error when opening the temperature file!'
687 
688 do m=1, 6; read(21,'(a)') ch_dummy; end do
689 
690 do n=1, 12 ! month counter
691  do m=1, 3; read(21,'(a)') ch_dummy; end do
692  do j=jmax, 0, -1
693  read(21,*) (temp_mm_present(j,i,n), i=0,imax)
694  end do
695 end do
696 
697 close(21, status='keep')
698 
699 !-------- Reading of LGM monthly-mean surface-temperature anomalies --------
700 
701 #if TSURFACE==5
702 
703 #if (GRID==0 || GRID==1)
704 
705 open(21, iostat=ios, &
706  file=inpath//'/scand/'//temp_mm_anom_file, &
707  recl=8192, status='old')
708 
709 #endif
710 
711 if (ios /= 0) stop ' sico_init: Error when opening the temperature anomaly file!'
712 
713 do m=1, 6; read(21,'(a)') ch_dummy; end do
714 
715 do n=1, 12 ! month counter
716  do m=1, 3; read(21,'(a)') ch_dummy; end do
717  do j=jmax, 0, -1
718  read(21,*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
719  end do
720 end do
721 
722 close(21, status='keep')
723 
724 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
725 
726 #endif
727 
728 !-------- Present reference elevation
729 ! (for precipitation and surface-temperature data) --------
730 
731 #if (GRID==0 || GRID==1)
732 
733 open(21, iostat=ios, &
734  file=inpath//'/scand/'//zs_present_file, &
735  recl=8192, status='old')
736 
737 #elif GRID==2
738 
739 stop ' sico_init: GRID==2 not allowed for the Scandinavia application!'
740 
741 #endif
742 
743 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
744 
745 do m=1, 6; read(21,'(a)') ch_dummy; end do
746 
747 do j=jmax, 0, -1
748  read(21,*) (zs_ref(j,i), i=0,imax)
749 end do
750 
751 close(21, status='keep')
752 
753 ! ------ Reset bathymetry data (sea floor elevation) to the
754 ! sea surface
755 
756 do i=0, imax
757 do j=0, jmax
758  if (maske_help(j,i).ge.2) zs_ref(j,i) = 0.0_dp
759 end do
760 end do
761 
762 !-------- Read data for delta_ts --------
763 
764 #if TSURFACE==4
765 
766 open(21, iostat=ios, &
767  file=inpath//'/general/'//grip_temp_file, &
768  status='old')
769 if (ios /= 0) &
770  stop ' sico_init: Error when opening the data file for delta_ts!'
771 
772 read(21,*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
773 
774 if (ch_dummy /= '#') then
775  write(6,*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
776  write(6,*) ' not defined in data file for delta_ts!'
777  stop
778 end if
779 
780 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
781 
782 allocate(griptemp(0:ndata_grip))
783 
784 do n=0, ndata_grip
785  read(21,*) d_dummy, griptemp(n)
786 end do
787 
788 close(21, status='keep')
789 
790 #endif
791 
792 !-------- Read data for the glacial index --------
793 
794 #if TSURFACE==5
795 
796 open(21, iostat=ios, &
797  file=inpath//'/general/'//glac_ind_file, status='old')
798 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
799 
800 read(21,*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
801 
802 if (ch_dummy /= '#') then
803  write(6,*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
804  write(6,*) ' not defined in glacial-index file!'
805  stop
806 end if
807 
808 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
809 
810 allocate(glacial_index(0:ndata_gi))
811 
812 do n=0, ndata_gi
813  read(21,*) d_dummy, glacial_index(n)
814 end do
815 
816 close(21, status='keep')
817 
818 #endif
819 
820 !-------- Read data for z_sl --------
821 
822 #if SEA_LEVEL==3
823 
824 open(21, iostat=ios, &
825  file=inpath//'/general/'//sea_level_file, &
826  status='old')
827 if (ios /= 0) &
828  stop ' sico_init: Error when opening the data file for z_sl!'
829 
830 read(21,*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
831 
832 if (ch_dummy /= '#') then
833  write(6,*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
834  write(6,*) ' not defined in data file for z_sl!'
835  stop
836 end if
837 
838 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
839 
840 allocate(specmap_zsl(0:ndata_specmap))
841 
842 do n=0, ndata_specmap
843  read(21,*) d_dummy, specmap_zsl(n)
844 end do
845 
846 close(21, status='keep')
847 
848 #endif
849 
850 !-------- Determination of the geothermal heat flux --------
851 
852 #if Q_GEO_MOD==1
853 
854 ! ------ Constant value
855 
856 do i=0, imax
857 do j=0, jmax
858  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
859 end do
860 end do
861 
862 #elif Q_GEO_MOD==2
863 
864 ! ------ Read data from file
865 
866 open(21, iostat=ios, &
867  file=inpath//'/scand/'//q_geo_file, &
868  recl=8192, status='old')
869 
870 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
871 
872 do m=1, 6; read(21,'(a)') ch_dummy; end do
873 
874 do j=jmax, 0, -1
875  read(21,*) (q_geo(j,i), i=0,imax)
876 end do
877 
878 close(21, status='keep')
879 
880 do i=0, imax
881 do j=0, jmax
882  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
883 end do
884 end do
885 
886 #endif
887 
888 !-------- Reading of tabulated kei function--------
889 
890 #if (REBOUND==0 || REBOUND==1)
891 
892 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
893  ! dummy values
894 #elif REBOUND==2
895 
896 call read_kei
897 
898 #endif
899 
900 !-------- Definition of initial values --------
901 
902 ! ------ Present topography
903 
904 #if ANF_DAT==1
905 
906 stop ' sico_init: ANF_DAT==1 not allowed for scand application!'
907 
908 ! ------ Ice-free, relaxed bedrock
909 
910 #elif ANF_DAT==2
911 
912 call topography2(dxi, deta)
913 
914 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
915 
916 call boundary(time_init, dtime, dxi, deta, &
917  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
918 
919 do i=0, imax
920 do j=0, jmax
921 
922  q_bm(j,i) = 0.0_dp
923  q_tld(j,i) = 0.0_dp
924  q_b_tot(j,i) = 0.0_dp
925 
926  p_b_w(j,i) = 0.0_dp
927  h_w(j,i) = 0.0_dp
928 
929  do kc=0, kcmax
930  temp_c(kc,j,i) = temp_s(j,i)
931  age_c(kc,j,i) = 15000.0_dp*year_sec
932  end do
933 
934  do kt=0, ktmax
935  omega_t(kt,j,i) = 0.0_dp
936  age_t(kt,j,i) = 15000.0_dp*year_sec
937  end do
938 
939  do kr=0, krmax
940  temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
941  *(1.0_dp-real(kr,dp)/real(krmax,dp))
942 ! (linear temperature distribution according to the
943 ! geothermal heat flux)
944  end do
945 
946 end do
947 end do
948 
949 call calc_enhance(time_init)
950 
951 ! ------ Read initial state from output of previous simulation
952 
953 #elif ANF_DAT==3
954 
955 #if NETCDF==1
956 call topography3(dxi, deta, z_sl, anfdatname)
957 #elif NETCDF==2
958 call topography3_nc(dxi, deta, z_sl, anfdatname)
959 #else
960 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
961 #endif
962 
963 call boundary(time_init, dtime, dxi, deta, &
964  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
965 
966 q_b_tot = q_bm + q_tld
967 
968 call calc_enhance(time_init)
969 
970 #endif
971 
972 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
973 
974 #if (GRID==0 || GRID==1)
975 
976 do ir=-imax, imax
977 do jr=-jmax, jmax
978  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
979  ! distortion due to stereographic projection not accounted for
980 end do
981 end do
982 
983 #endif
984 
985 !-------- General abbreviations --------
986 
987 ea = exp(deform)
988 
989 do kc=0, kcmax
990  zeta_c(kc) = kc*dzeta_c
991  eaz_c(kc) = exp(deform*zeta_c(kc))
992  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
993 end do
994 
995 do kt=0, ktmax
996  zeta_t(kt) = kt*dzeta_t
997 end do
998 
999 !-------- Initial velocities --------
1000 
1001 call calc_temp_melt
1002 
1003 call calc_vxy_b_sia(time, z_sl)
1004 call calc_vxy_sia(dzeta_c, dzeta_t)
1005 
1006 #if MARGIN==3
1007 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
1008 #endif
1009 
1010 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1011 
1012 #if MARGIN==3
1013 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1014 #endif
1015 
1016 !-------- Initialize time-series files --------
1017 
1018 ! ------ Time-series file for the ice sheet on the whole
1019 
1020 open(12, iostat=ios, &
1021  file=outpath//'/'//trim(runname)//'.ser', &
1022  status='new')
1023 if (ios /= 0) &
1024  stop ' sico_init: Error when opening the ser file!'
1025 
1026 if (forcing_flag == 1) then
1027 
1028  write(12,1102)
1029  write(12,1103)
1030 
1031 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1032 
1033  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1034  ' H_max(m) zs_max(m) V_g(m^3)', &
1035  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1036  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1037  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1038  1103 format('----------------------------------------------------', &
1039  '---------------------------------------')
1040 
1041 #elif OUTSER_V_A==2
1042 
1043  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1044  ' V(m^3) V_g(m^3) V_f(m^3)', &
1045  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1046  ' H_max(m) zs_max(m)', &
1047  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1048  ' A_t(m^2) V_bm(m^3/a)', &
1049  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1050  1103 format('----------------------------------------------------', &
1051  '---------------------------------------')
1052 
1053 #else
1054  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1055 #endif
1056 
1057 else if (forcing_flag == 2) then
1058 
1059  write(12,1112)
1060  write(12,1113)
1061 
1062 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1063 
1064  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1065  ' H_max(m) zs_max(m) V_g(m^3)', &
1066  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1067  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1068  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1069  1113 format('----------------------------------------------------', &
1070  '---------------------------------------')
1071 
1072 #elif OUTSER_V_A==2
1073 
1074  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1075  ' V(m^3) V_g(m^3) V_f(m^3)', &
1076  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1077  ' H_max(m) zs_max(m)', &
1078  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1079  ' A_t(m^2) V_bm(m^3/a)', &
1080  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1081  1113 format('----------------------------------------------------', &
1082  '---------------------------------------')
1083 
1084 #else
1085  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1086 #endif
1087 
1088 end if
1089 
1090 ! ------ Time-series file for the ice-thickness field
1091 
1092 #if OUTSER==2
1093 
1094 open(13, iostat=ios, &
1095  file=outpath//'/'//trim(runname)//'.thk', &
1096  status='new')
1097 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1098 
1099 if (forcing_flag == 1) then
1100 
1101  write(13,1104)
1102  write(13,1105)
1103 
1104  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1105  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1106  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1107  1105 format('----------------------------------------------------')
1108 
1109 else if (forcing_flag == 2) then
1110 
1111  write(13,1114)
1112  write(13,1115)
1113 
1114  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1115  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1116  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1117  1115 format('----------------------------------------------------')
1118 
1119 end if
1120 
1121 #endif
1122 
1123 !-------- Output of the initial state --------
1124 
1125 #if OUTPUT==1
1126 
1127 #if ERGDAT==0
1128  flag_3d_output = .false.
1129 #elif ERGDAT==1
1130  flag_3d_output = .true.
1131 #else
1132  stop ' sico_init: ERGDAT must be either 0 or 1!'
1133 #endif
1134 
1135 #if NETCDF==1
1136  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1137  flag_3d_output, ndat2d, ndat3d)
1138 #elif NETCDF==2
1139  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1140  flag_3d_output, ndat2d, ndat3d)
1141 #else
1142  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1143 #endif
1144 
1145 #elif OUTPUT==2
1146 
1147 if (time_output(1) <= time_init+eps) then
1148 
1149 #if ERGDAT==0
1150  flag_3d_output = .false.
1151 #elif ERGDAT==1
1152  flag_3d_output = .true.
1153 #else
1154  stop ' sico_init: ERGDAT must be either 0 or 1!'
1155 #endif
1156 
1157 #if NETCDF==1
1158  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1159  flag_3d_output, ndat2d, ndat3d)
1160 #elif NETCDF==2
1161  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1162  flag_3d_output, ndat2d, ndat3d)
1163 #else
1164  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1165 #endif
1166 
1167 end if
1168 
1169 #elif OUTPUT==3
1170 
1171  flag_3d_output = .false.
1172 
1173 #if NETCDF==1
1174  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1175  flag_3d_output, ndat2d, ndat3d)
1176 #elif NETCDF==2
1177  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1178  flag_3d_output, ndat2d, ndat3d)
1179 #else
1180  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1181 #endif
1182 
1183 if (time_output(1) <= time_init+eps) then
1184 
1185  flag_3d_output = .true.
1186 
1187 #if NETCDF==1
1188  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1189  flag_3d_output, ndat2d, ndat3d)
1190 #elif NETCDF==2
1191  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1192  flag_3d_output, ndat2d, ndat3d)
1193 #else
1194  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1195 #endif
1196 
1197 end if
1198 
1199 #else
1200  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
1201 #endif
1202 
1203 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1204 #if OUTSER==2
1205 call output3(time_init, delta_ts, glac_index, z_sl)
1206 #endif
1207 
1208 end subroutine sico_init
1209 !
Generated on Mon May 20 2013 10:18:36 for SICOPOLIS V3.0 by   doxygen 1.8.1