30
31
32
33 USE elbufdef_mod
34
35
36
37#include "implicit_f.inc"
38
39
40
41
42#include "param_c.inc"
43
44#include "com01_c.inc"
45
46
47
48 INTEGER, INTENT(IN) :: NEL, NG, MATID(NPROPMI, *), LIST(),NELG
49 my_real,
INTENT(IN) :: grav0, depth(*), pm(npropm, *), bufmat(*)
51 TYPE(ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET, INTENT(IN) :: ELBUF_TAB
52
53
54
55 INTEGER :: I,ISOLVER, K
56 my_real :: r1, c1, p0, pgrav, rho10, rho20, rho1, rho2, gam, rho0,
58 TYPE(G_BUFEL_), POINTER :: GBUF
59 TYPE(BUF_MAT_) ,POINTER :: MBUF
60
61
62
63
64
65
66
67
68
69 gbuf => elbuf_tab(ng)%GBUF
70
71 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
72
73 p0 = bufmat(9)
74 psh= bufmat(16)
75
76 rho10 = bufmat(11)
77 c1 = bufmat(4)
78 r1 = bufmat(6)
79
80 gam = bufmat(5)
81 rho20 = bufmat(12)
82 isolver = bufmat(17)
83
84 IF(psurf==zero .AND. isolver<=1)THEN
85 psurf=p0
86 print *, "**WARNING : INIGRAV CARD, PREF PARAMETER MUST BE A TOTAL PRESSURE WITH LAW37, SETTING PREF=P0"
87 ENDIF
88
89 DO k = 1, nel
90 i = list(k)
91 alpha1 = mbuf%VAR(i + (4 - 1) * nelg)
93 rho0 = alpha1 * rho10 +
alpha2 * rho20
94 pgrav = psurf - rho0 * grav0 * depth(k)
95 rho1 = (pgrav-p0)/r1 + rho10
96 rho2 = rho20 * (pgrav/p0) ** (one / gam)
97 gbuf%RHO(i) = alpha1 * rho1 +
alpha2 * rho2
98 mbuf%VAR(i + (4 - 1) * nelg) = alpha1
99 mbuf%VAR(i + (5 - 1) * nelg) = one - alpha1
100 mbuf%VAR(i + (2 - 1) * nelg) = rho2
101 mbuf%VAR(i + (3 - 1) * nelg) = rho1
102 mbuf%VAR(i + (1 - 1) * nelg) = alpha1 * rho1
103 gbuf%SIG(i) = - (pgrav-p0-psh)
104 gbuf%SIG(i + nelg) = - (pgrav-p0-psh)
105 gbuf%SIG(i + 2 * nelg) = - (pgrav-p0-psh)
106 ENDDO
107