36 mean_accum, mean_accum_inv, &
37 dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38 time, time_init, time_end, time_output, &
39 dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40 z_sl, dzsl_dtau, z_mar, &
42 ndat2d, ndat3d, n_output, &
50 integer(i4b),
intent(out) :: ii((imax+1)*(jmax+1)), &
51 jj((imax+1)*(jmax+1)), &
53 integer(i4b),
intent(out) :: ndat2d, ndat3d
54 integer(i4b),
intent(out) :: n_output
55 real(dp),
intent(out) :: delta_ts, glac_index
56 real(dp),
intent(out) :: mean_accum, mean_accum_inv
57 real(dp),
intent(out) :: dtime, dtime_temp, dtime_wss, &
59 real(dp),
intent(out) :: time, time_init, time_end, time_output(100)
60 real(dp),
intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
61 real(dp),
intent(out) :: z_sl, dzsl_dtau, z_mar
62 character(len=100),
intent(out) :: runname
64 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
67 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
68 real(dp) :: time_init0, time_end0, time_output0(100)
70 character (len=100) :: anfdatname
71 character (len=256) :: shell_command
73 logical :: flag_3d_output
75 character(len=64),
parameter :: fmt1 =
'(a)', &
79 character(len= 8) :: ch_imax
80 character(len=128) :: fmt4
82 write(ch_imax, fmt=
'(i8)') imax
83 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)'
103 #if (CALCZS==4 || MARGIN==3)
106 call lis_initialize(ierr)
116 if ( trim(version) /= trim(sico_version) ) &
117 stop
' sico_init: SICOPOLIS version not compatible with header file!'
122 #if (GRID==0 || GRID==1)
124 if (abs(dx-80.0_dp) < eps)
then
126 if ((imax /= 156).or.(jmax /= 156))
then
127 stop
' sico_init: IMAX and/or JMAX wrong!'
130 else if (abs(dx-40.0_dp) < eps)
then
132 if ((imax /= 312).or.(jmax /= 312))
then
133 stop
' sico_init: IMAX and/or JMAX wrong!'
136 else if (abs(dx-20.0_dp) < eps)
then
138 if ((imax /= 624).or.(jmax /= 624))
then
139 stop
' sico_init: IMAX and/or JMAX wrong!'
143 stop
' sico_init: DX wrong!'
148 stop
' sico_init: GRID==2 not allowed for the northern hemisphere application!'
156 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
157 stop
' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
160 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
161 stop
' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
168 stop
' sico_init: Option ADV_HOR==1 (central differences) not defined!'
177 #elif (TSURFACE == 5)
186 dtime_temp0 = dtime_temp0
188 dtime_wss0 = dtime_wss0
193 dzeta_c = 1.0_dp/
real(kcmax,dp)
194 dzeta_t = 1.0_dp/
real(ktmax,dp)
195 dzeta_r = 1.0_dp/
real(krmax,dp)
208 if ((i <= imax).and.(j <= jmax))
then
220 anfdatname = anfdatname
222 #if defined(YEAR_ZERO)
223 year_zero = year_zero
225 year_zero = 2000.0_dp
228 time_init0 = time_init0
229 time_end0 = time_end0
230 dtime_ser0 = dtime_ser0
232 #if ( OUTPUT==1 || OUTPUT==3 )
233 dtime_out0 = dtime_out0
236 #if ( OUTPUT==2 || OUTPUT==3 )
238 time_output0( 1) = time_out0_01
239 time_output0( 2) = time_out0_02
240 time_output0( 3) = time_out0_03
241 time_output0( 4) = time_out0_04
242 time_output0( 5) = time_out0_05
243 time_output0( 6) = time_out0_06
244 time_output0( 7) = time_out0_07
245 time_output0( 8) = time_out0_08
246 time_output0( 9) = time_out0_09
247 time_output0(10) = time_out0_10
248 time_output0(11) = time_out0_11
249 time_output0(12) = time_out0_12
250 time_output0(13) = time_out0_13
251 time_output0(14) = time_out0_14
252 time_output0(15) = time_out0_15
253 time_output0(16) = time_out0_16
254 time_output0(17) = time_out0_17
255 time_output0(18) = time_out0_18
256 time_output0(19) = time_out0_19
257 time_output0(20) = time_out0_20
262 shell_command =
'if [ ! -d'
263 shell_command = trim(shell_command)//
' '//outpath
264 shell_command = trim(shell_command)//
' '//
'] ; then mkdir'
265 shell_command = trim(shell_command)//
' '//outpath
266 shell_command = trim(shell_command)//
' '//
'; fi'
267 call system(trim(shell_command))
270 open(10, iostat=ios, &
271 file=outpath//
'/'//trim(runname)//
'.log', &
274 stop
' sico_init: Error when opening the log file!'
276 write(10, fmt=trim(fmt1))
'Computational domain:'
278 write(10, fmt=trim(fmt1))
'Antarctica'
280 write(10, fmt=trim(fmt1))
'Greenland'
282 write(10, fmt=trim(fmt1))
'Scandinavia and Eurasia'
284 write(10, fmt=trim(fmt1))
'Northern hemisphere'
285 #elif defined(EMTP2SGE)
286 write(10, fmt=trim(fmt1))
'EISMINT Phase 2 Simplified Geometry Experiment'
288 write(10, fmt=trim(fmt1))
'North polar cap of Mars'
290 write(10, fmt=trim(fmt1))
'South polar cap of Mars'
292 stop
' sico_init: No valid domain specified!'
294 write(10, fmt=trim(fmt1))
' '
296 write(10, fmt=trim(fmt2))
'imax =', imax
297 write(10, fmt=trim(fmt2))
'jmax =', jmax
298 write(10, fmt=trim(fmt2))
'kcmax =', kcmax
299 write(10, fmt=trim(fmt2))
'ktmax =', ktmax
300 write(10, fmt=trim(fmt2))
'krmax =', krmax
301 write(10, fmt=trim(fmt1))
' '
303 write(10, fmt=trim(fmt3))
'a =', deform
304 write(10, fmt=trim(fmt1))
' '
306 #if (GRID==0 || GRID==1)
307 write(10, fmt=trim(fmt3))
'x0 =', x0
308 write(10, fmt=trim(fmt3))
'y0 =', y0
309 write(10, fmt=trim(fmt3))
'dx =', dx
311 stop
' sico_init: GRID==2 not allowed for this application!'
313 write(10, fmt=trim(fmt1))
' '
315 write(10, fmt=trim(fmt3))
'year_zero =', year_zero
316 write(10, fmt=trim(fmt3))
'time_init =', time_init0
317 write(10, fmt=trim(fmt3))
'time_end =', time_end0
318 write(10, fmt=trim(fmt3))
'dtime =', dtime0
319 write(10, fmt=trim(fmt3))
'dtime_temp =', dtime_temp0
321 write(10, fmt=trim(fmt3))
'dtime_wss =', dtime_wss0
323 #if ( OUTPUT==1 || OUTPUT==3 )
324 write(10, fmt=trim(fmt3))
'dtime_out =', dtime_out0
326 write(10, fmt=trim(fmt3))
'dtime_ser =', dtime_ser0
327 write(10, fmt=trim(fmt1))
' '
329 write(10, fmt=trim(fmt1))
'zs_present file = '//zs_present_file
331 write(10, fmt=trim(fmt1))
'zl_present file = '//zl_present_file
333 write(10, fmt=trim(fmt1))
'zl0 file = '//zl0_file
334 write(10, fmt=trim(fmt1))
'mask_present file = '//mask_present_file
335 write(10, fmt=trim(fmt1))
' '
337 write(10, fmt=trim(fmt1))
'Physical-parameter file = '//phys_para_file
338 write(10, fmt=trim(fmt1))
' '
340 #if (CALCZS==3 || CALCZS==4)
341 write(10, fmt=trim(fmt3))
'ovi_weight =', ovi_weight
343 write(10, fmt=trim(fmt3))
'omega_sor =', omega_sor
345 write(10, fmt=trim(fmt1))
' '
348 write(10, fmt=trim(fmt1))
'temp_mm_present file = '//temp_mm_present_file
350 write(10, fmt=trim(fmt3))
'delta_ts0 =', delta_ts0
352 write(10, fmt=trim(fmt3))
'sine_amplit =', sine_amplit
353 write(10, fmt=trim(fmt3))
'sine_period =', sine_period
355 write(10, fmt=trim(fmt1))
'GRIP file = '//grip_temp_file
356 write(10, fmt=trim(fmt3))
'grip_temp_fact =', grip_temp_fact
358 write(10, fmt=trim(fmt1))
'Glacial-index file = '//glac_ind_file
359 write(10, fmt=trim(fmt1))
'temp_mm_anom file = '//temp_mm_anom_file
360 write(10, fmt=trim(fmt3))
'temp_mm_anom fact = ', temp_mm_anom_fact
363 write(10, fmt=trim(fmt1))
'precip_mm_present file = '//precip_mm_present_file
365 write(10, fmt=trim(fmt3))
'accfact =', accfact
366 #elif ( ACCSURFACE==2 || ACCSURFACE==3 )
367 write(10, fmt=trim(fmt3))
'gamma_s =', gamma_s
368 #elif ( ACCSURFACE==5 )
369 write(10, fmt=trim(fmt1))
'precip_mm_anom file = '//precip_mm_anom_file
370 write(10, fmt=trim(fmt3))
'precip_mm_anom fact = ', precip_mm_anom_fact
372 #if (ACCSURFACE <= 3)
373 write(10, fmt=trim(fmt2))
'ELEV_DESERT =', elev_desert
374 #if (ELEV_DESERT == 1)
375 write(10, fmt=trim(fmt3))
'gamma_p =', gamma_p
376 write(10, fmt=trim(fmt3))
'zs_thresh =', zs_thresh
381 write(10, fmt=trim(fmt3))
'lambda_lti =', lambda_lti
382 write(10, fmt=trim(fmt3))
'temp_lti =', temp_lti
386 write(10, fmt=trim(fmt3))
'z_sl0 =', z_sl0
388 write(10, fmt=trim(fmt1))
'sea-level file = '//sea_level_file
390 write(10, fmt=trim(fmt1))
' '
393 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
394 write(10, fmt=trim(fmt3))
'z_mar =', z_mar
395 write(10, fmt=trim(fmt1))
' '
396 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
397 || marine_ice_calving==6 || marine_ice_calving==7 )
398 write(10, fmt=trim(fmt3))
'fact_z_mar =', fact_z_mar
399 write(10, fmt=trim(fmt1))
' '
400 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
401 write(10, fmt=trim(fmt3))
'calv_uw_coeff =', calv_uw_coeff
402 write(10, fmt=trim(fmt3))
'r1_calv_uw =', r1_calv_uw
403 write(10, fmt=trim(fmt3))
'r2_calv_uw =', r2_calv_uw
404 write(10, fmt=trim(fmt1))
' '
407 #if ICE_SHELF_CALVING==2
408 write(10, fmt=trim(fmt3))
'H_calv =', h_calv
409 write(10, fmt=trim(fmt1))
' '
413 write(10, fmt=trim(fmt3))
'c_slide =', c_slide
414 write(10, fmt=trim(fmt3))
'gamma_slide =', gamma_slide
415 write(10, fmt=trim(fmt2))
'p_weert =', p_weert
416 write(10, fmt=trim(fmt2))
'q_weert =', q_weert
418 write(10, fmt=trim(fmt3))
'red_pres_limit_fact =', red_pres_limit_fact
420 write(10, fmt=trim(fmt1))
' '
423 write(10, fmt=trim(fmt3))
'q_geo =', q_geo
425 write(10, fmt=trim(fmt1))
'q_geo file = '//q_geo_file
427 write(10, fmt=trim(fmt1))
' '
429 #if defined(MARINE_ICE_BASAL_MELTING)
430 write(10, fmt=trim(fmt2))
'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
431 #if ( MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3 )
432 write(10, fmt=trim(fmt3))
'qbm_marine =', qbm_marine
434 write(10, fmt=trim(fmt1))
' '
438 write(10, fmt=trim(fmt3))
'frac_llra =', frac_llra
439 write(10, fmt=trim(fmt1))
' '
443 write(10, fmt=trim(fmt3))
'gr_size =', gr_size
444 write(10, fmt=trim(fmt1))
' '
447 write(10, fmt=trim(fmt3))
'sigma_res =', sigma_res
448 write(10, fmt=trim(fmt1))
' '
451 write(10, fmt=trim(fmt2))
'ENHMOD =', enhmod
452 write(10, fmt=trim(fmt3))
'enh_fact =', enh_fact
453 #if ( ENHMOD==2 || ENHMOD==3 )
454 write(10, fmt=trim(fmt3))
'enh_intg =', enh_intg
457 write(10, fmt=trim(fmt3))
'age_trans =', age_trans_0
460 write(10, fmt=trim(fmt3))
'date_trans1 =', date_trans1_0
461 write(10, fmt=trim(fmt3))
'date_trans2 =', date_trans2_0
462 write(10, fmt=trim(fmt3))
'date_trans3 =', date_trans3_0
465 write(10, fmt=trim(fmt3))
'enh_shelf =', enh_shelf
467 write(10, fmt=trim(fmt1))
' '
469 write(10, fmt=trim(fmt3))
'numdiff_H_t =', numdiff_h_t
470 write(10, fmt=trim(fmt3))
'tau_cts =', tau_cts
471 write(10, fmt=trim(fmt3))
'H_min =', h_min
472 write(10, fmt=trim(fmt3))
'vh_max =', vh_max
473 write(10, fmt=trim(fmt3))
'hd_min =', hd_min
474 write(10, fmt=trim(fmt3))
'hd_max =', hd_max
476 write(10, fmt=trim(fmt3))
'qbm_min =', qbm_min
477 #elif defined(QB_MIN) /* obsolete */
478 write(10, fmt=trim(fmt3))
'qb_min =', qb_min
481 write(10, fmt=trim(fmt3))
'qbm_max =', qbm_max
482 #elif defined(QB_MAX) /* obsolete */
483 write(10, fmt=trim(fmt3))
'qb_max =', qb_max
485 write(10, fmt=trim(fmt3))
'age_min =', age_min
486 write(10, fmt=trim(fmt3))
'age_max =', age_max
487 write(10, fmt=trim(fmt3))
'mean_accum =', mean_accum
489 write(10, fmt=trim(fmt3))
'age_diff =', agediff
491 write(10, fmt=trim(fmt1))
' '
493 write(10, fmt=trim(fmt2))
'GRID =', grid
494 write(10, fmt=trim(fmt2))
'CALCMOD =', calcmod
495 write(10, fmt=trim(fmt2))
'FLOW_LAW =', flow_law
496 write(10, fmt=trim(fmt2))
'FIN_VISC =', fin_visc
497 write(10, fmt=trim(fmt2))
'MARGIN =', margin
499 write(10, fmt=trim(fmt2))
'MARINE_ICE_FORMATION =', marine_ice_formation
500 write(10, fmt=trim(fmt2))
'MARINE_ICE_CALVING =', marine_ice_calving
502 write(10, fmt=trim(fmt2))
'ICE_SHELF_CALVING =', ice_shelf_calving
504 write(10, fmt=trim(fmt2))
'ANF_DAT =', anf_dat
505 write(10, fmt=trim(fmt2))
'REBOUND =', rebound
506 write(10, fmt=trim(fmt2))
'Q_LITHO =', q_litho
507 write(10, fmt=trim(fmt2))
'ZS_EVOL =', zs_evol
508 write(10, fmt=trim(fmt2))
'CALCZS =', calczs
509 write(10, fmt=trim(fmt2))
'ADV_HOR =', adv_hor
510 write(10, fmt=trim(fmt2))
'ADV_VERT =', adv_vert
511 write(10, fmt=trim(fmt2))
'TOPOGRAD =', topograd
512 write(10, fmt=trim(fmt2))
'TSURFACE =', tsurface
513 write(10, fmt=trim(fmt2))
'ACCSURFACE =', accsurface
515 write(10, fmt=trim(fmt2))
'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
517 write(10, fmt=trim(fmt2))
'SOLID_PRECIP =', solid_precip
518 write(10, fmt=trim(fmt2))
'ABLSURFACE =', ablsurface
519 write(10, fmt=trim(fmt2))
'SEA_LEVEL =', sea_level
520 write(10, fmt=trim(fmt2))
'SLIDE_LAW =', slide_law
521 write(10, fmt=trim(fmt2))
'Q_GEO_MOD =', q_geo_mod
522 write(10, fmt=trim(fmt1))
' '
524 #if defined(OUT_TIMES)
525 write(10, fmt=trim(fmt2))
'OUT_TIMES =', out_times
527 write(10, fmt=trim(fmt2))
'OUTPUT =', output
528 write(10, fmt=trim(fmt2))
'OUTSER =', outser
529 #if defined(OUTSER_V_A)
530 write(10, fmt=trim(fmt2))
'OUTSER_V_A =', outser_v_a
532 #if ( OUTPUT==1 || OUTPUT==2 )
533 write(10, fmt=trim(fmt2))
'ERGDAT =', ergdat
535 write(10, fmt=trim(fmt2))
'NETCDF =', netcdf
536 write(10, fmt=trim(fmt1))
' '
537 #if ( OUTPUT==2 || OUTPUT==3 )
538 write(10, fmt=trim(fmt2))
'n_output =', n_output
540 write(10, fmt=trim(fmt3))
'time_output =', time_output0(n)
542 write(10, fmt=trim(fmt1))
' '
545 write(10, fmt=trim(fmt1))
'Program version and date: '//version//
' / '//date
547 close(10, status=
'keep')
551 year_zero = year_zero*year_sec
552 time_init = time_init0*year_sec
553 time_end = time_end0*year_sec
554 dtime = dtime0*year_sec
555 dtime_temp = dtime_temp0*year_sec
557 dtime_wss = dtime_wss0*year_sec
559 dtime_ser = dtime_ser0*year_sec
560 #if ( OUTPUT==1 || OUTPUT==3 )
561 dtime_out = dtime_out0*year_sec
563 #if ( OUTPUT==2 || OUTPUT==3 )
565 time_output(n) = time_output0(n)*year_sec
569 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
570 stop
' sico_init: dtime_temp must be a multiple of dtime!'
573 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
574 stop
' sico_init: dtime_wss must be a multiple of dtime!'
581 #if (GRID==0 || GRID==1)
583 open(21, iostat=ios, &
584 file=inpath//
'/nhem/'//precip_mm_present_file, &
585 recl=16384, status=
'old')
589 stop
' sico_init: GRID==2 not allowed for the northern hemisphere application!'
593 if (ios /= 0) stop
' sico_init: Error when opening the precipitation file!'
595 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
598 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do
600 read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
604 close(21, status=
'keep')
608 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
615 #if (GRID==0 || GRID==1)
617 open(21, iostat=ios, &
618 file=inpath//
'/nhem/'//precip_mm_anom_file, &
619 recl=16384, status=
'old')
623 if (ios /= 0) stop
' sico_init: Error when opening the precip anomaly file!'
625 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
628 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do
630 read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
634 close(21, status=
'keep')
636 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
641 #if (PRECIP_ANOM_INTERPOL==1)
643 gamma_precip_lgm_anom(j,i,n) = 0.0_dp
645 #elif (PRECIP_ANOM_INTERPOL==2)
647 gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
650 stop
' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
660 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
663 mean_accum_inv = 1.0_dp/mean_accum
667 #if (GRID==0 || GRID==1)
669 open(24, iostat=ios, &
670 file=inpath//
'/nhem/'//mask_present_file, &
671 recl=1024, status=
'old')
675 stop
' sico_init: GRID==2 not allowed for the northern hemisphere application!'
680 stop
' sico_init: Error when opening the mask file!'
682 do m=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do
685 read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
688 close(24, status=
'keep')
692 #if (GRID==0 || GRID==1)
694 open(21, iostat=ios, &
695 file=inpath//
'/nhem/'//temp_mm_present_file, &
696 recl=16384, status=
'old')
700 stop
' sico_init: GRID==2 not allowed for the northern hemisphere application!'
704 if (ios /= 0) stop
' sico_init: Error when opening the temperature file!'
706 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
709 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do
711 read(21, fmt=*) (temp_mm_present(j,i,n), i=0,imax)
715 close(21, status=
'keep')
721 #if (GRID==0 || GRID==1)
723 open(21, iostat=ios, &
724 file=inpath//
'/nhem/'//temp_mm_anom_file, &
725 recl=16384, status=
'old')
729 if (ios /= 0) stop
' sico_init: Error when opening the temperature anomaly file!'
731 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
734 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do
736 read(21, fmt=*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
740 close(21, status=
'keep')
742 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
749 #if (GRID==0 || GRID==1)
751 open(21, iostat=ios, &
752 file=inpath//
'/nhem/'//zs_present_file, &
753 recl=16384, status=
'old')
757 stop
' sico_init: GRID==2 not allowed for the northern hemisphere application!'
761 if (ios /= 0) stop
' sico_init: Error when opening the zs file!'
763 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
766 read(21, fmt=*) (zs_ref(j,i), i=0,imax)
769 close(21, status=
'keep')
776 if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
784 open(21, iostat=ios, &
785 file=inpath//
'/general/'//grip_temp_file, &
788 stop
' sico_init: Error when opening the data file for delta_ts!'
790 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
792 if (ch_dummy /=
'#')
then
793 write(6, fmt=*)
' sico_init: grip_time_min, grip_time_stp, grip_time_max'
794 write(6, fmt=*)
' not defined in data file for delta_ts!'
798 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
800 allocate(griptemp(0:ndata_grip))
803 read(21, fmt=*) d_dummy, griptemp(n)
806 close(21, status=
'keep')
814 open(21, iostat=ios, &
815 file=inpath//
'/general/'//glac_ind_file, status=
'old')
816 if (ios /= 0) stop
' sico_init: Error when opening the glacial-index file!'
818 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
820 if (ch_dummy /=
'#')
then
821 write(6, fmt=*)
' sico_init: gi_time_min, gi_time_stp, gi_time_max'
822 write(6, fmt=*)
' not defined in glacial-index file!'
826 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
828 allocate(glacial_index(0:ndata_gi))
831 read(21, fmt=*) d_dummy, glacial_index(n)
834 close(21, status=
'keep')
842 open(21, iostat=ios, &
843 file=inpath//
'/general/'//sea_level_file, &
846 stop
' sico_init: Error when opening the data file for z_sl!'
848 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
850 if (ch_dummy /=
'#')
then
851 write(6, fmt=*)
' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
852 write(6, fmt=*)
' not defined in data file for z_sl!'
856 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
858 allocate(specmap_zsl(0:ndata_specmap))
860 do n=0, ndata_specmap
861 read(21, fmt=*) d_dummy, specmap_zsl(n)
864 close(21, status=
'keep')
876 q_geo(j,i) = q_geo *1.0e-03_dp
884 open(21, iostat=ios, &
885 file=inpath//
'/nhem/'//q_geo_file, &
886 recl=16384, status=
'old')
888 if (ios /= 0) stop
' sico_init: Error when opening the qgeo file!'
890 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
893 read(21, fmt=*) (q_geo(j,i), i=0,imax)
896 close(21, status=
'keep')
900 q_geo(j,i) = q_geo(j,i) *1.0e-03_dp
908 #if (REBOUND==0 || REBOUND==1)
910 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
924 stop
' sico_init: ANF_DAT==1 not allowed for northern hemisphere application!'
934 call
boundary(time_init, dtime, dxi, deta, &
935 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
942 q_b_tot(j,i) = 0.0_dp
948 temp_c(kc,j,i) = temp_s(j,i)
949 age_c(kc,j,i) = 15000.0_dp*year_sec
953 omega_t(kt,j,i) = 0.0_dp
954 age_t(kt,j,i) = 15000.0_dp*year_sec
958 temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
959 *(1.0_dp-
real(kr,dp)/
real(krmax,dp))
978 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
981 call
boundary(time_init, dtime, dxi, deta, &
982 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
984 q_b_tot = q_bm + q_tld
992 #if (GRID==0 || GRID==1)
996 dist_dxdy(jr,ir) = sqrt( (
real(ir,dp)*dxi)**2 + (
real(jr,dp)*deta)**2 )
1008 zeta_c(kc) = kc*dzeta_c
1009 eaz_c(kc) = exp(deform*zeta_c(kc))
1010 eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
1014 zeta_t(kt) = kt*dzeta_t
1025 call
calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
1038 open(12, iostat=ios, &
1039 file=outpath//
'/'//trim(runname)//
'.ser', &
1042 stop
' sico_init: Error when opening the ser file!'
1044 if (forcing_flag == 1)
then
1049 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1051 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1052 ' H_max(m) zs_max(m) V_g(m^3)', &
1053 ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1054 ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1055 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1056 1103
format(
'----------------------------------------------------', &
1057 '---------------------------------------')
1061 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1062 ' V(m^3) V_g(m^3) V_f(m^3)', &
1063 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1064 ' H_max(m) zs_max(m)', &
1065 ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1066 ' A_t(m^2) V_bm(m^3/a)', &
1067 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1068 1103
format(
'----------------------------------------------------', &
1069 '---------------------------------------')
1072 stop
' sico_init: OUTSER_V_A must be either 1 or 2!'
1075 else if (forcing_flag == 2)
then
1080 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1082 1112
format(
' t(a) glac_ind(1) z_sl(m)',/, &
1083 ' H_max(m) zs_max(m) V_g(m^3)', &
1084 ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1085 ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1086 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1087 1113
format(
'----------------------------------------------------', &
1088 '---------------------------------------')
1092 1112
format(
' t(a) glac_ind(1) z_sl(m)',/, &
1093 ' V(m^3) V_g(m^3) V_f(m^3)', &
1094 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1095 ' H_max(m) zs_max(m)', &
1096 ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1097 ' A_t(m^2) V_bm(m^3/a)', &
1098 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1099 1113
format(
'----------------------------------------------------', &
1100 '---------------------------------------')
1103 stop
' sico_init: OUTSER_V_A must be either 1 or 2!'
1112 open(13, iostat=ios, &
1113 file=outpath//
'/'//trim(runname)//
'.thk', &
1115 if (ios /= 0) stop
' sico_init: Error when opening the thk file!'
1117 if (forcing_flag == 1)
then
1122 1104
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1123 ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1124 ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1125 1105
format(
'----------------------------------------------------')
1127 else if (forcing_flag == 2)
then
1132 1114
format(
' t(a) glac_ind(1) z_sl(m)',/, &
1133 ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1134 ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1135 1115
format(
'----------------------------------------------------')
1146 flag_3d_output = .false.
1148 flag_3d_output = .true.
1150 stop
' sico_init: ERGDAT must be either 0 or 1!'
1154 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1155 flag_3d_output, ndat2d, ndat3d)
1157 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1158 flag_3d_output, ndat2d, ndat3d)
1160 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1165 if (time_output(1) <= time_init+eps)
then
1168 flag_3d_output = .false.
1170 flag_3d_output = .true.
1172 stop
' sico_init: ERGDAT must be either 0 or 1!'
1176 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1177 flag_3d_output, ndat2d, ndat3d)
1179 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1180 flag_3d_output, ndat2d, ndat3d)
1182 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1189 flag_3d_output = .false.
1192 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1193 flag_3d_output, ndat2d, ndat3d)
1195 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1196 flag_3d_output, ndat2d, ndat3d)
1198 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1201 if (time_output(1) <= time_init+eps)
then
1203 flag_3d_output = .true.
1206 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1207 flag_3d_output, ndat2d, ndat3d)
1209 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1210 flag_3d_output, ndat2d, ndat3d)
1212 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1218 stop
' sico_init: OUTPUT must be either 1, 2 or 3!'
1221 call
output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1223 call
output3(time_init, delta_ts, glac_index, z_sl)