SW_VAR_DENS_MODEL  0.9
Dept-averagedgas-particlemodel
inpout_2d.f90
Go to the documentation of this file.
1 !********************************************************************************
3 !
6 !
10 !
11 !********************************************************************************
12 
13 MODULE inpout_2d
14 
15  ! -- Variables for the namelist RUN_PARAMETERS
16  USE parameters_2d, ONLY : t_start , t_end , dt_output
17 
18  USE solver_2d, ONLY : verbose_level
19 
20  ! -- Variables for the namelist NEWRUN_PARAMETERS
21  USE geometry_2d, ONLY : x0 , y0 , comp_cells_x , comp_cells_y , cell_size
22  USE geometry_2d, ONLY : topography_profile , n_topography_profile_x , &
23  n_topography_profile_y
24  USE init_2d, ONLY : riemann_interface
25  USE parameters_2d, ONLY : riemann_flag , rheology_flag , energy_flag , &
26  topo_change_flag , radial_source_flag , collapsing_volume_flag , &
27  liquid_flag , gas_flag
28 
29  ! -- Variables for the namelist INITIAL_CONDITIONS
30  USE parameters_2d, ONLY : released_volume , x_release , y_release
31  USE parameters_2d, ONLY : velocity_mod_release , velocity_ang_release
32  USE parameters_2d, ONLY : alphas_init , sed_vol_perc
33  USE parameters_2d, ONLY : t_init
34 
35  ! -- Variables for the namelist LEFT_STATE
36  USE init_2d, ONLY : hb_w , u_w , v_w , alphas_w , t_w
37 
38  ! -- Variables for the namelist RIGHT_STATE
39  USE init_2d, ONLY : hb_e , u_e , v_e , alphas_e , t_e
40 
41  ! -- Variables for the namelists LEFT/RIGHT_BOUNDARY_CONDITIONS
42  USE parameters_2d, ONLY : bc
43 
44  ! -- Variables for the namelist NUMERIC_PARAMETERS
45  USE parameters_2d, ONLY : solver_scheme, dt0 , max_dt , cfl, limiter , theta, &
46  reconstr_coeff , interfaces_relaxation , n_rk
47 
48  ! -- Variables for the namelist EXPL_TERMS_PARAMETERS
49  USE constitutive_2d, ONLY : grav
50 
51  ! -- Variables for the namelist RADIAL_SOURCE_PARAMETERS
52  USE parameters_2d, ONLY : x_source , y_source , r_source , vel_source , &
53  t_source , h_source , alphas_source , alphal_source , time_param
54 
55  ! -- Variables for the namelist COLLAPSING_VOLUME_PARAMETERS
56  USE parameters_2d, ONLY : x_collapse , y_collapse , r_collapse , t_collapse , &
57  h_collapse , alphas_collapse
58 
59  ! -- Variables for the namelist TEMPERATURE_PARAMETERS
60  USE constitutive_2d, ONLY : emissivity , exp_area_fract , enne , emme , &
61  atm_heat_transf_coeff , thermal_conductivity , t_env , t_ground , c_p
62 
63  ! -- Variables for the namelist RHEOLOGY_PARAMETERS
64  USE parameters_2d, ONLY : rheology_model
65  USE constitutive_2d, ONLY : mu , xi , tau , nu_ref , visc_par , t_ref
66  USE constitutive_2d, ONLY : alpha2 , beta2 , alpha1_coeff , beta1 , kappa ,n_td
67  USE constitutive_2d, ONLY : friction_factor
68 
69  ! --- Variables for the namelist SOLID_TRANSPORT_PARAMETERS
70  USE parameters_2d, ONLY : n_solid
71  USE constitutive_2d, ONLY : rho_s , diam_s , sp_heat_s
72  USE constitutive_2d, ONLY : settling_flag , erosion_coeff
73  USE constitutive_2d, ONLY : t_s_substrate
74 
75  ! --- Variables for the namelist GAS_TRANSPORT_PARAMETERS
76  USE constitutive_2d, ONLY : sp_heat_a , sp_gas_const_a , kin_visc_a , pres , &
77  t_ambient , entrainment_flag
78 
79  ! --- Variables for the namelist LIQUID_TRANSPORT_PARAMETERS
80  USE constitutive_2d, ONLY : sp_heat_l , rho_l , kin_visc_l
81 
82 
83  IMPLICIT NONE
84 
85  CHARACTER(LEN=40) :: run_name
86  CHARACTER(LEN=40) :: bak_name
87  CHARACTER(LEN=40) :: input_file
88  CHARACTER(LEN=40) :: output_file
89  CHARACTER(LEN=40) :: restart_file
90  CHARACTER(LEN=40) :: probes_file
91  CHARACTER(LEN=40) :: output_file_2d
92  CHARACTER(LEN=40) :: output_esri_file
93  CHARACTER(LEN=40) :: output_max_file
94  CHARACTER(LEN=40) :: runout_file
95 
96  INTEGER, PARAMETER :: input_unit = 7
97  INTEGER, PARAMETER :: backup_unit = 8
98  INTEGER, PARAMETER :: output_unit = 9
99  INTEGER, PARAMETER :: restart_unit = 10
100  INTEGER, PARAMETER :: probes_unit = 11
101  INTEGER, PARAMETER :: output_unit_2d = 12
102  INTEGER, PARAMETER :: output_esri_unit = 13
103  INTEGER, PARAMETER :: output_max_unit = 14
104  INTEGER, PARAMETER :: dem_esri_unit = 15
105  INTEGER, PARAMETER :: runout_unit = 16
106  INTEGER, PARAMETER :: dakota_unit = 17
107 
109  INTEGER :: output_idx
110 
115  LOGICAL :: restart
116 
121  LOGICAL :: output_esri_flag
122 
127  LOGICAL :: output_phys_flag
128 
133  LOGICAL :: output_cons_flag
134 
139  LOGICAL :: output_runout_flag
140 
141  ! -- Variables for the namelists WEST_BOUNDARY_CONDITIONS
142  TYPE(bc) :: h_bcw , hu_bcw , hv_bcw , t_bcw
143  TYPE(bc), ALLOCATABLE :: alphas_bcw(:)
144 
145  ! -- Variables for the namelists EAST_BOUNDARY_CONDITIONS
146  TYPE(bc) :: h_bce , hu_bce , hv_bce , t_bce
147  TYPE(bc), ALLOCATABLE :: alphas_bce(:)
148 
149  ! -- Variables for the namelists SOUTH_BOUNDARY_CONDITIONS
150  TYPE(bc) :: h_bcs , hu_bcs , hv_bcs , t_bcs
151  TYPE(bc), ALLOCATABLE :: alphas_bcs(:)
152 
153  ! -- Variables for the namelists NORTH_BOUNDARY_CONDITIONS
154  TYPE(bc) :: h_bcn , hu_bcn , hv_bcn , t_bcn
155  TYPE(bc), ALLOCATABLE :: alphas_bcn(:)
156 
157 
158  ! parameters to read a dem file
159  INTEGER :: ncols, nrows, nodata_value
160 
161  REAL*8 :: xllcorner, yllcorner, cellsize
162 
163  LOGICAL :: write_first_q
164 
165  INTEGER :: n_probes
166 
167  REAL*8, ALLOCATABLE :: probes_coords(:,:)
168 
169  REAL*8 :: dt_runout
170 
171  REAL*8, ALLOCATABLE :: h_old(:,:)
172 
173  REAL*8 :: x0_runout, y0_runout , init_runout , eps_stop
174 
175  REAL*8 :: rho0_s(1000) , diam0_s(1000) , sp_heat0_s(1000), erosion_coeff0(1000)
176 
177  REAL*8 :: alpha1_ref
178 
179  namelist / run_parameters / run_name , restart , t_start , t_end , dt_output ,&
180  output_cons_flag , output_esri_flag , output_phys_flag , &
181  output_runout_flag , verbose_level
182 
183  namelist / restart_parameters / restart_file, t_init, t_ambient , sed_vol_perc
184 
185  namelist / newrun_parameters / x0 , y0 , comp_cells_x , comp_cells_y , &
186  cell_size , rheology_flag , riemann_flag , energy_flag , liquid_flag , &
187  radial_source_flag , collapsing_volume_flag , topo_change_flag , gas_flag
188 
189  namelist / initial_conditions / released_volume , x_release , y_release , &
190  velocity_mod_release , velocity_ang_release , t_init , t_ambient
191 
192  namelist / left_state / riemann_interface , hb_w , u_w , v_w , alphas_w , t_w
193 
194  namelist / right_state / hb_e , u_e , v_e , alphas_e , t_e
195 
196  namelist / west_boundary_conditions / h_bcw , hu_bcw , hv_bcw , alphas_bcw , &
197  t_bcw
198 
199  namelist / east_boundary_conditions / h_bce , hu_bce , hv_bce , alphas_bce , &
200  t_bce
201 
202  namelist / south_boundary_conditions / h_bcs , hu_bcs , hv_bcs , alphas_bcs , &
203  t_bcs
204 
205  namelist / north_boundary_conditions / h_bcn , hu_bcn , hv_bcn , alphas_bcn , &
206  t_bcn
207 
208  namelist / numeric_parameters / solver_scheme, dt0 , max_dt , cfl, limiter , &
209  theta , reconstr_coeff , interfaces_relaxation , n_rk
210 
211  namelist / expl_terms_parameters / grav
212 
213  namelist / radial_source_parameters / x_source , y_source , r_source , &
214  vel_source , t_source , h_source , alphas_source , alphal_source , time_param
215 
216  namelist / collapsing_volume_parameters / x_collapse , y_collapse , &
217  r_collapse , t_collapse , h_collapse , alphas_collapse
218 
219  namelist / temperature_parameters / emissivity , atm_heat_transf_coeff , &
220  thermal_conductivity , exp_area_fract , c_p , enne , emme , t_env , &
221  t_ground
222 
223  namelist / rheology_parameters / rheology_model , mu , xi , tau , nu_ref , &
224  visc_par , t_ref , alpha2 , beta2 , alpha1_ref , beta1 , kappa , n_td , &
225  friction_factor
226 
227  namelist / runout_parameters / x0_runout , y0_runout , dt_runout , &
228  eps_stop
229 
230  namelist / solid_transport_parameters / n_solid , rho0_s , diam0_s , &
231  sp_heat0_s , erosion_coeff0 , settling_flag , t_s_substrate
232 
233  namelist / gas_transport_parameters / sp_heat_a , sp_gas_const_a , kin_visc_a,&
234  pres , t_ambient , entrainment_flag
235 
236  namelist / liquid_transport_parameters / sp_heat_l , rho_l , kin_visc_l
237 
238 CONTAINS
239 
240  !******************************************************************************
242  !
246  !
250  !
251  !******************************************************************************
252 
253  SUBROUTINE init_param
254 
255  USE parameters_2d , ONLY : n_vars
256 
257  IMPLICIT none
258 
259  LOGICAL :: lexist
260 
261  INTEGER :: j , k
262 
263  n_vars = 3
264 
265  !-- Inizialization of the Variables for the namelist RUN_PARAMETERS
266  run_name = 'default'
267  restart = .false.
268  t_start = 0.0
269  t_end = 5.0d-2
270  dt_output = 5.d-3
271  output_cons_flag = .true.
272  output_esri_flag = .true.
273  output_phys_flag = .true.
274  output_runout_flag = .false.
275  verbose_level = 0
276 
277  !-- Inizialization of the Variables for the namelist restart parameters
278  restart_file = ''
279  t_init = 0.d0
280  t_ambient = 0.d0
281 
282  !-- Inizialization of the Variables for the namelist newrun_parameters
283  x0 = 0.d0
284  y0 = 0.d0
285  comp_cells_x = 1000
286  comp_cells_y = 1
287  cell_size = 1.0d-3
288  rheology_flag = .false.
289  riemann_flag =.true.
290  energy_flag = .false.
291  topo_change_flag = .false.
292  radial_source_flag = .false.
293  collapsing_volume_flag = .false.
294  liquid_flag = .false.
295  gas_flag = .true.
296 
297  !-- Inizialization of the Variables for the namelist left_state
298  riemann_interface = 0.5d0
299  hb_w = 2.0d0
300  u_w = 0.0d0
301  v_w = 0.d0
302  ! alphas_W = 0.5D0
303  t_w = -1.d0
304 
305  !-- Inizialization of the Variables for the namelist right_state
306  hb_e = 1.0d0
307  u_e = 0.0d0
308  u_e = 0.0d0
309  ! alphas_E = 0.5D0
310  t_e = -1.d0
311 
312  !-- Inizialization of the Variables for the namelist west boundary conditions
313  h_bcw%flag = -1
314  h_bcw%value = 0.d0
315 
316  hu_bcw%flag = 1
317  hu_bcw%value = 0.d0
318 
319  hv_bcw%flag = 1
320  hv_bcw%value = 0.d0
321 
322  !alphas_bcW%flag = 1
323  !alphas_bcW%value = 0.d0
324 
325  !-- Inizialization of the Variables for the namelist east boundary conditions
326  h_bce%flag = -1
327  h_bce%value = 0.d0
328 
329  hu_bce%flag = 1
330  hu_bce%value = 0.d0
331 
332  hv_bce%flag = 1
333  hv_bce%value = 0.d0
334 
335  !alphas_bcE%flag = 1
336  !alphas_bcE%value = 0.d0
337 
338  !-- Inizialization of the Variables for the namelist south boundary conditions
339  h_bcs%flag = -1
340  h_bcs%value = 0.d0
341 
342  hu_bcs%flag = 1
343  hu_bcs%value = 0.d0
344 
345  hv_bcs%flag = 1
346  hv_bcs%value = 0.d0
347 
348  !alphas_bcS%flag = 1
349  !alphas_bcS%value = 0.d0
350 
351  !-- Inizialization of the Variables for the namelist north boundary conditions
352  h_bcn%flag = -1
353  h_bcn%value = 0.d0
354 
355  hu_bcn%flag = 1
356  hu_bcn%value = 0.d0
357 
358  hv_bcn%flag = 1
359  hv_bcn%value = 0.d0
360 
361  !alphas_bcN%flag = 1
362  !alphas_bcN%value = 0.d0
363 
364  !-- Inizialization of the Variables for the namelist NUMERIC_PARAMETERS
365  dt0 = 1.d-4
366  max_dt = 1.d-3
367  solver_scheme = 'KT'
368  n_rk = 2
369  cfl = 0.24d0
370  limiter(1:n_vars+2) = 1
371  theta=1.0
372  reconstr_coeff = 1.0
373 
374  !-- Inizialization of the Variables for the namelist EXPL_TERMS_PARAMETERS
375  grav = 9.81d0
376 
377  !-- Inizialization of the Variables for the namelist TEMPERATURE_PARAMETERS
378  exp_area_fract = 0.5d0
379  emissivity = 0.0d0 ! no radiation to atmosphere
380  atm_heat_transf_coeff = 0.0d0 ! no convection to atmosphere
381  thermal_conductivity = 0.0d0 ! no conduction to ground
382  enne = 4.0d0
383  emme = 12.d0
384  t_env = 300.0d0
385  t_ground = 1200.0d0
386  c_p = 1200.d0
387 
388  !-- Inizialization of the Variables for the namelist RHEOLOGY_PARAMETERS
389  rheology_model = 0
390  nu_ref = 0.0d0
391  mu = 0.d0
392  xi = 0.d0
393  tau = 0.d0
394  t_ref = 0.d0
395  visc_par = 0.0d0
396 
397  !-- Inizialization of the Variables for the namelist RUNOUT_PARAMETERS
398  x0_runout = -1
399  y0_runout = -1
400  dt_runout = 60
401  eps_stop = 0.d0
402 
403  !-------------- Check if input file exists ----------------------------------
404  input_file = 'SW_VAR_DENS_MODEL.inp'
405 
406  INQUIRE (file=input_file,exist=lexist)
407 
408  IF (lexist .EQV. .false.) THEN
409 
410  OPEN(input_unit,file=input_file,status='NEW')
411 
412  WRITE(input_unit, run_parameters )
413  WRITE(input_unit, newrun_parameters )
414  WRITE(input_unit, left_state )
415  WRITE(input_unit, right_state )
416  WRITE(input_unit, numeric_parameters )
417  WRITE(input_unit, west_boundary_conditions )
418  WRITE(input_unit, east_boundary_conditions )
419  WRITE(input_unit, south_boundary_conditions )
420  WRITE(input_unit, north_boundary_conditions )
421  WRITE(input_unit, expl_terms_parameters )
422 
423  n_topography_profile_x = 2
424  n_topography_profile_y = 2
425 
426  ALLOCATE( topography_profile( 3 , n_topography_profile_x , &
427  n_topography_profile_y) )
428 
429  topography_profile(1,1,1) = 0.d0
430  topography_profile(2,1,1) = 0.d0
431  topography_profile(3,1,1) = 0.d0
432 
433  topography_profile(1,1,2) = 0.d0
434  topography_profile(2,1,2) = 1.d0
435  topography_profile(3,1,2) = 0.d0
436 
437  topography_profile(1,2,1) = 1.d0
438  topography_profile(2,2,1) = 0.d0
439  topography_profile(3,2,1) = 0.d0
440 
441  topography_profile(1,2,2) = 1.d0
442  topography_profile(2,2,2) = 1.d0
443  topography_profile(3,2,2) = 0.d0
444 
445 
446 
447  WRITE(input_unit,*) '''TOPOGRAPHY_PROFILE'''
448  WRITE(input_unit,*) n_topography_profile_x
449  WRITE(input_unit,*) n_topography_profile_y
450 
451  DO j = 1, n_topography_profile_x
452 
453  DO k = 1, n_topography_profile_y
454 
455  WRITE(*,108) topography_profile(1:3,j,k)
456  WRITE(input_unit,108) topography_profile(1:3,j,k)
457 
458 108 FORMAT(3(1x,e14.7))
459 
460  ENDDO
461 
462  END DO
463 
464  CLOSE(input_unit)
465 
466  WRITE(*,*) 'Input file SW_VAR_DENS_MODEL.inp not found'
467  WRITE(*,*) 'A new one with default values has been created'
468  stop
469 
470  ELSE
471 
472  END IF
473 
474  ! output file index
475  output_idx = 0
476 
477  ! -------------- Initialize values for checks during input reading ----------
478  h_bcw%flag = -1
479  hu_bcw%flag = -1
480  hv_bcw%flag = -1
481  !alphas_bcW%flag = -1
482  t_bcw%flag = -1
483 
484  h_bce%flag = -1
485  hu_bce%flag = -1
486  hv_bce%flag = -1
487  !alphas_bcE%flag = -1
488  t_bce%flag = -1
489 
490  h_bcs%flag = -1
491  hu_bcs%flag = -1
492  hv_bcs%flag = -1
493  !alphas_bcS%flag = -1
494  t_bcs%flag = -1
495 
496  h_bcn%flag = -1
497  hu_bcn%flag = -1
498  hv_bcn%flag = -1
499  !alphas_bcN%flag = -1
500  t_bcn%flag = -1
501 
502  ! sed_vol_perc = -1.D0
503 
504  rheology_model = -1
505  mu = -1
506  xi = -1
507  tau = -1
508  nu_ref = -1
509  visc_par = -1
510  t_ref = -1
511  friction_factor = -1
512 
513  alpha2 = -1
514  beta2 = -1
515  alpha1_ref = -1
516  beta1 = -1
517  kappa = -1
518  n_td = -1
519  rho0_s = -1
520 
521  exp_area_fract = -1.d0
522  emissivity = -1.d0
523  atm_heat_transf_coeff = -1.d0
524  thermal_conductivity = -1.d0
525  enne = -1.d0
526  emme = -1.d0
527  t_env = -1.d0
528  t_ground = -1.d0
529  c_p = -1.d0
530 
531  grav = -1.d0
532 
533  x0_runout = -1.d0
534  y0_runout = -1.d0
535 
536  !- Variables for the namelist SOLID_TRANSPORT_PARAMETERS
537  ! rho_s = -1.D0
538  settling_flag = .false.
539  erosion_coeff0 = -1.d0
540  n_solid = -1
541  t_s_substrate = -1.d0
542 
543  !- Variables for the namelist GAS_TRANSPORT_PARAMETERS
544  sp_heat_a = -1.d0
545  sp_gas_const_a = -1.d0
546  kin_visc_a = -1.d0
547  pres = -1.d0
548  t_ambient = -1.d0
549 
550  !- Variables for the namelist LIQUID_TRANSPORT_PARAMETERS
551  sp_heat_l = -1.d0
552  rho_l = -1.d0
553  kin_visc_l = -1.d0
554 
555  !- Variables for the namelist RADIAL_SOURCE_PARAMETERS
556  t_source = -1.d0
557  h_source = -1.d0
558  r_source = -1.d0
559  vel_source = -1.d0
560  time_param(1:4) = -1.d0
561 
562  !- Variables for the namelist COLLAPSING_VOLUME_PARAMETERS
563  t_collapse = -1.d0
564  h_collapse = -1.d0
565  r_collapse = -1.d0
566 
567  END SUBROUTINE init_param
568 
569  !******************************************************************************
571  !
576  !
580  !
581  !******************************************************************************
582 
583  SUBROUTINE read_param
584 
585  USE parameters_2d, ONLY : n_vars , n_eqns
586  USE parameters_2d, ONLY : limiter
587  USE parameters_2d, ONLY : bcw , bce , bcs , bcn
588 
589  USE geometry_2d, ONLY : deposit
590 
591  USE constitutive_2d, ONLY : rho_a_amb
592  USE constitutive_2d, ONLY : kin_visc_c
593 
594  IMPLICIT none
595 
596  REAL*8 :: max_cfl
597 
598  LOGICAL :: tend1
599  CHARACTER(LEN=80) :: card
600 
601  INTEGER :: j,k
602 
603  INTEGER :: dot_idx
604 
605  CHARACTER(LEN=3) :: check_file
606 
607  LOGICAL :: lexist
608 
609  CHARACTER(LEN=15) :: chara
610 
611  INTEGER :: ios
612 
613  REAL*8 :: expA , expB , Tc
614 
615 
616  OPEN(input_unit,file=input_file,status='old')
617 
618  ! ------- READ run_parameters NAMELIST -----------------------------------
619  READ(input_unit, run_parameters,iostat=ios )
620 
621  IF ( ios .NE. 0 ) THEN
622 
623  WRITE(*,*) 'IOSTAT=',ios
624  WRITE(*,*) 'ERROR: problem with namelist RUN_PARAMETERS'
625  WRITE(*,*) 'Please check the input file'
626  stop
627 
628  ELSE
629 
630  rewind(input_unit)
631 
632  END IF
633 
634  ! ------- READ newrun_parameters NAMELIST --------------------------------
635  READ(input_unit,newrun_parameters,iostat=ios)
636 
637  IF ( ios .NE. 0 ) THEN
638 
639  WRITE(*,*) 'IOSTAT=',ios
640  WRITE(*,*) 'ERROR: problem with namelist NEWRUN_PARAMETERS'
641  WRITE(*,*) 'Please check the input file'
642  stop
643 
644  ELSE
645 
646  rewind(input_unit)
647 
648  END IF
649 
650  IF ( ( comp_cells_x .EQ. 1 ) .OR. ( comp_cells_y .EQ. 1 ) ) THEN
651 
652  WRITE(*,*) '----- 1D SIMULATION -----'
653 
654  ELSE
655 
656  WRITE(*,*) '----- 2D SIMULATION -----'
657 
658  END IF
659 
660  IF ( ( .NOT. liquid_flag ) .AND. ( .NOT. gas_flag ) ) THEN
661 
662  WRITE(*,*) 'IOSTAT=',ios
663  WRITE(*,*) 'ERROR: problem with namelist NEWRUN_PARAMETERS'
664  WRITE(*,*) 'One of these parameters must be set to .TRUE.'
665  WRITE(*,*) 'LIQUID_FLAG',liquid_flag
666  WRITE(*,*) 'GAS_FLAG',liquid_flag
667  WRITE(*,*) 'Please check the input file'
668  stop
669 
670  END IF
671 
672  ! ------- READ gas_transport_parameters NAMELIST --------------------------
673 
674  READ(input_unit, gas_transport_parameters,iostat=ios)
675 
676  IF ( ios .NE. 0 ) THEN
677 
678  WRITE(*,*) 'IOSTAT=',ios
679  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
680  WRITE(*,*) 'Please check the input file'
681  stop
682 
683  ELSE
684 
685  rewind(input_unit)
686 
687  END IF
688 
689  IF ( sp_heat_a .EQ. -1.d0 ) THEN
690 
691  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
692  WRITE(*,*) 'SP_HEAT_a =' , sp_heat_a
693  WRITE(*,*) 'Please check the input file'
694  stop
695 
696  END IF
697 
698  IF ( sp_gas_const_a .EQ. -1.d0 ) THEN
699 
700  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
701  WRITE(*,*) 'SP_GAS_CONST_a =' , sp_gas_const_a
702  WRITE(*,*) 'Please check the input file'
703  stop
704 
705  END IF
706 
707  IF ( kin_visc_a .EQ. -1.d0 ) THEN
708 
709  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
710  WRITE(*,*) 'KIN_VISC_CONST_a =' , kin_visc_a
711  WRITE(*,*) 'Please check the input file'
712  stop
713 
714  ELSE
715 
716  IF ( gas_flag ) THEN
717 
718  WRITE(*,*) 'CARRIER PHASE: gas'
719  WRITE(*,*) 'Carrier phase kinematic viscosity:',kin_visc_a
720  kin_visc_c = kin_visc_a
721 
722  END IF
723 
724  END IF
725 
726  IF ( pres .EQ. -1.d0 ) THEN
727 
728  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
729  WRITE(*,*) 'pres =' , pres
730  WRITE(*,*) 'Please check the input file'
731  stop
732 
733  END IF
734 
735  IF ( t_ambient .EQ. -1.d0 ) THEN
736 
737  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
738  WRITE(*,*) 'T_ambient =' , t_ambient
739  WRITE(*,*) 'Please check the input file'
740  stop
741 
742  END IF
743 
744  IF ( ( .NOT. gas_flag ) .AND. ( liquid_flag .AND. entrainment_flag ) ) THEN
745 
746  WRITE(*,*) 'ERROR: problem with namelist GAS_TRANSPORT_PARAMETERS'
747  WRITE(*,*) 'LIQUID_FLAG',liquid_flag
748  WRITE(*,*) 'ENTRAINMENT_FLAG =' , entrainment_flag
749  WRITE(*,*) 'Please check the input file'
750  stop
751 
752  END IF
753 
754  rho_a_amb = pres / ( sp_gas_const_a * t_ambient )
755  WRITE(*,*) 'Ambient density = ',rho_a_amb,' (kg/m3)'
756 
757  ! ------- READ liquid_transport_parameters NAMELIST -------------------------
758 
759  n_vars = 4
760 
761  IF ( liquid_flag ) THEN
762 
763  IF ( gas_flag ) n_vars = n_vars + 1
764 
765  READ(input_unit, liquid_transport_parameters,iostat=ios)
766 
767  IF ( ios .NE. 0 ) THEN
768 
769  WRITE(*,*) 'IOSTAT=',ios
770  WRITE(*,*) 'ERROR: problem with namelist LIQUID_TRANSPORT_PARAMETERS'
771  WRITE(*,*) 'Please check the input file'
772  stop
773 
774  ELSE
775 
776  rewind(input_unit)
777 
778  END IF
779 
780  IF ( sp_heat_l .EQ. -1.d0 ) THEN
781 
782  WRITE(*,*) 'ERROR: problem with namelist LIQUID_TRANSPORT_PARAMETERS'
783  WRITE(*,*) 'SP_HEAT_l =' , sp_heat_l
784  WRITE(*,*) 'Please check the input file'
785  stop
786 
787  END IF
788 
789  IF ( rho_l .EQ. -1.d0 ) THEN
790 
791  WRITE(*,*) 'ERROR: problem with namelist LIQUID_TRANSPORT_PARAMETERS'
792  WRITE(*,*) 'RHO_L =' , rho_l
793  WRITE(*,*) 'Please check the input file'
794  stop
795 
796  END IF
797 
798  IF ( kin_visc_l .EQ. -1.d0 ) THEN
799 
800  WRITE(*,*) 'ERROR: problem with namelist LIQUID_TRANSPORT_PARAMETERS'
801  WRITE(*,*) 'KIN_VISC_L =' , kin_visc_l
802  WRITE(*,*) 'Please check the input file'
803  stop
804 
805  ELSE
806 
807  IF ( .NOT. gas_flag ) THEN
808 
809  WRITE(*,*) 'CARRIER PHASE: liquid'
810  WRITE(*,*) 'Carrier phase kinematic viscosity:',kin_visc_l
811  kin_visc_c = kin_visc_l
812 
813  END IF
814 
815 
816  END IF
817 
818  END IF
819 
820  ! ------- READ solid_transport_parameters NAMELIST --------------------------
821 
822  READ(input_unit, solid_transport_parameters,iostat=ios)
823 
824  IF ( ios .NE. 0 ) THEN
825 
826  WRITE(*,*) 'IOSTAT=',ios
827  WRITE(*,*) 'ERROR: problem with namelist SOLID_TRANSPORT_PARAMETERS'
828  WRITE(*,*) 'Please check the input file'
829  stop
830 
831  ELSE
832 
833  rewind(input_unit)
834 
835  END IF
836 
837  IF ( n_solid .LT. 1 ) THEN
838 
839  WRITE(*,*) 'ERROR: problem with namelist SOLID_TRANSPORT_PARAMETERS'
840  WRITE(*,*) 'n_solid =' , n_solid
841  WRITE(*,*) 'Please check the input file'
842  stop
843 
844  END IF
845 
846  IF ( any(rho0_s(1:n_solid) .EQ. -1.d0 ) ) THEN
847 
848  WRITE(*,*) 'ERROR: problem with namelist SOLID_TRANSPORT_PARAMETERS'
849  WRITE(*,*) 'RHO_s =' , rho0_s(1:n_solid)
850  WRITE(*,*) 'Please check the input file'
851  stop
852 
853  END IF
854 
855  IF ( any(sp_heat0_s(1:n_solid) .EQ. -1.d0 ) ) THEN
856 
857  WRITE(*,*) 'ERROR: problem with namelist SOLID_TRANSPORT_PARAMETERS'
858  WRITE(*,*) 'RHO_s =' , rho0_s(1:n_solid)
859  WRITE(*,*) 'Please check the input file'
860  stop
861 
862  END IF
863 
864  IF ( any(erosion_coeff0(1:n_solid) .LT. 0.d0 ) ) THEN
865 
866  WRITE(*,*) 'ERROR: problem with namelist SOLID_TRANSPORT_PARAMETERS'
867  WRITE(*,*) 'EROSION_COEFF =' , erosion_coeff
868  WRITE(*,*) 'Please check the input file'
869  stop
870 
871  END IF
872 
873  IF ( t_s_substrate .LT. 0.d0 ) THEN
874 
875  WRITE(*,*) 'ERROR: problem with namelist SOLID_TRANSPORT_PARAMETERS'
876  WRITE(*,*) 'T_s_substrate =' , t_s_substrate
877  WRITE(*,*) 'Please check the input file'
878  stop
879 
880  END IF
881 
882  n_vars = n_vars + n_solid
883  n_eqns = n_vars
884 
885  ALLOCATE( alphas_bcw(n_solid) , alphas_bce(n_solid) , alphas_bcs(n_solid) , &
886  alphas_bcn(n_solid) )
887 
888  alphas_bcw(1:n_solid)%flag = -1
889  alphas_bce(1:n_solid)%flag = -1
890  alphas_bcs(1:n_solid)%flag = -1
891  alphas_bcn(1:n_solid)%flag = -1
892 
893 
894  ALLOCATE( bcw(n_vars) , bce(n_vars) , bcs(n_vars) , bcn(n_vars) )
895 
896  bcw(1:n_vars)%flag = -1
897  bce(1:n_vars)%flag = -1
898  bcs(1:n_vars)%flag = -1
899  bcn(1:n_vars)%flag = -1
900 
901  ALLOCATE( rho_s(n_solid) , diam_s(n_solid) , sp_heat_s(n_solid) )
902 
903  ALLOCATE( alphas_init(n_solid) , sed_vol_perc(n_solid) )
904 
905  ALLOCATE( erosion_coeff(n_solid) )
906 
907  rho_s(1:n_solid) = rho0_s(1:n_solid)
908  diam_s(1:n_solid) = diam0_s(1:n_solid)
909  sp_heat_s(1:n_solid) = sp_heat0_s(1:n_solid)
910  erosion_coeff(1:n_solid) = erosion_coeff0(1:n_solid)
911 
912  ALLOCATE( deposit( comp_cells_x , comp_cells_y , n_solid ) )
913 
914  deposit(1:comp_cells_x,1:comp_cells_y,1:n_solid ) = 0.d0
915 
916  IF ( restart ) THEN
917 
918  ! ------- READ restart_parameters NAMELIST --------------------------
919  READ(input_unit,restart_parameters,iostat=ios)
920 
921  IF ( ios .NE. 0 ) THEN
922 
923  WRITE(*,*) 'IOSTAT=',ios
924  WRITE(*,*) 'ERROR: problem with namelist RESTART_PARAMETERS'
925  WRITE(*,*) 'Please check the input file'
926  stop
927 
928  ELSE
929 
930  dot_idx = scan(restart_file, ".", .true.)
931 
932  check_file = restart_file(dot_idx+1:dot_idx+3)
933 
934  IF ( check_file .EQ. 'asc' ) THEN
935 
936  IF ( ( any(sed_vol_perc .LT. 0.d0 ) ) .OR. ( any(sed_vol_perc .GT. 100.d0 ) ) ) &
937  THEN
938 
939  WRITE(*,*) 'ERROR: problem with namelist RESTART_PARAMETERS'
940  WRITE(*,*) 'SED_VOL_PERC =' , sed_vol_perc
941  stop
942 
943  END IF
944 
945  alphas_init = 1.d-2 * sed_vol_perc
946  WRITE(*,*) 'INITIAL VOLUME FRACTION OF SOLIDS:', alphas_init
947  rewind(input_unit)
948 
949 
950  IF ( t_init*t_ambient .EQ. 0.d0 ) THEN
951 
952  WRITE(*,*) 'T_init=',t_init
953  WRITE(*,*) 'T_ambient=',t_ambient
954  WRITE(*,*) 'Add the variables to the namelist RESTART_PARAMETERS'
955  stop
956 
957  END IF
958 
959  END IF
960 
961  END IF
962 
963  ELSE
964 
965  IF ( riemann_flag ) THEN
966 
967  READ(input_unit,left_state,iostat=ios)
968 
969  IF ( ios .NE. 0 ) THEN
970 
971  WRITE(*,*) 'IOSTAT=',ios
972  WRITE(*,*) 'ERROR: problem with namelist LEFT_PARAMETERS'
973  WRITE(*,*) 'Please check the input file'
974  stop
975 
976  ELSE
977 
978  rewind(input_unit)
979 
980  IF ( t_w .EQ. -1.d0 ) THEN
981 
982  WRITE(*,*) 'ERROR: problem with namelist LEFT_PARAMETERS'
983  WRITE(*,*) 'Initial temperature T_W not defined'
984  stop
985 
986  END IF
987 
988  END IF
989 
990  READ(input_unit,right_state,iostat=ios)
991 
992  IF ( ios .NE. 0 ) THEN
993 
994  WRITE(*,*) 'IOSTAT=',ios
995  WRITE(*,*) 'ERROR: problem with namelist RIGHT_PARAMETERS'
996  WRITE(*,*) 'Please check the input file'
997  stop
998 
999  ELSE
1000 
1001  rewind(input_unit)
1002 
1003  IF ( t_e .EQ. -1.d0 ) THEN
1004 
1005  WRITE(*,*) 'ERROR: problem with namelist RIGHT_PARAMETERS'
1006  WRITE(*,*) 'Initial temperature T_E not defined'
1007  stop
1008 
1009  END IF
1010 
1011 
1012  END IF
1013 
1014  ELSE
1015 
1016 
1017  END IF
1018 
1019  END IF
1020 
1021  ! ------- READ numeric_parameters NAMELIST ----------------------------------
1022 
1023  READ(input_unit,numeric_parameters)
1024 
1025  IF ( ios .NE. 0 ) THEN
1026 
1027  WRITE(*,*) 'IOSTAT=',ios
1028  WRITE(*,*) 'ERROR: problem with namelist NUMERIC_PARAMETERS'
1029  WRITE(*,*) 'Please check the input file'
1030  stop
1031 
1032  ELSE
1033 
1034  rewind(input_unit)
1035 
1036  END IF
1037 
1038  IF ( ( solver_scheme .NE. 'LxF' ) .AND. ( solver_scheme .NE. 'KT' ) .AND. &
1039  ( solver_scheme .NE. 'GFORCE' ) .AND. ( solver_scheme .NE. 'UP' ) ) THEN
1040 
1041  WRITE(*,*) 'WARNING: no correct solver scheme selected',solver_scheme
1042  WRITE(*,*) 'Chose between: LxF, GFORCE or KT'
1043  stop
1044 
1045  END IF
1046 
1047  IF ( ( solver_scheme.EQ.'LxF' ) .OR. ( solver_scheme.EQ.'GFORCE' ) ) THEN
1048 
1049  max_cfl = 1.0
1050 
1051  ELSE
1052 
1053  IF ( ( comp_cells_x .EQ. 1 ) .OR. ( comp_cells_y .EQ. 1 ) ) THEN
1054 
1055  max_cfl = 0.50d0
1056 
1057  ELSE
1058 
1059  max_cfl = 0.25d0
1060 
1061  END IF
1062 
1063  END IF
1064 
1065 
1066  IF ( ( cfl .GT. max_cfl ) .OR. ( cfl .LT. 0.d0 ) ) THEN
1067 
1068  WRITE(*,*) 'WARNING: wrong value of cfl ',cfl
1069  WRITE(*,*) 'Choose a value between 0.0 and ',max_cfl
1070  READ(*,*)
1071 
1072  END IF
1073 
1074  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'Limiters',limiter(1:n_vars)
1075 
1076  limiter(n_vars+1) = limiter(2)
1077  limiter(n_vars+2) = limiter(3)
1078 
1079  IF ( ( maxval(limiter(1:n_vars)) .GT. 3 ) .OR. &
1080  ( minval(limiter(1:n_vars)) .LT. 0 ) ) THEN
1081 
1082  WRITE(*,*) 'WARNING: wrong limiter ',limiter(1:n_vars)
1083  WRITE(*,*) 'Choose among: none, minmod,superbee,van_leer'
1084  stop
1085 
1086  END IF
1087 
1088 
1089  WRITE(*,*) 'Linear reconstruction and b. c. applied to variables:'
1090  WRITE(*,*) 'h,hu,hv,T,alphas'
1091 
1092  IF ( ( reconstr_coeff .GT. 1.0d0 ) .OR. ( reconstr_coeff .LT. 0.d0 ) ) THEN
1093 
1094  WRITE(*,*) 'WARNING: wrong value of reconstr_coeff ',reconstr_coeff
1095  WRITE(*,*) 'Change the value between 0.0 and 1.0 in the input file'
1096  READ(*,*)
1097 
1098  END IF
1099 
1100  ! ------- READ boundary_conditions NAMELISTS --------------------------------
1101 
1102  IF ( comp_cells_x .GT. 1 ) THEN
1103 
1104  ! --------- West boundary conditions -------------------------------------
1105 
1106  READ(input_unit,west_boundary_conditions,iostat=ios)
1107 
1108  IF ( ios .NE. 0 ) THEN
1109 
1110  WRITE(*,*) 'IOSTAT=',ios
1111  WRITE(*,*) 'ERROR: problem with namelist WEST_BOUNDARY_CONDITIONS'
1112  WRITE(*,*) 'Please check the input file'
1113  stop
1114 
1115  ELSE
1116 
1117  rewind(input_unit)
1118 
1119  END IF
1120 
1121  IF ( ( h_bcw%flag .EQ. -1 ) ) THEN
1122 
1123  WRITE(*,*) 'ERROR: problem with namelist WEST_BOUNDARY_CONDITIONS'
1124  WRITE(*,*) 'B.C. for h not set properly'
1125  WRITE(*,*) 'Please check the input file'
1126  stop
1127 
1128  END IF
1129 
1130  IF ( comp_cells_x .GT. 1 ) THEN
1131 
1132  IF ( hu_bcw%flag .EQ. -1 ) THEN
1133 
1134  WRITE(*,*) 'ERROR: problem with namelist WEST_BOUNDARY_CONDITIONS'
1135  WRITE(*,*) 'B.C. for hu not set properly'
1136  WRITE(*,*) 'Please check the input file'
1137  stop
1138 
1139  ELSE
1140 
1141  ! hu_bcW%flag = 1
1142  ! hu_bcW%value = 0.D0
1143 
1144  END IF
1145 
1146  END IF
1147 
1148  IF ( comp_cells_y .GT. 1 ) THEN
1149 
1150  IF ( hv_bcw%flag .EQ. -1 ) THEN
1151 
1152  WRITE(*,*) 'ERROR: problem with namelist WEST_BOUNDARY_CONDITIONS'
1153  WRITE(*,*) 'B.C. for hv not set properly'
1154  WRITE(*,*) 'Please check the input file'
1155  stop
1156 
1157  ELSE
1158 
1159  hv_bcw%flag = 1
1160  hv_bcw%value = 0.d0
1161 
1162  END IF
1163 
1164  END IF
1165 
1166  IF ( any(alphas_bcw%flag .EQ. -1 ) ) THEN
1167 
1168  WRITE(*,*) 'ERROR: problem with namelist WEST_BOUNDARY_CONDITIONS'
1169  WRITE(*,*) 'B.C. for sediment conentration not set properly'
1170  WRITE(*,*) 'Please check the input file'
1171  WRITE(*,*) 'alphas_bcW'
1172  WRITE(*,*) alphas_bcw
1173  stop
1174 
1175  END IF
1176 
1177  IF ( t_bcw%flag .EQ. -1 ) THEN
1178 
1179  WRITE(*,*) 'ERROR: problem with namelist WEST_BOUNDARY_CONDITIONS'
1180  WRITE(*,*) 'B.C. for temperature not set properly'
1181  WRITE(*,*) 'Please check the input file'
1182  stop
1183 
1184  END IF
1185 
1186  ! set the approriate boundary conditions
1187 
1188  bcw(1) = h_bcw
1189  bcw(2) = hu_bcw
1190  bcw(3) = hv_bcw
1191 
1192 
1193  ! ------------- East boundary conditions --------------------------------
1194 
1195  READ(input_unit,east_boundary_conditions,iostat=ios)
1196 
1197  IF ( ios .NE. 0 ) THEN
1198 
1199  WRITE(*,*) 'IOSTAT=',ios
1200  WRITE(*,*) 'ERROR: problem with namelist EAST_BOUNDARY_CONDITIONS'
1201  WRITE(*,*) 'Please check the input file'
1202  stop
1203 
1204  ELSE
1205 
1206  rewind(input_unit)
1207 
1208  END IF
1209 
1210  IF ( ( h_bce%flag .EQ. -1 ) ) THEN
1211 
1212  WRITE(*,*) 'ERROR: problem with namelist EAST_BOUNDARY_CONDITIONS'
1213  WRITE(*,*) 'B.C. for h not set properly'
1214  WRITE(*,*) 'Please check the input file'
1215  stop
1216 
1217  END IF
1218 
1219  IF ( comp_cells_x .GT. 1 ) THEN
1220 
1221  IF ( hu_bce%flag .EQ. -1 ) THEN
1222 
1223  WRITE(*,*) 'ERROR: problem with namelist EAST_BOUNDARY_CONDITIONS'
1224  WRITE(*,*) 'B.C. for hu not set properly'
1225  WRITE(*,*) 'Please check the input file'
1226  stop
1227 
1228  END IF
1229 
1230  ELSE
1231 
1232  hu_bce%flag = 1
1233  hu_bce%value = 0.d0
1234 
1235  END IF
1236 
1237  IF ( comp_cells_y .GT. 1 ) THEN
1238 
1239  IF ( hv_bce%flag .EQ. -1 ) THEN
1240 
1241  WRITE(*,*) 'ERROR: problem with namelist EAST_BOUNDARY_CONDITIONS'
1242  WRITE(*,*) 'B.C. for hv not set properly'
1243  WRITE(*,*) 'Please check the input file'
1244  stop
1245 
1246  END IF
1247 
1248  ELSE
1249 
1250  hv_bce%flag = 1
1251  hv_bce%value = 0.d0
1252 
1253 
1254  END IF
1255 
1256  IF ( any(alphas_bce%flag .EQ. -1 ) ) THEN
1257 
1258  WRITE(*,*) 'ERROR: problem with namelist EAST_BOUNDARY_CONDITIONS'
1259  WRITE(*,*) 'B.C. for sediment concentration not set properly'
1260  WRITE(*,*) 'Please check the input file'
1261  WRITE(*,*) 'alphas_bcE'
1262  WRITE(*,*) alphas_bce
1263  stop
1264 
1265  END IF
1266 
1267  IF ( t_bce%flag .EQ. -1 ) THEN
1268 
1269  WRITE(*,*) 'ERROR: problem with namelist EAST_BOUNDARY_CONDITIONS'
1270  WRITE(*,*) 'B.C. for temperature not set properly'
1271  WRITE(*,*) 'Please check the input file'
1272  stop
1273 
1274  END IF
1275 
1276  bce(1) = h_bce
1277  bce(2) = hu_bce
1278  bce(3) = hv_bce
1279 
1280  END IF
1281 
1282  IF ( comp_cells_y .GT. 1 ) THEN
1283 
1284  ! --------------- South boundary conditions ------------------------------
1285 
1286  READ(input_unit,south_boundary_conditions,iostat=ios)
1287 
1288  IF ( ios .NE. 0 ) THEN
1289 
1290  WRITE(*,*) 'IOSTAT=',ios
1291  WRITE(*,*) 'ERROR: problem with namelist SOUTH_BOUNDARY_CONDITIONS'
1292  WRITE(*,*) 'Please check the input file'
1293  stop
1294 
1295  ELSE
1296 
1297  rewind(input_unit)
1298 
1299  END IF
1300 
1301  IF ( ( h_bcs%flag .EQ. -1 ) ) THEN
1302 
1303  WRITE(*,*) 'ERROR: problem with namelist SOUTH_BOUNDARY_CONDITIONS'
1304  WRITE(*,*) 'B.C. for h not set properly'
1305  WRITE(*,*) 'Please check the input file'
1306  stop
1307 
1308  END IF
1309 
1310  IF ( comp_cells_x .GT. 1 ) THEN
1311 
1312  IF ( hu_bcs%flag .EQ. -1 ) THEN
1313 
1314  WRITE(*,*) 'ERROR: problem with namelist SOUTH_BOUNDARY_CONDITIONS'
1315  WRITE(*,*) 'B.C. for hu not set properly'
1316  WRITE(*,*) 'Please check the input file'
1317  stop
1318 
1319  END IF
1320 
1321  ELSE
1322 
1323  hu_bcs%flag = 1
1324  hu_bcs%value = 0.d0
1325 
1326  END IF
1327 
1328  IF ( comp_cells_y .GT. 1 ) THEN
1329 
1330  IF ( hv_bcs%flag .EQ. -1 ) THEN
1331 
1332  WRITE(*,*) 'ERROR: problem with namelist SOUTH_BOUNDARY_CONDITIONS'
1333  WRITE(*,*) 'B.C. for hv not set properly'
1334  WRITE(*,*) 'Please check the input file'
1335  stop
1336 
1337  END IF
1338 
1339  ELSE
1340 
1341  hv_bcs%flag = 1
1342  hv_bcs%value = 0.d0
1343 
1344  END IF
1345 
1346  IF ( any(alphas_bcs%flag .EQ. -1 ) ) THEN
1347 
1348  WRITE(*,*) 'ERROR: problem with namelist SOUTH_BOUNDARY_CONDITIONS'
1349  WRITE(*,*) 'B.C. for sediment concentrations not set properly'
1350  WRITE(*,*) 'Please check the input file'
1351  WRITE(*,*) 'alphas_bcS'
1352  WRITE(*,*) alphas_bcs
1353  stop
1354 
1355  END IF
1356 
1357  IF ( t_bcs%flag .EQ. -1 ) THEN
1358 
1359  WRITE(*,*) 'ERROR: problem with namelist SOUTH_BOUNDARY_CONDITIONS'
1360  WRITE(*,*) 'B.C. for temperature not set properly'
1361  WRITE(*,*) 'Please check the input file'
1362  stop
1363 
1364  END IF
1365 
1366  bcs(1) = h_bcs
1367  bcs(2) = hu_bcs
1368  bcs(3) = hv_bcs
1369 
1370  ! ---------------- North boundary conditions ----------------------------
1371 
1372  READ(input_unit,north_boundary_conditions,iostat=ios)
1373 
1374  IF ( ios .NE. 0 ) THEN
1375 
1376  WRITE(*,*) 'IOSTAT=',ios
1377  WRITE(*,*) 'ERROR: problem with namelist NORTH_BOUNDARY_CONDITIONS'
1378  WRITE(*,*) 'Please check the input file'
1379  stop
1380 
1381  ELSE
1382 
1383  rewind(input_unit)
1384 
1385  END IF
1386 
1387 
1388  IF ( ( h_bcn%flag .EQ. -1 ) ) THEN
1389 
1390  WRITE(*,*) 'ERROR: problem with namelist NORTH_BOUNDARY_CONDITIONS'
1391  WRITE(*,*) 'B.C. for h not set properly'
1392  WRITE(*,*) 'Please check the input file'
1393  stop
1394 
1395  END IF
1396 
1397 
1398  IF ( comp_cells_x .GT. 1 ) THEN
1399 
1400  IF ( hu_bcn%flag .EQ. -1 ) THEN
1401 
1402  WRITE(*,*) 'ERROR: problem with namelist NORTH_BOUNDARY_CONDITIONS'
1403  WRITE(*,*) 'B.C. for hu not set properly'
1404  WRITE(*,*) 'Please check the input file'
1405  stop
1406 
1407  END IF
1408 
1409  ELSE
1410 
1411  hu_bcn%flag = 1
1412  hu_bcn%value = 0.d0
1413 
1414  END IF
1415 
1416  IF ( comp_cells_y .GT. 1 ) THEN
1417 
1418  IF ( hv_bcn%flag .EQ. -1 ) THEN
1419 
1420  WRITE(*,*) 'ERROR: problem with namelist NORTH_BOUNDARY_CONDITIONS'
1421  WRITE(*,*) 'B.C. for hv not set properly'
1422  WRITE(*,*) 'Please check the input file'
1423  stop
1424 
1425  END IF
1426 
1427  ELSE
1428 
1429  hv_bcn%flag = 1
1430  hv_bcn%value = 0.d0
1431 
1432  END IF
1433 
1434  IF ( any(alphas_bcn%flag .EQ. -1 ) ) THEN
1435 
1436  WRITE(*,*) 'ERROR: problem with namelist NORTH_BOUNDARY_CONDITIONS'
1437  WRITE(*,*) 'B.C. for sediment concentrations not set properly'
1438  WRITE(*,*) 'Please check the input file'
1439  WRITE(*,*) 'alphas_bcN'
1440  WRITE(*,*) alphas_bcn
1441  stop
1442 
1443  END IF
1444 
1445  IF ( t_bcn%flag .EQ. -1 ) THEN
1446 
1447  WRITE(*,*) 'ERROR: problem with namelist NORTH_BOUNDARY_CONDITIONS'
1448  WRITE(*,*) 'B.C. for temperature not set properly'
1449  WRITE(*,*) 'Please check the input file'
1450  stop
1451 
1452  END IF
1453 
1454  bcn(1) = h_bcn
1455  bcn(2) = hu_bcn
1456  bcn(3) = hv_bcn
1457 
1458  END IF
1459 
1460  bcw(4) = t_bcw
1461  bce(4) = t_bce
1462  bcs(4) = t_bcs
1463  bcn(4) = t_bcn
1464 
1465  bcw(5:4+n_solid) = alphas_bcw(1:n_solid)
1466  bce(5:4+n_solid) = alphas_bce(1:n_solid)
1467  bcs(5:4+n_solid) = alphas_bcs(1:n_solid)
1468  bcn(5:4+n_solid) = alphas_bcn(1:n_solid)
1469 
1470  WRITE(*,*) 'bcE',bce
1471 
1472  ! ------- READ expl_terms_parameters NAMELIST -------------------------------
1473 
1474  READ(input_unit, expl_terms_parameters,iostat=ios)
1475 
1476  IF ( ios .NE. 0 ) THEN
1477 
1478  WRITE(*,*) 'IOSTAT=',ios
1479  WRITE(*,*) 'ERROR: problem with namelist EXPL_TERMS_PARAMETERS'
1480  WRITE(*,*) 'Please check the input file'
1481  stop
1482 
1483  ELSE
1484 
1485  rewind(input_unit)
1486 
1487  END IF
1488 
1489  IF ( grav .EQ. -1.d0 ) THEN
1490 
1491  WRITE(*,*) 'ERROR: problem with namelist EXPL_TERMS_PARAMETERS'
1492  WRITE(*,*) 'GRAV not set properly'
1493  WRITE(*,*) 'Please check the input file'
1494  stop
1495 
1496  END IF
1497 
1498  ! ------- READ radial_source_parameters NAMELIST ----------------------------
1499 
1500  IF ( radial_source_flag ) THEN
1501 
1502  alphal_source = -1.d0
1503 
1504  ALLOCATE( alphas_source(n_solid) )
1505 
1506  READ(input_unit,radial_source_parameters,iostat=ios)
1507 
1508  IF ( ios .NE. 0 ) THEN
1509 
1510  WRITE(*,*) 'IOSTAT=',ios
1511  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1512  WRITE(*,*) 'Please check the input file'
1513  stop
1514 
1515  ELSE
1516 
1517  rewind(input_unit)
1518 
1519  IF ( t_source .EQ. -1.d0 ) THEN
1520 
1521  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1522  WRITE(*,*) 'PLEASE CHEC VALUE OF T_SOURCE',t_source
1523  stop
1524 
1525  END IF
1526 
1527  IF ( h_source .EQ. -1.d0 ) THEN
1528 
1529  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1530  WRITE(*,*) 'PLEASE CHEC VALUE OF H_SOURCE',h_source
1531  stop
1532 
1533  END IF
1534 
1535  IF ( r_source .EQ. -1.d0 ) THEN
1536 
1537  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1538  WRITE(*,*) 'PLEASE CHEC VALUE OF R_SOURCE',r_source
1539  stop
1540 
1541  END IF
1542 
1543  IF ( vel_source .EQ. -1.d0 ) THEN
1544 
1545  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1546  WRITE(*,*) 'PLEASE CHEC VALUE OF VEL_SOURCE',vel_source
1547  stop
1548 
1549  END IF
1550 
1551  IF ( ( x_source - r_source ) .LE. x0 + cell_size ) THEN
1552 
1553  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1554  WRITE(*,*) 'SOURCE TOO LARGE'
1555  WRITE(*,*) ' x_source - radius ',x_source-r_source
1556  stop
1557 
1558  END IF
1559 
1560  IF ( ( x_source + r_source ) .GE. x0+(comp_cells_x-1)*cell_size ) THEN
1561 
1562  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1563  WRITE(*,*) 'SOURCE TOO LARGE'
1564  WRITE(*,*) ' x_source + radius ',x_source+r_source
1565  stop
1566 
1567  END IF
1568 
1569  IF ( ( y_source - r_source ) .LE. y0 + cell_size ) THEN
1570 
1571  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1572  WRITE(*,*) 'SOURCE TOO LARGE'
1573  WRITE(*,*) ' y_source - radius ',y_source-r_source
1574  stop
1575 
1576  END IF
1577 
1578  IF ( ( y_source + r_source ) .GE. y0 +(comp_cells_y-1)*cell_size ) THEN
1579 
1580  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1581  WRITE(*,*) 'SOURCE TOO LARGE'
1582  WRITE(*,*) ' y_source + radius ',y_source+r_source
1583  stop
1584 
1585  END IF
1586 
1587  IF ( gas_flag .AND. liquid_flag ) THEN
1588 
1589  IF ( alphal_source .LT. 0.d0 ) THEN
1590 
1591  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1592  WRITE(*,*) 'PLEASE CHECK VALUE OF ALPHAL_SOURCE',alphal_source
1593  stop
1594 
1595  END IF
1596 
1597  END IF
1598 
1599  IF ( any(time_param .LT. 0.d0 ) ) THEN
1600 
1601  WRITE(*,*)
1602  WRITE(*,*) 'WARNING: problem with namelist RADIAL_SOURCEPARAMETERS'
1603  WRITE(*,*) 'time_param =' , time_param
1604  time_param(1) = t_end
1605  time_param(2) = t_end
1606  time_param(3) = 0.d0
1607  time_param(4) = t_end
1608  WRITE(*,*) 'CHANGED TO time_param =',time_param
1609  WRITE(*,*) 'Radial source now constant in time'
1610  WRITE(*,*)
1611 
1612  ELSE
1613 
1614  IF ( time_param(2) .GT. time_param(1) ) THEN
1615 
1616  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCEPARAMETERS'
1617  WRITE(*,*) 'time_param(1),time_param(2) =' , time_param(1:2)
1618  WRITE(*,*) 'time_param(1) must be larger than time_param(2)'
1619  stop
1620 
1621  END IF
1622 
1623  IF ( time_param(3) .GT. ( 0.5d0*time_param(2) ) ) THEN
1624 
1625  WRITE(*,*) 'ERROR: problem with namelist RADIAL_SOURCE_PARAMETERS'
1626  WRITE(*,*) 'time_param(3) =', time_param(3)
1627  WRITE(*,*) 'time_param(3) must be smaller than 0.5*time_param(2)'
1628  stop
1629 
1630  END IF
1631 
1632 
1633  END IF
1634 
1635  END IF
1636 
1637  END IF
1638 
1639 
1640 
1641 
1642  ! ------- READ collapsing_volume_parameters NAMELIST ------------------------
1643 
1644  IF ( collapsing_volume_flag ) THEN
1645 
1646  ALLOCATE( alphas_collapse(n_solid) )
1647 
1648  READ(input_unit,collapsing_volume_parameters,iostat=ios)
1649 
1650  IF ( ios .NE. 0 ) THEN
1651 
1652  WRITE(*,*) 'IOSTAT=',ios
1653  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1654  WRITE(*,*) 'Please check the input file'
1655  stop
1656 
1657  ELSE
1658 
1659  rewind(input_unit)
1660 
1661  IF ( t_collapse .EQ. -1.d0 ) THEN
1662 
1663  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1664  WRITE(*,*) 'PLEASE CHECK VALUE OF T_COLLAPSE',t_collapse
1665  stop
1666 
1667  END IF
1668 
1669  IF ( h_collapse .EQ. -1.d0 ) THEN
1670 
1671  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1672  WRITE(*,*) 'PLEASE CHECK VALUE OF H_COLLAPSE',h_collapse
1673  stop
1674 
1675  END IF
1676 
1677  IF ( r_collapse .EQ. -1.d0 ) THEN
1678 
1679  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1680  WRITE(*,*) 'PLEASE CHECK VALUE OF R_COLLAPSE',r_collapse
1681  stop
1682 
1683  END IF
1684 
1685  IF ( ( x_collapse - r_collapse ) .LE. x0 + cell_size ) THEN
1686 
1687  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1688  WRITE(*,*) 'COLLAPSING VOLUME TOO LARGE'
1689  WRITE(*,*) ' x_collapse - radius ',x_collapse-r_collapse
1690  stop
1691 
1692  END IF
1693 
1694  IF ( ( x_collapse + r_collapse ) .GE. x0+(comp_cells_x-1)*cell_size ) THEN
1695 
1696  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1697  WRITE(*,*) 'COLLAPSING VOLUME TOO LARGE'
1698  WRITE(*,*) ' x_collapse + radius ',x_collapse+r_collapse
1699  stop
1700 
1701  END IF
1702 
1703  IF ( ( y_collapse - r_collapse ) .LE. y0 + cell_size ) THEN
1704 
1705  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1706  WRITE(*,*) 'COLLAPSING VOLUME TOO LARGE'
1707  WRITE(*,*) ' y_collapse - radius ',y_collapse-r_collapse
1708  stop
1709 
1710  END IF
1711 
1712  IF ( ( y_collapse + r_collapse ) .GE. y0 +(comp_cells_y-1)*cell_size ) THEN
1713 
1714  WRITE(*,*) 'ERROR: problem with namelist COLLAPSING_VOLUME_PARAMETERS'
1715  WRITE(*,*) 'COLLAPSING VOLUME TOO LARGE'
1716  WRITE(*,*) ' y_collapse + radius ',y_collapse+r_collapse
1717  stop
1718 
1719  END IF
1720 
1721  END IF
1722 
1723  END IF
1724 
1725 
1726  ! ------- READ temperature_parameters NAMELIST ------------------------------
1727 
1728  READ(input_unit, temperature_parameters,iostat=ios)
1729 
1730  IF ( ios .NE. 0 ) THEN
1731 
1732  WRITE(*,*) 'IOSTAT=',ios
1733  WRITE(*,*) 'ERROR: problem with namelist TEMPERATURE_PARAMETERS'
1734  WRITE(*,*) 'Please check the input file'
1735  stop
1736 
1737  ELSE
1738 
1739  rewind(input_unit)
1740 
1741  END IF
1742 
1743  ! ------- READ rheology_parameters NAMELIST ---------------------------------
1744 
1745  IF ( rheology_flag ) THEN
1746 
1747  READ(input_unit, rheology_parameters,iostat=ios)
1748 
1749  IF ( ios .NE. 0 ) THEN
1750 
1751  WRITE(*,*) 'IOSTAT=',ios
1752  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1753  WRITE(*,*) 'Please check the input file'
1754  stop
1755 
1756  ELSE
1757 
1758  rewind(input_unit)
1759 
1760  END IF
1761 
1762  IF ( rheology_model .EQ. 0 ) THEN
1763 
1764  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1765  WRITE(*,*) 'RHEOLOGY_FLAG' , rheology_flag , 'RHEOLOGY_MODEL =' , &
1766  rheology_model
1767  WRITE(*,*) 'Please check the input file'
1768  stop
1769 
1770  ELSEIF ( rheology_model .EQ. 1 ) THEN
1771 
1772  IF ( ( mu .EQ. -1.d0 ) .AND. ( xi .EQ. -1.d0 ) ) THEN
1773 
1774  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1775  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
1776  WRITE(*,*) 'MU =' , mu ,' XI =' , xi
1777  WRITE(*,*) 'Please check the input file'
1778  stop
1779 
1780  END IF
1781 
1782  IF ( ( t_ref .NE. -1.d0 ) .OR. ( nu_ref .NE. -1.d0 ) .OR. &
1783  ( visc_par .NE. -1.d0 ) .OR. ( tau .NE. -1.d0 ) ) THEN
1784 
1785  WRITE(*,*) 'WARNING: parameters not used in RHEOLOGY_PARAMETERS'
1786  IF ( t_ref .NE. -1.d0 ) WRITE(*,*) 'T_ref =',t_ref
1787  IF ( nu_ref .NE. -1.d0 ) WRITE(*,*) 'nu_ref =',nu_ref
1788  IF ( visc_par .NE. -1.d0 ) WRITE(*,*) 'visc_par =',visc_par
1789  IF ( tau .NE. -1.d0 ) WRITE(*,*) 'tau =',tau
1790  WRITE(*,*) 'Press ENTER to continue'
1791  READ(*,*)
1792 
1793  END IF
1794 
1795  ELSEIF ( rheology_model .EQ. 2 ) THEN
1796 
1797  IF ( tau .EQ. -1.d0 ) THEN
1798 
1799  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1800  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
1801  WRITE(*,*) 'TAU =' , tau
1802  WRITE(*,*) 'Please check the input file'
1803  stop
1804 
1805  END IF
1806 
1807  IF ( ( t_ref .NE. -1.d0 ) .OR. ( nu_ref .NE. -1.d0 ) .OR. &
1808  ( visc_par .NE. -1.d0 ) .OR. ( mu .NE. -1.d0 ) .OR. &
1809  ( xi .NE. -1.d0 ) ) THEN
1810 
1811  WRITE(*,*) 'WARNING: parameters not used in RHEOLOGY_PARAMETERS'
1812  IF ( t_ref .NE. -1.d0 ) WRITE(*,*) 'T_ref =',t_ref
1813  IF ( nu_ref .NE. -1.d0 ) WRITE(*,*) 'nu_ref =',nu_ref
1814  IF ( visc_par .NE. -1.d0 ) WRITE(*,*) 'visc_par =',visc_par
1815  IF ( mu .NE. -1.d0 ) WRITE(*,*) 'mu =',mu
1816  IF ( xi .NE. -1.d0 ) WRITE(*,*) 'xi =',xi
1817  WRITE(*,*) 'Press ENTER to continue'
1818  READ(*,*)
1819 
1820 
1821  END IF
1822 
1823  ELSEIF ( rheology_model .EQ. 3 ) THEN
1824 
1825  IF ( nu_ref .EQ. -1.d0 ) THEN
1826 
1827  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1828  WRITE(*,*) 'NU_REF =' , nu_ref
1829  WRITE(*,*) 'Please check the input file'
1830  stop
1831 
1832  END IF
1833 
1834  IF ( visc_par .EQ. -1.d0 ) THEN
1835 
1836  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1837  WRITE(*,*) 'VISC_PAR =' , visc_par
1838  WRITE(*,*) 'Please check the input file'
1839  stop
1840 
1841  ELSEIF ( visc_par .EQ. 0.d0 ) THEN
1842 
1843  WRITE(*,*) 'WARNING: temperature and momentum uncoupled'
1844  WRITE(*,*) 'VISC_PAR =' , visc_par
1845  WRITE(*,*) 'Press ENTER to continue'
1846  READ(*,*)
1847 
1848  ELSE
1849 
1850  IF ( t_ref .EQ. -1.d0 ) THEN
1851 
1852  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1853  WRITE(*,*) 'T_REF =' , t_ref
1854  WRITE(*,*) 'Please check the input file'
1855  stop
1856 
1857  END IF
1858 
1859  END IF
1860 
1861  IF ( ( mu .NE. -1.d0 ) .OR. ( xi .NE. -1.d0 ) .OR. ( tau .NE. -1.d0 ) ) &
1862  THEN
1863 
1864  WRITE(*,*) 'WARNING: parameters not used in RHEOLOGY_PARAMETERS'
1865  IF ( mu .NE. -1.d0 ) WRITE(*,*) 'mu =',mu
1866  IF ( xi .NE. -1.d0 ) WRITE(*,*) 'xi =',xi
1867  IF ( tau .NE. -1.d0 ) WRITE(*,*) 'tau =',tau
1868  WRITE(*,*) 'Press ENTER to continue'
1869  READ(*,*)
1870 
1871  END IF
1872 
1873  ELSEIF ( rheology_model .EQ. 4 ) THEN
1874 
1875  IF ( gas_flag .OR. ( .NOT. liquid_flag ) ) THEN
1876 
1877  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1878  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
1879  WRITE(*,*) 'GAS FLAG = ' , gas_flag
1880  WRITE(*,*) 'LIQUID FLAG = ' , liquid_flag
1881  stop
1882 
1883  END IF
1884 
1885  IF ( any(sed_vol_perc .EQ. -1.d0 ) ) THEN
1886 
1887  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1888  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
1889  WRITE(*,*) 'SED_VOL_PERC = ' , sed_vol_perc
1890  stop
1891 
1892  END IF
1893 
1894  IF ( alpha2 .EQ. -1.d0 ) THEN
1895 
1896  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1897  WRITE(*,*) 'ALPHA2 =' , alpha2
1898  WRITE(*,*) 'Please check the input file'
1899  stop
1900 
1901  END IF
1902 
1903  IF ( beta2 .EQ. -1.d0 ) THEN
1904 
1905  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1906  WRITE(*,*) 'BETA2 =' , beta2
1907  WRITE(*,*) 'Please check the input file'
1908  stop
1909 
1910  END IF
1911 
1912  IF ( t_ref .LE. 273.15d0 ) THEN
1913 
1914  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1915  WRITE(*,*) 'T_REF =' , t_ref
1916  WRITE(*,*) 'Please check the input file'
1917  stop
1918 
1919  END IF
1920 
1921  IF ( alpha1_ref .EQ. -1.d0 ) THEN
1922 
1923  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1924  WRITE(*,*) 'ALPHA1 =' , alpha1_ref
1925  WRITE(*,*) 'Please check the input file'
1926  stop
1927 
1928  ELSE
1929 
1930  tc = t_ref - 273.15d0
1931 
1932  IF ( tc .LT. 20.d0 ) THEN
1933 
1934  expa = 1301.d0 / ( 998.333d0 + 8.1855d0 * ( tc - 20.d0 ) &
1935  + 0.00585d0 * ( tc - 20.d0 )**2 ) - 1.30223d0
1936 
1937  alpha1_coeff = alpha1_ref / ( 1.d-3 * 10.d0**expa )
1938 
1939  ELSE
1940 
1941  expb = ( 1.3272d0 * ( 20.d0 - tc ) - 0.001053d0 * &
1942  ( tc - 20.d0 )**2 ) / ( tc + 105.0d0 )
1943 
1944  alpha1_coeff = alpha1_ref / ( 1.002d-3 * 10.d0**expb )
1945 
1946  END IF
1947 
1948  WRITE(*,*) 'alpha1 coefficient:',alpha1_coeff
1949 
1950  END IF
1951 
1952  IF ( beta1 .EQ. -1.d0 ) THEN
1953 
1954  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1955  WRITE(*,*) 'BETA1 =' , beta1
1956  WRITE(*,*) 'Please check the input file'
1957  stop
1958 
1959  END IF
1960 
1961  IF ( kappa .EQ. -1.d0 ) THEN
1962 
1963  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1964  WRITE(*,*) 'KAPPA =' , kappa
1965  WRITE(*,*) 'Please check the input file'
1966  stop
1967 
1968  END IF
1969 
1970  IF ( n_td .EQ. -1.d0 ) THEN
1971 
1972  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1973  WRITE(*,*) 'N_TD =' , n_td
1974  WRITE(*,*) 'Please check the input file'
1975  stop
1976 
1977  END IF
1978 
1979  ELSEIF ( rheology_model .EQ. 5 ) THEN
1980 
1981  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
1982  WRITE(*,*) 'Kurganov & Petrova Example 5'
1983 
1984  ELSEIF ( rheology_model .EQ. 6 ) THEN
1985 
1986  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
1987  WRITE(*,*) 'Bursik & Woods'
1988 
1989  IF ( friction_factor .LT. 0.d0 ) THEN
1990 
1991  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
1992  WRITE(*,*) 'FRICTION_FACTOR =' , friction_factor
1993  WRITE(*,*) 'Please check the input file'
1994  stop
1995 
1996  END IF
1997 
1998  ELSE
1999 
2000  WRITE(*,*) 'ERROR: problem with namelist RHEOLOGY_PARAMETERS'
2001  WRITE(*,*) 'RHEOLOGY_MODEL =' , rheology_model
2002  WRITE(*,*) 'Please check the input file'
2003  stop
2004 
2005  END IF
2006 
2007 
2008  END IF
2009 
2010  ! ---------------------------------------------------------------------------
2011 
2012  IF ( verbose_level .GE. 1 ) WRITE(*,*)
2013 
2014  WRITE(*,*) 'Searching for DEM file'
2015 
2016  INQUIRE(file='topography_dem.asc',exist=lexist)
2017 
2018  IF(lexist)THEN
2019 
2020  OPEN(2001, file='topography_dem.asc', status='old', action='read')
2021 
2022  ELSE
2023 
2024  WRITE(*,*) 'no dem file'
2025  stop
2026 
2027  ENDIF
2028 
2029  READ(2001,*) chara, ncols
2030  READ(2001,*) chara, nrows
2031  READ(2001,*) chara, xllcorner
2032  READ(2001,*) chara, yllcorner
2033  READ(2001,*) chara, cellsize
2034  READ(2001,*) chara, nodata_value
2035 
2036  ! The values read from the DEM files are associated to the center of the
2037  ! pixels. x0 is the left margin of the computational domain and has to be
2038  ! greater than the center of the first pixel.
2039  IF ( x0 - ( xllcorner + 0.5d0 * cellsize ) .LT. -1.d-10 ) THEN
2040 
2041  WRITE(*,*) 'Computational domain problem'
2042  WRITE(*,*) 'x0 < xllcorner+0.5*cellsize',x0,xllcorner+0.5d0*cellsize
2043  stop
2044 
2045  END IF
2046 
2047  ! The right margin of the computational domain should be smaller then the
2048  ! center of the last pixel
2049  IF ( x0 + ( comp_cells_x ) * cell_size .GT. &
2050  xllcorner + ( 0.5d0 + ncols ) * cellsize ) THEN
2051 
2052  WRITE(*,*) 'Computational domain problem'
2053  WRITE(*,*) 'right edge > xllcorner+ncols*cellsize', &
2054  x0+comp_cells_x*cell_size , xllcorner+(0.5d0+ncols)*cellsize
2055  stop
2056 
2057  END IF
2058 
2059  IF ( y0 - ( yllcorner+0.5d0*cellsize ) .LT. -1.d-10 ) THEN
2060 
2061  WRITE(*,*) 'Computational domain problem'
2062  WRITE(*,*) 'y0 < yllcorner+0.5*cellsize',y0,yllcorner+0.5d0*cellsize
2063  stop
2064 
2065  END IF
2066 
2067  IF ( abs( ( y0 + comp_cells_y * cell_size ) - ( yllcorner + 0.5d0 + &
2068  nrows * cellsize ) ) .LT. 1.d-10 ) THEN
2069 
2070  WRITE(*,*) 'Computational domain problem'
2071  WRITE(*,*) 'top edge > yllcorner+nrows*cellsize', &
2072  y0+comp_cells_y*cell_size , yllcorner+(0.5d0+nrows)*cellsize
2073  stop
2074 
2075  END IF
2076 
2077 
2078  WRITE(*,*) 'Reading DEM file'
2079  WRITE(*,*) 'ncols',ncols
2080  WRITE(*,*) 'nrows',nrows
2081 
2082  n_topography_profile_x = ncols
2083 
2084  n_topography_profile_y = nrows
2085 
2086  ALLOCATE( topography_profile( 3 , n_topography_profile_x , &
2087  n_topography_profile_y) )
2088 
2089  DO j=1,n_topography_profile_x
2090 
2091  topography_profile(1,j,:) = xllcorner + ( j - 0.5d0 ) * cellsize
2092 
2093  ENDDO
2094 
2095  DO k=1,n_topography_profile_y
2096 
2097  topography_profile(2,:,k) = yllcorner + ( k - 0.5d0 ) * cellsize
2098 
2099  ENDDO
2100 
2101  ! Read topography values (starting at the upper-left corner)
2102 
2103  DO k=1,n_topography_profile_y
2104 
2105  WRITE(*,fmt="(A1,A,t21,F6.2,A)",advance="NO") achar(13), &
2106  & " Percent Complete: " , &
2107  ( REAL(k) / REAL(n_topography_profile_y))*100.0, "%"
2108 
2109  READ(2001,*) topography_profile(3,:,n_topography_profile_y-k+1)
2110 
2111  ENDDO
2112 
2113  topography_profile(3,:,:) = max(0.d0,topography_profile(3,:,:))
2114 
2115  WRITE(*,*) ''
2116 
2117  CLOSE(2001)
2118 
2119 
2120  ! ------- READ runout_parameters NAMELIST -----------------------------------
2121 
2122  IF ( output_runout_flag ) THEN
2123 
2124  READ(input_unit, runout_parameters,iostat=ios)
2125 
2126  IF ( ios .NE. 0 ) THEN
2127 
2128  WRITE(*,*) 'IOSTAT=',ios
2129  WRITE(*,*) 'ERROR: problem with namelist RUNOUT_PARAMETERS'
2130  WRITE(*,*) 'Please check the input file'
2131  stop
2132 
2133  ELSE
2134 
2135  rewind(input_unit)
2136 
2137  END IF
2138 
2139  IF ( ( x0_runout .EQ. -1.d0 ) .AND. ( y0_runout .EQ. -1.d0 ) ) THEN
2140 
2141  WRITE(*,*) 'Runout reference location not defined'
2142 
2143  ELSE
2144 
2145  IF ( x0_runout .LT. x0 ) THEN
2146 
2147  WRITE(*,*) 'Computational domain problem'
2148  WRITE(*,*) 'x0_runout < x0',x0,x0_runout
2149  stop
2150 
2151  END IF
2152 
2153  IF ( x0 .GT. x0+comp_cells_x*cell_size ) THEN
2154 
2155  WRITE(*,*) 'Computational domain problem'
2156  WRITE(*,*) 'x0_runout > x0+comp_cells_x*cell_size' , x0 , &
2157  x0_runout+comp_cells_x*cell_size
2158  stop
2159 
2160  END IF
2161 
2162  IF ( y0_runout .LT. y0 ) THEN
2163 
2164  WRITE(*,*) 'Computational domain problem'
2165  WRITE(*,*) 'y0_runout < y0',y0,y0_runout
2166  stop
2167 
2168  END IF
2169 
2170  IF ( y0 .GT. y0+comp_cells_y*cell_size ) THEN
2171 
2172  WRITE(*,*) 'Computational domain problem'
2173  WRITE(*,*) 'y0_runout > y0+comp_cells_y*cell_size' , y0 , &
2174  y0_runout+comp_cells_y*cell_size
2175  stop
2176 
2177  END IF
2178 
2179  END IF
2180 
2181  runout_file = trim(run_name)//'_runout'//'.txt'
2182 
2183  OPEN(runout_unit,file=runout_file,status='unknown',form='formatted')
2184 
2185  END IF
2186 
2187 
2188  !------ search for check points --------------------------------------------
2189 
2190  tend1 = .false.
2191 
2192  WRITE(*,*) 'Searching for topography_profile'
2193 
2194  n_probes = 0
2195 
2196  probes_search: DO
2197 
2198  READ(input_unit,*, end = 300 ) card
2199 
2200  IF( trim(card) == 'PROBES_COORDS' ) THEN
2201 
2202  EXIT probes_search
2203 
2204  END IF
2205 
2206  END DO probes_search
2207 
2208 
2209  READ(input_unit,*) n_probes
2210 
2211  WRITE(*,*) 'n_probes ',n_probes
2212 
2213  ALLOCATE( probes_coords( 2 , n_probes ) )
2214 
2215  DO k = 1, n_probes
2216 
2217  READ(input_unit,*) probes_coords( 1:2 , k )
2218 
2219  IF ( verbose_level.GE.0 ) WRITE(*,*) k , probes_coords( 1:2 , k )
2220 
2221  END DO
2222 
2223  GOTO 310
2224 300 tend1 = .true.
2225 310 CONTINUE
2226 
2227 
2228 
2229 
2230  ! ----- end search for check points -----------------------------------------
2231 
2232 
2233  CLOSE( input_unit )
2234 
2235  bak_name = trim(run_name)//'.bak'
2236 
2237  OPEN(backup_unit,file=bak_name,status='unknown')
2238 
2239  WRITE(backup_unit, run_parameters )
2240 
2241  IF ( restart ) THEN
2242 
2243  WRITE(backup_unit,newrun_parameters)
2244  WRITE(backup_unit,restart_parameters)
2245 
2246  ELSE
2247 
2248  WRITE(backup_unit,newrun_parameters)
2249 
2250  IF ( riemann_flag) THEN
2251 
2252  WRITE(backup_unit,left_state)
2253  WRITE(backup_unit,right_state)
2254 
2255  ELSE
2256 
2257  WRITE(backup_unit,initial_conditions)
2258 
2259  END IF
2260 
2261  END IF
2262 
2263  WRITE(backup_unit, numeric_parameters )
2264 
2265  IF ( comp_cells_x .GT. 1 ) THEN
2266 
2267  WRITE(backup_unit,west_boundary_conditions)
2268  WRITE(backup_unit,east_boundary_conditions)
2269 
2270  END IF
2271 
2272  IF ( comp_cells_y .GT. 1 ) THEN
2273 
2274  WRITE(backup_unit,north_boundary_conditions)
2275  WRITE(backup_unit,south_boundary_conditions)
2276 
2277  END IF
2278 
2279  WRITE(backup_unit, expl_terms_parameters )
2280 
2281  WRITE(backup_unit,temperature_parameters)
2282 
2283  WRITE(backup_unit,solid_transport_parameters)
2284  WRITE(backup_unit,gas_transport_parameters)
2285  WRITE(backup_unit,liquid_transport_parameters)
2286 
2287  IF ( rheology_flag ) WRITE(backup_unit,rheology_parameters)
2288 
2289  IF ( output_runout_flag ) WRITE(backup_unit, runout_parameters)
2290 
2291  IF ( n_probes .GT. 0 ) THEN
2292 
2293  WRITE(backup_unit,*) '''PROBES_COORDS'''
2294  WRITE(backup_unit,*) n_probes
2295 
2296  DO k = 1,n_probes
2297 
2298  WRITE(backup_unit,109) probes_coords(1:2,k)
2299 
2300 109 FORMAT(2(1x,e14.7))
2301 
2302  END DO
2303 
2304  END IF
2305 
2306  CLOSE(backup_unit)
2307 
2308  END SUBROUTINE read_param
2309 
2310  !******************************************************************************
2312  !
2317  !
2321  !
2322  !******************************************************************************
2323 
2324  SUBROUTINE update_param
2325 
2326  IMPLICIT none
2327 
2328  INTEGER :: ios
2329 
2330  CHARACTER(LEN=40) :: run_name_org
2331  LOGICAL :: restart_org
2332  REAL*8 :: t_start_org
2333  REAL*8 :: t_end_org
2334  REAL*8 :: dt_output_org
2335  LOGICAL :: output_cons_flag_org
2336  LOGICAL :: output_phys_flag_org
2337  LOGICAL :: output_esri_flag_org
2338  INTEGER :: verbose_level_org
2339 
2340  run_name_org = run_name
2341  restart_org = restart
2342  t_start_org = t_start
2343  t_end_org = t_end
2344  dt_output_org = dt_output
2345  output_cons_flag_org = output_cons_flag
2346  output_phys_flag_org = output_phys_flag
2347  output_esri_flag_org = output_esri_flag
2348  verbose_level_org = verbose_level
2349 
2350 
2351  OPEN(input_unit,file=input_file,status='old')
2352 
2353  ! ------- READ run_parameters NAMELIST -----------------------------------
2354  READ(input_unit, run_parameters,iostat=ios )
2355 
2356  IF ( ios .NE. 0 ) THEN
2357 
2358  WRITE(*,*) 'IOSTAT=',ios
2359  WRITE(*,*) 'ERROR: problem with namelist RUN_PARAMETERS'
2360  WRITE(*,*) 'Please check the input file'
2361  stop
2362 
2363  ELSE
2364 
2365  rewind(input_unit)
2366 
2367  END IF
2368 
2369  CLOSE(input_unit)
2370 
2371  IF ( t_end_org .NE. t_end ) THEN
2372 
2373  WRITE(*,*) 'Modified input file: t_end =',t_end
2374 
2375  END IF
2376 
2377  IF ( dt_output_org .NE. dt_output ) THEN
2378 
2379  WRITE(*,*) 'Modified input file: dt_output =',dt_output
2380 
2381  END IF
2382 
2383  IF ( output_cons_flag_org .NEQV. output_cons_flag ) THEN
2384 
2385  WRITE(*,*) 'Modified input file: output_cons_flag =',output_cons_flag
2386 
2387  END IF
2388 
2389  IF ( output_phys_flag_org .NEQV. output_phys_flag ) THEN
2390 
2391  WRITE(*,*) 'Modified input file: output_phys_flag =',output_phys_flag
2392 
2393  END IF
2394 
2395  IF ( output_esri_flag_org .NEQV. output_esri_flag ) THEN
2396 
2397  WRITE(*,*) 'Modified input file: output_esri_flag =',output_esri_flag
2398 
2399  END IF
2400 
2401  IF ( verbose_level_org .NE. verbose_level ) THEN
2402 
2403  WRITE(*,*) 'Modified input file: verbose_level =',verbose_level
2404 
2405  END IF
2406 
2407  run_name_org = run_name
2408  restart_org = restart
2409  t_start_org = t_start
2410 
2411 
2412  END SUBROUTINE update_param
2413 
2414 
2415  !******************************************************************************
2417  !
2420  !
2424  !
2425  !******************************************************************************
2426 
2427  SUBROUTINE read_solution
2428 
2429  ! External procedures
2430  USE geometry_2d, ONLY : interp_2d_scalarb , regrid_scalar
2431  USE solver_2d, ONLY : allocate_solver_variables
2432 
2433  ! External variables
2434  USE geometry_2d, ONLY : comp_cells_x , x0 , comp_cells_y , y0 , dx , dy
2435  USE geometry_2d, ONLY : b_cent
2436  USE init_2d, ONLY : thickness_init
2437  USE parameters_2d, ONLY : n_vars
2438  USE solver_2d, ONLY : q
2439 
2440  IMPLICIT none
2441 
2442  CHARACTER(LEN=15) :: chara
2443 
2444  INTEGER :: j,k
2445 
2446  INTEGER :: dot_idx
2447 
2448  LOGICAL :: lexist
2449 
2450  CHARACTER(LEN=30) :: string
2451 
2452  CHARACTER(LEN=3) :: check_file
2453 
2454  INTEGER :: ncols , nrows , nodata_value
2455 
2456  REAL*8 :: xllcorner , yllcorner , cellsize
2457 
2458  REAL*8 :: xj , yk
2459 
2460  REAL*8, ALLOCATABLE :: thickness_input(:,:)
2461 
2462  REAL*8, ALLOCATABLE :: x1(:) , y1(:)
2463 
2464  REAL*8 :: x2 , y2
2465 
2466  REAL*8 :: xl , xr , yl , yr
2467 
2468  REAL*8 :: rho_c , rho_m , mass_fract(n_solid)
2469 
2470  REAL*8 :: sp_heat_c
2471 
2472  INTEGER :: solid_idx
2473 
2474  INTEGER :: i_vars
2475 
2476  INQUIRE (file=restart_file,exist=lexist)
2477 
2478  WRITE(*,*)
2479  ! WRITE(*,*) 'READ INIT',restart_file,lexist,restart_unit
2480 
2481  IF ( lexist .EQV. .false.) THEN
2482 
2483  WRITE(*,*) 'Restart: ',trim(restart_file),' not found'
2484  stop
2485 
2486  ELSE
2487 
2488  OPEN(restart_unit,file=restart_file,status='old')
2489 
2490  WRITE(*,*) 'Restart: ',trim(restart_file), ' found'
2491 
2492  END IF
2493 
2494  dot_idx = scan(restart_file, ".", .true.)
2495 
2496  check_file = restart_file(dot_idx+1:dot_idx+3)
2497 
2498  IF ( check_file .EQ. 'asc' ) THEN
2499 
2500  IF ( liquid_flag .AND. gas_flag ) THEN
2501 
2502  WRITE(*,*) 'ERROR: problem with namelist NEWRUN_PARAMETERS'
2503  WRITE(*,*) 'When restarting from .asc file only'
2504  WRITE(*,*) 'one of these parameters must be set to .TRUE.'
2505  WRITE(*,*) 'LIQUID_FLAG',liquid_flag
2506  WRITE(*,*) 'GAS_FLAG',liquid_flag
2507  WRITE(*,*) 'Please check the input file'
2508  CLOSE(restart_unit)
2509  stop
2510 
2511  ELSEIF ( ( .NOT.liquid_flag ) .AND. ( .NOT. gas_flag ) ) THEN
2512 
2513  WRITE(*,*) 'ERROR: problem with namelist NEWRUN_PARAMETERS'
2514  WRITE(*,*) 'When restarting from .asc file only one'
2515  WRITE(*,*) 'of these parameters must be set to .TRUE.'
2516  WRITE(*,*) 'LIQUID_FLAG',liquid_flag
2517  WRITE(*,*) 'GAS_FLAG',liquid_flag
2518  WRITE(*,*) 'Please check the input file'
2519  CLOSE(restart_unit)
2520  stop
2521 
2522  ELSEIF ( gas_flag ) THEN
2523 
2524  WRITE(*,*) 'Carrier phase: gas'
2525 
2526  ELSEIF ( liquid_flag ) THEN
2527 
2528  WRITE(*,*) 'Carrier phase: liquid'
2529 
2530  END IF
2531 
2532  READ(restart_unit,*) chara, ncols
2533  READ(restart_unit,*) chara, nrows
2534  READ(restart_unit,*) chara, xllcorner
2535  READ(restart_unit,*) chara, yllcorner
2536  READ(restart_unit,*) chara, cellsize
2537  READ(restart_unit,*) chara, nodata_value
2538 
2539  ALLOCATE( thickness_init(comp_cells_x,comp_cells_y) )
2540  ALLOCATE( thickness_input(ncols,nrows) )
2541 
2542  IF ( ( xllcorner - x0 ) .GT. 1.d-5*cellsize ) THEN
2543 
2544  WRITE(*,*)
2545  WRITE(*,*) 'WARNING: initial solution and domain extent'
2546  WRITE(*,*) 'xllcorner greater than x0', xllcorner , x0
2547 
2548  END IF
2549 
2550  IF ( ( yllcorner - y0 ) .GT. 1.d-5*cellsize ) THEN
2551 
2552  WRITE(*,*)
2553  WRITE(*,*) 'WARNING: initial solution and domain extent'
2554  WRITE(*,*) 'yllcorner greater then y0', yllcorner , y0
2555 
2556  END IF
2557 
2558  IF ( x0+cell_size*(comp_cells_x+1) - ( xllcorner+cellsize*(ncols+1) ) &
2559  .GT. 1.d-5*cellsize ) THEN
2560 
2561  WRITE(*,*)
2562  WRITE(*,*) 'WARNING: initial solution and domain extent'
2563  WRITE(*,*) 'xrrcorner greater than ', xllcorner , x0
2564 
2565  END IF
2566 
2567  IF ( x0+cell_size*(comp_cells_x+1) - ( xllcorner+cellsize*(ncols+1) ) &
2568  .GT. 1.d-5*cellsize ) THEN
2569 
2570  WRITE(*,*)
2571  WRITE(*,*) 'WARNING: initial solution and domain extent'
2572  WRITE(*,*) 'yllcorner greater then y0', yllcorner , y0
2573 
2574  END IF
2575 
2576 
2577  IF ( cellsize .NE. cell_size ) THEN
2578 
2579  WRITE(*,*)
2580  WRITE(*,*) 'WARNING: changing resolution of restart'
2581  WRITE(*,*) 'cellsize not equal to cell_size', cellsize , cell_size
2582  WRITE(*,*)
2583 
2584  END IF
2585 
2586  DO k=1,nrows
2587 
2588  WRITE(*,fmt="(A1,A,t21,F6.2,A)",advance="NO") achar(13), &
2589  & " Percent Complete: ",( REAL(k) / REAL(nrows))*100.0, "%"
2590 
2591  READ(restart_unit,*) thickness_input(1:ncols,nrows-k+1)
2592 
2593  ENDDO
2594 
2595  WRITE(*,*)
2596  WRITE(*,*) 'Total volume from restart =',cellsize**2*sum(thickness_input)
2597 
2598  WHERE ( thickness_init .EQ. nodata_value )
2599 
2600  thickness_init = 0.d0
2601 
2602  END WHERE
2603 
2604 
2605 
2606  ALLOCATE( x1(ncols) , y1(nrows) )
2607 
2608  DO j=1,ncols
2609 
2610  x1(j) = xllcorner + (j-1)*cellsize
2611 
2612  END DO
2613 
2614  DO k=1,nrows
2615 
2616  y1(k) = yllcorner + (k-1)*cellsize
2617 
2618  END DO
2619 
2620  DO j=1,comp_cells_x
2621 
2622  x2 = x0 + (j-1)*cell_size
2623 
2624  DO k=1,comp_cells_y
2625 
2626  y2 = y0 + (k-1)*cell_size
2627 
2628  CALL interp_2d_scalarb( x1 , y1 , thickness_input , x2 , y2 , &
2629  thickness_init(j,k) )
2630 
2631  END DO
2632 
2633  END DO
2634 
2635  !----- NEW INITIALIZATION OF THICKNESS FROM RESTART
2636  DEALLOCATE( x1 , y1 )
2637 
2638  ALLOCATE( x1(ncols+1) , y1(nrows+1) )
2639 
2640  DO j=1,ncols+1
2641 
2642  x1(j) = xllcorner + (j-1)*cellsize
2643 
2644  END DO
2645 
2646  DO k=1,nrows+1
2647 
2648  y1(k) = yllcorner + (k-1)*cellsize
2649 
2650  END DO
2651 
2652  DO j=1,comp_cells_x
2653 
2654  xl = x0 + (j-1)*cell_size
2655  xr = x0 + (j)*cell_size
2656 
2657  DO k=1,comp_cells_y
2658 
2659  yl = y0 + (k-1)*cell_size
2660  yr = y0 + (k)*cell_size
2661 
2662  CALL regrid_scalar( x1 , y1 , thickness_input , xl , xr , yl , &
2663  yr , thickness_init(j,k) )
2664 
2665  END DO
2666 
2667  END DO
2668 
2669 
2670  !----- END NEW INITIALIZATION OF THICKNESS FROM RESTART
2671 
2672  IF ( gas_flag ) THEN
2673 
2674  rho_c = pres / ( sp_gas_const_a * t_init )
2675  sp_heat_c = sp_heat_a
2676 
2677  ELSE
2678 
2679  rho_c = rho_l
2680  sp_heat_c = sp_heat_l
2681 
2682  END IF
2683 
2684  rho_m = sum( rho_s(1:n_solid)*alphas_init(1:n_solid) ) + ( 1.d0 - &
2685  sum( alphas_init(1:n_solid) ) ) * rho_c
2686 
2687  mass_fract = rho_s * alphas_init / rho_m
2688 
2689  q(1,:,:) = thickness_init(:,:) * rho_m
2690 
2691  WRITE(*,*) 'Total volume on computational grid =',cell_size**2 * &
2692  sum( thickness_init(:,:) )
2693  WRITE(*,*) 'Total mass on computational grid =',cell_size**2 * &
2694  sum( q(1,:,:) )
2695 
2696  ! rhom*h*u
2697  q(2,:,:) = 0.d0
2698  ! rhom*h*v
2699  q(3,:,:) = 0.d0
2700 
2701  ! energy (total or internal)
2702  q(4,:,:) = 0.d0
2703 
2704  WHERE ( thickness_init .GT. 0.d0 )
2705 
2706  q(4,:,:) = q(1,:,:) * t_init * ( sum( mass_fract(1:n_solid) * &
2707  sp_heat_s(1:n_solid) ) + &
2708  ( 1.d0 - sum( mass_fract ) ) * sp_heat_l )
2709 
2710  END WHERE
2711 
2712  DO solid_idx=5,4+n_solid
2713 
2714  ! rhos*h*alphas
2715  q(solid_idx,:,:) = 0.d0
2716 
2717  WHERE ( thickness_init .GT. 0.d0 )
2718 
2719  q(solid_idx,:,:) = thickness_init(:,:) * alphas_init(solid_idx-4) *&
2720  rho_s(solid_idx-4)
2721 
2722  END WHERE
2723 
2724  END DO
2725 
2726  WRITE(*,*) 'MAXVAL(q(5,:,:))',maxval(q(5:4+n_solid,:,:))
2727 
2728  WRITE(*,*) 'Total sediment volume =',cell_size**2*sum( thickness_init* &
2729  sum(alphas_init) )
2730 
2731  output_idx = 0
2732 
2733  IF ( verbose_level .GE. 1 ) THEN
2734 
2735  WRITE(*,*) 'Min q(1,:,:) =',minval(q(1,:,:))
2736  WRITE(*,*) 'Max q(1,:,:) =',maxval(q(1,:,:))
2737  WRITE(*,*) 'SUM(q(1,:,:)) =',sum(q(1,:,:))
2738 
2739  DO k=1,nrows
2740 
2741  WRITE(*,*) k,b_cent(:,k)
2742  READ(*,*)
2743 
2744  END DO
2745 
2746  WRITE(*,*) 'SUM(B_cent(:,:)) =',sum(b_cent(:,:))
2747  READ(*,*)
2748 
2749  END IF
2750 
2751  DEALLOCATE( thickness_input )
2752 
2753  WRITE(*,*) 'n_vars',n_vars
2754 
2755  ELSEIF ( check_file .EQ. 'q_2' ) THEN
2756 
2757  DO k=1,comp_cells_y
2758 
2759  DO j=1,comp_cells_x
2760 
2761  READ(restart_unit,'(2e20.12,*(e20.12))') xj , yk , &
2762  (q(i_vars,j,k),i_vars=1,n_vars)
2763 
2764  IF ( q(1,j,k) .LE. 0.d0 ) q(1:n_vars,j,k) = 0.d0
2765 
2766  ENDDO
2767 
2768  READ(restart_unit,*)
2769 
2770  END DO
2771 
2772  WRITE(*,*) 'Total mass =',dx*dy* sum( q(1,:,:) )
2773 
2774  DO solid_idx=5,4+n_solid
2775 
2776  WRITE(*,*) 'Total sediment mass =',dx*dy* sum( q(solid_idx,:,:) )
2777 
2778  END DO
2779 
2780  j = scan(restart_file, '.' , .true. )
2781 
2782  string = trim(restart_file(j-4:j-1))
2783  READ( string,* ) output_idx
2784 
2785  WRITE(*,*)
2786  WRITE(*,*) 'Starting from output index ',output_idx
2787 
2788  ! Set this flag to 0 to not overwrite the initial condition
2789 
2790  ELSE
2791 
2792  WRITE(*,*) 'Restart file not in the right format (*.asc or *)'
2793  stop
2794 
2795  END IF
2796 
2797  CLOSE(restart_unit)
2798 
2799 
2800  END SUBROUTINE read_solution
2801 
2802  !******************************************************************************
2804  !
2809  !
2811  !
2815  !
2816  !******************************************************************************
2817 
2818  SUBROUTINE output_solution(time)
2819 
2820  ! external procedures
2821  USE constitutive_2d, ONLY : qc_to_qp, mixt_var
2822 
2823  USE geometry_2d, ONLY : comp_cells_x , b_cent , comp_cells_y , x_comp, &
2824  y_comp , deposit
2825 
2826  USE parameters_2d, ONLY : n_vars
2827  USE parameters_2d, ONLY : t_output , dt_output
2828  USE parameters_2d, ONLY : t_steady
2829 
2830  USE solver_2d, ONLY : q
2831 
2832  IMPLICIT none
2833 
2834  REAL*8, INTENT(IN) :: time
2835 
2836  CHARACTER(LEN=4) :: idx_string
2837 
2838  REAL*8 :: qp(n_vars+2)
2839 
2840  REAL*8 :: B_out
2841 
2842  REAL*8 :: r_u , r_v , r_h , r_alphas(n_solid) , r_T , r_Ri , r_rho_m
2843  REAL*8 :: r_rho_c
2844  REAL*8 :: r_red_grav
2845 
2846 
2847  INTEGER :: j,k
2848  INTEGER :: i
2849  INTEGER :: i_vars
2850 
2851  output_idx = output_idx + 1
2852 
2853  idx_string = lettera(output_idx-1)
2854 
2855  IF ( output_cons_flag ) THEN
2856 
2857  output_file_2d = trim(run_name)//'_'//idx_string//'.q_2d'
2858 
2859  WRITE(*,*) 'WRITING ',output_file_2d
2860 
2861  OPEN(output_unit_2d,file=output_file_2d,status='unknown',form='formatted')
2862 
2863  !WRITE(output_unit_2d,1002) x0,dx,comp_cells_x,y0,dy,comp_cells_y,t
2864 
2865  DO k = 1,comp_cells_y
2866 
2867  DO j = 1,comp_cells_x
2868 
2869  DO i = 1,n_vars
2870 
2871  ! Exponents with more than 2 digits cause problems reading
2872  ! into matlab... reset tiny values to zero:
2873  IF ( dabs(q(i,j,k)) .LT. 1d-99) q(i,j,k) = 0.d0
2874 
2875  ENDDO
2876 
2877  WRITE(output_unit_2d,'(2e20.12,*(e20.12))') x_comp(j), y_comp(k), &
2878  (q(i_vars,j,k),i_vars=1,n_vars)
2879 
2880  ENDDO
2881 
2882  WRITE(output_unit_2d,*) ' '
2883 
2884  END DO
2885 
2886  WRITE(output_unit_2d,*) ' '
2887  WRITE(output_unit_2d,*) ' '
2888 
2889  CLOSE(output_unit_2d)
2890 
2891  END IF
2892 
2893  IF ( output_phys_flag ) THEN
2894 
2895  output_file_2d = trim(run_name)//'_'//idx_string//'.p_2d'
2896 
2897  WRITE(*,*) 'WRITING ',output_file_2d
2898 
2899  OPEN(output_unit_2d,file=output_file_2d,status='unknown',form='formatted')
2900 
2901  DO k = 1,comp_cells_y
2902 
2903  DO j = 1,comp_cells_x
2904 
2905  CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars+2))
2906 
2907  CALL mixt_var(qp(1:n_vars+2),r_ri,r_rho_m,r_rho_c,r_red_grav)
2908 
2909  r_h = qp(1)
2910  r_u = qp(n_vars+1)
2911  r_v = qp(n_vars+2)
2912  r_t = qp(4)
2913  r_alphas(1:n_solid) = qp(5:4+n_solid)
2914 
2915  IF ( dabs( r_h ) .LT. 1d-99) r_h = 0.d0
2916  IF ( dabs( r_u ) .LT. 1d-99) r_u = 0.d0
2917  IF ( dabs( r_v ) .LT. 1d-99) r_v = 0.d0
2918  IF ( dabs(b_cent(j,k)) .LT. 1d-99) THEN
2919 
2920  b_out = 0.d0
2921 
2922  ELSE
2923 
2924  b_out = b_cent(j,k)
2925 
2926  END IF
2927 
2928  DO i=1,n_solid
2929 
2930  IF ( dabs( r_alphas(i) ) .LT. 1d-99) r_alphas(i) = 0.d0
2931  IF ( dabs( deposit(j,k,i) ) .LT. 1d-99) deposit(j,k,i) = 0.d0
2932 
2933  END DO
2934 
2935  IF ( dabs( r_t ) .LT. 1d-99) r_t = 0.d0
2936  IF ( dabs( r_rho_m ) .LT. 1d-99) r_rho_m = 0.d0
2937  IF ( dabs( r_red_grav ) .LT. 1d-99) r_red_grav = 0.d0
2938 
2939  WRITE(output_unit_2d,1010) x_comp(j), y_comp(k), r_h , r_u , r_v , &
2940  b_out , r_h + b_out , r_alphas , r_t , r_rho_m , r_red_grav , &
2941  deposit(j,k,:)
2942 
2943  END DO
2944 
2945  WRITE(output_unit_2d,*) ' '
2946 
2947  ENDDO
2948 
2949  WRITE(output_unit_2d,*) ' '
2950  WRITE(output_unit_2d,*) ' '
2951 
2952  CLOSE(output_unit_2d)
2953 
2954  END IF
2955 
2956 1010 FORMAT(100es13.5e2)
2957 
2958  t_output = time + dt_output
2959 
2960  IF ( output_esri_flag ) THEN
2961 
2962  CALL output_esri(output_idx)
2963 
2964  IF ( ( time .GE. t_end ) .OR. ( time .GE. t_steady ) ) THEN
2965 
2966  CALL output_max
2967 
2968  END IF
2969 
2970  END IF
2971 
2972  IF ( n_probes .GT. 0 ) CALL output_probes(output_idx)
2973 
2974  END SUBROUTINE output_solution
2975 
2976 
2977  !******************************************************************************
2979  !
2982  !
2986  !
2987  !******************************************************************************
2988 
2989  SUBROUTINE output_max
2990 
2991  USE geometry_2d, ONLY : grid_output
2992  USE solver_2d, ONLY : q1max
2993 
2994  INTEGER :: j
2995 
2996 
2997  !Save max thickness
2998  output_max_file = trim(run_name)//'_max.asc'
2999 
3000  WRITE(*,*) 'WRITING ',output_max_file
3001 
3002  OPEN(output_max_unit,file=output_max_file,status='unknown',form='formatted')
3003 
3004  grid_output = -9999
3005 
3006  WHERE ( q1max(:,:).GE. 1.d-5 )
3007 
3008  grid_output = q1max(:,:)
3009 
3010  END WHERE
3011 
3012  WRITE(output_max_unit,'(A,I5)') 'ncols ', comp_cells_x
3013  WRITE(output_max_unit,'(A,I5)') 'nrows ', comp_cells_y
3014  WRITE(output_max_unit,'(A,F15.3)') 'xllcorner ', x0
3015  WRITE(output_max_unit,'(A,F15.3)') 'yllcorner ', y0
3016  WRITE(output_max_unit,'(A,F15.3)') 'cellsize ', cell_size
3017  WRITE(output_max_unit,'(A,I5)') 'NODATA_value ', -9999
3018 
3019  DO j = comp_cells_y,1,-1
3020 
3021  WRITE(output_max_unit,*) grid_output(1:comp_cells_x,j)
3022 
3023  ENDDO
3024 
3025  CLOSE(output_max_unit)
3026 
3027  RETURN
3028 
3029  END SUBROUTINE output_max
3030 
3031  !******************************************************************************
3033  !
3036  !
3038  !
3042  !
3043  !******************************************************************************
3044 
3045  SUBROUTINE output_esri(output_idx)
3046 
3047  USE geometry_2d, ONLY : b_cent , grid_output , deposit
3048  ! USE geometry_2d, ONLY : comp_interfaces_x , comp_interfaces_y
3049  USE solver_2d, ONLY : qp
3050 
3051  IMPLICIT NONE
3052 
3053  INTEGER, INTENT(IN) :: output_idx
3054 
3055  CHARACTER(LEN=4) :: idx_string
3056  CHARACTER(LEN=4) :: isolid_string
3057  INTEGER :: j
3058  INTEGER :: i_solid
3059 
3060  IF ( output_idx .EQ. 1 ) THEN
3061 
3062  OPEN(dem_esri_unit,file='dem_esri.asc',status='unknown',form='formatted')
3063 
3064  WRITE(dem_esri_unit,'(A,I5)') 'ncols ', comp_cells_x
3065  WRITE(dem_esri_unit,'(A,I5)') 'nrows ', comp_cells_y
3066  WRITE(dem_esri_unit,'(A,F15.3)') 'xllcorner ', x0
3067  WRITE(dem_esri_unit,'(A,F15.3)') 'yllcorner ', y0
3068  WRITE(dem_esri_unit,'(A,F15.3)') 'cellsize ', cell_size
3069  WRITE(dem_esri_unit,'(A,I5)') 'NODATA_value ', -9999
3070 
3071  DO j = comp_cells_y,1,-1
3072 
3073  WRITE(dem_esri_unit,*) b_cent(1:comp_cells_x,j)
3074 
3075  ENDDO
3076 
3077  CLOSE(dem_esri_unit)
3078 
3079  END IF
3080 
3081  idx_string = lettera(output_idx-1)
3082 
3083  !Save thickness
3084  output_esri_file = trim(run_name)//'_'//idx_string//'.asc'
3085 
3086  WRITE(*,*) 'WRITING ',output_esri_file
3087 
3088  OPEN(output_esri_unit,file=output_esri_file,status='unknown',form='formatted')
3089 
3090  grid_output = -9999
3091 
3092  WHERE ( qp(1,:,:).GE. 1.d-5 )
3093 
3094  grid_output = qp(1,:,:)
3095 
3096  END WHERE
3097 
3098  WRITE(output_esri_unit,'(A,I5)') 'ncols ', comp_cells_x
3099  WRITE(output_esri_unit,'(A,I5)') 'nrows ', comp_cells_y
3100  WRITE(output_esri_unit,'(A,F15.3)') 'xllcorner ', x0
3101  WRITE(output_esri_unit,'(A,F15.3)') 'yllcorner ', y0
3102  WRITE(output_esri_unit,'(A,F15.3)') 'cellsize ', cell_size
3103  WRITE(output_esri_unit,'(A,I5)') 'NODATA_value ', -9999
3104 
3105  DO j = comp_cells_y,1,-1
3106 
3107  WRITE(output_esri_unit,*) grid_output(1:comp_cells_x,j)
3108 
3109  ENDDO
3110 
3111  CLOSE(output_esri_unit)
3112 
3113  !Save temperature
3114  output_esri_file = trim(run_name)//'_T_'//idx_string//'.asc'
3115 
3116  WRITE(*,*) 'WRITING ',output_esri_file
3117 
3118  OPEN(output_esri_unit,file=output_esri_file,status='unknown',form='formatted')
3119 
3120  grid_output = -9999
3121 
3122  WHERE ( qp(1,:,:) .GE. 1.d-5 )
3123 
3124  grid_output = qp(4,:,:)
3125 
3126  END WHERE
3127 
3128  WRITE(output_esri_unit,'(A,I5)') 'ncols ', comp_cells_x
3129  WRITE(output_esri_unit,'(A,I5)') 'nrows ', comp_cells_y
3130  WRITE(output_esri_unit,'(A,F15.3)') 'xllcorner ', x0
3131  WRITE(output_esri_unit,'(A,F15.3)') 'yllcorner ', y0
3132  WRITE(output_esri_unit,'(A,F15.3)') 'cellsize ', cell_size
3133  WRITE(output_esri_unit,'(A,I5)') 'NODATA_value ', -9999
3134 
3135  DO j = comp_cells_y,1,-1
3136 
3137  WRITE(output_esri_unit,*) grid_output(1:comp_cells_x,j)
3138 
3139  ENDDO
3140 
3141  CLOSE(output_esri_unit)
3142 
3143  !Save deposit
3144  DO i_solid=1,n_solid
3145 
3146  isolid_string = lettera(i_solid)
3147 
3148  output_esri_file = trim(run_name)//'_dep_'//isolid_string//'_'//idx_string//'.asc'
3149 
3150  WRITE(*,*) 'WRITING ',output_esri_file
3151 
3152  OPEN(output_esri_unit,file=output_esri_file,status='unknown',form='formatted')
3153 
3154  grid_output = -9999
3155 
3156  WHERE ( deposit(:,:,i_solid) .GT. 0.d0 )
3157 
3158  grid_output = deposit(:,:,i_solid)
3159 
3160  END WHERE
3161 
3162  WRITE(output_esri_unit,'(A,I5)') 'ncols ', comp_cells_x
3163  WRITE(output_esri_unit,'(A,I5)') 'nrows ', comp_cells_y
3164  WRITE(output_esri_unit,'(A,F15.3)') 'xllcorner ', x0
3165  WRITE(output_esri_unit,'(A,F15.3)') 'yllcorner ', y0
3166  WRITE(output_esri_unit,'(A,F15.3)') 'cellsize ', cell_size
3167  WRITE(output_esri_unit,'(A,I5)') 'NODATA_value ', -9999
3168 
3169  DO j = comp_cells_y,1,-1
3170 
3171  WRITE(output_esri_unit,'(2000ES12.3E3)') grid_output(1:comp_cells_x,j)
3172 
3173  ENDDO
3174 
3175  CLOSE(output_esri_unit)
3176 
3177  END DO
3178 
3179  RETURN
3180 
3181  END SUBROUTINE output_esri
3182 
3183  SUBROUTINE close_units
3184 
3185  IMPLICIT NONE
3186 
3187  IF ( output_runout_flag) CLOSE(runout_unit)
3188 
3189  END SUBROUTINE close_units
3190 
3191  !******************************************************************************
3193  !
3196  !
3198  !
3202  !
3203  !******************************************************************************
3204 
3205  CHARACTER*4 FUNCTION lettera(k)
3206  IMPLICIT NONE
3207  CHARACTER ones,tens,hund,thou
3208  !
3209  INTEGER :: k
3210  !
3211  INTEGER :: iten, ione, ihund, ithou
3212  !
3213  ithou=int(k/1000)
3214  ihund=int((k-(ithou*1000))/100)
3215  iten=int((k-(ithou*1000)-(ihund*100))/10)
3216  ione=k-ithou*1000-ihund*100-iten*10
3217  ones=char(ione+48)
3218  tens=char(iten+48)
3219  hund=char(ihund+48)
3220  thou=char(ithou+48)
3221  lettera=thou//hund//tens//ones
3222  !
3223  RETURN
3224  END FUNCTION lettera
3225 
3226  !******************************************************************************
3228  !
3236  !******************************************************************************
3237 
3238  SUBROUTINE output_probes(output_idx)
3239 
3240  USE geometry_2d, ONLY : x_comp , y_comp
3241  USE solver_2d, ONLY : q
3242 
3243  USE geometry_2d, ONLY : interp_2d_scalarb
3244 
3245  IMPLICIT NONE
3246 
3247  INTEGER, INTENT(IN) :: output_idx
3248 
3249  CHARACTER(LEN=4) :: idx_string
3250 
3251  REAL*8 :: f2
3252 
3253  INTEGER :: k
3254 
3255  idx_string = lettera(output_idx-1)
3256 
3257  !Save thickness
3258  probes_file = trim(run_name)//'_'//idx_string//'.prb'
3259 
3260  OPEN(probes_unit,file=probes_file,status='unknown',form='formatted')
3261 
3262  DO k=1,n_probes
3263 
3264  CALL interp_2d_scalarb( x_comp , y_comp , q(1,:,:) , &
3265  probes_coords(1,k) , probes_coords(2,k) , f2 )
3266 
3267  WRITE(probes_unit,'(3e20.12)') probes_coords(1,k) , probes_coords(2,k) ,f2
3268 
3269  END DO
3270 
3271  CLOSE(probes_unit)
3272 
3273  END SUBROUTINE output_probes
3274 
3275  !******************************************************************************
3277  !
3286  !******************************************************************************
3287 
3288  SUBROUTINE output_runout(time,stop_flag)
3289 
3290  USE geometry_2d, ONLY : x_comp , y_comp , b_cent , dx , dy
3291  USE parameters_2d, ONLY : t_runout
3292  USE solver_2d, ONLY : q , q0 , dt
3293 
3294 
3295  IMPLICIT NONE
3296 
3297  REAL*8, INTENT(IN) :: time
3298  LOGICAL, INTENT(INOUT) :: stop_flag
3299 
3300  REAL*8, ALLOCATABLE :: X(:,:), Y(:,:) , dist(:,:)
3301  INTEGER :: sX, sY
3302  INTEGER :: imax(2) , imin(2)
3303  REAL*8 :: max_mom
3304  REAL*8 :: area , area0 , dareaRel_dt
3305  REAL*8 :: dhRel_dt
3306 
3307  sx = size(x_comp)
3308  sy = size(y_comp)
3309 
3310  ALLOCATE( x(sx,sy) , y(sx,sy) , dist(sx,sy) )
3311 
3312  x(:,:) = spread( x_comp, 2, sy )
3313  y(:,:) = spread( y_comp, 1, sx )
3314 
3315  dist(:,:) = 0.d0
3316 
3317  IF ( time .EQ. t_start ) THEN
3318 
3319  ALLOCATE( h_old(sx,sy) )
3320 
3321  h_old(:,:) = q(1,:,:)
3322 
3323  IF ( ( x0_runout .EQ. -1 ) .AND. ( y0_runout .EQ. -1 ) ) THEN
3324 
3325  WHERE( q(1,:,:) > 1.d-5 ) dist = b_cent
3326  imin = maxloc( dist )
3327 
3328  x0_runout = x(imin(1),imin(2))
3329  y0_runout = y(imin(1),imin(2))
3330 
3331  WRITE(*,*) 'Runout calculated as linear distance from: (' , &
3332  x0_runout ,',',y0_runout,')'
3333 
3334  dist(:,:) = 0.d0
3335 
3336  WHERE( q(1,:,:) >1.d-5 ) dist = dsqrt( (x-x0_runout)**2 &
3337  + ( y - y0_runout )**2 )
3338 
3339  imax = maxloc( dist )
3340 
3341  init_runout = dist(imax(1),imax(2))
3342 
3343  END IF
3344 
3345  max_mom = 0.d0
3346 
3347  ELSE
3348 
3349  WHERE( h_old(:,:) > 1.d-5 ) dist = dsqrt( q(2,:,:)**2 + q(3,:,:)**2 )
3350 
3351  max_mom = maxval( dist )
3352 
3353  h_old(:,:) = q(1,:,:)
3354 
3355  END IF
3356 
3357  dist(:,:) = 0.d0
3358 
3359  WHERE( q(1,:,:) > 1.d-5 ) dist = dsqrt( ( x - x0_runout )**2 &
3360  + ( y - y0_runout )**2 )
3361 
3362  imax = maxloc( dist )
3363 
3364  OPEN(dakota_unit,file='dakota.txt',status='replace',form='formatted')
3365 
3366  WRITE(dakota_unit,*) 'final runout =', dist(imax(1),imax(2)) - init_runout
3367 
3368  CLOSE(dakota_unit)
3369 
3370  area0 = dx*dy*count(q0(1,:,:).GT.1.d-7)
3371  area = dx*dy*count(q(1,:,:).GT.1.d-7)
3372 
3373  WRITE(runout_unit,'(A,F12.3,A,F12.3,A,F14.3)') 'Time (s) = ',time , &
3374  ' Runout (m) = ',dist(imax(1),imax(2)) - init_runout,' Area (m^2) = ', &
3375  area
3376 
3377  CALL flush(runout_unit)
3378 
3379  dist(:,:) = 0.d0
3380 
3381  WHERE( q(1,:,:) > 1.d-5 ) dist = abs(q(1,:,:)-q0(1,:,:)) / &
3382  max(q(1,:,:),q0(1,:,:))
3383 
3384  IF ( time .GT. t_start ) THEN
3385 
3386  darearel_dt = abs( area - area0 ) / ( area * dt )
3387 
3388  dhrel_dt = sum( dist / dt ) / count(q(1,:,:).GT.1.d-5)
3389 
3390  IF ( ( max(darearel_dt,dhrel_dt) .LT. eps_stop ) .AND. &
3391  (.NOT.stop_flag) ) THEN
3392 
3393  WRITE(*,*) 'Steady solution reached'
3394  WRITE(*,*) 'dareaRel_dt',darearel_dt
3395  WRITE(*,*) 'dhRel_dt',dhrel_dt
3396  stop_flag = .true.
3397 
3398  END IF
3399 
3400  END IF
3401 
3402  DEALLOCATE( x , y , dist )
3403 
3404  t_runout = time + dt_runout
3405 
3406  END SUBROUTINE output_runout
3407 
3408 END MODULE inpout_2d
3409 
constitutive_2d::kin_visc_a
real *8 kin_visc_a
Kinematic viscosity of air (units: m2 s-1)
Definition: constitutive_2d.f90:112
constitutive_2d::emissivity
real *8 emissivity
emissivity (eps in Costa & Macedonio, 2005)
Definition: constitutive_2d.f90:66
parameters_2d::max_dt
real *8 max_dt
Largest time step allowed.
Definition: parameters_2d.f90:16
inpout_2d::xllcorner
real *8 xllcorner
Definition: inpout_2d.f90:161
inpout_2d::probes_unit
integer, parameter probes_unit
Probes data unit.
Definition: inpout_2d.f90:100
inpout_2d::probes_coords
real *8, dimension(:,:), allocatable probes_coords
Definition: inpout_2d.f90:167
geometry_2d::dy
real *8 dy
Control volumes size.
Definition: geometry_2d.f90:93
parameters_2d::t_collapse
real *8 t_collapse
Definition: parameters_2d.f90:91
parameters_2d::n_rk
integer n_rk
Runge-Kutta order.
Definition: parameters_2d.f90:131
inpout_2d::run_name
character(len=40) run_name
Name of the run.
Definition: inpout_2d.f90:85
constitutive_2d::sp_heat_l
real *8 sp_heat_l
Sepcific heat of liquid (units: J K-1 kg-1)
Definition: constitutive_2d.f90:151
parameters_2d::t_steady
real *8 t_steady
end time when reached steady solution
Definition: parameters_2d.f90:159
inpout_2d::eps_stop
real *8 eps_stop
Definition: inpout_2d.f90:173
init_2d::alphas_w
real *8, dimension(:), allocatable alphas_w
Left sediment concentration.
Definition: init_2d.f90:26
geometry_2d::x_comp
real *8, dimension(:), allocatable x_comp
Location of the centers (x) of the control volume of the domain.
Definition: geometry_2d.f90:14
geometry_2d::y0
real *8 y0
Bottom of the physical domain.
Definition: geometry_2d.f90:94
constitutive_2d::nu_ref
real *8 nu_ref
reference kinematic viscosity [m2/s]
Definition: constitutive_2d.f90:45
constitutive_2d::t_ground
real *8 t_ground
temperature of lava-ground interface
Definition: constitutive_2d.f90:72
solver_2d::q0
real *8, dimension(:,:,:), allocatable q0
Conservative variables.
Definition: solver_2d.f90:44
inpout_2d::x0_runout
real *8 x0_runout
Definition: inpout_2d.f90:173
inpout_2d::output_probes
subroutine output_probes(output_idx)
Write solution at selected points on file.
Definition: inpout_2d.f90:3239
geometry_2d::b_cent
real *8, dimension(:,:), allocatable b_cent
Topography at the centers of the control volumes.
Definition: geometry_2d.f90:41
geometry_2d::comp_cells_x
integer comp_cells_x
Number of control volumes x in the comp. domain.
Definition: geometry_2d.f90:98
parameters_2d::vel_source
real *8 vel_source
Definition: parameters_2d.f90:78
constitutive_2d::beta2
real *8 beta2
2nd param for yield strenght empirical relationship (O&#39;Brian et al, 1993)
Definition: constitutive_2d.f90:83
parameters_2d::alphas_init
real *8, dimension(:), allocatable alphas_init
Initial sediment concentration in the pile of material.
Definition: parameters_2d.f90:120
init_2d::thickness_init
real *8, dimension(:,:), allocatable thickness_init
Definition: init_2d.f90:17
geometry_2d::y_comp
real *8, dimension(:), allocatable y_comp
Location of the centers (y) of the control volume of the domain.
Definition: geometry_2d.f90:20
parameters_2d::rheology_flag
logical rheology_flag
Flag to choose if we add the rheology.
Definition: parameters_2d.f90:59
inpout_2d::output_esri_unit
integer, parameter output_esri_unit
Esri Output unit.
Definition: inpout_2d.f90:102
geometry_2d::n_topography_profile_y
integer n_topography_profile_y
Definition: geometry_2d.f90:88
constitutive_2d::mixt_var
subroutine mixt_var(qpj, r_Ri, r_rho_m, r_rho_c, r_red_grav)
Physical variables.
Definition: constitutive_2d.f90:569
inpout_2d::alphas_bcs
type(bc), dimension(:), allocatable alphas_bcs
Definition: inpout_2d.f90:151
constitutive_2d::entrainment_flag
logical entrainment_flag
flag to activate air entrainment
Definition: constitutive_2d.f90:17
inpout_2d::read_solution
subroutine read_solution
Read the solution from the restart unit.
Definition: inpout_2d.f90:2428
parameters_2d::t_init
real *8 t_init
Initial temperature of the pile of material.
Definition: parameters_2d.f90:117
parameters_2d::liquid_flag
logical liquid_flag
Definition: parameters_2d.f90:67
inpout_2d::alpha1_ref
real *8 alpha1_ref
Definition: inpout_2d.f90:177
inpout_2d::y0_runout
real *8 y0_runout
Definition: inpout_2d.f90:173
init_2d::u_e
real *8 u_e
Right velocity x.
Definition: init_2d.f90:30
constitutive_2d::c_p
real *8 c_p
specific heat [J kg-1 K-1]
Definition: constitutive_2d.f90:54
init_2d::v_w
real *8 v_w
Left velocity y.
Definition: init_2d.f90:25
inpout_2d::nodata_value
integer nodata_value
Definition: inpout_2d.f90:159
parameters_2d
Parameters.
Definition: parameters_2d.f90:7
geometry_2d::dx
real *8 dx
Control volumes size.
Definition: geometry_2d.f90:90
inpout_2d::output_esri_flag
logical output_esri_flag
Flag to save the output in esri ascii format *.asc.
Definition: inpout_2d.f90:121
parameters_2d::y_release
real *8 y_release
Initial y-coordinate of the pile.
Definition: parameters_2d.f90:104
parameters_2d::t_end
real *8 t_end
end time for the run
Definition: parameters_2d.f90:155
geometry_2d::comp_cells_y
integer comp_cells_y
Number of control volumes y in the comp. domain.
Definition: geometry_2d.f90:100
solver_2d
Numerical solver.
Definition: solver_2d.f90:12
inpout_2d::ncols
integer ncols
Definition: inpout_2d.f90:159
inpout_2d::restart_file
character(len=40) restart_file
Name of the restart file.
Definition: inpout_2d.f90:89
inpout_2d::hv_bcw
type(bc) hv_bcw
Definition: inpout_2d.f90:142
inpout_2d::h_bcs
type(bc) h_bcs
Definition: inpout_2d.f90:150
constitutive_2d::xi
real *8 xi
Definition: constitutive_2d.f90:24
constitutive_2d::kin_visc_l
real *8 kin_visc_l
Kinematic viscosity of liquid (units: m2 s-1)
Definition: constitutive_2d.f90:115
parameters_2d::rheology_model
integer rheology_model
choice of the rheology model
Definition: parameters_2d.f90:65
inpout_2d::probes_file
character(len=40) probes_file
Name of the probes file.
Definition: inpout_2d.f90:90
parameters_2d::bcw
type(bc), dimension(:), allocatable bcw
bcW&flag defines the west boundary condition:
Definition: parameters_2d.f90:178
inpout_2d
Input/Output module.
Definition: inpout_2d.f90:13
inpout_2d::hv_bce
type(bc) hv_bce
Definition: inpout_2d.f90:146
init_2d::riemann_interface
real *8 riemann_interface
Riemann problem interface relative position. It is a value between 0 and 1.
Definition: init_2d.f90:21
constitutive_2d::sp_gas_const_a
real *8 sp_gas_const_a
Specific gas constant of air (units: J kg-1 K-1)
Definition: constitutive_2d.f90:109
init_2d
Initial solution.
Definition: init_2d.f90:8
inpout_2d::update_param
subroutine update_param
Read the input file.
Definition: inpout_2d.f90:2325
inpout_2d::hv_bcn
type(bc) hv_bcn
Definition: inpout_2d.f90:154
init_2d::u_w
real *8 u_w
Left velocity x.
Definition: init_2d.f90:24
constitutive_2d::thermal_conductivity
real *8 thermal_conductivity
thermal conductivity [W m-1 K-1] (k in Costa & Macedonio, 2005)
Definition: constitutive_2d.f90:75
parameters_2d::velocity_ang_release
real *8 velocity_ang_release
Initial velocity direction (angle in degree): .
Definition: parameters_2d.f90:114
constitutive_2d::t_s_substrate
real *8 t_s_substrate
temperature of solid substrate (units: K)
Definition: constitutive_2d.f90:142
geometry_2d::regrid_scalar
subroutine regrid_scalar(xin, yin, fin, xl, xr, yl, yr, fout)
Scalar regrid (2D)
Definition: geometry_2d.f90:882
parameters_2d::velocity_mod_release
real *8 velocity_mod_release
Initial velocity module of the pile.
Definition: parameters_2d.f90:107
constitutive_2d::sp_heat_a
real *8 sp_heat_a
Specific heat of air (units: J K-1 kg-1)
Definition: constitutive_2d.f90:106
geometry_2d::n_topography_profile_x
integer n_topography_profile_x
Definition: geometry_2d.f90:88
constitutive_2d::beta1
real *8 beta1
2nd param for fluid viscosity empirical relationship (O&#39;Brian et al, 1993)
Definition: constitutive_2d.f90:89
parameters_2d::released_volume
real *8 released_volume
Initial volume of the flow.
Definition: parameters_2d.f90:98
inpout_2d::output_esri_file
character(len=40) output_esri_file
Name of the esri output files.
Definition: inpout_2d.f90:92
solver_2d::dt
real *8 dt
Time step.
Definition: solver_2d.f90:121
constitutive_2d::kin_visc_c
real *8 kin_visc_c
Kinematic viscosity of carrier phase (units: m2 s-1)
Definition: constitutive_2d.f90:118
parameters_2d::time_param
real *8, dimension(4) time_param
Definition: parameters_2d.f90:83
parameters_2d::n_vars
integer n_vars
Number of conservative variables.
Definition: parameters_2d.f90:124
parameters_2d::topo_change_flag
logical topo_change_flag
Definition: parameters_2d.f90:71
inpout_2d::output_max
subroutine output_max
Write the maximum thickness in ESRI format.
Definition: inpout_2d.f90:2990
inpout_2d::dakota_unit
integer, parameter dakota_unit
Definition: inpout_2d.f90:106
parameters_2d::alphas_collapse
real *8, dimension(:), allocatable alphas_collapse
Definition: parameters_2d.f90:93
parameters_2d::cfl
real *8 cfl
Courant-Friedrichs-Lewy parameter.
Definition: parameters_2d.f90:18
parameters_2d::alphas_source
real *8, dimension(:), allocatable alphas_source
Definition: parameters_2d.f90:81
constitutive_2d::enne
real *8 enne
thermal boundary layer fraction of total thickness
Definition: constitutive_2d.f90:69
parameters_2d::dt0
real *8 dt0
Initial time step.
Definition: parameters_2d.f90:14
constitutive_2d
Constitutive equations.
Definition: constitutive_2d.f90:4
init_2d::alphas_e
real *8, dimension(:), allocatable alphas_e
Right sediment concentration.
Definition: init_2d.f90:32
inpout_2d::h_bcw
type(bc) h_bcw
Definition: inpout_2d.f90:142
inpout_2d::sp_heat0_s
real *8, dimension(1000) sp_heat0_s
Definition: inpout_2d.f90:175
constitutive_2d::grav
real *8 grav
gravitational acceleration
Definition: constitutive_2d.f90:20
inpout_2d::init_runout
real *8 init_runout
Definition: inpout_2d.f90:173
inpout_2d::t_bcw
type(bc) t_bcw
Definition: inpout_2d.f90:142
constitutive_2d::mu
real *8 mu
drag coefficients (Voellmy-Salm model)
Definition: constitutive_2d.f90:23
constitutive_2d::visc_par
real *8 visc_par
viscosity parameter [K-1] (b in Table 1 Costa & Macedonio, 2005)
Definition: constitutive_2d.f90:48
constitutive_2d::alpha1_coeff
real *8 alpha1_coeff
ratio between reference value from input and computed values from eq.
Definition: constitutive_2d.f90:86
parameters_2d::alphal_source
real *8 alphal_source
Definition: parameters_2d.f90:82
parameters_2d::x_release
real *8 x_release
Initial x-coordiante of the pile.
Definition: parameters_2d.f90:101
inpout_2d::nrows
integer nrows
Definition: inpout_2d.f90:159
parameters_2d::r_collapse
real *8 r_collapse
Definition: parameters_2d.f90:90
constitutive_2d::diam_s
real *8, dimension(:), allocatable diam_s
Diameter of sediments ( units: m )
Definition: constitutive_2d.f90:127
constitutive_2d::n_td
real *8 n_td
Mannings roughness coefficient ( units: T L^(-1/3) )
Definition: constitutive_2d.f90:95
inpout_2d::restart
logical restart
Flag to start a run from a previous output: .
Definition: inpout_2d.f90:115
constitutive_2d::erosion_coeff
real *8, dimension(:), allocatable erosion_coeff
erosion model coefficient (units: m-1 )
Definition: constitutive_2d.f90:139
geometry_2d::x0
real *8 x0
Left of the physical domain.
Definition: geometry_2d.f90:91
constitutive_2d::settling_flag
logical settling_flag
Flag to determine if sedimentation is active.
Definition: constitutive_2d.f90:133
inpout_2d::output_runout
subroutine output_runout(time, stop_flag)
Write runout on file.
Definition: inpout_2d.f90:3289
inpout_2d::output_solution
subroutine output_solution(time)
Write the solution on the output unit.
Definition: inpout_2d.f90:2819
inpout_2d::h_bce
type(bc) h_bce
Definition: inpout_2d.f90:146
inpout_2d::backup_unit
integer, parameter backup_unit
Backup input data unit.
Definition: inpout_2d.f90:97
constitutive_2d::rho_l
real *8 rho_l
liquid density (units: kg m-3)
Definition: constitutive_2d.f90:148
inpout_2d::lettera
character *4 function lettera(k)
Numeric to String conversion.
Definition: inpout_2d.f90:3206
parameters_2d::solver_scheme
character(len=20) solver_scheme
Finite volume method: .
Definition: parameters_2d.f90:151
inpout_2d::output_phys_flag
logical output_phys_flag
Flag to save the physical variables on file *.p_2d.
Definition: inpout_2d.f90:127
inpout_2d::h_bcn
type(bc) h_bcn
Definition: inpout_2d.f90:154
geometry_2d::grid_output
real *8, dimension(:,:), allocatable grid_output
Solution in ascii grid format (ESRI)
Definition: geometry_2d.f90:50
inpout_2d::output_unit_2d
integer, parameter output_unit_2d
Output data 2D unit.
Definition: inpout_2d.f90:101
constitutive_2d::t_env
real *8 t_env
evironment temperature [K]
Definition: constitutive_2d.f90:33
constitutive_2d::alpha2
real *8 alpha2
1st param for yield strenght empirical relationship (O&#39;Brian et al, 1993)
Definition: constitutive_2d.f90:80
parameters_2d::x_source
real *8 x_source
Definition: parameters_2d.f90:75
inpout_2d::output_unit
integer, parameter output_unit
Output data unit.
Definition: inpout_2d.f90:98
inpout_2d::output_max_file
character(len=40) output_max_file
Name of the esri max. thick. file.
Definition: inpout_2d.f90:93
geometry_2d
Grid module.
Definition: geometry_2d.f90:7
inpout_2d::hu_bcs
type(bc) hu_bcs
Definition: inpout_2d.f90:150
geometry_2d::cell_size
real *8 cell_size
Definition: geometry_2d.f90:102
parameters_2d::gas_flag
logical gas_flag
Definition: parameters_2d.f90:69
inpout_2d::runout_unit
integer, parameter runout_unit
Definition: inpout_2d.f90:105
inpout_2d::output_cons_flag
logical output_cons_flag
Flag to save the conservative variables on file *.q_2d.
Definition: inpout_2d.f90:133
parameters_2d::theta
real *8 theta
Van Leer limiter parameter.
Definition: parameters_2d.f90:153
parameters_2d::riemann_flag
logical riemann_flag
Flag to choose the sort of problem to solve.
Definition: parameters_2d.f90:53
inpout_2d::output_file
character(len=40) output_file
Name of the output files.
Definition: inpout_2d.f90:88
inpout_2d::rho0_s
real *8, dimension(1000) rho0_s
Definition: inpout_2d.f90:175
init_2d::t_w
real *8 t_w
Left temperature.
Definition: init_2d.f90:27
constitutive_2d::atm_heat_transf_coeff
real *8 atm_heat_transf_coeff
atmospheric heat trasnfer coefficient [W m-2 K-1] (lambda in C&M, 2005)
Definition: constitutive_2d.f90:57
init_2d::v_e
real *8 v_e
Right velocity y.
Definition: init_2d.f90:31
parameters_2d::t_output
real *8 t_output
time of the next output
Definition: parameters_2d.f90:156
constitutive_2d::rho_a_amb
real *8 rho_a_amb
Ambient density of air ( units: kg m-3 )
Definition: constitutive_2d.f90:103
inpout_2d::read_param
subroutine read_param
Read the input file.
Definition: inpout_2d.f90:584
inpout_2d::t_bcn
type(bc) t_bcn
Definition: inpout_2d.f90:154
constitutive_2d::pres
real *8 pres
ambient pressure (units: Pa)
Definition: constitutive_2d.f90:145
geometry_2d::topography_profile
real *8, dimension(:,:,:), allocatable topography_profile
Definition: geometry_2d.f90:61
constitutive_2d::t_ref
real *8 t_ref
reference temperature [K]
Definition: constitutive_2d.f90:42
constitutive_2d::qc_to_qp
subroutine qc_to_qp(qc, qp)
Conservative to physical variables.
Definition: constitutive_2d.f90:678
inpout_2d::input_file
character(len=40) input_file
File with the run parameters.
Definition: inpout_2d.f90:87
inpout_2d::output_max_unit
integer, parameter output_max_unit
Esri max thick. output unit.
Definition: inpout_2d.f90:103
inpout_2d::hu_bce
type(bc) hu_bce
Definition: inpout_2d.f90:146
parameters_2d::y_source
real *8 y_source
Definition: parameters_2d.f90:76
geometry_2d::interp_2d_scalarb
subroutine interp_2d_scalarb(x1, y1, f1, x2, y2, f2)
Scalar interpolation (2D)
Definition: geometry_2d.f90:570
inpout_2d::output_esri
subroutine output_esri(output_idx)
Write the thickness in ESRI format.
Definition: inpout_2d.f90:3046
inpout_2d::alphas_bce
type(bc), dimension(:), allocatable alphas_bce
Definition: inpout_2d.f90:147
parameters_2d::h_collapse
real *8 h_collapse
Definition: parameters_2d.f90:92
inpout_2d::dt_runout
real *8 dt_runout
Definition: inpout_2d.f90:169
parameters_2d::energy_flag
logical energy_flag
Flag to choose the equation for temperature to solve.
Definition: parameters_2d.f90:47
inpout_2d::hu_bcw
type(bc) hu_bcw
Definition: inpout_2d.f90:142
constitutive_2d::exp_area_fract
real *8 exp_area_fract
fractional area of the exposed inner core (f in C&M, 2005)
Definition: constitutive_2d.f90:60
parameters_2d::limiter
integer, dimension(10) limiter
Limiter for the slope in the linear reconstruction: .
Definition: parameters_2d.f90:144
constitutive_2d::sp_heat_s
real *8, dimension(:), allocatable sp_heat_s
Specific heat of solids (units: J K-1 kg-1)
Definition: constitutive_2d.f90:130
inpout_2d::hv_bcs
type(bc) hv_bcs
Definition: inpout_2d.f90:150
inpout_2d::alphas_bcw
type(bc), dimension(:), allocatable alphas_bcw
Definition: inpout_2d.f90:143
inpout_2d::output_idx
integer output_idx
Counter for the output files.
Definition: inpout_2d.f90:109
inpout_2d::close_units
subroutine close_units
Definition: inpout_2d.f90:3184
inpout_2d::cellsize
real *8 cellsize
Definition: inpout_2d.f90:161
parameters_2d::bcn
type(bc), dimension(:), allocatable bcn
bcN&flag defines the north boundary condition:
Definition: parameters_2d.f90:208
inpout_2d::t_bce
type(bc) t_bce
Definition: inpout_2d.f90:146
parameters_2d::t_runout
real *8 t_runout
time of the next runout output
Definition: parameters_2d.f90:158
init_2d::t_e
real *8 t_e
Right temperature.
Definition: init_2d.f90:33
inpout_2d::hu_bcn
type(bc) hu_bcn
Definition: inpout_2d.f90:154
inpout_2d::erosion_coeff0
real *8, dimension(1000) erosion_coeff0
Definition: inpout_2d.f90:175
parameters_2d::n_eqns
integer n_eqns
Number of equations.
Definition: parameters_2d.f90:125
geometry_2d::deposit
real *8, dimension(:,:,:), allocatable deposit
deposit for the different classes
Definition: geometry_2d.f90:59
inpout_2d::output_file_2d
character(len=40) output_file_2d
Name of the output files.
Definition: inpout_2d.f90:91
inpout_2d::dem_esri_unit
integer, parameter dem_esri_unit
Computational grid Esri fmt unit.
Definition: inpout_2d.f90:104
constitutive_2d::t_ambient
real *8 t_ambient
Temperature of ambient air (units: K)
Definition: constitutive_2d.f90:121
inpout_2d::bak_name
character(len=40) bak_name
Backup file for the parameters.
Definition: inpout_2d.f90:86
inpout_2d::init_param
subroutine init_param
Initialization of the variables read from the input file.
Definition: inpout_2d.f90:254
inpout_2d::h_old
real *8, dimension(:,:), allocatable h_old
Definition: inpout_2d.f90:171
parameters_2d::h_source
real *8 h_source
Definition: parameters_2d.f90:80
parameters_2d::reconstr_coeff
real *8 reconstr_coeff
Slope coefficient in the linear reconstruction.
Definition: parameters_2d.f90:22
parameters_2d::bcs
type(bc), dimension(:), allocatable bcs
bcS&flag defines the south boundary condition:
Definition: parameters_2d.f90:198
constitutive_2d::emme
real *8 emme
velocity boundary layer fraction of total thickness
Definition: constitutive_2d.f90:51
inpout_2d::write_first_q
logical write_first_q
Definition: inpout_2d.f90:163
inpout_2d::n_probes
integer n_probes
Definition: inpout_2d.f90:165
inpout_2d::t_bcs
type(bc) t_bcs
Definition: inpout_2d.f90:150
inpout_2d::alphas_bcn
type(bc), dimension(:), allocatable alphas_bcn
Definition: inpout_2d.f90:155
constitutive_2d::tau
real *8 tau
drag coefficients (plastic model)
Definition: constitutive_2d.f90:30
parameters_2d::bce
type(bc), dimension(:), allocatable bce
bcE&flag defines the east boundary condition:
Definition: parameters_2d.f90:188
constitutive_2d::friction_factor
real *8 friction_factor
drag coefficients (B&W model)
Definition: constitutive_2d.f90:27
solver_2d::q1max
real *8, dimension(:,:), allocatable q1max
Maximum over time of thickness.
Definition: solver_2d.f90:86
parameters_2d::y_collapse
real *8 y_collapse
Definition: parameters_2d.f90:89
inpout_2d::yllcorner
real *8 yllcorner
Definition: inpout_2d.f90:161
parameters_2d::dt_output
real *8 dt_output
time interval for the output of the solution
Definition: parameters_2d.f90:157
parameters_2d::radial_source_flag
logical radial_source_flag
Definition: parameters_2d.f90:73
inpout_2d::input_unit
integer, parameter input_unit
Input data unit.
Definition: inpout_2d.f90:96
constitutive_2d::kappa
real *8 kappa
Empirical resistance parameter (dimensionless, input parameter)
Definition: constitutive_2d.f90:92
inpout_2d::restart_unit
integer, parameter restart_unit
Restart data unit.
Definition: inpout_2d.f90:99
inpout_2d::runout_file
character(len=40) runout_file
Name of the runout file.
Definition: inpout_2d.f90:94
solver_2d::allocate_solver_variables
subroutine allocate_solver_variables
Memory allocation.
Definition: solver_2d.f90:216
init_2d::hb_w
real *8 hb_w
Left height.
Definition: init_2d.f90:23
parameters_2d::x_collapse
real *8 x_collapse
Definition: parameters_2d.f90:88
parameters_2d::collapsing_volume_flag
logical collapsing_volume_flag
Definition: parameters_2d.f90:86
inpout_2d::diam0_s
real *8, dimension(1000) diam0_s
Definition: inpout_2d.f90:175
init_2d::hb_e
real *8 hb_e
Right height.
Definition: init_2d.f90:29
parameters_2d::n_solid
integer n_solid
Number of solid classes.
Definition: parameters_2d.f90:127
constitutive_2d::rho_s
real *8, dimension(:), allocatable rho_s
Density of sediments ( units: kg m-3 )
Definition: constitutive_2d.f90:124
parameters_2d::interfaces_relaxation
logical interfaces_relaxation
Flag to add the relaxation terms after the linear reconstruction: .
Definition: parameters_2d.f90:28
solver_2d::qp
real *8, dimension(:,:,:), allocatable qp
Physical variables ( )
Definition: solver_2d.f90:107
inpout_2d::output_runout_flag
logical output_runout_flag
Flag to save the max runout at ouput times.
Definition: inpout_2d.f90:139
parameters_2d::t_start
real *8 t_start
initial time for the run
Definition: parameters_2d.f90:154
parameters_2d::sed_vol_perc
real *8, dimension(:), allocatable sed_vol_perc
Definition: parameters_2d.f90:122
parameters_2d::r_source
real *8 r_source
Definition: parameters_2d.f90:77
parameters_2d::bc
Definition: parameters_2d.f90:163
parameters_2d::t_source
real *8 t_source
Definition: parameters_2d.f90:79
solver_2d::q
real *8, dimension(:,:,:), allocatable q
Conservative variables.
Definition: solver_2d.f90:42
Generated on Fri Dec 13 2019 21:38:13 for SW_VAR_DENS_MODEL by   doxygen 1.8.11