32
33
34
35
36
37
38
39#include "implicit_f.inc"
40
41
42
43#include "mvsiz_p.inc"
44
45
46
47#include "vect01_c.inc"
48#include "param_c.inc"
49
50
51
52 my_real,
INTENT(IN) :: pm(npropm, *), geo(npropg, *), aire(*), vol(*)
53 my_real,
INTENT(INOUT) :: dtx(*)
54 my_real,
INTENT(IN),
DIMENSION(:),
TARGET :: bufmat(*), deltax(*)
55 INTEGER, INTENT(IN) :: PID(*),MAT(*),IPM(NPROPMI, *)
56
57
58
59 INTEGER :: I, IADBUF,IFLG(MVSIZ), BIJ(4)
60 my_real,
DIMENSION(:),
POINTER :: uparam
61 my_real :: ssp(mvsiz) , rho0(mvsiz) , bulk(mvsiz) ,
62 . av1(mvsiz) , av2(mvsiz) , av3(mvsiz) , av4(mvsiz) ,
63 . rho10(mvsiz), rho20(mvsiz), rho30(mvsiz), rho40(mvsiz),
64 . ssp1(mvsiz) , ssp2(mvsiz) , ssp3(mvsiz) , ssp4(mvsiz),
65 . vis(mvsiz) , rho0_bak(mvsiz)
66 INTEGER :: IAV(4),IRHO(4),ISSP(4),NUPARAM
67
68
69 iav(1:4) = (/ 4,5,6,46 /)
70 irho(1:4) = (/ 9,10,11,47 /)
71 issp(1:4) = (/ 124,174,224,273 /)
72
73
74
75 iflg(1) = 0
76 DO i=1,llt
77 iadbuf = ipm(7,mat(i))
78 nuparam = ipm(9,mat(i))
79 uparam =>bufmat(iadbuf:iadbuf+nuparam-1)
80 bij(1:4) = nint(uparam(277:280))
81 av1(i) = uparam(iav(bij(1)))
82 av2(i) = uparam(iav(bij(2)))
83 av3(i) = uparam(iav(bij(3)))
84 av4(i) = uparam(iav(bij(4)))
85 rho10(i) = uparam(irho(bij(1)))
86 rho20(i) = uparam(irho(bij(2)))
87 rho30(i) = uparam(irho(bij(3)))
88 rho40(i) = uparam(irho(bij(4)))
89 iflg(i) = nint(uparam(31))
90 vis(i) = zero
91 ssp1(i) = uparam(issp(bij(1)))
92 ssp2(i) = uparam(issp(bij(2)))
93 ssp3(i) = uparam(issp(bij(3)))
94 ssp4(i) = uparam(issp(bij(4)))
95 enddo
96
97 IF (iflg(1) == 6) THEN
98 DO i=lft,llt
99 dtx(i)=ep20
100 ENDDO
101 RETURN
102 ENDIF
103
104 DO i=lft,llt
105 rho0(i) = av1(i)*rho10(i) + av2(i)*rho20(i) + av3(i)*rho30(i) + av4(i)*rho40(i)
106 ENDDO
107
108 DO i=lft,llt
109 bulk(i) = av1(i)*rho10(i)*ssp1(i)**2
110 . + av2(i)*rho20(i)*ssp2(i)**2
111 . + av3(i)*rho30(i)*ssp3(i)**2
112 . + av4(i)*rho40(i)*ssp4(i)**2
113 IF (rho0(i) > zero) THEN
114 ssp(i) = sqrt(bulk(i)/rho0(i))
115 rho0_bak(i)=rho0(i)
116 ELSE
117 ssp(i) = em20
118 rho0_bak(i)=em20
119 ENDIF
120 enddo
121
122
123
124
125 IF(jsph==0)THEN
126 CALL dtel(ssp,pm,geo,pid,mat, rho0_bak, vis, deltax, aire, vol, dtx)
127 ELSE
128 CALL dtsph(ssp,pm,geo,pid,mat, rho0_bak, vis, deltax, vol, dtx)
129 ENDIF
130
131 RETURN
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
subroutine dtsph(ssp, pm, geo, pid, mat, rho0, vis, deltax, vol, dtx)