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