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