SW_VAR_DENS_MODEL  0.9
Dept-averagedgas-particlemodel
constitutive_2d.f90
Go to the documentation of this file.
1 !********************************************************************************
3 !********************************************************************************
4 MODULE constitutive_2d
5 
6  USE parameters_2d, ONLY : n_eqns , n_vars , n_solid
7  USE parameters_2d, ONLY : rheology_flag , rheology_model , energy_flag , &
8  liquid_flag , gas_flag
9  USE parameters_2d, ONLY : sed_vol_perc
10 
11  IMPLICIT none
12 
14  LOGICAL, ALLOCATABLE :: implicit_flag(:)
15 
17  LOGICAL :: entrainment_flag
18 
20  REAL*8 :: grav
21 
23  REAL*8 :: mu
24  REAL*8 :: xi
25 
27  REAL*8 :: friction_factor
28 
30  REAL*8 :: tau
31 
33  REAL*8 :: t_env
34 
36  REAL*8 :: rad_coeff
37 
39  REAL*8 :: frict_coeff
40 
42  REAL*8 :: t_ref
43 
45  REAL*8 :: nu_ref
46 
48  REAL*8 :: visc_par
49 
51  REAL*8 :: emme
52 
54  REAL*8 :: c_p
55 
57  REAL*8 :: atm_heat_transf_coeff
58 
60  REAL*8 :: exp_area_fract
61 
63  REAL*8, PARAMETER :: sbconst = 5.67d-8
64 
66  REAL*8 :: emissivity
67 
69  REAL*8 :: enne
70 
72  REAL*8 :: t_ground
73 
75  REAL*8 :: thermal_conductivity
76 
77  !--- START Lahars rheology model parameters
78 
80  REAL*8 :: alpha2 ! (units: kg m-1 s-2)
81 
83  REAL*8 :: beta2 ! (units: nondimensional)
84 
86  REAL*8 :: alpha1_coeff ! (units: nondimensional )
87 
89  REAL*8 :: beta1 ! (units: nondimensional, input parameter)
90 
92  REAL*8 :: kappa
93 
95  REAL*8 :: n_td
96 
97  !--- END Lahars rheology model parameters
98 
100  REAL*8 :: sp_heat_c ! ( initialized from input)
101 
103  REAL*8 :: rho_a_amb
104 
106  REAL*8 :: sp_heat_a
107 
109  REAL*8 :: sp_gas_const_a
110 
112  REAL*8 :: kin_visc_a
113 
115  REAL*8 :: kin_visc_l
116 
118  REAL*8 :: kin_visc_c
119 
121  REAL*8 :: t_ambient
122 
124  REAL*8, ALLOCATABLE :: rho_s(:)
125 
127  REAL*8, ALLOCATABLE :: diam_s(:)
128 
130  REAL*8, ALLOCATABLE :: sp_heat_s(:)
131 
133  LOGICAL :: settling_flag
134 
136  REAL*8 :: settling_vel
137 
139  REAL*8, ALLOCATABLE :: erosion_coeff(:)
140 
142  REAL*8 :: t_s_substrate
143 
145  REAL*8 :: pres
146 
148  REAL*8 :: rho_l
149 
151  REAL*8 :: sp_heat_l
152 
153 CONTAINS
154 
155  !******************************************************************************
157  !
161  !******************************************************************************
162 
163  SUBROUTINE init_problem_param
164 
165  USE parameters_2d, ONLY : n_nh
166 
167  ALLOCATE( implicit_flag(n_eqns) )
168 
169  implicit_flag(1:n_eqns) = .false.
170  implicit_flag(2) = .true.
171  implicit_flag(3) = .true.
172 
173  ! Temperature
174  implicit_flag(4) = .true.
175 
176  ! Solid volume fraction
177  implicit_flag(5:4+n_solid) = .true.
178 
179  n_nh = count( implicit_flag )
180 
181  RETURN
182 
183  END SUBROUTINE init_problem_param
184 
185  !******************************************************************************
187  !
198  !
201  !
203  !******************************************************************************
204 
205  SUBROUTINE r_phys_var(r_qj,r_h,r_u,r_v,r_alphas,r_rho_m,r_T,r_alphal)
206 
207  USE parameters_2d, ONLY : eps_sing
208  IMPLICIT none
209 
210  REAL*8, INTENT(IN) :: r_qj(n_vars)
211  REAL*8, INTENT(OUT) :: r_h
212  REAL*8, INTENT(OUT) :: r_u
213  REAL*8, INTENT(OUT) :: r_v
214  REAL*8, INTENT(OUT) :: r_alphas(n_solid)
215  REAL*8, INTENT(OUT) :: r_rho_m
216  REAL*8, INTENT(OUT) :: r_T
217  REAL*8, INTENT(OUT) :: r_alphal
218 
219  REAL*8 :: r_inv_rhom
220  REAL*8 :: r_xs(n_solid)
221  REAL*8 :: r_xs_tot
222 
223  REAL*8 :: r_Ri
224  REAL*8 :: r_xl
225  REAL*8 :: r_xc
226  REAL*8 :: r_alphac
227  REAL*8 :: r_rho_c
228  REAL*8 :: r_red_grav
229  REAL*8 :: r_sp_heat_mix
230 
231  ! compute solid mass fractions
232  IF ( r_qj(1) .GT. 1.d-25 ) THEN
233 
234  r_xs(1:n_solid) = r_qj(5:4+n_solid) / r_qj(1)
235 
236  ELSE
237 
238  r_xs(1:n_solid) = 0.d0
239 
240  END IF
241 
242  r_xs_tot = sum(r_xs)
243 
244  IF ( gas_flag .AND. liquid_flag ) THEN
245 
246  ! compute liquid mass fraction
247  IF ( r_qj(1) .GT. 1.d-25 ) THEN
248 
249  r_xl = r_qj(n_vars) / r_qj(1)
250 
251  ELSE
252 
253  r_xl = 0.d0
254 
255  END IF
256 
257  ! compute carrier phase (gas) mass fraction
258  r_xc = 1.d0 - r_xs_tot - r_xl
259 
260  ! specific heat of the mixutre: mass average of sp. heat pf phases
261  r_sp_heat_mix = sum( r_xs(1:n_solid) * sp_heat_s(1:n_solid) ) &
262  + r_xl * sp_heat_l + r_xc * sp_heat_c
263 
264  ELSE
265 
266  ! compute carrier phase (gas or liquid) mass fraction
267  r_xc = 1.d0 - r_xs_tot
268 
269  ! specific heaf of the mixutre: mass average of sp. heat pf phases
270  r_sp_heat_mix = sum( r_xs(1:n_solid) * sp_heat_s(1:n_solid) ) &
271  + r_xc * sp_heat_c
272 
273  END IF
274 
275  ! compute temperature from energy
276  IF ( r_qj(1) .GT. 1.d-25 ) THEN
277 
278  IF ( energy_flag ) THEN
279 
280  r_t = ( r_qj(4) - ( r_qj(2)**2 + r_qj(3)**2 ) / ( 2.d0*r_qj(1) ) ) / &
281  ( r_qj(1) * r_sp_heat_mix )
282 
283  ELSE
284 
285  r_t = r_qj(4) / ( r_qj(1) * r_sp_heat_mix )
286 
287  END IF
288 
289  IF ( r_t .LE. 0.d0 ) r_t = t_ambient
290 
291  ELSE
292 
293  r_t = t_ambient
294 
295  END IF
296 
297  IF ( gas_flag ) THEN
298 
299  ! carrier phase is gas
300  r_rho_c = pres / ( sp_gas_const_a * r_t )
301  sp_heat_c = sp_heat_a
302 
303  ELSE
304 
305  ! carrier phase is liquid
306  r_rho_c = rho_l
307  sp_heat_c = sp_heat_l
308 
309  END IF
310 
311  IF ( gas_flag .AND. liquid_flag ) THEN
312 
313  r_inv_rhom = ( sum(r_xs(1:n_solid) / rho_s(1:n_solid)) + r_xl / rho_l &
314  + r_xc / r_rho_c )
315 
316  r_rho_m = 1.d0 / r_inv_rhom
317 
318  r_alphal = r_xl * r_rho_m / rho_l
319 
320  ELSE
321 
322  r_inv_rhom = ( sum(r_xs(1:n_solid) / rho_s(1:n_solid)) + r_xc / r_rho_c )
323 
324  r_rho_m = 1.d0 / r_inv_rhom
325 
326  END IF
327 
328  ! convert from mass fraction to volume fraction
329  r_alphas(1:n_solid) = r_xs(1:n_solid) * r_rho_m / rho_s(1:n_solid)
330 
331  ! convert from mass fraction to volume fraction
332  r_alphac = r_xc * r_rho_m / r_rho_c
333 
334  r_h = r_qj(1) / r_rho_m
335 
336  ! reduced gravity
337  r_red_grav = ( r_rho_m - rho_a_amb ) / r_rho_m * grav
338 
339  ! velocity components
340  IF ( r_qj(1) .GT. eps_sing ) THEN
341 
342  r_u = r_qj(2) / r_qj(1)
343  r_v = r_qj(3) / r_qj(1)
344 
345  ELSE
346 
347  r_u = dsqrt(2.d0) * r_qj(1) * r_qj(2) / dsqrt( r_qj(1)**4 + eps_sing**4 )
348  r_v = dsqrt(2.d0) * r_qj(1) * r_qj(3) / dsqrt( r_qj(1)**4 + eps_sing**4 )
349 
350  END IF
351 
352  ! Richardson number
353  IF ( ( r_u**2 + r_v**2 ) .GT. 0.d0 ) THEN
354 
355  r_ri = r_red_grav * r_h / ( r_u**2 + r_v**2 )
356 
357  ELSE
358 
359  r_ri = 0.d0
360 
361  END IF
362 
363  RETURN
364 
365  END SUBROUTINE r_phys_var
366 
367  !******************************************************************************
369  !
380  !
383  !
385  !******************************************************************************
386 
387  SUBROUTINE c_phys_var(c_qj,h,u,v,T,rho_m,red_grav,alphas)
388 
389  USE complexify
390  USE parameters_2d, ONLY : eps_sing
391  IMPLICIT none
392 
393  COMPLEX*16, INTENT(IN) :: c_qj(n_vars)
394  COMPLEX*16, INTENT(OUT) :: h
395  COMPLEX*16, INTENT(OUT) :: u
396  COMPLEX*16, INTENT(OUT) :: v
397  COMPLEX*16, INTENT(OUT) :: T
398  COMPLEX*16, INTENT(OUT) :: rho_m
399  COMPLEX*16, INTENT(OUT) :: red_grav
400  COMPLEX*16, INTENT(OUT) :: alphas(n_solid)
401 
402  COMPLEX*16 :: inv_rhom
403  COMPLEX*16 :: xs(n_solid)
404  COMPLEX*16 :: xs_tot
405  COMPLEX*16 :: Ri
406  COMPLEX*16 :: xl
407  COMPLEX*16 :: xc
408  COMPLEX*16 :: alphal
409  COMPLEX*16 :: alphac
410  COMPLEX*16 :: sp_heat_mix
411  COMPLEX*16 :: rho_c
412 
413 
414  ! compute solid mass fractions
415  IF ( dble(c_qj(1)) .GT. 1.d-25 ) THEN
416 
417  xs(1:n_solid) = c_qj(5:4+n_solid) / c_qj(1)
418 
419  ELSE
420 
421  xs(1:n_solid) = dcmplx(0.d0,0.d0)
422 
423  END IF
424 
425  xs_tot = sum(xs)
426 
427  IF ( gas_flag .AND. liquid_flag ) THEN
428 
429  ! compute liquid mass fraction
430  IF ( dble(c_qj(1)) .GT. 1.d-25 ) THEN
431 
432  xl = c_qj(n_vars) / c_qj(1)
433 
434  ELSE
435 
436  xl = dcmplx(0.d0,0.d0)
437 
438  END IF
439 
440  ! compute carrier phase (gas) mass fraction
441  xc = dcmplx(1.d0,0.d0) - xs_tot - xl
442 
443  ! specific heaf of the mixutre: mass average of sp. heat pf phases
444  sp_heat_mix = sum( xs(1:n_solid) * sp_heat_s(1:n_solid) ) &
445  + xl * sp_heat_l + xc * sp_heat_c
446 
447  ELSE
448 
449  ! compute carrier phase (gas or liquid) mass fraction
450  xc = dcmplx(1.d0,0.d0) - xs_tot
451 
452  ! specific heaf of the mixutre: mass average of sp. heat pf phases
453  sp_heat_mix = sum( xs(1:n_solid) * sp_heat_s(1:n_solid) ) + xc * sp_heat_c
454 
455  END IF
456 
457  ! compute temperature from energy
458  IF ( dble(c_qj(1)) .GT. 1.d-25 ) THEN
459 
460  IF ( energy_flag ) THEN
461 
462  t = ( c_qj(4) - ( c_qj(2)**2 + c_qj(3)**2 ) / ( 2.d0*c_qj(1) ) ) / &
463  ( c_qj(1) * sp_heat_mix )
464 
465  ELSE
466 
467  t = c_qj(4) / ( c_qj(1) * sp_heat_mix )
468 
469  END IF
470 
471  IF ( dble(t) .LE. 0.d0 ) t = dcmplx(t_ambient,0.d0)
472 
473  ELSE
474 
475  t = dcmplx(t_ambient,0.d0)
476 
477  END IF
478 
479  IF ( gas_flag ) THEN
480 
481  ! carrier phase is gas
482  rho_c = pres / ( sp_gas_const_a * t )
483  sp_heat_c = sp_heat_a
484 
485  ELSE
486 
487  ! carrier phase is liquid
488  rho_c = rho_l
489  sp_heat_c = sp_heat_l
490 
491  END IF
492 
493  IF ( gas_flag .AND. liquid_flag ) THEN
494 
495  inv_rhom = ( sum(xs(1:n_solid) / rho_s(1:n_solid)) + xl / rho_l &
496  + xc / rho_c )
497 
498  rho_m = 1.d0 / inv_rhom
499 
500  alphal = xl * rho_m / rho_l
501 
502  ELSE
503 
504  inv_rhom = ( sum(xs(1:n_solid) / rho_s(1:n_solid)) + xc / rho_c )
505 
506  rho_m = 1.d0 / inv_rhom
507 
508  END IF
509 
510  ! convert from mass fraction to volume fraction
511  alphas(1:n_solid) = xs(1:n_solid) * rho_m / rho_s(1:n_solid)
512 
513  ! convert from mass fraction to volume fraction
514  alphac = xc * rho_m / rho_c
515 
516  h = c_qj(1) / rho_m
517 
518  ! reduced gravity
519  red_grav = ( rho_m - rho_a_amb ) / rho_m * grav
520 
521  ! velocity components
522  IF ( dble( c_qj(1) ) .GT. eps_sing ) THEN
523 
524  u = c_qj(2) / c_qj(1)
525  v = c_qj(3) / c_qj(1)
526 
527  ELSE
528 
529  u = dsqrt(2.d0) * c_qj(1) * c_qj(2) / cdsqrt( c_qj(1)**4 + eps_sing**4 )
530  v = dsqrt(2.d0) * c_qj(1) * c_qj(3) / cdsqrt( c_qj(1)**4 + eps_sing**4 )
531 
532  END IF
533 
534  ! Richardson number
535  IF ( dble( u**2 + v**2 ) .GT. 0.d0 ) THEN
536 
537  ri = red_grav * h / ( u**2 + v**2 )
538 
539  ELSE
540 
541  ri = dcmplx(0.d0,0.d0)
542 
543  END IF
544 
545  RETURN
546 
547  END SUBROUTINE c_phys_var
548 
549 
550  !******************************************************************************
552  !
561  !
564  !
566  !******************************************************************************
567 
568  SUBROUTINE mixt_var(qpj,r_Ri,r_rho_m,r_rho_c,r_red_grav)
569 
570  IMPLICIT none
571 
572  REAL*8, INTENT(IN) :: qpj(n_vars+2)
573  REAL*8, INTENT(OUT) :: r_Ri
574  REAL*8, INTENT(OUT) :: r_rho_m
575  REAL*8, INTENT(OUT) :: r_rho_c
576  REAL*8, INTENT(OUT) :: r_red_grav
577  REAL*8 :: r_u
578  REAL*8 :: r_v
579  REAL*8 :: r_h
580  REAL*8 :: r_alphas(n_solid)
581  REAL*8 :: r_T
582  REAL*8 :: r_alphal
583 
584  r_h = qpj(1)
585 
586  IF ( qpj(1) .LE. 0.d0 ) THEN
587 
588  r_u = 0.d0
589  r_v = 0.d0
590  r_t = t_ambient
591  r_alphas(1:n_solid) = 0.d0
592  r_red_grav = 0.d0
593  r_ri = 0.d0
594  r_rho_m = rho_a_amb
595  IF ( gas_flag .AND. liquid_flag ) r_alphal = 0.d0
596 
597  RETURN
598 
599  END IF
600 
601  r_u = qpj(n_vars+1)
602  r_v = qpj(n_vars+2)
603  r_t = qpj(4)
604  r_alphas(1:n_solid) = qpj(5:4+n_solid)
605 
606  IF ( gas_flag ) THEN
607 
608  ! continuous phase is air
609  r_rho_c = pres / ( sp_gas_const_a * r_t )
610 
611  ELSE
612 
613  ! continuous phase is liquid
614  r_rho_c = rho_l
615 
616  END IF
617 
618  IF ( gas_flag .AND. liquid_flag ) THEN
619 
620  r_alphal = qpj(n_vars)
621 
622  ! density of mixture of carrier (gas), liquid and solids
623  r_rho_m = ( 1.d0 - sum(r_alphas) - r_alphal ) * r_rho_c &
624  + sum( r_alphas * rho_s ) + r_alphal * rho_l
625 
626  ELSE
627 
628  ! density of mixture of carrier phase and solids
629  r_rho_m = ( 1.d0 - sum(r_alphas) ) * r_rho_c + sum( r_alphas * rho_s )
630 
631  END IF
632 
633  ! reduced gravity
634  r_red_grav = ( r_rho_m - rho_a_amb ) / r_rho_m * grav
635 
636  ! Richardson number
637  IF ( ( r_u**2 + r_v**2 ) .GT. 0.d0 ) THEN
638 
639  r_ri = r_red_grav * r_h / ( r_u**2 + r_v**2 )
640 
641  ELSE
642 
643  r_ri = 0.d0
644 
645  END IF
646 
647  RETURN
648 
649  END SUBROUTINE mixt_var
650 
651  !******************************************************************************
653  !
669  !
671  !
674  !
675  !******************************************************************************
676 
677  SUBROUTINE qc_to_qp(qc,qp)
678 
679  IMPLICIT none
680 
681  REAL*8, INTENT(IN) :: qc(n_vars)
682  REAL*8, INTENT(OUT) :: qp(n_vars+2)
683 
684  REAL*8 :: r_h
685  REAL*8 :: r_u
686  REAL*8 :: r_v
687  REAL*8 :: r_alphas(n_solid)
688  REAL*8 :: r_rho_m
689  REAL*8 :: r_T
690  REAL*8 :: r_alphal
691 
692 
693  IF ( qc(1) .GT. 0.d0 ) THEN
694 
695  CALL r_phys_var(qc,r_h,r_u,r_v,r_alphas,r_rho_m,r_t,r_alphal)
696 
697  qp(1) = r_h
698 
699  qp(2) = r_h*r_u
700  qp(3) = r_h*r_v
701 
702  qp(4) = r_t
703  qp(5:4+n_solid) = r_alphas(1:n_solid)
704 
705  IF ( gas_flag .AND. liquid_flag ) qp(n_vars) = r_alphal
706 
707  qp(n_vars+1) = r_u
708  qp(n_vars+2) = r_v
709 
710  ELSE
711 
712  qp(1) = 0.d0 ! h
713  qp(2) = 0.d0 ! hu
714  qp(3) = 0.d0 ! hv
715  qp(4) = t_ambient ! T
716  qp(5:n_vars) = 0.d0 ! alphas
717  qp(n_vars+1) = 0.d0 ! u
718  qp(n_vars+2) = 0.d0 ! v
719 
720  END IF
721 
722  RETURN
723 
724  END SUBROUTINE qc_to_qp
725 
726  !******************************************************************************
728  !
743  !
745  !
748  !
749  !******************************************************************************
750 
751  SUBROUTINE qp_to_qc(qp,B,qc)
752 
753  IMPLICIT none
754 
755  REAL*8, INTENT(IN) :: qp(n_vars+2)
756  REAL*8, INTENT(IN) :: B
757  REAL*8, INTENT(OUT) :: qc(n_vars)
758 
759  REAL*8 :: r_sp_heat_mix
760  REAL*8 :: sum_sl
761 
762  REAL*8 :: r_u
763  REAL*8 :: r_v
764  REAL*8 :: r_h
765  REAL*8 :: r_hu
766  REAL*8 :: r_hv
767  REAL*8 :: r_alphas(n_solid)
768  REAL*8 :: r_xl
769  REAL*8 :: r_xc
770  REAL*8 :: r_T
771  REAL*8 :: r_alphal
772  REAL*8 :: r_alphac
773  REAL*8 :: r_rho_m
774  REAL*8 :: r_rho_c
775  REAL*8 :: r_xs(n_solid)
776 
777  r_h = qp(1)
778 
779  IF ( r_h .GT. 0.d0 ) THEN
780 
781  r_hu = qp(2)
782  r_hv = qp(3)
783 
784  r_u = qp(n_vars+1)
785  r_v = qp(n_vars+2)
786 
787  ELSE
788 
789  r_hu = 0.d0
790  r_hv = 0.d0
791 
792  r_u = 0.d0
793  r_v = 0.d0
794 
795  qc(1:n_vars) = 0.d0
796  RETURN
797 
798  END IF
799 
800  r_t = qp(4)
801 
802  r_alphas(1:n_solid) = qp(5:4+n_solid)
803 
804  IF ( gas_flag ) THEN
805 
806  ! carrier phase is gas
807  r_rho_c = pres / ( sp_gas_const_a * r_t )
808  sp_heat_c = sp_heat_a
809 
810  ELSE
811 
812  ! carrier phase is liquid
813  r_rho_c = rho_l
814  sp_heat_c = sp_heat_l
815 
816  END IF
817 
818  IF ( gas_flag .AND. liquid_flag ) THEN
819 
820  ! mixture of gas, liquid and solid
821  r_alphal = qp(n_vars)
822 
823  ! check and correction on dispersed phases volume fractions
824  IF ( ( sum(r_alphas) + r_alphal ) .GT. 1.d0 ) THEN
825 
826  sum_sl = sum(r_alphas) + r_alphal
827  r_alphas(1:n_solid) = r_alphas(1:n_solid) / sum_sl
828  r_alphal = r_alphal / sum_sl
829 
830  ELSEIF ( ( sum(r_alphas) + r_alphal ) .LT. 0.d0 ) THEN
831 
832  r_alphas(1:n_solid) = 0.d0
833  r_alphal = 0.d0
834 
835  END IF
836 
837  ! carrier phase volume fraction
838  r_alphac = 1.d0 - sum(r_alphas) - r_alphal
839 
840  ! volume averaged mixture density: carrier (gas) + solids + liquid
841  r_rho_m = r_alphac * r_rho_c + sum( r_alphas * rho_s ) + r_alphal * rho_l
842 
843  ! liquid mass fraction
844  r_xl = r_alphal * rho_l / r_rho_m
845 
846  ! solid mass fractions
847  r_xs(1:n_solid) = r_alphas(1:n_solid) * rho_s(1:n_solid) / r_rho_m
848 
849  ! carrier (gas) mass fraction
850  r_xc = r_alphac * r_rho_c / r_rho_m
851 
852  ! mass averaged mixture specific heat
853  r_sp_heat_mix = sum( r_xs*sp_heat_s ) + r_xl*sp_heat_l + r_xc*sp_heat_c
854 
855  ELSE
856 
857  ! mixture of carrier phase ( gas or liquid ) and solid
858 
859  ! check and corrections on dispersed phases
860  IF ( sum(r_alphas) .GT. 1.d0 ) THEN
861 
862  r_alphas(1:n_solid) = r_alphas(1:n_solid) / sum(r_alphas)
863 
864  ELSEIF ( sum(r_alphas).LT. 0.d0 ) THEN
865 
866  r_alphas(1:n_solid) = 0.d0
867 
868  END IF
869 
870  ! carrier (gas or liquid) volume fraction
871  r_alphac = 1.d0 - sum(r_alphas)
872 
873  ! volume averaged mixture density: carrier (gas or liquid) + solids
874  r_rho_m = r_alphac * r_rho_c + sum( r_alphas * rho_s )
875 
876  ! solid mass fractions
877  r_xs(1:n_solid) = r_alphas(1:n_solid) * rho_s(1:n_solid) / r_rho_m
878 
879  ! carrier (gas or liquid) mass fraction
880  r_xc = r_alphac * r_rho_c / r_rho_m
881 
882  ! mass averaged mixture specific heat
883  r_sp_heat_mix = sum( r_xs * sp_heat_s ) + r_xc * sp_heat_c
884 
885  END IF
886 
887  qc(1) = r_rho_m * r_h
888  qc(2) = r_rho_m * r_hu
889  qc(3) = r_rho_m * r_hv
890 
891  IF ( energy_flag ) THEN
892 
893  IF ( r_h .GT. 0.d0 ) THEN
894 
895  ! total energy (internal and kinetic)
896  qc(4) = r_h * r_rho_m * ( r_sp_heat_mix * r_t &
897  + 0.5d0 * ( r_hu**2 + r_hv**2 ) / r_h**2 )
898 
899  ELSE
900 
901  qc(4) = 0.d0
902 
903  END IF
904 
905  ELSE
906 
907  ! internal energy
908  qc(4) = r_h * r_rho_m * r_sp_heat_mix * r_t
909 
910  END IF
911 
912  qc(5:4+n_solid) = r_h * r_alphas(1:n_solid) * rho_s(1:n_solid)
913 
914  IF ( gas_flag .AND. liquid_flag ) qc(n_vars) = r_h * r_alphal * rho_l
915 
916  RETURN
917 
918  END SUBROUTINE qp_to_qc
919 
920  !******************************************************************************
922  !
931  !******************************************************************************
932 
933  SUBROUTINE qp_to_qp2(qpj,Bj,qp2j)
934 
935  IMPLICIT none
936 
937  REAL*8, INTENT(IN) :: qpj(n_vars)
938  REAL*8, INTENT(IN) :: Bj
939  REAL*8, INTENT(OUT) :: qp2j(3)
940 
941  qp2j(1) = qpj(1) + bj
942 
943  IF ( qpj(1) .LE. 0.d0 ) THEN
944 
945  qp2j(2) = 0.d0
946  qp2j(3) = 0.d0
947 
948  ELSE
949 
950  qp2j(2) = qpj(2)/qpj(1)
951  qp2j(3) = qpj(3)/qpj(1)
952 
953  END IF
954 
955  RETURN
956 
957  END SUBROUTINE qp_to_qp2
958 
959  !******************************************************************************
961  !
970  !******************************************************************************
971 
972  SUBROUTINE eval_local_speeds_x(qpj,vel_min,vel_max)
973 
974  IMPLICIT none
975 
976  REAL*8, INTENT(IN) :: qpj(n_vars+2)
977 
978  REAL*8, INTENT(OUT) :: vel_min(n_vars) , vel_max(n_vars)
979 
980  REAL*8 :: r_h
981  REAL*8 :: r_u
982  REAL*8 :: r_v
983  REAL*8 :: r_Ri
984  REAL*8 :: r_rho_m
985  REAL*8 :: r_rho_c
986  REAL*8 :: r_red_grav
987 
988  CALL mixt_var(qpj,r_ri,r_rho_m,r_rho_c,r_red_grav)
989 
990  r_h = qpj(1)
991  r_u = qpj(n_vars+1)
992  r_v = qpj(n_vars+2)
993 
994  IF ( r_red_grav * r_h .LT. 0.d0 ) THEN
995 
996  vel_min(1:n_eqns) = r_u
997  vel_max(1:n_eqns) = r_u
998 
999  ELSE
1000 
1001  vel_min(1:n_eqns) = r_u - dsqrt( r_red_grav * r_h )
1002  vel_max(1:n_eqns) = r_u + dsqrt( r_red_grav * r_h )
1003 
1004  END IF
1005 
1006  RETURN
1007 
1008  END SUBROUTINE eval_local_speeds_x
1009 
1010  !******************************************************************************
1012  !
1021  !******************************************************************************
1022 
1023  SUBROUTINE eval_local_speeds_y(qpj,vel_min,vel_max)
1024 
1025  IMPLICIT none
1026 
1027  REAL*8, INTENT(IN) :: qpj(n_vars+2)
1028  REAL*8, INTENT(OUT) :: vel_min(n_vars) , vel_max(n_vars)
1029 
1030  REAL*8 :: r_h
1031  REAL*8 :: r_u
1032  REAL*8 :: r_v
1033  REAL*8 :: r_Ri
1034  REAL*8 :: r_rho_m
1035  REAL*8 :: r_rho_c
1036  REAL*8 :: r_red_grav
1037 
1038  CALL mixt_var(qpj,r_ri,r_rho_m,r_rho_c,r_red_grav)
1039 
1040  r_h = qpj(1)
1041  r_u = qpj(n_vars+1)
1042  r_v = qpj(n_vars+2)
1043 
1044  IF ( r_red_grav * r_h .LT. 0.d0 ) THEN
1045 
1046  vel_min(1:n_eqns) = r_v
1047  vel_max(1:n_eqns) = r_v
1048 
1049  ELSE
1050 
1051  vel_min(1:n_eqns) = r_v - dsqrt( r_red_grav * r_h )
1052  vel_max(1:n_eqns) = r_v + dsqrt( r_red_grav * r_h )
1053 
1054  END IF
1055 
1056  RETURN
1057 
1058  END SUBROUTINE eval_local_speeds_y
1059 
1060  !******************************************************************************
1062  !
1070  !
1073  !
1074  !******************************************************************************
1075 
1076  SUBROUTINE eval_fluxes(qcj,qpj,dir,flux)
1077 
1078  IMPLICIT none
1079 
1080  REAL*8, INTENT(IN) :: qcj(n_vars)
1081  REAL*8, INTENT(IN) :: qpj(n_vars+2)
1082  INTEGER, INTENT(IN) :: dir
1083 
1084  REAL*8, INTENT(OUT) :: flux(n_eqns)
1085 
1086  REAL*8 :: r_h
1087  REAL*8 :: r_u
1088  REAL*8 :: r_v
1089  REAL*8 :: r_Ri
1090  REAL*8 :: r_rho_m
1091  REAL*8 :: r_rho_c
1092  REAL*8 :: r_red_grav
1093 
1094  pos_thick:IF ( qcj(1) .GT. 0.d0 ) THEN
1095 
1096  r_h = qpj(1)
1097  r_u = qpj(n_vars+1)
1098  r_v = qpj(n_vars+2)
1099 
1100  CALL mixt_var(qpj,r_ri,r_rho_m,r_rho_c,r_red_grav)
1101 
1102  IF ( dir .EQ. 1 ) THEN
1103 
1104  ! Mass flux in x-direction: u * ( rhom * h )
1105  flux(1) = r_u * qcj(1)
1106 
1107  ! x-momentum flux in x-direction + hydrostatic pressure term
1108  flux(2) = r_u * qcj(2) + 0.5d0 * r_rho_m * r_red_grav * r_h**2
1109 
1110  ! y-momentum flux in x-direction: u * ( rho * h * v )
1111  flux(3) = r_u * qcj(3)
1112 
1113  IF ( energy_flag ) THEN
1114 
1115  ! ENERGY flux in x-direction
1116  flux(4) = r_u * ( qcj(4) + 0.5d0 * r_rho_m * r_red_grav * r_h**2 )
1117 
1118  ELSE
1119 
1120  ! Temperature flux in x-direction: u * ( h * T )
1121  flux(4) = r_u * qcj(4)
1122 
1123  END IF
1124 
1125  ! Mass flux of solid in x-direction: u * ( h * alphas * rhos )
1126  flux(5:4+n_solid) = r_u * qcj(5:4+n_solid)
1127 
1128  ! Solid flux can't be larger than total flux
1129  IF ( ( flux(1) .GT. 0.d0 ) .AND. ( sum(flux(5:4+n_solid)) / flux(1) &
1130  .GT. 1.d0 ) ) THEN
1131 
1132  flux(5:4+n_solid) = &
1133  flux(5:4+n_solid) / sum(flux(5:4+n_solid)) * flux(1)
1134 
1135  END IF
1136 
1137  ! Mass flux of liquid in x-direction: u * ( h * alphal * rhol )
1138  IF ( gas_flag .AND. liquid_flag ) flux(n_vars) = r_u * qcj(n_vars)
1139 
1140  ELSEIF ( dir .EQ. 2 ) THEN
1141 
1142  ! flux G (derivated wrt y in the equations)
1143  flux(1) = r_v * qcj(1)
1144 
1145  flux(2) = r_v * qcj(2)
1146 
1147  flux(3) = r_v * qcj(3) + 0.5d0 * r_rho_m * r_red_grav * r_h**2
1148 
1149  IF ( energy_flag ) THEN
1150 
1151  ! ENERGY flux in x-direction
1152  flux(4) = r_v * ( qcj(4) + 0.5d0 * r_rho_m * r_red_grav * r_h**2 )
1153 
1154  ELSE
1155 
1156  ! Temperature flux in y-direction
1157  flux(4) = r_v * qcj(4)
1158 
1159  END IF
1160 
1161  ! Mass flux of solid in y-direction: v * ( h * alphas * rhos )
1162  flux(5:4+n_solid) = r_v * qcj(5:4+n_solid)
1163 
1164  ! Solid flux can't be larger than total flux
1165  IF ( ( flux(1) .GT. 0.d0 ) .AND. ( sum(flux(5:4+n_solid)) / flux(1) &
1166  .GT. 1.d0 ) ) THEN
1167 
1168  flux(5:4+n_solid) = &
1169  flux(5:4+n_solid) / sum(flux(5:4+n_solid)) * flux(1)
1170 
1171  END IF
1172 
1173  ! Mass flux of liquid in x-direction: u * ( h * alphal * rhol )
1174  IF ( gas_flag .AND. liquid_flag ) flux(n_vars) = r_v * qcj(n_vars)
1175 
1176  END IF
1177 
1178  ELSE
1179 
1180  flux(1:n_eqns) = 0.d0
1181 
1182  ENDIF pos_thick
1183 
1184  RETURN
1185 
1186  END SUBROUTINE eval_fluxes
1187 
1188  !******************************************************************************
1190  !
1200  !
1203  !
1204  !******************************************************************************
1205 
1206  SUBROUTINE eval_nonhyperbolic_terms( c_qj , c_nh_term_impl , r_qj , &
1207  r_nh_term_impl )
1208 
1209  USE complexify
1210 
1211  IMPLICIT NONE
1212 
1213  COMPLEX*16, INTENT(IN), OPTIONAL :: c_qj(n_vars)
1214  COMPLEX*16, INTENT(OUT), OPTIONAL :: c_nh_term_impl(n_eqns)
1215  REAL*8, INTENT(IN), OPTIONAL :: r_qj(n_vars)
1216  REAL*8, INTENT(OUT), OPTIONAL :: r_nh_term_impl(n_eqns)
1217 
1218  COMPLEX*16 :: h
1219  COMPLEX*16 :: u
1220  COMPLEX*16 :: v
1221  COMPLEX*16 :: T
1222  COMPLEX*16 :: rho_m
1223  COMPLEX*16 :: red_grav
1224  COMPLEX*16 :: alphas(n_solid)
1225 
1226  COMPLEX*16 :: qj(n_vars)
1227  COMPLEX*16 :: nh_term(n_eqns)
1228  COMPLEX*16 :: forces_term(n_eqns)
1229 
1230  COMPLEX*16 :: mod_vel
1231  COMPLEX*16 :: gamma
1232  REAL*8 :: h_threshold
1233 
1234  INTEGER :: i
1235 
1236  !--- Lahars rheology model variables
1237 
1239  COMPLEX*16 :: Tc
1240 
1241  COMPLEX*16 :: expA , expB
1242 
1244  COMPLEX*16 :: alpha1 ! (units: kg m-1 s-1 )
1245 
1247  COMPLEX*16 :: fluid_visc
1248 
1250  COMPLEX*16 :: s_f
1251 
1253  COMPLEX*16 :: s_v
1254 
1256  COMPLEX*16 :: s_td
1257 
1258  IF ( ( thermal_conductivity .GT. 0.d0 ) .OR. ( emme .GT. 0.d0 ) ) THEN
1259 
1260  h_threshold = 1.d-10
1261 
1262  ELSE
1263 
1264  h_threshold = 0.d0
1265 
1266  END IF
1267 
1268  IF ( present(c_qj) .AND. present(c_nh_term_impl) ) THEN
1269 
1270  qj = c_qj
1271 
1272  ELSEIF ( present(r_qj) .AND. present(r_nh_term_impl) ) THEN
1273 
1274  DO i = 1,n_vars
1275 
1276  qj(i) = dcmplx( r_qj(i) )
1277 
1278  END DO
1279 
1280  ELSE
1281 
1282  WRITE(*,*) 'Constitutive, eval_fluxes: problem with arguments'
1283  stop
1284 
1285  END IF
1286 
1287  ! initialize and evaluate the forces terms
1288  forces_term(1:n_eqns) = dcmplx(0.d0,0.d0)
1289 
1290  IF (rheology_flag) THEN
1291 
1292  CALL c_phys_var(qj,h,u,v,t,rho_m,red_grav,alphas)
1293 
1294  mod_vel = cdsqrt( u**2 + v**2 )
1295 
1296  ! Voellmy Salm rheology
1297  IF ( rheology_model .EQ. 1 ) THEN
1298 
1299  IF ( dble(mod_vel) .NE. 0.d0 ) THEN
1300 
1301  ! IMPORTANT: grav3_surf is always negative
1302  forces_term(2) = forces_term(2) - rho_m * ( u / mod_vel ) * &
1303  ( grav / xi ) * mod_vel ** 2
1304 
1305  forces_term(3) = forces_term(3) - rho_m * ( v / mod_vel ) * &
1306  ( grav / xi ) * mod_vel ** 2
1307 
1308  ENDIF
1309 
1310  ! Plastic rheology
1311  ELSEIF ( rheology_model .EQ. 2 ) THEN
1312 
1313  IF ( dble(mod_vel) .NE. 0.d0 ) THEN
1314 
1315  forces_term(2) = forces_term(2) - rho_m * tau * (u/mod_vel)
1316 
1317  forces_term(3) = forces_term(3) - rho_m * tau * (v/mod_vel)
1318 
1319  ENDIF
1320 
1321  ! Temperature dependent rheology
1322  ELSEIF ( rheology_model .EQ. 3 ) THEN
1323 
1324  IF ( dble(h) .GT. h_threshold ) THEN
1325 
1326  ! Equation 6 from Costa & Macedonio, 2005
1327  gamma = 3.d0 * nu_ref / h * cdexp( - visc_par * ( t - t_ref ) )
1328 
1329  ELSE
1330 
1331  ! Equation 6 from Costa & Macedonio, 2005
1332  gamma = 3.d0 * nu_ref / h_threshold * cdexp( - visc_par &
1333  * ( t - t_ref ) )
1334 
1335  END IF
1336 
1337  IF ( dble(mod_vel) .NE. 0.d0 ) THEN
1338 
1339  ! Last R.H.S. term in equation 2 from Costa & Macedonio, 2005
1340  forces_term(2) = forces_term(2) - rho_m * gamma * u
1341 
1342  ! Last R.H.S. term in equation 3 from Costa & Macedonio, 2005
1343  forces_term(3) = forces_term(3) - rho_m * gamma * v
1344 
1345  ENDIF
1346 
1347  ! Lahars rheology (O'Brien 1993, FLO2D)
1348  ELSEIF ( rheology_model .EQ. 4 ) THEN
1349 
1350  ! alpha1 here has units: kg m-1 s-1
1351  ! in Table 2 from O'Brien 1988, the values reported have different
1352  ! units ( poises). 1poises = 0.1 kg m-1 s-1
1353 
1354  h_threshold = 1.d-20
1355 
1356  ! convert from Kelvin to Celsius
1357  tc = t - 273.15d0
1358 
1359  ! the dependance of viscosity on temperature is modeled with the
1360  ! equation presented at:
1361  ! https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118131473.app3
1362  !
1363  ! In addition, we use a reference value provided in input at a
1364  ! reference temperature. This value is used to scale the equation
1365  IF ( tc .LT. 20.d0 ) THEN
1366 
1367  expa = 1301.d0 / ( 998.333d0 + 8.1855d0 * ( tc - 20.d0 ) &
1368  + 0.00585d0 * ( tc - 20.d0 )**2 ) - 1.30223d0
1369 
1370  alpha1 = alpha1_coeff * 1.d-3 * 10.d0**expa
1371 
1372  ELSE
1373 
1374  expb = ( 1.3272d0 * ( 20.d0 - tc ) - 0.001053d0 * &
1375  ( tc - 20.d0 )**2 ) / ( tc + 105.0d0 )
1376 
1377  alpha1 = alpha1_coeff * 1.002d-3 * 10.d0**expb
1378 
1379  END IF
1380 
1381  ! Fluid viscosity
1382  fluid_visc = alpha1 * cdexp( beta1 * sum(alphas) )
1383 
1384  IF ( dble(h) .GT. h_threshold ) THEN
1385 
1386  ! Viscous slope component (dimensionless)
1387  s_v = kappa * fluid_visc * mod_vel / ( 8.d0 * rho_m * grav *h**2 )
1388 
1389  ! Turbulent dispersive component (dimensionless)
1390  s_td = n_td**2 * mod_vel**2 / ( h**(4.d0/3.d0) )
1391 
1392  ELSE
1393 
1394  ! Viscous slope component (dimensionless)
1395  s_v = kappa * fluid_visc * mod_vel / ( 8.d0 * h_threshold**2 )
1396 
1397  ! Turbulent dispersive components (dimensionless)
1398  s_td = n_td**2 * (mod_vel**2) / ( h_threshold**(4.d0/3.d0) )
1399 
1400  END IF
1401 
1402  ! Total implicit friction slope (dimensionless)
1403  s_f = s_v + s_td
1404 
1405  IF ( mod_vel .GT. 0.d0 ) THEN
1406 
1407  ! same units of dqc(2)/dt: kg m-1 s-2
1408  forces_term(2) = forces_term(2) - grav * rho_m * h * &
1409  ( u / mod_vel ) * s_f
1410 
1411  ! same units of dqc(3)/dt: kg m-1 s-2
1412  forces_term(3) = forces_term(3) - grav * rho_m * h * &
1413  ( v / mod_vel ) * s_f
1414 
1415  END IF
1416 
1417  ELSEIF ( rheology_model .EQ. 5 ) THEN
1418 
1419  tau = 1.d-3 / ( 1.d0 + 10.d0 * h ) * mod_vel
1420 
1421  IF ( dble(mod_vel) .NE. 0.d0 ) THEN
1422 
1423  forces_term(2) = forces_term(2) - rho_m * tau * ( u / mod_vel )
1424  forces_term(3) = forces_term(3) - rho_m * tau * ( v / mod_vel )
1425 
1426  END IF
1427 
1428 
1429  ELSEIF ( rheology_model .EQ. 6 ) THEN
1430 
1431  IF ( dble(mod_vel) .NE. 0.d0 ) THEN
1432 
1433  forces_term(2) = forces_term(2) - rho_m * ( u / mod_vel ) * &
1434  friction_factor * mod_vel ** 2
1435 
1436  forces_term(3) = forces_term(3) - rho_m * ( v / mod_vel ) * &
1437  friction_factor * mod_vel ** 2
1438 
1439  ENDIF
1440 
1441  ENDIF
1442 
1443  ENDIF
1444 
1445  nh_term = forces_term
1446 
1447  IF ( present(c_qj) .AND. present(c_nh_term_impl) ) THEN
1448 
1449  c_nh_term_impl = nh_term
1450 
1451  ELSEIF ( present(r_qj) .AND. present(r_nh_term_impl) ) THEN
1452 
1453  r_nh_term_impl = dble( nh_term )
1454 
1455  END IF
1456 
1457  RETURN
1458 
1459  END SUBROUTINE eval_nonhyperbolic_terms
1460 
1461  !******************************************************************************
1463  !
1471  !
1474  !
1475  !******************************************************************************
1476 
1477  SUBROUTINE eval_nh_semi_impl_terms( grav3_surf , qcj , nh_semi_impl_term )
1478 
1479  IMPLICIT NONE
1480 
1481  REAL*8, INTENT(IN) :: grav3_surf
1482 
1483  REAL*8, INTENT(IN) :: qcj(n_vars)
1484  REAL*8, INTENT(OUT) :: nh_semi_impl_term(n_eqns)
1485 
1486  REAL*8 :: forces_term(n_eqns)
1487 
1488  REAL*8 :: mod_vel
1489 
1490  REAL*8 :: h_threshold
1491 
1492  !--- Lahars rheology model variables
1494  REAL*8 :: tau_y
1496  REAL*8 :: s_y
1497 
1498  REAL*8 :: r_h
1499  REAL*8 :: r_u
1500  REAL*8 :: r_v
1501  REAL*8 :: r_alphas(n_solid)
1502  REAL*8 :: r_rho_m
1503  REAL*8 :: r_T
1504  REAL*8 :: r_alphal
1505 
1506 
1507  ! initialize and evaluate the forces terms
1508  forces_term(1:n_eqns) = 0.d0
1509 
1510  IF (rheology_flag) THEN
1511 
1512  CALL r_phys_var(qcj,r_h,r_u,r_v,r_alphas,r_rho_m,r_t,r_alphal)
1513 
1514  mod_vel = dsqrt( r_u**2 + r_v**2 )
1515 
1516  ! Voellmy Salm rheology
1517  IF ( rheology_model .EQ. 1 ) THEN
1518 
1519  IF ( mod_vel .GT. 0.d0 ) THEN
1520 
1521  ! units of dqc(2)/dt=d(rho h v)/dt (kg m−1 s−2)
1522  forces_term(2) = forces_term(2) - r_rho_m * ( r_u / mod_vel ) * &
1523  mu * r_h * ( - grav * grav3_surf )
1524 
1525  ! units of dqc(3)/dt=d(rho h v)/dt (kg m−1 s−2)
1526  forces_term(3) = forces_term(3) - r_rho_m * ( r_v / mod_vel ) * &
1527  mu * r_h * ( - grav * grav3_surf )
1528 
1529  END IF
1530 
1531  ! Plastic rheology
1532  ELSEIF ( rheology_model .EQ. 2 ) THEN
1533 
1534 
1535  ! Temperature dependent rheology
1536  ELSEIF ( rheology_model .EQ. 3 ) THEN
1537 
1538 
1539  ! Lahars rheology (O'Brien 1993, FLO2D)
1540  ELSEIF ( rheology_model .EQ. 4 ) THEN
1541 
1542  h_threshold = 1.d-20
1543 
1544  ! Yield strength (units: kg m−1 s−2)
1545  tau_y = alpha2 * dexp( beta2 * sum(r_alphas) )
1546 
1547  IF ( r_h .GT. h_threshold ) THEN
1548 
1549  ! Yield slope component (dimensionless)
1550  s_y = tau_y / ( grav * r_rho_m * r_h )
1551 
1552  ELSE
1553 
1554  ! Yield slope component (dimensionless)
1555  s_y = tau_y / ( grav * r_rho_m * h_threshold )
1556 
1557  END IF
1558 
1559  IF ( mod_vel .GT. 0.d0 ) THEN
1560 
1561  ! units of dqc(2)/dt (kg m−1 s−2)
1562  forces_term(2) = forces_term(2) - grav * r_rho_m * r_h * &
1563  ( r_u / mod_vel ) * s_y
1564 
1565  ! units of dqc(3)/dt (kg m−1 s−2)
1566  forces_term(3) = forces_term(3) - grav * r_rho_m * r_h * &
1567  ( r_v / mod_vel ) * s_y
1568 
1569  END IF
1570 
1571  ELSEIF ( rheology_model .EQ. 5 ) THEN
1572 
1573  ENDIF
1574 
1575  ENDIF
1576 
1577  nh_semi_impl_term = forces_term
1578 
1579  RETURN
1580 
1581  END SUBROUTINE eval_nh_semi_impl_terms
1582 
1583  !******************************************************************************
1585  !
1596  !
1599  !
1600  !******************************************************************************
1601 
1602  SUBROUTINE eval_expl_terms( Bprimej_x , Bprimej_y , source_xy , qpj , qcj , &
1603  expl_term )
1604 
1605  IMPLICIT NONE
1606 
1607  REAL*8, INTENT(IN) :: Bprimej_x
1608  REAL*8, INTENT(IN) :: Bprimej_y
1609  REAL*8, INTENT(IN) :: source_xy
1610 
1611  REAL*8, INTENT(IN) :: qpj(n_vars+2)
1612  REAL*8, INTENT(IN) :: qcj(n_vars)
1613  REAL*8, INTENT(OUT) :: expl_term(n_eqns)
1614 
1615  REAL*8 :: r_h
1616  REAL*8 :: r_u
1617  REAL*8 :: r_v
1618  REAL*8 :: r_Ri
1619  REAL*8 :: r_rho_m
1620  REAL*8 :: r_rho_c
1621  REAL*8 :: r_red_grav
1622 
1623  expl_term(1:n_eqns) = 0.d0
1624 
1625  IF ( qpj(1) .LE. 0.d0 ) RETURN
1626 
1627  r_h = qpj(1)
1628  r_u = qpj(n_vars+1)
1629  r_v = qpj(n_vars+2)
1630 
1631  CALL mixt_var(qpj,r_ri,r_rho_m,r_rho_c,r_red_grav)
1632 
1633  expl_term(2) = r_red_grav * r_rho_m * r_h * bprimej_x
1634 
1635  expl_term(3) = r_red_grav * r_rho_m * r_h * bprimej_y
1636 
1637  IF ( energy_flag ) THEN
1638 
1639  expl_term(4) = r_red_grav * r_rho_m * r_h * ( r_u * bprimej_x &
1640  + r_v * bprimej_y )
1641 
1642  ELSE
1643 
1644  expl_term(4) = 0.d0
1645 
1646  END IF
1647 
1648  RETURN
1649 
1650  END SUBROUTINE eval_expl_terms
1651 
1652  !******************************************************************************
1654  !
1661  !
1664  !
1665  !******************************************************************************
1666 
1667  SUBROUTINE eval_erosion_dep_term( qpj , dt , erosion_term , deposition_term )
1668 
1669  IMPLICIT NONE
1670 
1671  REAL*8, INTENT(IN) :: qpj(n_vars+2)
1672  REAL*8, INTENT(IN) :: dt
1673 
1674  REAL*8, INTENT(OUT) :: erosion_term(n_solid)
1675  REAL*8, INTENT(OUT) :: deposition_term(n_solid)
1676 
1677  REAL*8 :: mod_vel
1678 
1679  REAL*8 :: hind_exp
1680  REAL*8 :: alpha_max
1681 
1682  INTEGER :: i_solid
1683 
1684  REAL*8 :: r_h
1685  REAL*8 :: r_u
1686  REAL*8 :: r_v
1687  REAL*8 :: r_alphas(n_solid)
1688  REAL*8 :: r_Ri
1689  REAL*8 :: r_rho_m
1690  REAL*8 :: r_rho_c
1691  REAL*8 :: r_red_grav
1692 
1693 
1694  ! parameters for Michaels and Bolger (1962) sedimentation correction
1695  alpha_max = 0.6d0
1696  hind_exp = 4.65d0
1697 
1698  deposition_term(1:n_solid) = 0.d0
1699  erosion_term(1:n_solid) = 0.d0
1700 
1701  IF ( qpj(1) .LE. 0.d0 ) RETURN
1702 
1703  r_h = qpj(1)
1704  r_u = qpj(n_vars+1)
1705  r_v = qpj(n_vars+2)
1706  r_alphas(1:n_solid) = qpj(5:4+n_solid)
1707 
1708  CALL mixt_var(qpj,r_ri,r_rho_m,r_rho_c,r_red_grav)
1709 
1710  DO i_solid=1,n_solid
1711 
1712  IF ( ( r_alphas(i_solid) .GT. 0.d0 ) .AND. ( settling_flag ) ) THEN
1713 
1714  settling_vel = settling_velocity( diam_s(i_solid) , rho_s(i_solid) , &
1715  r_rho_c , i_solid )
1716 
1717  deposition_term(i_solid) = r_alphas(i_solid) * settling_vel * &
1718  ( 1.d0 - min( 1.d0 , sum( r_alphas ) / alpha_max ) )**hind_exp
1719 
1720  deposition_term(i_solid) = min( deposition_term(i_solid) , &
1721  r_h * r_alphas(i_solid) / dt )
1722 
1723  IF ( deposition_term(i_solid) .LT. 0.d0 ) THEN
1724 
1725  WRITE(*,*) 'eval_erosion_dep_term'
1726  WRITE(*,*) 'deposition_term(i_solid)',deposition_term(i_solid)
1727  READ(*,*)
1728 
1729  END IF
1730 
1731  END IF
1732 
1733  mod_vel = dsqrt( r_u**2 + r_v**2 )
1734 
1735  IF ( r_h .GT. 1.d-2) THEN
1736 
1737  ! empirical formulation (see Fagents & Baloga 2006, Eq. 5)
1738  ! here we use the solid volume fraction instead of relative density
1739  ! This term has units: m s-1
1740  erosion_term(i_solid) = erosion_coeff(i_solid) * mod_vel * r_h &
1741  * ( 1.d0 - sum(r_alphas) )
1742 
1743  ELSE
1744 
1745  erosion_term(i_solid) = 0.d0
1746 
1747  END IF
1748 
1749  END DO
1750 
1751  RETURN
1752 
1753  END SUBROUTINE eval_erosion_dep_term
1754 
1755 
1756  !******************************************************************************
1758  !
1766  !
1769  !
1770  !******************************************************************************
1771 
1772  SUBROUTINE eval_topo_term( qpj , deposition_avg_term , erosion_avg_term , &
1773  eqns_term, deposit_term )
1774 
1775  IMPLICIT NONE
1776 
1777  REAL*8, INTENT(IN) :: qpj(n_vars+2)
1778  REAL*8, INTENT(IN) :: deposition_avg_term(n_solid)
1779  REAL*8, INTENT(IN) :: erosion_avg_term(n_solid)
1780 
1781  REAL*8, INTENT(OUT):: eqns_term(n_eqns)
1782  REAL*8, INTENT(OUT):: deposit_term(n_solid)
1783 
1784  REAL*8 :: entr_coeff
1785  REAL*8 :: air_entr
1786  REAL*8 :: mag_vel
1787 
1788  REAL*8 :: r_h
1789  REAL*8 :: r_u
1790  REAL*8 :: r_v
1791  REAL*8 :: r_T
1792  REAL*8 :: r_Ri
1793  REAL*8 :: r_rho_m
1794  REAL*8 :: r_rho_c
1795  REAL*8 :: r_red_grav
1796 
1797 
1798  IF ( qpj(1) .LE. 0.d0 ) THEN
1799 
1800  eqns_term(1:n_eqns) = 0.d0
1801  deposit_term(1:n_solid) = 0.d0
1802 
1803  RETURN
1804 
1805  END IF
1806 
1807  CALL mixt_var(qpj,r_ri,r_rho_m,r_rho_c,r_red_grav)
1808 
1809  r_h = qpj(1)
1810  r_u = qpj(n_vars+1)
1811  r_v = qpj(n_vars+2)
1812  r_t = qpj(4)
1813 
1814  mag_vel = dsqrt( r_u**2.d0 + r_v**2.d0 )
1815 
1816  IF ( entrainment_flag .AND. ( mag_vel**2 .GT. 0.d0 ) .AND. &
1817  ( r_h .GT. 0.d0 ) ) THEN
1818 
1819  entr_coeff = 0.075d0 / dsqrt( 1.d0 + 718.d0 * max(0.d0,r_ri)**2.4 )
1820 
1821  air_entr = entr_coeff * mag_vel
1822 
1823  ELSE
1824 
1825  air_entr = 0.d0
1826 
1827  END IF
1828 
1829  eqns_term(1:n_eqns) = 0.d0
1830 
1831  ! free surface (topography+flow) equation
1832  eqns_term(1) = sum( rho_s * ( erosion_avg_term - deposition_avg_term ) ) + &
1833  rho_a_amb * air_entr
1834 
1835  ! x-momenutm equation
1836  eqns_term(2) = - r_u * sum( rho_s * deposition_avg_term )
1837 
1838  ! y-momentum equation
1839  eqns_term(3) = - r_v * sum( rho_s * deposition_avg_term )
1840 
1841  ! Temperature/Energy equation
1842  IF ( energy_flag ) THEN
1843 
1844  eqns_term(4) = - r_t * sum( rho_s * sp_heat_s * deposition_avg_term ) &
1845  - 0.5d0 * mag_vel**2 * sum( rho_s * deposition_avg_term ) &
1846  + t_s_substrate * sum( rho_s * sp_heat_s * erosion_avg_term ) &
1847  + t_ambient * sp_heat_a * rho_a_amb * air_entr
1848 
1849  ELSE
1850 
1851  eqns_term(4) = - r_t * sum( rho_s * sp_heat_s * deposition_avg_term ) &
1852  + t_s_substrate * sum( rho_s * sp_heat_s * erosion_avg_term ) &
1853  + t_ambient * sp_heat_a * rho_a_amb * air_entr
1854 
1855  END IF
1856 
1857  ! solid phase thickness equation
1858  eqns_term(5:4+n_solid) = rho_s(1:n_solid) * ( erosion_avg_term(1:n_solid) &
1859  - deposition_avg_term(1:n_solid) )
1860 
1861  ! solid deposit rate terms
1862  deposit_term(1:n_solid) = deposition_avg_term(1:n_solid) &
1863  - erosion_avg_term(1:n_solid)
1864 
1865  RETURN
1866 
1867  END SUBROUTINE eval_topo_term
1868 
1869  !******************************************************************************
1871  !
1881  !
1884  !
1885  !******************************************************************************
1886 
1887  SUBROUTINE eval_source_bdry( time, vect_x , vect_y , source_bdry )
1888 
1889  USE parameters_2d, ONLY : h_source , vel_source , t_source , alphas_source ,&
1890  alphal_source , time_param
1891 
1892  IMPLICIT NONE
1893 
1894  REAL*8, INTENT(IN) :: time
1895  REAL*8, INTENT(IN) :: vect_x
1896  REAL*8, INTENT(IN) :: vect_y
1897  REAL*8, INTENT(OUT) :: source_bdry(n_vars)
1898 
1899  REAL*8 :: t_rem
1900  REAL*8 :: t_coeff
1901  REAL*8 :: pi_g
1902 
1903  IF ( time .GE. time_param(4) ) THEN
1904 
1905  ! The exponents of t_coeff are such that Ri does not depend on t_coeff
1906  source_bdry(1) = 0.d0
1907  source_bdry(2) = 0.d0
1908  source_bdry(3) = 0.d0
1909  source_bdry(4) = t_source
1910  source_bdry(5:4+n_solid) = alphas_source(1:n_solid)
1911 
1912  IF ( gas_flag .AND. liquid_flag ) THEN
1913 
1914  source_bdry(n_vars) = alphal_source
1915 
1916  END IF
1917 
1918  source_bdry(n_vars+1) = 0.d0
1919  source_bdry(n_vars+2) = 0.d0
1920 
1921  RETURN
1922 
1923  END IF
1924 
1925  t_rem = dmod( time , time_param(1) )
1926 
1927  pi_g = 4.d0 * datan(1.d0)
1928 
1929  t_coeff = 0.d0
1930 
1931  IF ( time_param(3) .EQ. 0.d0 ) THEN
1932 
1933  IF ( t_rem .LE. time_param(2) ) t_coeff = 1.d0
1934 
1935  ELSE
1936 
1937  IF ( t_rem .LT. time_param(3) ) THEN
1938 
1939  t_coeff = 0.5d0 * ( 1.d0 - dcos( pi_g * t_rem / time_param(3) ) )
1940 
1941  ELSEIF ( t_rem .LE. ( time_param(2) - time_param(3) ) ) THEN
1942 
1943  t_coeff = 1.d0
1944 
1945  ELSEIF ( t_rem .LE. time_param(2) ) THEN
1946 
1947  t_coeff = 0.5d0 * ( 1.d0 + dcos( pi_g * ( ( t_rem - time_param(2) ) &
1948  / time_param(3) + 1.d0 ) ) )
1949 
1950  END IF
1951 
1952  END IF
1953 
1954  ! The exponents of t_coeff are such that Ri does not depend on t_coeff
1955  source_bdry(1) = t_coeff * h_source
1956  source_bdry(2) = t_coeff**1.5d0 * h_source * vel_source * vect_x
1957  source_bdry(3) = t_coeff**1.5d0 * h_source * vel_source * vect_y
1958  source_bdry(4) = t_source
1959  source_bdry(5:4+n_solid) = alphas_source(1:n_solid)
1960 
1961  IF ( gas_flag .AND. liquid_flag ) THEN
1962 
1963  source_bdry(n_vars) = alphal_source
1964 
1965  END IF
1966 
1967  source_bdry(n_vars+1) = t_coeff**0.5d0 * vel_source * vect_x
1968  source_bdry(n_vars+2) = t_coeff**0.5d0 * vel_source * vect_y
1969 
1970  RETURN
1971 
1972  END SUBROUTINE eval_source_bdry
1973 
1974  !------------------------------------------------------------------------------
1976  !
1984  !
1987  !
1988  !------------------------------------------------------------------------------
1989 
1990  REAL*8 FUNCTION settling_velocity(diam,rhos,rhoc,i_solid)
1991 
1992  IMPLICIT NONE
1993 
1994  REAL*8, INTENT(IN) :: diam
1995  REAL*8, INTENT(IN) :: rhos
1996  REAL*8, INTENT(IN) :: rhoc
1997  INTEGER, INTENT(IN) :: i_solid
1998 
1999  REAL*8 :: Rey
2000  REAL*8 :: C_D
2001 
2002  ! loop variables
2003  INTEGER :: i
2004  REAL*8 :: const_part
2005  REAL*8 :: C_D_old
2006  REAL*8 :: set_vel_old
2007 
2008  INTEGER :: dig
2009 
2010  c_d = 1.d0
2011 
2012  const_part = dsqrt( 0.75d0 * ( rhos / rhoc - 1.d0 ) * diam * grav )
2013 
2014  settling_velocity = const_part / dsqrt( c_d )
2015 
2016  rey = diam * settling_velocity / kin_visc_c
2017 
2018  IF ( rey .LE. 1000.d0 ) THEN
2019 
2020  c_d_loop:DO i=1,20
2021 
2022  set_vel_old = settling_velocity
2023  c_d_old = c_d
2024  c_d = 24.d0 / rey * ( 1.d0 + 0.15d0 * rey**(0.687) )
2025 
2026  settling_velocity = const_part / dsqrt( c_d )
2027 
2028  IF ( dabs( set_vel_old - settling_velocity ) / set_vel_old &
2029  .LT. 1.d-6 ) THEN
2030 
2031  ! round to first three significative digits
2032  dig = floor(log10(set_vel_old))
2033  settling_velocity = 10.d0**(dig-3) &
2034  * floor( 10.0**(-dig+3)*set_vel_old )
2035 
2036  EXIT c_d_loop
2037 
2038  END IF
2039 
2040  rey = diam * settling_velocity / kin_visc_c
2041 
2042  END DO c_d_loop
2043 
2044  END IF
2045 
2046  RETURN
2047 
2048  END FUNCTION settling_velocity
2049 
2050 END MODULE constitutive_2d
2051 
2052 
constitutive_2d::kin_visc_a
real *8 kin_visc_a
Kinematic viscosity of air (units: m2 s-1)
Definition: constitutive_2d.f90:112
constitutive_2d::emissivity
real *8 emissivity
emissivity (eps in Costa & Macedonio, 2005)
Definition: constitutive_2d.f90:66
constitutive_2d::r_phys_var
subroutine r_phys_var(r_qj, r_h, r_u, r_v, r_alphas, r_rho_m, r_T, r_alphal)
Physical variables.
Definition: constitutive_2d.f90:206
constitutive_2d::sp_heat_l
real *8 sp_heat_l
Sepcific heat of liquid (units: J K-1 kg-1)
Definition: constitutive_2d.f90:151
constitutive_2d::nu_ref
real *8 nu_ref
reference kinematic viscosity [m2/s]
Definition: constitutive_2d.f90:45
constitutive_2d::t_ground
real *8 t_ground
temperature of lava-ground interface
Definition: constitutive_2d.f90:72
constitutive_2d::eval_local_speeds_x
subroutine eval_local_speeds_x(qpj, vel_min, vel_max)
Local Characteristic speeds x direction.
Definition: constitutive_2d.f90:973
parameters_2d::vel_source
real *8 vel_source
Definition: parameters_2d.f90:78
constitutive_2d::beta2
real *8 beta2
2nd param for yield strenght empirical relationship (O'Brian et al, 1993)
Definition: constitutive_2d.f90:83
constitutive_2d::sbconst
real *8, parameter sbconst
Stephan-Boltzmann constant [W m-2 K-4].
Definition: constitutive_2d.f90:63
parameters_2d::rheology_flag
logical rheology_flag
Flag to choose if we add the rheology.
Definition: parameters_2d.f90:59
constitutive_2d::mixt_var
subroutine mixt_var(qpj, r_Ri, r_rho_m, r_rho_c, r_red_grav)
Physical variables.
Definition: constitutive_2d.f90:569
constitutive_2d::entrainment_flag
logical entrainment_flag
flag to activate air entrainment
Definition: constitutive_2d.f90:17
constitutive_2d::rad_coeff
real *8 rad_coeff
radiative coefficient
Definition: constitutive_2d.f90:36
parameters_2d::liquid_flag
logical liquid_flag
Definition: parameters_2d.f90:67
constitutive_2d::c_p
real *8 c_p
specific heat [J kg-1 K-1]
Definition: constitutive_2d.f90:54
constitutive_2d::c_phys_var
subroutine c_phys_var(c_qj, h, u, v, T, rho_m, red_grav, alphas)
Physical variables.
Definition: constitutive_2d.f90:388
parameters_2d
Parameters.
Definition: parameters_2d.f90:7
constitutive_2d::qp_to_qp2
subroutine qp_to_qp2(qpj, Bj, qp2j)
Additional Physical variables.
Definition: constitutive_2d.f90:934
constitutive_2d::xi
real *8 xi
Definition: constitutive_2d.f90:24
constitutive_2d::kin_visc_l
real *8 kin_visc_l
Kinematic viscosity of liquid (units: m2 s-1)
Definition: constitutive_2d.f90:115
parameters_2d::rheology_model
integer rheology_model
choice of the rheology model
Definition: parameters_2d.f90:65
constitutive_2d::sp_gas_const_a
real *8 sp_gas_const_a
Specific gas constant of air (units: J kg-1 K-1)
Definition: constitutive_2d.f90:109
constitutive_2d::thermal_conductivity
real *8 thermal_conductivity
thermal conductivity [W m-1 K-1] (k in Costa & Macedonio, 2005)
Definition: constitutive_2d.f90:75
constitutive_2d::t_s_substrate
real *8 t_s_substrate
temperature of solid substrate (units: K)
Definition: constitutive_2d.f90:142
constitutive_2d::sp_heat_a
real *8 sp_heat_a
Specific heat of air (units: J K-1 kg-1)
Definition: constitutive_2d.f90:106
constitutive_2d::eval_nh_semi_impl_terms
subroutine eval_nh_semi_impl_terms(grav3_surf, qcj, nh_semi_impl_term)
Non-Hyperbolic semi-implicit terms.
Definition: constitutive_2d.f90:1478
constitutive_2d::beta1
real *8 beta1
2nd param for fluid viscosity empirical relationship (O'Brian et al, 1993)
Definition: constitutive_2d.f90:89
constitutive_2d::kin_visc_c
real *8 kin_visc_c
Kinematic viscosity of carrier phase (units: m2 s-1)
Definition: constitutive_2d.f90:118
parameters_2d::time_param
real *8, dimension(4) time_param
Definition: parameters_2d.f90:83
parameters_2d::n_vars
integer n_vars
Number of conservative variables.
Definition: parameters_2d.f90:124
parameters_2d::alphas_source
real *8, dimension(:), allocatable alphas_source
Definition: parameters_2d.f90:81
constitutive_2d::sp_heat_c
real *8 sp_heat_c
Specific heat of carrier phase (gas or liquid)
Definition: constitutive_2d.f90:100
constitutive_2d::eval_expl_terms
subroutine eval_expl_terms(Bprimej_x, Bprimej_y, source_xy, qpj, qcj, expl_term)
Explicit Forces term.
Definition: constitutive_2d.f90:1604
constitutive_2d::enne
real *8 enne
thermal boundary layer fraction of total thickness
Definition: constitutive_2d.f90:69
constitutive_2d::eval_fluxes
subroutine eval_fluxes(qcj, qpj, dir, flux)
Hyperbolic Fluxes.
Definition: constitutive_2d.f90:1077
constitutive_2d
Constitutive equations.
Definition: constitutive_2d.f90:4
constitutive_2d::grav
real *8 grav
gravitational acceleration
Definition: constitutive_2d.f90:20
constitutive_2d::mu
real *8 mu
drag coefficients (Voellmy-Salm model)
Definition: constitutive_2d.f90:23
constitutive_2d::visc_par
real *8 visc_par
viscosity parameter [K-1] (b in Table 1 Costa & Macedonio, 2005)
Definition: constitutive_2d.f90:48
constitutive_2d::alpha1_coeff
real *8 alpha1_coeff
ratio between reference value from input and computed values from eq.
Definition: constitutive_2d.f90:86
parameters_2d::alphal_source
real *8 alphal_source
Definition: parameters_2d.f90:82
constitutive_2d::implicit_flag
logical, dimension(:), allocatable implicit_flag
flag used for size of implicit non linear-system
Definition: constitutive_2d.f90:14
constitutive_2d::diam_s
real *8, dimension(:), allocatable diam_s
Diameter of sediments ( units: m )
Definition: constitutive_2d.f90:127
constitutive_2d::n_td
real *8 n_td
Mannings roughness coefficient ( units: T L^(-1/3) )
Definition: constitutive_2d.f90:95
constitutive_2d::erosion_coeff
real *8, dimension(:), allocatable erosion_coeff
erosion model coefficient (units: m-1 )
Definition: constitutive_2d.f90:139
constitutive_2d::settling_flag
logical settling_flag
Flag to determine if sedimentation is active.
Definition: constitutive_2d.f90:133
constitutive_2d::eval_source_bdry
subroutine eval_source_bdry(time, vect_x, vect_y, source_bdry)
Internal boundary source fluxes.
Definition: constitutive_2d.f90:1888
constitutive_2d::rho_l
real *8 rho_l
liquid density (units: kg m-3)
Definition: constitutive_2d.f90:148
constitutive_2d::eval_local_speeds_y
subroutine eval_local_speeds_y(qpj, vel_min, vel_max)
Local Characteristic speeds y direction.
Definition: constitutive_2d.f90:1024
constitutive_2d::init_problem_param
subroutine init_problem_param
Initialization of relaxation flags.
Definition: constitutive_2d.f90:164
constitutive_2d::t_env
real *8 t_env
evironment temperature [K]
Definition: constitutive_2d.f90:33
constitutive_2d::alpha2
real *8 alpha2
1st param for yield strenght empirical relationship (O'Brian et al, 1993)
Definition: constitutive_2d.f90:80
constitutive_2d::eval_erosion_dep_term
subroutine eval_erosion_dep_term(qpj, dt, erosion_term, deposition_term)
Erosion/Deposition term.
Definition: constitutive_2d.f90:1668
parameters_2d::gas_flag
logical gas_flag
Definition: parameters_2d.f90:69
constitutive_2d::atm_heat_transf_coeff
real *8 atm_heat_transf_coeff
atmospheric heat trasnfer coefficient [W m-2 K-1] (lambda in C&M, 2005)
Definition: constitutive_2d.f90:57
constitutive_2d::settling_vel
real *8 settling_vel
Hindered settling velocity (units: m s-1 )
Definition: constitutive_2d.f90:136
constitutive_2d::rho_a_amb
real *8 rho_a_amb
Ambient density of air ( units: kg m-3 )
Definition: constitutive_2d.f90:103
constitutive_2d::pres
real *8 pres
ambient pressure (units: Pa)
Definition: constitutive_2d.f90:145
constitutive_2d::t_ref
real *8 t_ref
reference temperature [K]
Definition: constitutive_2d.f90:42
constitutive_2d::qc_to_qp
subroutine qc_to_qp(qc, qp)
Conservative to physical variables.
Definition: constitutive_2d.f90:678
parameters_2d::energy_flag
logical energy_flag
Flag to choose the equation for temperature to solve.
Definition: parameters_2d.f90:47
constitutive_2d::exp_area_fract
real *8 exp_area_fract
fractional area of the exposed inner core (f in C&M, 2005)
Definition: constitutive_2d.f90:60
constitutive_2d::sp_heat_s
real *8, dimension(:), allocatable sp_heat_s
Specific heat of solids (units: J K-1 kg-1)
Definition: constitutive_2d.f90:130
parameters_2d::n_eqns
integer n_eqns
Number of equations.
Definition: parameters_2d.f90:125
constitutive_2d::qp_to_qc
subroutine qp_to_qc(qp, B, qc)
Physical to conservative variables.
Definition: constitutive_2d.f90:752
constitutive_2d::t_ambient
real *8 t_ambient
Temperature of ambient air (units: K)
Definition: constitutive_2d.f90:121
constitutive_2d::eval_nonhyperbolic_terms
subroutine eval_nonhyperbolic_terms(c_qj, c_nh_term_impl, r_qj, r_nh_term_impl)
Non-Hyperbolic terms.
Definition: constitutive_2d.f90:1208
parameters_2d::h_source
real *8 h_source
Definition: parameters_2d.f90:80
constitutive_2d::emme
real *8 emme
velocity boundary layer fraction of total thickness
Definition: constitutive_2d.f90:51
parameters_2d::eps_sing
real *8 eps_sing
parameter for desingularization
Definition: parameters_2d.f90:20
constitutive_2d::tau
real *8 tau
drag coefficients (plastic model)
Definition: constitutive_2d.f90:30
constitutive_2d::friction_factor
real *8 friction_factor
drag coefficients (B&W model)
Definition: constitutive_2d.f90:27
constitutive_2d::frict_coeff
real *8 frict_coeff
friction coefficient
Definition: constitutive_2d.f90:39
constitutive_2d::kappa
real *8 kappa
Empirical resistance parameter (dimensionless, input parameter)
Definition: constitutive_2d.f90:92
parameters_2d::n_solid
integer n_solid
Number of solid classes.
Definition: parameters_2d.f90:127
constitutive_2d::rho_s
real *8, dimension(:), allocatable rho_s
Density of sediments ( units: kg m-3 )
Definition: constitutive_2d.f90:124
constitutive_2d::eval_topo_term
subroutine eval_topo_term(qpj, deposition_avg_term, erosion_avg_term, eqns_term, deposit_term)
Topography modification related term.
Definition: constitutive_2d.f90:1774
parameters_2d::n_nh
integer n_nh
Number of non-hyperbolic terms.
Definition: parameters_2d.f90:129
constitutive_2d::settling_velocity
real *8 function settling_velocity(diam, rhos, rhoc, i_solid)
Settling velocity function.
Definition: constitutive_2d.f90:1991
parameters_2d::sed_vol_perc
real *8, dimension(:), allocatable sed_vol_perc
Definition: parameters_2d.f90:122
parameters_2d::t_source
real *8 t_source
Definition: parameters_2d.f90:79
Generated on Fri Dec 13 2019 21:38:13 for SW_VAR_DENS_MODEL by   doxygen 1.8.11