33
34
35
37
38
39
40#include "implicit_f.inc"
41
42
43
44#include "com01_c.inc"
45#include "com04_c.inc"
46#include "com08_c.inc"
47#include "scr11_c.inc"
48#include "statr_c.inc"
49#include "stati_c.inc"
50#include "task_c.inc"
51#include "sms_c.inc"
52
53
54
55 INTEGER WEIGHT_MD(*)
56 my_real v(3,*), vr(3,*), a(3,*), ar(3,*),ms(*),in(*)
57 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
58
59 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
60
61
62
63 INTEGER J, I,LAST,N,NGR2USR,,ISTAT3
64 my_real encino, omega, uomega, domega, omega2, encint, encinn,encgrp,encgrpn
66
67 DATA encino /0.0/
68 DATA last /0/
69 SAVE encino,last
70
71 IF(istatg<0)THEN
72 istatg=
ngr2usr(-istatg,igrnod,ngrnod)
73 ENDIF
74 istat2=0
75 istat3=0
76 IF (ncycle==0) THEN
77 IF ((istat==2.OR.istat==3).AND.tst_stop==zero) tst_stop=tstop
78 END IF
79 IF (istat==2.AND.tt>=tst_start.AND.tt<=tst_stop) istat2=1
80 IF (istat==3.AND.tt>=tst_start.AND.tt<=tst_stop) istat3=1
81 IF((istat==1.OR.istat3==1).AND.istatg==0)THEN
82
83 omega = betate * dt12
84 uomega = one - omega
85 domega = two*betate
86 omega2 = (one-two*omega)**2
87 encint = encin
88
89 encinn = encint * omega2
90
91 IF(idtmins==0) THEN
92 DO i = 1, numnod
93 a(1,i) = -domega*v(1,i) + uomega*a(1,i)
94 a(2,i) = -domega*v(2,i) + uomega*a(2,i)
95 a(3,i) = -domega*v(3,i) + uomega*a(3,i)
96 ENDDO
97 IF (ispmd==0) wfext = wfext - encint + encinn
98 END IF
99
100 IF(iroddl/=0)THEN
101 encint = enrot
102 encinn = encint * omega2
103 DO i = 1, numnod
104 ar(1,i) = -domega*vr(1,i) + uomega*ar(1,i)
105 ar(2,i) = -domega*vr(2,i) + uomega*ar(2,i)
106 ar(3,i) = -domega*vr(3,i) + uomega*ar(3,i)
107 ENDDO
108 IF (ispmd==0) wfext = wfext - encint + encinn
109 ENDIF
110
111 ELSEIF(istat2==1.AND.istatg==0)THEN
112
113 last = last+1
114 encint = encin+enrot
115 IF(encint<encino)THEN
116 IF (ispmd==0)wfext = wfext - encint
117 encino=zero
118 last=0
119
120 DO j=1,3
121 DO i = 1, numnod
122 v(j,i) = zero
123 ENDDO
124 IF(iroddl/=0)THEN
125 DO i = 1, numnod
126 vr(j,i) = zero
127 ENDDO
128 ENDIF
129 ENDDO
130 ELSE
131 encino = encint
132 ENDIF
133 ELSEIF((istat==1.OR.istat3==1).AND.istatg/=0)THEN
134
135 omega = betate * dt12
136 uomega = one - omega
137 domega = two*betate
138 omega2 = (one-two*omega)**2
139 encint = encin+enrot
140 encinn = encint * omega2
141 encgrp = zero
142 encgrpr = zero
143
144 IF(idtmins==0) THEN
145 DO j=1,3
146#include "vectorize.inc"
147 DO n=1,igrnod(istatg)%NENTITY
148 i=igrnod(istatg)%ENTITY(n)
149 encgrp = encgrp + ms(i)*weight_md(i)*
150 . v(j,i)*v(j,i)
151 a(j,i) = -domega*v(j,i) + uomega*a(j,i)
152 ENDDO
153 ENDDO
154 END IF
155
156 IF(iroddl/=0)THEN
157 DO j=1,3
158#include "vectorize.inc"
159 DO n=1,igrnod(istatg)%NENTITY
160 i=igrnod(istatg)%ENTITY(n)
161 encgrpr = encgrpr + ms(i)*weight_md(i)*
162 . vr(j,i)*vr(j,i)
163 ar(j,i) = -domega*vr(j,i) + uomega*ar(j,i)
164 ENDDO
165 ENDDO
166 ENDIF
167
168 encgrp = half*encgrp
169 encgrpr = half*encgrpr
170 encgrpn = encgrp * omega2
171 encgrprn = encgrpr * omega2
172 IF (ispmd==0)wfext = wfext + encgrpn - encgrp + encgrprn - encgrpr
173
174 ELSEIF(istat2==1.AND.istatg/=0)THEN
175
176 last = last+1
177 encint = encin+enrot
178 IF(encint<encino)THEN
179 IF (ispmd==0)wfext = wfext - encint
180 encino=zero
181 last=0
182
183 DO j=1,3
184#include "vectorize.inc"
185 DO n=1,igrnod(istatg)%NENTITY
186 i=igrnod(istatg)%ENTITY(n)
187 v(j,i) = zero
188 ENDDO
189 IF(iroddl/=0)THEN
190#include "vectorize.inc"
191 DO n=1,igrnod(istatg)%NENTITY
192 i=igrnod(istatg)%ENTITY(n)
193 vr(j,i) = zero
194 ENDDO
195 ENDIF
196 ENDDO
197 ELSE
198 encino = encint
199 ENDIF
200 ENDIF
201
202 RETURN
integer function ngr2usr(iu, igr, ngr)