FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
heat_LIB_DFLUX.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
7 module m_heat_lib_dflux
8 contains
9  !C************************************************************************
10  !C* DFLUX_111(NN,XX,YY,ZZ,ASECT,LTYPE,VAL,VECT)
11  !C* DFLUX_231(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
12  !C* DFLUX_232(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
13  !C* DFLUX_241(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
14  !C* DFLUX_242(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
15  !C* DFLUX_341(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
16  !C* DFLUX_342(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
17  !C* DFLUX_351(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
18  !C* DFLUX_352(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
19  !C* DFLUX_361(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
20  !C* DFLUX_362(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
21  !C* DFLUX_731(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
22  !C* DFLUX_741(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
23  !C*--------------------------------------------------------------------*
24  subroutine heat_dflux_111(NN,XX,YY,ZZ,ASECT,LTYPE,val,VECT)
25  !C*--------------------------------------------------------------------*
26  !C***
27  !C*** SET 111 DFLUX
28  !C***
29  !C* LTYPE=0 : BODY FLUX
30  use hecmw
31  implicit none
32  !C* I/F VARIABLES
33  integer(kind=kint) nn, ltype
34  real(kind=kreal) xx(nn),yy(nn),zz(nn),asect,val,vect(nn)
35  !C* LOCAL VARIABLES
36  integer(kind=kint) i
37  real(kind=kreal) dx, dy, dz, al, vv
38  !C*
39  !C*
40  if( ltype.EQ.0 ) then
41  asect = 1.0d0 * asect
42  else
43  asect = 0.0d0
44  endif
45 
46  do i=1,nn
47  vect(i)=0.0
48  enddo
49  !C*
50  dx = xx(2) - xx(1)
51  dy = yy(2) - yy(1)
52  dz = zz(2) - zz(1)
53  al = dsqrt( dx*dx + dy*dy + dz*dz )
54  vv = asect * al
55  !
56  vect(1)= -val*vv/2.0d0
57  vect(2)= -val*vv/2.0d0
58  !
59  return
60 
61  end subroutine heat_dflux_111
62  !C*--------------------------------------------------------------------*
63  subroutine heat_dflux_231(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
64  !C*--------------------------------------------------------------------*
65  !C***
66  !C*** SET 231 DFLUX
67  !C***
68  !C* LTYPE=0 : BODY FLUX
69  use hecmw
70  implicit none
71  !C* I/F VARIABLES
72  integer(kind=kint) :: nn, ltype
73  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
74  !C* LOCAL VARIABLES
75  integer(kind=kint) :: ivol, isuf, i, lx
76  real(kind=kreal) :: ri, gx, gy, xsum
77  real(kind=kreal) :: v1x, v1y, v1z, v2x, v2y, v2z, v3x, v3y, v3z, aa, vv
78  real(kind=kreal) :: xg(2), wgt(2), h(4), hr(4), hs(4)
79  integer(kind=kint) :: nod(2)
80  !*************************
81  ! GAUSS INTEGRATION POINT
82  !*************************
83  data xg/-0.5773502691896,0.5773502691896/
84  data wgt/1.0,1.0/
85  !C*
86  ivol=0
87  isuf=0
88  if( ltype.EQ.0 ) then
89  ivol=1
90  else if( ltype.GE.1 ) then
91  isuf=1
92  if(ltype.EQ.1) then
93  nod(1)=1
94  nod(2)=2
95  else if(ltype.EQ.2) then
96  nod(1)=2
97  nod(2)=3
98  else if(ltype.EQ.3) then
99  nod(1)=3
100  nod(2)=1
101  endif
102  endif
103  !** SURFACE LOAD
104  if( isuf.EQ.1 ) then
105  do i=1,nn
106  vect(i)=0.0
107  enddo
108  !** INTEGRATION OVER SURFACE
109  do lx=1,2
110  ri=xg(lx)
111  h(1)=0.5*(1.0-ri)
112  h(2)=0.5*(1.0+ri)
113  hr(1)=-0.5
114  hr(2)= 0.5
115  gx=0.0
116  gy=0.0
117  do i=1,2
118  gx=gx+hr(i)*xx(nod(i))
119  gy=gy+hr(i)*yy(nod(i))
120  enddo
121  xsum=dsqrt(gx*gx+gy*gy)*thick
122  do i=1,2
123  vect(nod(i))=vect(nod(i))-xsum*wgt(lx)*h(i)*val
124  enddo
125  enddo
126  endif
127  !** VOLUME LOAD
128  if( ivol.EQ.1 ) then
129  do i=1,nn
130  vect(i)=0.0
131  enddo
132 
133  v1x=xx(2)-xx(1)
134  v1y=yy(2)-yy(1)
135  v1z=zz(2)-zz(1)
136  v2x=xx(3)-xx(1)
137  v2y=yy(3)-yy(1)
138  v2z=zz(3)-zz(1)
139  v3x= v1y*v2z-v1z*v2y
140  v3y= v1z*v2x-v1x*v2z
141  v3z= v1x*v2y-v1y*v2x
142 
143  aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
144  vv=aa*thick
145  vect(1)= -val*vv/3.0
146  vect(2)= -val*vv/3.0
147  vect(3)= -val*vv/3.0
148 
149  endif
150  return
151 
152  end subroutine heat_dflux_231
153  !C*--------------------------------------------------------------------*
154  subroutine heat_dflux_232(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
155  !C*--------------------------------------------------------------------*
156  !C***
157  !C*** SET 232 DFLUX
158  !C***
159  !C* LTYPE=0 : BODY FLUX
160  use hecmw
161  implicit none
162  !C* I/F VARIABLES
163  integer(kind=kint) :: nn, ltype
164  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
165  !C* LOCAL VARIABLES
166  integer(kind=kint) :: ivol, isuf, i, l1, l2, lx
167  real(kind=kreal) :: ri, gx, gy, xsum, x1, x2, x3, xl1, xl2
168  real(kind=kreal) :: xj11, xj12, xj21, xj22, det, wg
169  real(kind=kreal) :: xg(3), wgt(3), h(6), hr(3)
170  real(kind=kreal) :: hl1(6), hl2(6), hl3(6)
171  integer(kind=kint) :: nod(3)
172  !*************************
173  ! GAUSS INTEGRATION POINT
174  !*************************
175  xg(1) = -0.7745966692
176  xg(2) = 0.0
177  xg(3) = 0.7745966692
178  wgt(1) = 0.5555555555
179  wgt(2) = 0.8888888888
180  wgt(3) = 0.5555555555
181  !C*
182  ivol=0
183  isuf=0
184  if( ltype.EQ.0 ) then
185  ivol=1
186  else if( ltype.GE.1 ) then
187  isuf=1
188  if(ltype.EQ.1) then
189  nod(1)=1
190  nod(2)=2
191  nod(3)=4
192  else if(ltype.EQ.2) then
193  nod(1)=2
194  nod(2)=3
195  nod(3)=5
196  else if(ltype.EQ.3) then
197  nod(1)=3
198  nod(2)=1
199  nod(3)=6
200  endif
201  endif
202  do i=1,nn
203  vect(i)=0.0
204  enddo
205  !** SURFACE LOAD
206  if( isuf.EQ.1 ) then
207  !** INTEGRATION OVER SURFACE
208  do lx=1,3
209  ri=xg(lx)
210  h(1)=-0.5*(1.0-ri)*ri
211  h(2)= 0.5*(1.0+ri)*ri
212  h(3)=1.0-ri*ri
213  hr(1)= ri-0.5
214  hr(2)= ri+0.5
215  hr(3)=-2.0*ri
216  gx=0.0
217  gy=0.0
218  do i=1,3
219  gx=gx+hr(i)*xx(nod(i))
220  gy=gy+hr(i)*yy(nod(i))
221  enddo
222  xsum=dsqrt(gx*gx+gy*gy)
223  wg=wgt(lx)*val*thick
224  vect(nod(1))=vect(nod(1))-h(1)*xsum*wg
225  vect(nod(2))=vect(nod(2))-h(2)*xsum*wg
226  vect(nod(3))=vect(nod(3))-h(3)*xsum*wg
227  enddo
228  endif
229  !** VOLUME LOAD
230  if( ivol.EQ.1 ) then
231 
232  !* LOOP OVER ALL INTEGRATION POINTS
233  do l2=1,3
234  xl2=xg(l2)
235  x2 =(xl2+1.0)*0.5
236  do l1=1,3
237  xl1=xg(l1)
238  x1=0.5*(1.0-x2)*(xl1+1.0)
239  ! INTERPOLATION FUNCTION
240  x3=1.0-x1-x2
241  h(1)= x1*(2.0*x1-1.)
242  h(2)= x2*(2.0*x2-1.)
243  h(3)= x3*(2.0*x3-1.)
244  h(4)= 4.0*x1*x2
245  h(5)= 4.0*x2*x3
246  h(6)= 4.0*x1*x3
247  ! DERIVATIVE OF INTERPOLATION FUNCTION
248  ! FOR L1-COORDINATE
249  hl1(1)=4.0*x1-1.0
250  hl1(2)= 0.0
251  hl1(3)= 0.0
252  hl1(4)= 4.0*x2
253  hl1(5)= 0.0
254  hl1(6)= 4.0*x3
255  ! FOR L2-COORDINATE
256  hl2(1)= 0.0
257  hl2(2)= 4.0*x2-1.0
258  hl2(3)= 0.0
259  hl2(4)= 4.0*x1
260  hl2(5)= 4.0*x3
261  hl2(6)= 0.0
262  ! FOR L3-COORDINATE
263  hl3(1)= 0.0
264  hl3(2)= 0.0
265  hl3(3)= 4.0*x3-1.0
266  hl3(4)= 0.0
267  hl3(5)= 4.0*x2
268  hl3(6)= 4.0*x1
269  ! JACOBI MATRIX
270  xj11=0.0
271  xj21=0.0
272  xj12=0.0
273  xj22=0.0
274  do i=1,nn
275  xj11=xj11+(hl1(i)-hl3(i))*xx(i)
276  xj21=xj21+(hl2(i)-hl3(i))*xx(i)
277  xj12=xj12+(hl1(i)-hl3(i))*yy(i)
278  xj22=xj22+(hl2(i)-hl3(i))*yy(i)
279  enddo
280  det=xj11*xj22-xj21*xj12
281  wg=wgt(l1)*wgt(l2)*det*thick*(1.0-x2)*0.25
282  do i = 1, nn
283  vect(i) = vect(i) - h(i)*wg*val
284  enddo
285  enddo
286  enddo
287  endif
288  return
289 
290  end subroutine heat_dflux_232
291  !C*--------------------------------------------------------------------*
292  subroutine heat_dflux_241(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
293  !C*--------------------------------------------------------------------*
294  !C***
295  !C*** SET 241 DFLUX
296  !C***
297  !C* LTYPE=0 : BODY FLUX
298  use hecmw
299  implicit none
300  !C* I/F VARIABLES
301  integer(kind=kint) :: nn, ltype
302  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
303  !C* LOCAL VARIABLES
304  integer(kind=kint) :: ivol, isuf, i, lx, ly
305  real(kind=kreal) :: gx, gy, xsum
306  real(kind=kreal) :: ri, si, rp, sp, rm, sm
307  real(kind=kreal) :: xj11, xj12, xj21, xj22, det, wg
308  real(kind=kreal) :: xg(2), wgt(2), h(4), hr(4), hs(4)
309  integer(kind=kint) :: nod(2)
310  !*************************
311  ! GAUSS INTEGRATION POINT
312  !*************************
313  data xg/-0.5773502691896,0.5773502691896/
314  data wgt/1.0,1.0/
315  !C*
316  ivol=0
317  isuf=0
318  if( ltype.EQ.0 ) then
319  ivol=1
320  else if( ltype.GE.1 ) then
321  isuf=1
322  if(ltype.EQ.1) then
323  nod(1)=1
324  nod(2)=2
325  else if(ltype.EQ.2) then
326  nod(1)=2
327  nod(2)=3
328  else if(ltype.EQ.3) then
329  nod(1)=3
330  nod(2)=4
331  else if(ltype.EQ.4) then
332  nod(1)=4
333  nod(2)=1
334  endif
335  endif
336  do i=1,nn
337  vect(i)=0.0
338  enddo
339  !** SURFACE LOAD
340  if( isuf.EQ.1 ) then
341  !** INTEGRATION OVER SURFACE
342  do lx=1,2
343  ri=xg(lx)
344  h(1)=0.5*(1.0-ri)
345  h(2)=0.5*(1.0+ri)
346  hr(1)=-0.5
347  hr(2)= 0.5
348  gx=0.0
349  gy=0.0
350  do i=1,2
351  gx=gx+hr(i)*xx(nod(i))
352  gy=gy+hr(i)*yy(nod(i))
353  enddo
354  xsum=dsqrt(gx*gx+gy*gy)*thick
355  do i=1,2
356  vect(nod(i))=vect(nod(i))-xsum*wgt(lx)*h(i)*val
357  enddo
358  enddo
359  endif
360  !** VOLUME LOAD
361  if( ivol.EQ.1 ) then
362  do lx=1,2
363  ri=xg(lx)
364  do ly=1,2
365  si=xg(ly)
366  rp=1.0+ri
367  sp=1.0+si
368  rm=1.0-ri
369  sm=1.0-si
370  h(1)=0.25*rm*sm
371  h(2)=0.25*rp*sm
372  h(3)=0.25*rp*sp
373  h(4)=0.25*rm*sp
374  hr(1)=-.25*sm
375  hr(2)= .25*sm
376  hr(3)= .25*sp
377  hr(4)=-.25*sp
378  hs(1)=-.25*rm
379  hs(2)=-.25*rp
380  hs(3)= .25*rp
381  hs(4)= .25*rm
382  xj11=0.0
383  xj21=0.0
384  xj12=0.0
385  xj22=0.0
386  do i=1,nn
387  xj11=xj11+hr(i)*xx(i)
388  xj21=xj21+hs(i)*xx(i)
389  xj12=xj12+hr(i)*yy(i)
390  xj22=xj22+hs(i)*yy(i)
391  enddo
392  det=xj11*xj22-xj21*xj12
393  wg=wgt(lx)*wgt(ly)*det*thick
394  do i =1,nn
395  vect(i)=vect(i)-wg*h(i)*val
396  enddo
397  enddo
398  enddo
399  endif
400  return
401 
402  end subroutine heat_dflux_241
403  !C*--------------------------------------------------------------------*
404  subroutine heat_dflux_242(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
405  !C*--------------------------------------------------------------------*
406  !C***
407  !C*** SET 242 DFLUX
408  !C***
409  !C* LTYPE=0 : BODY FLUX
410  use hecmw
411  implicit none
412  !C* I/F VARIABLES
413  integer(kind=kint) :: nn, ltype
414  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
415  !C* LOCAL VARIABLES
416  integer(kind=kint) :: ivol, isuf, i, lx, ly
417  real(kind=kreal) :: gx, gy, xsum
418  real(kind=kreal) :: ri, si, rp, sp, rm, sm
419  real(kind=kreal) :: xj11, xj12, xj21, xj22, det, wg
420  real(kind=kreal) :: xg(3), wgt(3), h(8), hr(8), hs(8)
421  integer(kind=kint) :: nod(3)
422  !*************************
423  ! GAUSS INTEGRATION POINT
424  !*************************
425  xg(1) = -0.7745966692
426  xg(2) = 0.0
427  xg(3) = 0.7745966692
428  wgt(1) = 0.5555555555
429  wgt(2) = 0.8888888888
430  wgt(3) = 0.5555555555
431  !C*
432  ivol=0
433  isuf=0
434  if( ltype.EQ.0 ) then
435  ivol=1
436  else if( ltype.GE.1 ) then
437  isuf=1
438  if(ltype.EQ.1) then
439  nod(1)=1
440  nod(2)=2
441  nod(3)=5
442  else if(ltype.EQ.2) then
443  nod(1)=2
444  nod(2)=3
445  nod(3)=6
446  else if(ltype.EQ.3) then
447  nod(1)=3
448  nod(2)=4
449  nod(3)=7
450  else if(ltype.EQ.4) then
451  nod(1)=4
452  nod(2)=1
453  nod(3)=8
454  endif
455  endif
456  do i=1,nn
457  vect(i)=0.0
458  enddo
459  !** SURFACE LOAD
460  if( isuf.EQ.1 ) then
461  !** INTEGRATION OVER SURFACE
462  do lx=1,3
463  ri=xg(lx)
464  h(1)=-0.5*(1.0-ri)*ri
465  h(2)= 0.5*(1.0+ri)*ri
466  h(3)=1.0-ri*ri
467  hr(1)= ri-0.5
468  hr(2)= ri+0.5
469  hr(3)=-2.0*ri
470  gx=0.0
471  gy=0.0
472  do i=1,3
473  gx=gx+hr(i)*xx(nod(i))
474  gy=gy+hr(i)*yy(nod(i))
475  enddo
476  xsum=dsqrt(gx*gx+gy*gy)*thick
477  vect(nod(1))=vect(nod(1))-h(1)*wgt(lx)*val*xsum
478  vect(nod(2))=vect(nod(2))-h(2)*wgt(lx)*val*xsum
479  vect(nod(3))=vect(nod(3))-h(3)*wgt(lx)*val*xsum
480  enddo
481  endif
482  !** VOLUME LOAD
483  if( ivol.EQ.1 ) then
484  do lx=1,3
485  ri=xg(lx)
486  do ly=1,3
487  si=xg(ly)
488  rp=1.0+ri
489  sp=1.0+si
490  rm=1.0-ri
491  sm=1.0-si
492  h(1)=0.25*rm*sm*(-1.0-ri-si)
493  h(2)=0.25*rp*sm*(-1.0+ri-si)
494  h(3)=0.25*rp*sp*(-1.0+ri+si)
495  h(4)=0.25*rm*sp*(-1.0-ri+si)
496  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
497  h(6)=0.5*(1.0-si*si)*(1.0+ri)
498  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
499  h(8)=0.5*(1.0-si*si)*(1.0-ri)
500  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
501  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
502  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
503  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
504  hr(5)=-ri*(1.0-si)
505  hr(6)= 0.5*(1.0-si*si)
506  hr(7)=-ri*(1.0+si)
507  hr(8)=-0.5*(1.0-si*si)
508  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
509  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
510  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
511  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
512  hs(5)=-0.5*(1.0-ri*ri)
513  hs(6)=-si *(1.0+ri)
514  hs(7)= 0.5*(1.0-ri*ri)
515  hs(8)= -si*(1.0-ri)
516  xj11=0.0
517  xj21=0.0
518  xj12=0.0
519  xj22=0.0
520  do i=1,nn
521  xj11=xj11+hr(i)*xx(i)
522  xj21=xj21+hs(i)*xx(i)
523  xj12=xj12+hr(i)*yy(i)
524  xj22=xj22+hs(i)*yy(i)
525  enddo
526  det=xj11*xj22-xj21*xj12
527  wg=wgt(lx)*wgt(ly)*det*thick
528  do i =1,nn
529  vect(i)=vect(i)-wg*h(i)*val
530  enddo
531  enddo
532  enddo
533  endif
534  return
535 
536  end subroutine heat_dflux_242
537  !----------------------------------------------------------------------*
538  subroutine heat_dflux_341(NN,XX,YY,ZZ,LTYPE,val,VECT)
539  !----------------------------------------------------------------------*
540  !**
541  !** SET 341 DFLUX
542  !**
543  ! BF LTYPE=0 :BODY FLUX
544  ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
545  ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
546  ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
547  ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
548  use hecmw
549  implicit none
550  ! I/F VARIABLES
551  integer(kind=kint) :: nn, ltype
552  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), vect(nn)
553  ! LOCAL VARIABLES
554  real(kind=kreal) :: h(4), hr(4), hs(4), ht(4)
555  real(kind=kreal) :: xg(2), wgt(2)
556  real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
557  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33, det, wg
558  integer(kind=kint) :: ivol, isuf
559  integer(kind=kint) :: nod(4)
560  integer(kind=kint) :: L1, L2, L3, I
561  real(kind=kreal) :: val, aa
562  real(kind=kreal) :: v1x, v1y, v1z
563  real(kind=kreal) :: v2x, v2y, v2z
564  real(kind=kreal) :: v3x, v3y, v3z
565  real(kind=kreal) :: xl1, xl2, xl3
566  real(kind=kreal) :: x1, x2, x3
567  !*************************
568  ! GAUSS INTEGRATION POINT
569  !*************************
570  data xg/-0.5773502691896,0.5773502691896/
571  data wgt/1.0,1.0/
572  !
573  ivol=0
574  isuf=0
575  nod = 0
576  if( ltype.EQ.0 ) then
577  ivol=1
578  else
579  isuf=1
580  if(ltype.EQ.1) then
581  nod(1)=1
582  nod(2)=2
583  nod(3)=3
584  else if(ltype.EQ.2) then
585  nod(1)=4
586  nod(2)=2
587  nod(3)=1
588  else if(ltype.EQ.3) then
589  nod(1)=4
590  nod(2)=3
591  nod(3)=2
592  else if(ltype.EQ.4) then
593  nod(1)=4
594  nod(2)=1
595  nod(3)=3
596  endif
597  endif
598  ! CLEAR VECT
599  do i=1,nn
600  vect(i)=0.0
601  enddo
602  !** SURFACE LOAD
603  if( isuf.EQ.1 ) then
604  v1x=xx(nod(2))-xx(nod(1))
605  v1y=yy(nod(2))-yy(nod(1))
606  v1z=zz(nod(2))-zz(nod(1))
607  v2x=xx(nod(3))-xx(nod(1))
608  v2y=yy(nod(3))-yy(nod(1))
609  v2z=zz(nod(3))-zz(nod(1))
610  v3x= v1y*v2z-v1z*v2y
611  v3y= v1z*v2x-v1x*v2z
612  v3z= v1x*v2y-v1y*v2x
613  aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
614  vect(nod(1))= -val*aa/3.0
615  vect(nod(2))= -val*aa/3.0
616  vect(nod(3))= -val*aa/3.0
617  endif
618  !** VOLUME LOAD
619  if( ivol.EQ.1 ) then
620  ! LOOP FOR INTEGRATION POINTS
621  do l3=1,2
622  xl3=xg(l3)
623  x3 =(xl3+1.0)*0.5
624  do l2=1,2
625  xl2=xg(l2)
626  x2 =(1.0-x3)*(xl2+1.0)*0.5
627  do l1=1,2
628  xl1=xg(l1)
629  x1=(1.0-x2-x3)*(xl1+1.0)*0.5
630  ! INTERPOLATION FUNCTION
631  h(1)=x1
632  h(2)=x2
633  h(3)=x3
634  h(4)=1.0-x1-x2-x3
635  ! DERIVATIVE OF INTERPOLATION FUNCTION
636  ! FOR L1-COORDINATE
637  hr(1)= 1.0
638  hr(2)= 0.0
639  hr(3)= 0.0
640  hr(4)=-1.0
641  ! FOR L2-COORDINATE
642  hs(1)= 0.0
643  hs(2)= 1.0
644  hs(3)= 0.0
645  hs(4)=-1.0
646  ! FOR ZETA-COORDINATE
647  ht(1)= 0.0
648  ht(2)= 0.0
649  ht(3)= 1.0
650  ht(4)=-1.0
651  ! JACOBI MATRIX
652  xj11=hr(1)*xx(1)+hr(2)*xx(2)+hr(3)*xx(3)+hr(4)*xx(4)
653  xj21=hs(1)*xx(1)+hs(2)*xx(2)+hs(3)*xx(3)+hs(4)*xx(4)
654  xj31=ht(1)*xx(1)+ht(2)*xx(2)+ht(3)*xx(3)+ht(4)*xx(4)
655  xj12=hr(1)*yy(1)+hr(2)*yy(2)+hr(3)*yy(3)+hr(4)*yy(4)
656  xj22=hs(1)*yy(1)+hs(2)*yy(2)+hs(3)*yy(3)+hs(4)*yy(4)
657  xj32=ht(1)*yy(1)+ht(2)*yy(2)+ht(3)*yy(3)+ht(4)*yy(4)
658  xj13=hr(1)*zz(1)+hr(2)*zz(2)+hr(3)*zz(3)+hr(4)*zz(4)
659  xj23=hs(1)*zz(1)+hs(2)*zz(2)+hs(3)*zz(3)+hs(4)*zz(4)
660  xj33=ht(1)*zz(1)+ht(2)*zz(2)+ht(3)*zz(3)+ht(4)*zz(4)
661  ! DETERMINANT OF JACOBIAN
662  det=xj11*xj22*xj33 &
663  +xj12*xj23*xj31 &
664  +xj13*xj21*xj32 &
665  -xj13*xj22*xj31 &
666  -xj12*xj21*xj33 &
667  -xj11*xj23*xj32
668  !
669  do i = 1, nn
670  vect(i) = vect(i) + val*h(i)*det*(1.0-x3)*(1.0-x2-x3)*0.125
671  enddo
672 
673  enddo
674  enddo
675  enddo
676  endif
677  !
678  return
679 
680  end subroutine heat_dflux_341
681  !----------------------------------------------------------------------*
682  subroutine heat_dflux_342(NN,XX,YY,ZZ,LTYPE,val,VECT)
683  !----------------------------------------------------------------------*
684  !**
685  !** SET 342 DFLUX
686  !**
687  ! BF LTYPE=0 :BODY FLUX
688  ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
689  ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
690  ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
691  ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
692  use hecmw
693  implicit none
694  ! I/F VARIABLES
695  integer(kind=kint) :: NN, LTYPE
696  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), vect(nn)
697  ! LOCAL VARIABLES
698  real(kind=kreal) :: h(10)
699  real(kind=kreal) :: hl1(10), hl2(10), hl3(10), hl4(10)
700  real(kind=kreal) :: xg(3), wgt(3)
701  real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
702  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
703  real(kind=kreal) :: det, wg
704  integer(kind=kint) :: IVOL, ISUF
705  integer(kind=kint) :: NOD(10)
706  integer(kind=kint) :: LX, LY, LZ, I, IG1, IG2, L1, L2, L3
707  real(kind=kreal) :: val, xsum
708  real(kind=kreal) :: g1x, g1y, g1z
709  real(kind=kreal) :: g2x, g2y, g2z
710  real(kind=kreal) :: g3x, g3y, g3z
711  real(kind=kreal) :: xl1, xl2, xl3
712  real(kind=kreal) :: x1, x2, x3, x4
713  !*************************
714  ! GAUSS INTEGRATION POINT
715  !*************************
716  xg(1) = -0.7745966692
717  xg(2) = 0.0
718  xg(3) = 0.7745966692
719  wgt(1) = 0.5555555555
720  wgt(2) = 0.8888888888
721  wgt(3) = 0.5555555555
722  !
723  ivol=0
724  isuf=0
725  if( ltype.EQ.0 ) then
726  ivol=1
727  else
728  isuf=1
729  if(ltype.EQ.1) then
730  nod(1)=1
731  nod(2)=2
732  nod(3)=3
733  nod(4)=5
734  nod(5)=6
735  nod(6)=7
736  else if(ltype.EQ.2) then
737  nod(1)=4
738  nod(2)=2
739  nod(3)=1
740  nod(4)=9
741  nod(5)=5
742  nod(6)=8
743  else if(ltype.EQ.3) then
744  nod(1)=4
745  nod(2)=3
746  nod(3)=2
747  nod(4)=10
748  nod(5)=6
749  nod(6)=9
750  else if(ltype.EQ.4) then
751  nod(1)=4
752  nod(2)=1
753  nod(3)=3
754  nod(4)=8
755  nod(5)=7
756  nod(6)=10
757  endif
758  endif
759  ! CLEAR VECT
760  do i=1,10
761  vect(i)=0.0
762  enddo
763  !** SURFACE LOAD
764  if( isuf.EQ.1 ) then
765  ! INTEGRATION OVER SURFACE
766  do l2=1,3
767  xl2=xg(l2)
768  x2 =(xl2+1.0)*0.5
769  do l1=1,3
770  xl1=xg(l1)
771  x1=0.5*(1.0-x2)*(xl1+1.0)
772  ! INTERPOLATION FUNCTION
773  x3=1.0-x1-x2
774  h(1)= x1*(2.0*x1-1.)
775  h(2)= x2*(2.0*x2-1.)
776  h(3)= x3*(2.0*x3-1.)
777  h(4)= 4.0*x1*x2
778  h(5)= 4.0*x2*x3
779  h(6)= 4.0*x1*x3
780  ! DERIVATIVE OF INTERPOLATION FUNCTION
781  ! FOR L1-COORDINATE
782  hl1(1)=4.0*x1-1.0
783  hl1(2)= 0.0
784  hl1(3)= 0.0
785  hl1(4)= 4.0*x2
786  hl1(5)= 0.0
787  hl1(6)= 4.0*x3
788  ! FOR L2-COORDINATE
789  hl2(1)= 0.0
790  hl2(2)= 4.0*x2-1.0
791  hl2(3)= 0.0
792  hl2(4)= 4.0*x1
793  hl2(5)= 4.0*x3
794  hl2(6)= 0.0
795  ! FOR L3-COORDINATE
796  hl3(1)= 0.0
797  hl3(2)= 0.0
798  hl3(3)= 4.0*x3-1.0
799  hl3(4)= 0.0
800  hl3(5)= 4.0*x2
801  hl3(6)= 4.0*x1
802  ! JACOBI MATRIX
803  g1x=0.0
804  g1y=0.0
805  g1z=0.0
806  g2x=0.0
807  g2y=0.0
808  g2z=0.0
809  do i=1,6
810  g1x=g1x+(hl1(i)-hl3(i))*xx(nod(i))
811  g1y=g1y+(hl1(i)-hl3(i))*yy(nod(i))
812  g1z=g1z+(hl1(i)-hl3(i))*zz(nod(i))
813  g2x=g2x+(hl2(i)-hl3(i))*xx(nod(i))
814  g2y=g2y+(hl2(i)-hl3(i))*yy(nod(i))
815  g2z=g2z+(hl2(i)-hl3(i))*zz(nod(i))
816  enddo
817  g3x=g1y*g2z-g1z*g2y
818  g3y=g1z*g2x-g1x*g2z
819  g3z=g1x*g2y-g1y*g2x
820  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
821  g3x=g3x/xsum
822  g3y=g3y/xsum
823  g3z=g3z/xsum
824  ! JACOBI MATRIX
825  xj11=g1x
826  xj12=g1y
827  xj13=g1z
828  xj21=g2x
829  xj22=g2y
830  xj23=g2z
831  xj31=g3x
832  xj32=g3y
833  xj33=g3z
834  !DETERMINANT OF JACOBIAN
835  det=xj11*xj22*xj33 &
836  +xj12*xj23*xj31 &
837  +xj13*xj21*xj32 &
838  -xj13*xj22*xj31 &
839  -xj12*xj21*xj33 &
840  -xj11*xj23*xj32
841  wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25
842  do i=1,6
843  vect(nod(i))=vect(nod(i))-val*wg*h(i)
844  enddo
845  !
846  enddo
847  enddo
848  endif
849  !** VOLUME LOAD
850  if( ivol.EQ.1 ) then
851  ! LOOP FOR INTEGRATION POINTS
852  !DETERMINANT OF JACOBIAN
853  do l3=1,3
854  xl3=xg(l3)
855  x3 =(xl3+1.0)*0.5
856  do l2=1,3
857  xl2=xg(l2)
858  x2 =(1.0-x3)*(xl2+1.0)*0.5
859  do l1=1,3
860  xl1=xg(l1)
861  x1=(1.0-x2-x3)*(xl1+1.0)*0.5
862  !C* INTERPOLATION FUNCTION
863  x4=1.0-x1-x2-x3
864  h(1)=x1*(2.0*x1-1.0)
865  h(2)=x2*(2.0*x2-1.0)
866  h(3)=x3*(2.0*x3-1.0)
867  h(4)=x4*(2.0*x4-1.0)
868  h(5)=4.0*x1*x2
869  h(6)=4.0*x2*x3
870  h(7)=4.0*x1*x3
871  h(8)=4.0*x1*x4
872  h(9)=4.0*x2*x4
873  h(10)=4.0*x3*x4
874  !C* DERIVATIVE OF INTERPOLATION FUNCTION
875  !C* FOR L1-COORDINATE
876  hl1(1)= 4.0*x1-1.0
877  hl1(2)= 0.0
878  hl1(3)= 0.0
879  hl1(4)= 0.0
880  hl1(5)= 4.0*x2
881  hl1(6)= 0.0
882  hl1(7)= 4.0*x3
883  hl1(8)= 4.0*x4
884  hl1(9)= 0.0
885  hl1(10)= 0.0
886  !C* FOR L2-COORDINATE
887  hl2(1)= 0.0
888  hl2(2)= 4.0*x2-1.0
889  hl2(3)= 0.0
890  hl2(4)= 0.0
891  hl2(5)= 4.0*x1
892  hl2(6)= 4.0*x3
893  hl2(7)= 0.0
894  hl2(8)= 0.0
895  hl2(9)= 4.0*x4
896  hl2(10)= 0.0
897  !C* FOR L3-COORDINATE
898  hl3(1)= 0.0
899  hl3(2)= 0.0
900  hl3(3)= 4.0*x3-1.0
901  hl3(4)= 0.0
902  hl3(5)= 0.0
903  hl3(6)= 4.0*x2
904  hl3(7)= 4.0*x1
905  hl3(8)= 0.0
906  hl3(9)= 0.0
907  hl3(10)= 4.0*x4
908  !C* FOR L4-COORDINATE
909  hl4(1)= 0.0
910  hl4(2)= 0.0
911  hl4(3)= 0.0
912  hl4(4)= 4.0*x4-1.0
913  hl4(5)= 0.0
914  hl4(6)= 0.0
915  hl4(7)= 0.0
916  hl4(8)= 4.0*x1
917  hl4(9)= 4.0*x2
918  hl4(10)= 4.0*x3
919  !C* JACOBI MATRIX
920  xj11=0.0
921  xj21=0.0
922  xj31=0.0
923  xj12=0.0
924  xj22=0.0
925  xj32=0.0
926  xj13=0.0
927  xj23=0.0
928  xj33=0.0
929  do i=1,nn
930  xj11=xj11+(hl1(i)-hl4(i))*xx(i)
931  xj21=xj21+(hl2(i)-hl4(i))*xx(i)
932  xj31=xj31+(hl3(i)-hl4(i))*xx(i)
933  xj12=xj12+(hl1(i)-hl4(i))*yy(i)
934  xj22=xj22+(hl2(i)-hl4(i))*yy(i)
935  xj32=xj32+(hl3(i)-hl4(i))*yy(i)
936  xj13=xj13+(hl1(i)-hl4(i))*zz(i)
937  xj23=xj23+(hl2(i)-hl4(i))*zz(i)
938  xj33=xj33+(hl3(i)-hl4(i))*zz(i)
939  enddo
940  det=xj11*xj22*xj33 &
941  +xj12*xj23*xj31 &
942  +xj13*xj21*xj32 &
943  -xj13*xj22*xj31 &
944  -xj12*xj21*xj33 &
945  -xj11*xj23*xj32
946  det=-det
947  wg=det*wgt(l1)*wgt(l2)*wgt(l3)*(1.-x3)*(1.0-x2-x3)*0.125
948  do i = 1, nn
949  vect(i) = vect(i) - val*wg*h(i)
950  enddo
951  !
952  enddo
953  enddo
954  enddo
955  endif
956  !
957  return
958 
959  end subroutine heat_dflux_342
960  !----------------------------------------------------------------------*
961  subroutine heat_dflux_351(NN,XX,YY,ZZ,LTYPE,val,VECT)
962  !----------------------------------------------------------------------*
963  !**
964  !** SET 351 DFLUX
965  !**
966  ! BF LTYPE=0 :BODY FLUX
967  ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
968  ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
969  ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
970  ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
971  ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
972  use hecmw
973  implicit none
974  ! I/F VARIABLES
975  integer(kind=kint) :: NN, LTYPE
976  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
977  ! LOCAL VARIABLES
978  integer(kind=kint) :: IVOL, ISUF, I, IG1, IG2, L12, LZ
979  real(kind=kreal) :: ri, si, xsum
980  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
981  real(kind=kreal) :: v1x, v1y, v1z, v2x, v2y, v2z, v3x, v3y, v3z, aa
982  real(kind=kreal) :: zi, x1, x2
983  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
984  real(kind=kreal) :: det, wg
985  real(kind=kreal) :: h(6), hr(6), hs(6), ht(6), pl(6)
986  real(kind=kreal) :: xg(2), wgt(2), xg1(3), xg2(3), wgt1(3)
987  integer(kind=kint) :: NOD(4)
988  !*************************
989  ! GAUSS INTEGRATION POINT
990  !*************************
991  data xg / -0.5773502691896,0.5773502691896 /
992  data wgt / 1.0, 1.0 /
993  data xg1 / 0.6666666667, 0.16666666667, 0.16666666667 /
994  data xg2 / 0.1666666667, 0.66666666667, 0.16666666667 /
995  data wgt1/ 0.1666666667, 0.16666666667, 0.16666666667 /
996  !
997  ! SELECTION OF LOAD TYPE
998  ivol=0
999  isuf=0
1000  if( ltype.EQ.0 ) then
1001  ivol=1
1002  else
1003  isuf=1
1004  if( ltype.EQ.1 ) then
1005  nod(1) = 1
1006  nod(2) = 2
1007  nod(3) = 3
1008  isuf=2
1009  else if( ltype.EQ.2 ) then
1010  nod(1) = 6
1011  nod(2) = 5
1012  nod(3) = 4
1013  isuf=2
1014  else if( ltype.EQ.3 ) then
1015  nod(1) = 4
1016  nod(2) = 5
1017  nod(3) = 2
1018  nod(4) = 1
1019  else if( ltype.EQ.4 ) then
1020  nod(1) = 5
1021  nod(2) = 6
1022  nod(3) = 3
1023  nod(4) = 2
1024  else if( ltype.EQ.5 ) then
1025  nod(1) = 6
1026  nod(2) = 4
1027  nod(3) = 1
1028  nod(4) = 3
1029  endif
1030  endif
1031  do i=1,nn
1032  vect(i)=0.0
1033  enddo
1034  !
1035  !** SURFACE LOAD
1036  !
1037  if( isuf.EQ.1 ) then
1038  do ig2=1,2
1039  si=xg(ig2)
1040  do ig1=1,2
1041  ri=xg(ig1)
1042  !
1043  h(1)=0.25*(1.0-ri)*(1.0-si)
1044  h(2)=0.25*(1.0+ri)*(1.0-si)
1045  h(3)=0.25*(1.0+ri)*(1.0+si)
1046  h(4)=0.25*(1.0-ri)*(1.0+si)
1047  hr(1)=-.25*(1.0-si)
1048  hr(2)= .25*(1.0-si)
1049  hr(3)= .25*(1.0+si)
1050  hr(4)=-.25*(1.0+si)
1051  hs(1)=-.25*(1.0-ri)
1052  hs(2)=-.25*(1.0+ri)
1053  hs(3)= .25*(1.0+ri)
1054  hs(4)= .25*(1.0-ri)
1055  !
1056  g1x=0.0
1057  g1y=0.0
1058  g1z=0.0
1059  g2x=0.0
1060  g2y=0.0
1061  g2z=0.0
1062  !
1063  do i=1,4
1064  g1x=g1x+hr(i)*xx(nod(i))
1065  g1y=g1y+hr(i)*yy(nod(i))
1066  g1z=g1z+hr(i)*zz(nod(i))
1067  g2x=g2x+hs(i)*xx(nod(i))
1068  g2y=g2y+hs(i)*yy(nod(i))
1069  g2z=g2z+hs(i)*zz(nod(i))
1070  enddo
1071  !
1072  g3x=g1y*g2z-g1z*g2y
1073  g3y=g1z*g2x-g1x*g2z
1074  g3z=g1x*g2y-g1y*g2x
1075  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1076  do i=1,4
1077  vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1078  enddo
1079  !
1080  enddo
1081  enddo
1082  elseif( isuf.EQ.2 ) then
1083  v1x=xx(nod(2))-xx(nod(1))
1084  v1y=yy(nod(2))-yy(nod(1))
1085  v1z=zz(nod(2))-zz(nod(1))
1086  v2x=xx(nod(3))-xx(nod(1))
1087  v2y=yy(nod(3))-yy(nod(1))
1088  v2z=zz(nod(3))-zz(nod(1))
1089  v3x= v1y*v2z-v1z*v2y
1090  v3y= v1z*v2x-v1x*v2z
1091  v3z= v1x*v2y-v1y*v2x
1092  aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
1093  vect(nod(1))= -val*aa/3.0
1094  vect(nod(2))= -val*aa/3.0
1095  vect(nod(3))= -val*aa/3.0
1096  endif
1097  !
1098  !** VOLUME LOAD
1099  !
1100  if( ivol.EQ.1 ) then
1101  do i=1,6
1102  vect(i)=0.0
1103  enddo
1104  !
1105  do lz=1,2
1106  zi=xg(lz)
1107  do l12=1,3
1108  x1=xg1(l12)
1109  x2=xg2(l12)
1110 
1111  h(1)=x1*(1.0-zi)*0.5
1112  h(2)=x2*(1.0-zi)*0.5
1113  h(3)=(1.0-x1-x2)*(1.0-zi)*0.5
1114  h(4)=x1*(1.0+zi)*0.5
1115  h(5)=x2*(1.0+zi)*0.5
1116  h(6)=(1.0-x1-x2)*(1.0+zi)*0.5
1117  hr(1)= (1.0-zi)*0.5
1118  hr(2)= 0.0
1119  hr(3)=-(1.0-zi)*0.5
1120  hr(4)= (1.0+zi)*0.5
1121  hr(5)= 0.0
1122  hr(6)=-(1.0+zi)*0.5
1123  hs(1)= 0.0
1124  hs(2)= (1.0-zi)*0.5
1125  hs(3)= -(1.0-zi)*0.5
1126  hs(4)= 0.0
1127  hs(5)= (1.0+zi)*0.5
1128  hs(6)= -(1.0+zi)*0.5
1129  ht(1)= -0.5*x1
1130  ht(2)= -0.5*x2
1131  ht(3)= -0.5*(1.0-x1-x2)
1132  ht(4)= 0.5*x1
1133  ht(5)= 0.5*x2
1134  ht(6)= 0.5*(1.0-x1-x2)
1135  !
1136  xj11=hr(1)*xx(1)+hr(2)*xx(2)+hr(3)*xx(3)+hr(4)*xx(4) &
1137  +hr(5)*xx(5)+hr(6)*xx(6)
1138  xj21=hs(1)*xx(1)+hs(2)*xx(2)+hs(3)*xx(3)+hs(4)*xx(4) &
1139  +hs(5)*xx(5)+hs(6)*xx(6)
1140  xj31=ht(1)*xx(1)+ht(2)*xx(2)+ht(3)*xx(3)+ht(4)*xx(4) &
1141  +ht(5)*xx(5)+ht(6)*xx(6)
1142  xj12=hr(1)*yy(1)+hr(2)*yy(2)+hr(3)*yy(3)+hr(4)*yy(4) &
1143  +hr(5)*yy(5)+hr(6)*yy(6)
1144  xj22=hs(1)*yy(1)+hs(2)*yy(2)+hs(3)*yy(3)+hs(4)*yy(4) &
1145  +hs(5)*yy(5)+hs(6)*yy(6)
1146  xj32=ht(1)*yy(1)+ht(2)*yy(2)+ht(3)*yy(3)+ht(4)*yy(4) &
1147  +ht(5)*yy(5)+ht(6)*yy(6)
1148  xj13=hr(1)*zz(1)+hr(2)*zz(2)+hr(3)*zz(3)+hr(4)*zz(4) &
1149  +hr(5)*zz(5)+hr(6)*zz(6)
1150  xj23=hs(1)*zz(1)+hs(2)*zz(2)+hs(3)*zz(3)+hs(4)*zz(4) &
1151  +hs(5)*zz(5)+hs(6)*zz(6)
1152  xj33=ht(1)*zz(1)+ht(2)*zz(2)+ht(3)*zz(3)+ht(4)*zz(4) &
1153  +ht(5)*zz(5)+ht(6)*zz(6)
1154  !*DETERMINANT OF JACOBIAN
1155  det=xj11*xj22*xj33 &
1156  +xj12*xj23*xj31 &
1157  +xj13*xj21*xj32 &
1158  -xj13*xj22*xj31 &
1159  -xj12*xj21*xj33 &
1160  -xj11*xj23*xj32
1161  wg=wgt1(l12)*wgt(lz)*det
1162  do i=1,nn
1163  vect(i)=vect(i)-val*h(i)*wg
1164  enddo
1165  !
1166  enddo
1167  enddo
1168  endif
1169  !
1170  return
1171 
1172  end subroutine heat_dflux_351
1173  !----------------------------------------------------------------------*
1174  subroutine heat_dflux_352(NN,XX,YY,ZZ,LTYPE,val,VECT)
1175  !----------------------------------------------------------------------*
1176  !**
1177  !** SET 352 DFLUX
1178  !**
1179  ! BF LTYPE=0 :BODY FLUX
1180  ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
1181  ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
1182  ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
1183  ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
1184  ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
1185  use hecmw
1186  implicit none
1187  ! I/F VARIABLES
1188  integer(kind=kint) :: NN, LTYPE
1189  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
1190  ! LOCAL VARIABLES
1191  integer(kind=kint) :: IVOL, ISUF, I, IG1, IG2, LX, LY, LZ
1192  real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm, xsum
1193  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
1194  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
1195  real(kind=kreal) :: det, wg
1196  real(kind=kreal) :: h(15), hr(15), hs(15), ht(15), pl(15)
1197  real(kind=kreal) :: xg(3), wgt(3)
1198  integer(kind=kint) :: NOD(8)
1199  !*************************
1200  ! GAUSS INTEGRATION POINT
1201  !*************************
1202  xg(1) = -0.7745966692
1203  xg(2) = 0.0
1204  xg(3) = 0.7745966692
1205  wgt(1) = 0.5555555555
1206  wgt(2) = 0.8888888888
1207  wgt(3) = 0.5555555555
1208  !
1209  ! SELECTION OF LOAD TYPE
1210  !
1211  ivol=0
1212  isuf=0
1213  if( ltype.EQ.0 ) then
1214  ivol=1
1215  else
1216  isuf=1
1217  if( ltype.EQ.1 ) then
1218  nod(1) = 1
1219  nod(2) = 2
1220  nod(3) = 3
1221  nod(4) = 7
1222  nod(5) = 8
1223  nod(6) = 9
1224  isuf=2
1225  else if( ltype.EQ.2 ) then
1226  nod(1) = 6
1227  nod(2) = 5
1228  nod(3) = 4
1229  nod(4) = 11
1230  nod(5) = 10
1231  nod(6) = 12
1232  isuf=2
1233  else if( ltype.EQ.3 ) then
1234  nod(1) = 4
1235  nod(2) = 5
1236  nod(3) = 2
1237  nod(4) = 1
1238  nod(5) = 10
1239  nod(6) = 14
1240  nod(7) = 7
1241  nod(8) = 13
1242  else if( ltype.EQ.4 ) then
1243  nod(1) = 5
1244  nod(2) = 6
1245  nod(3) = 3
1246  nod(4) = 2
1247  nod(5) = 11
1248  nod(6) = 15
1249  nod(7) = 8
1250  nod(8) = 14
1251  else if( ltype.EQ.5 ) then
1252  nod(1) = 6
1253  nod(2) = 4
1254  nod(3) = 1
1255  nod(4) = 3
1256  nod(5) = 12
1257  nod(6) = 13
1258  nod(7) = 9
1259  nod(8) = 15
1260  endif
1261  endif
1262  do i=1,nn
1263  vect(i)=0.0
1264  enddo
1265  !
1266  !** SURFACE LOAD
1267  !
1268  if( isuf.EQ.1 ) then
1269  do ig2=1,3
1270  si=xg(ig2)
1271  do ig1=1,3
1272  ri=xg(ig1)
1273  rp=1.0+ri
1274  sp=1.0+si
1275  rm=1.0-ri
1276  sm=1.0-si
1277  h(1)=0.25*rm*sm*(-1.0-ri-si)
1278  h(2)=0.25*rp*sm*(-1.0+ri-si)
1279  h(3)=0.25*rp*sp*(-1.0+ri+si)
1280  h(4)=0.25*rm*sp*(-1.0-ri+si)
1281  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
1282  h(6)=0.5*(1.0-si*si)*(1.0+ri)
1283  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
1284  h(8)=0.5*(1.0-si*si)*(1.0-ri)
1285  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1286  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1287  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
1288  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
1289  hr(5)=-ri*(1.0-si)
1290  hr(6)= .5*(1.0-si*si)
1291  hr(7)=-ri*(1.0+si)
1292  hr(8)=-.5*(1.0-si*si)
1293  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1294  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1295  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
1296  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
1297  hs(5)=-.5*(1.0-ri*ri)
1298  hs(6)=-si*(1.0+ri)
1299  hs(7)= .5*(1.0-ri*ri)
1300  hs(8)=-si*(1.0-ri)
1301  g1x=0.0
1302  g1y=0.0
1303  g1z=0.0
1304  g2x=0.0
1305  g2y=0.0
1306  g2z=0.0
1307  g3x=0.0
1308  g3y=0.0
1309  g3z=0.0
1310  do i=1,8
1311  g1x=g1x+hr(i)*xx(nod(i))
1312  g1y=g1y+hr(i)*yy(nod(i))
1313  g1z=g1z+hr(i)*zz(nod(i))
1314  g2x=g2x+hs(i)*xx(nod(i))
1315  g2y=g2y+hs(i)*yy(nod(i))
1316  g2z=g2z+hs(i)*zz(nod(i))
1317  g3x=g3x+ht(i)*xx(nod(i))
1318  g3y=g3y+ht(i)*yy(nod(i))
1319  g3z=g3z+ht(i)*zz(nod(i))
1320  enddo
1321  g3x=g1y*g2z-g1z*g2y
1322  g3y=g1z*g2x-g1x*g2z
1323  g3z=g1x*g2y-g1y*g2x
1324  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1325  do i=1,8
1326  vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1327  enddo
1328  enddo
1329  enddo
1330  elseif( isuf.EQ.2 ) then
1331  do ig2=1,3
1332  si=xg(ig2)
1333  do ig1=1,3
1334  ri=xg(ig1)
1335  rp=1.0+ri
1336  sp=1.0+si
1337  rm=1.0-ri
1338  sm=1.0-si
1339  h(1)=0.25*rm*sm*(-1.0-ri-si)
1340  h(2)=0.25*rp*sm*(-1.0+ri-si)
1341  h(3)=0.5*sp*si
1342  h(4)=0.5*(1.0-ri*ri)*(1.0-si)
1343  h(5)=0.5*(1.0-si*si)*(1.0+ri)
1344  h(6)=0.5*(1.0-si*si)*(1.0-ri)
1345  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1346  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1347  hr(3)= 0.0
1348  hr(4)=-ri*(1.0-si)
1349  hr(5)= 0.5*(1.0-si*si)
1350  hr(6)=-0.5*(1.0-si*si)
1351  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1352  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1353  hs(3)=0.5*(1.0+2.0*si)
1354  hs(4)=-0.5*(1.0-ri*ri)
1355  hs(5)=-si *(1.0+ri)
1356  hs(6)=-si*(1.0-ri)
1357  g1x=0.0
1358  g1y=0.0
1359  g1z=0.0
1360  g2x=0.0
1361  g2y=0.0
1362  g2z=0.0
1363  g3x=0.0
1364  g3y=0.0
1365  g3z=0.0
1366  do i=1, 6
1367  g1x=g1x+hr(i)*xx(nod(i))
1368  g1y=g1y+hr(i)*yy(nod(i))
1369  g1z=g1z+hr(i)*zz(nod(i))
1370  g2x=g2x+hs(i)*xx(nod(i))
1371  g2y=g2y+hs(i)*yy(nod(i))
1372  g2z=g2z+hs(i)*zz(nod(i))
1373  g3x=g3x+ht(i)*xx(nod(i))
1374  g3y=g3y+ht(i)*yy(nod(i))
1375  g3z=g3z+ht(i)*zz(nod(i))
1376  enddo
1377  g3x=g1y*g2z-g1z*g2y
1378  g3y=g1z*g2x-g1x*g2z
1379  g3z=g1x*g2y-g1y*g2x
1380  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1381  do i=1,6
1382  vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1383  enddo
1384  enddo
1385  enddo
1386  endif
1387  !
1388  !** VOLUME LOAD
1389  !
1390  if( ivol.EQ.1 ) then
1391  do lx=1,3
1392  ri=xg(lx)
1393  do ly=1,3
1394  si=xg(ly)
1395  do lz=1,3
1396  ti=xg(lz)
1397  rp=1.0+ri
1398  sp=1.0+si
1399  tp=1.0+ti
1400  rm=1.0-ri
1401  sm=1.0-si
1402  tm=1.0-ti
1403  ! INTERPOLATION FUNCTION
1404  h(1)=-0.125*rm*sm*tm*(2.0+ri+si+ti)
1405  h(2)=-0.125*rp*sm*tm*(2.0-ri+si+ti)
1406  h(3)=-0.25*sp*tm*(1.0-si+ti)
1407  h(4)=-0.125*rm*sm*tp*(2.0+ri+si-ti)
1408  h(5)=-0.125*rp*sm*tp*(2.0-ri+si-ti)
1409  h(6)=0.25*sp*tp*(si+ti-1.0)
1410  h(7)=0.25*(1.0-ri**2)*sm*tm
1411  h(8)=0.25*rp*(1.0-si**2)*tm
1412  h(9)=0.25*rm*(1.0-si**2)*tm
1413  h(10)=0.25*(1.0-ri**2)*sm*tp
1414  h(11)=0.25*rp*(1.0-si**2)*tp
1415  h(12)=0.25*rm*(1.0-si**2)*tp
1416  h(13)=0.25*rm*sm*(1.0-ti**2)
1417  h(14)=0.25*rp*sm*(1.0-ti**2)
1418  h(15)=0.5*sp*(1.0-ti**2)
1419  ! DERIVATIVE OF INTERPOLATION FUNCTION
1420  ! FOR R-COORDINATE
1421  hr(1)=-0.125*rm*sm*tm+0.125*sm*tm*(2.0+ri+si+ti)
1422  hr(2)=+0.125*rp*sm*tm-0.125*sm*tm*(2.0-ri+si+ti)
1423  hr(3)=0.0
1424  hr(4)=-0.125*rm*sm*tp+0.125*sm*tp*(2.0+ri+si-ti)
1425  hr(5)=+0.125*rp*sm*tp-0.125*sm*tp*(2.0-ri+si-ti)
1426  hr(6)=0.0
1427  hr(7)=-0.50*ri*sm*tm
1428  hr(8)=+0.25*(1.0-si**2)*tm
1429  hr(9)=-0.25*(1.0-si**2)*tm
1430  hr(10)=-0.50*ri*sm*tp
1431  hr(11)=+0.25*(1.0-si**2)*tp
1432  hr(12)=-0.25*(1.0-si**2)*tp
1433  hr(13)=-0.25*sm*(1.0-ti**2)
1434  hr(14)=+0.25*sm*(1.0-ti**2)
1435  hr(15)=0.0
1436  ! FOR S-COORDINATE
1437  hs(1)=-0.125*rm*sm*tm+0.125*rm*tm*(2.0+ri+si+ti)
1438  hs(2)=-0.125*rp*sm*tm+0.125*rp*tm*(2.0-ri+si+ti)
1439  hs(3)=+0.25*tm*(2.0*si-ti)
1440  hs(4)=-0.125*rm*sm*tp+0.125*rm*tp*(2.0+ri+si-ti)
1441  hs(5)=-0.125*rp*sm*tp+0.125*rp*tp*(2.0-ri+si-ti)
1442  hs(6)=+0.25*tp*(2.0*si+ti)
1443  hs(7)=-0.25*(1.0-ri**2)*tm
1444  hs(8)=-0.50*rp*si*tm
1445  hs(9)=-0.50*rm*si*tm
1446  hs(10)=-0.25*(1.0-ri**2)*tp
1447  hs(11)=-0.50*rp*si*tp
1448  hs(12)=-0.50*rm*si*tp
1449  hs(13)=-0.25*rm*(1.0-ti**2)
1450  hs(14)=-0.25*rp*(1.0-ti**2)
1451  hs(15)=+0.50*(1.0-ti**2)
1452  ! FOR T-COORDINATE
1453  ht(1)=-0.125*rm*sm*tm+0.125*rm*sm*(2.0+ri+si+ti)
1454  ht(2)=-0.125*rp*sm*tm+0.125*rp*sm*(2.0-ri+si+ti)
1455  ht(3)=+0.25*sp*(2.0*ti-si)
1456  ht(4)=+0.125*rm*sm*tp-0.125*rm*sm*(2.0+ri+si-ti)
1457  ht(5)=+0.125*rp*sm*tp-0.125*rp*sm*(2.0-ri+si-ti)
1458  ht(6)=+0.25*sp*(si+2.0*ti)
1459  ht(7)=-0.25*(1.0-ri**2)*sm
1460  ht(8)=-0.25*rp*(1.0-si**2)
1461  ht(9)=-0.25*rm*(1.0-si**2)
1462  ht(10)=0.25*(1.0-ri**2)*sm
1463  ht(11)=0.25*rp*(1.0-si**2)
1464  ht(12)=0.25*rm*(1.0-si**2)
1465  ht(13)=-0.5*rm*sm*ti
1466  ht(14)=-0.5*rp*sm*ti
1467  ht(15)=-sp*ti
1468  ! JACOBI MATRIX
1469  xj11=0.0
1470  xj21=0.0
1471  xj31=0.0
1472  xj12=0.0
1473  xj22=0.0
1474  xj32=0.0
1475  xj13=0.0
1476  xj23=0.0
1477  xj33=0.0
1478  do i=1,nn
1479  xj11=xj11+hr(i)*xx(i)
1480  xj21=xj21+hs(i)*xx(i)
1481  xj31=xj31+ht(i)*xx(i)
1482  xj12=xj12+hr(i)*yy(i)
1483  xj22=xj22+hs(i)*yy(i)
1484  xj32=xj32+ht(i)*yy(i)
1485  xj13=xj13+hr(i)*zz(i)
1486  xj23=xj23+hs(i)*zz(i)
1487  xj33=xj33+ht(i)*zz(i)
1488  enddo
1489  !DETERMINANT OF JACOBIAN
1490  det=xj11*xj22*xj33 &
1491  +xj12*xj23*xj31 &
1492  +xj13*xj21*xj32 &
1493  -xj13*xj22*xj31 &
1494  -xj12*xj21*xj33 &
1495  -xj11*xj23*xj32
1496  wg=det*wgt(lx)*wgt(ly)*wgt(lz)
1497  do i=1,nn
1498  vect(i)=vect(i) - val*wg*h(i)
1499  enddo
1500  enddo
1501  enddo
1502  enddo
1503  endif
1504  return
1505 
1506  end subroutine heat_dflux_352
1507  !----------------------------------------------------------------------*
1508  subroutine heat_dflux_361(NN,XX,YY,ZZ,LTYPE,val,VECT)
1509  !----------------------------------------------------------------------*
1510  !**
1511  !** SET 361 DFLUX
1512  !**
1513  ! BF LTYPE=0 :BODY FLUX
1514  ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
1515  ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
1516  ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
1517  ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
1518  ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
1519  ! S6 LTYPE=6 :FLUX IN NORMAL-DIRECTION FOR FACE-6
1520  use hecmw
1521  implicit none
1522  ! I/F VARIABLES
1523  integer(kind=kint) :: NN, LTYPE
1524  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
1525  ! LOCAL VARIABLES
1526  real(kind=kreal) :: h(8), hr(8), hs(8), ht(8), pl(8)
1527  real(kind=kreal) :: xg(2), wgt(2)
1528  real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
1529  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33, det, wg
1530  integer(kind=kint) :: IVOL, ISUF
1531  integer(kind=kint) :: NOD(4)
1532  integer(kind=kint) :: IG1, IG2, LX, LY, LZ, I
1533  real(kind=kreal) :: vx, vy, vz
1534  real(kind=kreal) :: g1x, g1y, g1z
1535  real(kind=kreal) :: g2x, g2y, g2z
1536  real(kind=kreal) :: g3x, g3y, g3z
1537  real(kind=kreal) :: xsum
1538  !*************************
1539  ! GAUSS INTEGRATION POINT
1540  !*************************
1541  data xg/-0.5773502691896,0.5773502691896/
1542  data wgt/1.0,1.0/
1543  !
1544  ! SELECTION OF LOAD TYPE
1545  !
1546  ivol=0
1547  isuf=0
1548  if( ltype.EQ.0 ) then
1549  ivol=1
1550  else
1551  isuf=1
1552  if( ltype.EQ.1 ) then
1553  nod(1) = 1
1554  nod(2) = 2
1555  nod(3) = 3
1556  nod(4) = 4
1557  else if( ltype.EQ.2 ) then
1558  nod(1) = 8
1559  nod(2) = 7
1560  nod(3) = 6
1561  nod(4) = 5
1562  else if( ltype.EQ.3 ) then
1563  nod(1) = 5
1564  nod(2) = 6
1565  nod(3) = 2
1566  nod(4) = 1
1567  else if( ltype.EQ.4 ) then
1568  nod(1) = 6
1569  nod(2) = 7
1570  nod(3) = 3
1571  nod(4) = 2
1572  else if( ltype.EQ.5 ) then
1573  nod(1) = 7
1574  nod(2) = 8
1575  nod(3) = 4
1576  nod(4) = 3
1577  else if( ltype.EQ.6 ) then
1578  nod(1) = 8
1579  nod(2) = 5
1580  nod(3) = 1
1581  nod(4) = 4
1582  endif
1583  endif
1584  ! CLEAR VECT
1585  do i=1,nn
1586  vect(i)=0.0
1587  enddo
1588  !
1589  !** SURFACE LOAD
1590  !
1591  if( isuf.EQ.1 ) then
1592  ! INTEGRATION OVER SURFACE
1593  do ig2=1,2
1594  si=xg(ig2)
1595  do ig1=1,2
1596  ri=xg(ig1)
1597  h(1)=0.25*(1.0-ri)*(1.0-si)
1598  h(2)=0.25*(1.0+ri)*(1.0-si)
1599  h(3)=0.25*(1.0+ri)*(1.0+si)
1600  h(4)=0.25*(1.0-ri)*(1.0+si)
1601  hr(1)=-.25*(1.0-si)
1602  hr(2)= .25*(1.0-si)
1603  hr(3)= .25*(1.0+si)
1604  hr(4)=-.25*(1.0+si)
1605  hs(1)=-.25*(1.0-ri)
1606  hs(2)=-.25*(1.0+ri)
1607  hs(3)= .25*(1.0+ri)
1608  hs(4)= .25*(1.0-ri)
1609  g1x=0.0
1610  g1y=0.0
1611  g1z=0.0
1612  g2x=0.0
1613  g2y=0.0
1614  g2z=0.0
1615  do i=1,4
1616  g1x=g1x+hr(i)*xx(nod(i))
1617  g1y=g1y+hr(i)*yy(nod(i))
1618  g1z=g1z+hr(i)*zz(nod(i))
1619  g2x=g2x+hs(i)*xx(nod(i))
1620  g2y=g2y+hs(i)*yy(nod(i))
1621  g2z=g2z+hs(i)*zz(nod(i))
1622  enddo
1623  g3x=g1y*g2z-g1z*g2y
1624  g3y=g1z*g2x-g1x*g2z
1625  g3z=g1x*g2y-g1y*g2x
1626  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1627  do i=1,4
1628  vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1629  enddo
1630  enddo
1631  enddo
1632  endif
1633  !
1634  !** VOLUME LOAD
1635  !
1636  if( ivol.EQ.1 ) then
1637  do i=1,nn
1638  pl(i)=0.0
1639  enddo
1640  ! LOOP FOR INTEGRATION POINTS
1641  do lx=1,2
1642  ri=xg(lx)
1643  do ly=1,2
1644  si=xg(ly)
1645  do lz=1,2
1646  ti=xg(lz)
1647  rp=1.0+ri
1648  sp=1.0+si
1649  tp=1.0+ti
1650  rm=1.0-ri
1651  sm=1.0-si
1652  tm=1.0-ti
1653  ! INTERPOLATION FUNCTION
1654  h(1)=0.125*rm*sm*tm
1655  h(2)=0.125*rp*sm*tm
1656  h(3)=0.125*rp*sp*tm
1657  h(4)=0.125*rm*sp*tm
1658  h(5)=0.125*rm*sm*tp
1659  h(6)=0.125*rp*sm*tp
1660  h(7)=0.125*rp*sp*tp
1661  h(8)=0.125*rm*sp*tp
1662  ! DERIVATIVE OF INTERPOLATION FUNCTION
1663  ! FOR R-COORDINATE
1664  hr(1)=-.125*sm*tm
1665  hr(2)= .125*sm*tm
1666  hr(3)= .125*sp*tm
1667  hr(4)=-.125*sp*tm
1668  hr(5)=-.125*sm*tp
1669  hr(6)= .125*sm*tp
1670  hr(7)= .125*sp*tp
1671  hr(8)=-.125*sp*tp
1672  ! FOR S-COORDINATE
1673  hs(1)=-.125*rm*tm
1674  hs(2)=-.125*rp*tm
1675  hs(3)= .125*rp*tm
1676  hs(4)= .125*rm*tm
1677  hs(5)=-.125*rm*tp
1678  hs(6)=-.125*rp*tp
1679  hs(7)= .125*rp*tp
1680  hs(8)= .125*rm*tp
1681  ! FOR T-COORDINATE
1682  ht(1)=-.125*rm*sm
1683  ht(2)=-.125*rp*sm
1684  ht(3)=-.125*rp*sp
1685  ht(4)=-.125*rm*sp
1686  ht(5)= .125*rm*sm
1687  ht(6)= .125*rp*sm
1688  ht(7)= .125*rp*sp
1689  ht(8)= .125*rm*sp
1690  ! JACOBI MATRIX
1691  xj11=0.0
1692  xj21=0.0
1693  xj31=0.0
1694  xj12=0.0
1695  xj22=0.0
1696  xj32=0.0
1697  xj13=0.0
1698  xj23=0.0
1699  xj33=0.0
1700  do i=1,nn
1701  xj11=xj11+hr(i)*xx(i)
1702  xj21=xj21+hs(i)*xx(i)
1703  xj31=xj31+ht(i)*xx(i)
1704  xj12=xj12+hr(i)*yy(i)
1705  xj22=xj22+hs(i)*yy(i)
1706  xj32=xj32+ht(i)*yy(i)
1707  xj13=xj13+hr(i)*zz(i)
1708  xj23=xj23+hs(i)*zz(i)
1709  xj33=xj33+ht(i)*zz(i)
1710  enddo
1711  !DETERMINANT OF JACOBIAN
1712  det=xj11*xj22*xj33 &
1713  +xj12*xj23*xj31 &
1714  +xj13*xj21*xj32 &
1715  -xj13*xj22*xj31 &
1716  -xj12*xj21*xj33 &
1717  -xj11*xj23*xj32
1718  wg=wgt(lx)*wgt(ly)*wgt(lz)*det
1719  do i=1,nn
1720  vect(i)=vect(i)-h(i)*wg*val
1721  enddo
1722  enddo
1723  enddo
1724  enddo
1725  endif
1726  !*
1727  return
1728 
1729  end subroutine heat_dflux_361
1730  !----------------------------------------------------------------------*
1731  subroutine heat_dflux_362(NN,XX,YY,ZZ,LTYPE,val,VECT)
1732  !----------------------------------------------------------------------*
1733  !**
1734  !** SET 362 DFLUX
1735  !**
1736  ! BF LTYPE=0 :BODY FLUX
1737  ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
1738  ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
1739  ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
1740  ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
1741  ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
1742  ! S6 LTYPE=6 :FLUX IN NORMAL-DIRECTION FOR FACE-6
1743  use hecmw
1744  implicit none
1745  ! I/F VARIABLES
1746  integer(kind=kint) :: NN, LTYPE
1747  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
1748  ! LOCAL VARIABLES
1749  real(kind=kreal) :: h(20), hr(20), hs(20), ht(20)
1750  real(kind=kreal) :: xg(3), wgt(3)
1751  real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
1752  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
1753  real(kind=kreal) :: det, wg
1754  integer(kind=kint) :: IVOL, ISUF
1755  integer(kind=kint) :: NOD(8)
1756  integer(kind=kint) :: IG1, IG2, LX, LY, LZ, I
1757  real(kind=kreal) :: vx, vy, vz
1758  real(kind=kreal) :: g1x, g1y, g1z
1759  real(kind=kreal) :: g2x, g2y, g2z
1760  real(kind=kreal) :: g3x, g3y, g3z
1761  real(kind=kreal) :: xsum
1762  !*************************
1763  ! GAUSS INTEGRATION POINT
1764  !*************************
1765  xg(1) = -0.7745966692
1766  xg(2) = 0.0
1767  xg(3) = 0.7745966692
1768  wgt(1) = 0.5555555555
1769  wgt(2) = 0.8888888888
1770  wgt(3) = 0.5555555555
1771  !
1772  ! SELECTION OF LOAD TYPE
1773  !
1774  ivol=0
1775  isuf=0
1776  if( ltype.EQ.0 ) then
1777  ivol=1
1778  else
1779  isuf=1
1780  if( ltype.EQ.1 ) then
1781  nod(1) = 1
1782  nod(2) = 2
1783  nod(3) = 3
1784  nod(4) = 4
1785  nod(5) = 9
1786  nod(6) = 10
1787  nod(7) = 11
1788  nod(8) = 12
1789  else if( ltype.EQ.2 ) then
1790  nod(1) = 8
1791  nod(2) = 7
1792  nod(3) = 6
1793  nod(4) = 5
1794  nod(5) = 15
1795  nod(6) = 14
1796  nod(7) = 13
1797  nod(8) = 16
1798  else if( ltype.EQ.3 ) then
1799  nod(1) = 5
1800  nod(2) = 6
1801  nod(3) = 2
1802  nod(4) = 1
1803  nod(5) = 13
1804  nod(6) = 18
1805  nod(7) = 9
1806  nod(8) = 17
1807  else if( ltype.EQ.4 ) then
1808  nod(1) = 6
1809  nod(2) = 7
1810  nod(3) = 3
1811  nod(4) = 2
1812  nod(5) = 14
1813  nod(6) = 19
1814  nod(7) = 10
1815  nod(8) = 18
1816  else if( ltype.EQ.5 ) then
1817  nod(1) = 7
1818  nod(2) = 8
1819  nod(3) = 4
1820  nod(4) = 3
1821  nod(5) = 15
1822  nod(6) = 20
1823  nod(7) = 11
1824  nod(8) = 19
1825  else if( ltype.EQ.6 ) then
1826  nod(1) = 8
1827  nod(2) = 5
1828  nod(3) = 1
1829  nod(4) = 4
1830  nod(5) = 16
1831  nod(6) = 17
1832  nod(7) = 12
1833  nod(8) = 20
1834  endif
1835  endif
1836  ! CLEAR VECT
1837  do i=1,nn
1838  vect(i)=0.0
1839  enddo
1840  !
1841  !** SURFACE LOAD
1842  !
1843  if( isuf.EQ.1 ) then
1844  ! INTEGRATION OVER SURFACE
1845  do ig2=1,3
1846  si=xg(ig2)
1847  do ig1=1,3
1848  ri=xg(ig1)
1849  rp=1.0+ri
1850  sp=1.0+si
1851  rm=1.0-ri
1852  sm=1.0-si
1853  h(1)=0.25*rm*sm*(-1.0-ri-si)
1854  h(2)=0.25*rp*sm*(-1.0+ri-si)
1855  h(3)=0.25*rp*sp*(-1.0+ri+si)
1856  h(4)=0.25*rm*sp*(-1.0-ri+si)
1857  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
1858  h(6)=0.5*(1.0-si*si)*(1.0+ri)
1859  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
1860  h(8)=0.5*(1.0-si*si)*(1.0-ri)
1861  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1862  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1863  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
1864  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
1865  hr(5)=-ri*(1.0-si)
1866  hr(6)= .5*(1.0-si*si)
1867  hr(7)=-ri*(1.0+si)
1868  hr(8)=-.5*(1.0-si*si)
1869  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1870  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1871  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
1872  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
1873  hs(5)=-.5*(1.0-ri*ri)
1874  hs(6)=-si*(1.0+ri)
1875  hs(7)= .5*(1.0-ri*ri)
1876  hs(8)=-si*(1.0-ri)
1877  g1x=0.0
1878  g1y=0.0
1879  g1z=0.0
1880  g2x=0.0
1881  g2y=0.0
1882  g2z=0.0
1883  do i=1,8
1884  g1x=g1x+hr(i)*xx(nod(i))
1885  g1y=g1y+hr(i)*yy(nod(i))
1886  g1z=g1z+hr(i)*zz(nod(i))
1887  g2x=g2x+hs(i)*xx(nod(i))
1888  g2y=g2y+hs(i)*yy(nod(i))
1889  g2z=g2z+hs(i)*zz(nod(i))
1890  enddo
1891  g3x=g1y*g2z-g1z*g2y
1892  g3y=g1z*g2x-g1x*g2z
1893  g3z=g1x*g2y-g1y*g2x
1894  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1895  do i=1,8
1896  vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1897  enddo
1898  enddo
1899  enddo
1900  endif
1901  !
1902  !** VOLUME LOAD
1903  !
1904  if( ivol.EQ.1 ) then
1905  ! LOOP FOR INTEGRATION POINTS
1906  do lx=1,3
1907  ri=xg(lx)
1908  do ly=1,3
1909  si=xg(ly)
1910  do lz=1,3
1911  ti=xg(lz)
1912  rp=1.0+ri
1913  sp=1.0+si
1914  tp=1.0+ti
1915  rm=1.0-ri
1916  sm=1.0-si
1917  tm=1.0-ti
1918  ! INTERPOLATION FUNCTION
1919  h(1)=-0.125*rm*sm*tm*(2.0+ri+si+ti)
1920  h(2)=-0.125*rp*sm*tm*(2.0-ri+si+ti)
1921  h(3)=-0.125*rp*sp*tm*(2.0-ri-si+ti)
1922  h(4)=-0.125*rm*sp*tm*(2.0+ri-si+ti)
1923  h(5)=-0.125*rm*sm*tp*(2.0+ri+si-ti)
1924  h(6)=-0.125*rp*sm*tp*(2.0-ri+si-ti)
1925  h(7)=-0.125*rp*sp*tp*(2.0-ri-si-ti)
1926  h(8)=-0.125*rm*sp*tp*(2.0+ri-si-ti)
1927  h(9)=0.25*(1.0-ri**2)*sm*tm
1928  h(10)=0.25*rp*(1.0-si**2)*tm
1929  h(11)=0.25*(1.0-ri**2)*sp*tm
1930  h(12)=0.25*rm*(1.0-si**2)*tm
1931  h(13)=0.25*(1.0-ri**2)*sm*tp
1932  h(14)=0.25*rp*(1.0-si**2)*tp
1933  h(15)=0.25*(1.0-ri**2)*sp*tp
1934  h(16)=0.25*rm*(1.0-si**2)*tp
1935  h(17)=0.25*rm*sm*(1.0-ti**2)
1936  h(18)=0.25*rp*sm*(1.0-ti**2)
1937  h(19)=0.25*rp*sp*(1.0-ti**2)
1938  h(20)=0.25*rm*sp*(1.0-ti**2)
1939  ! DERIVATIVE OF INTERPOLATION FUNCTION
1940  ! FOR R-COORDINATE
1941  hr(1)=-0.125*rm*sm*tm+0.125*sm*tm*(2.0+ri+si+ti)
1942  hr(2)=+0.125*rp*sm*tm-0.125*sm*tm*(2.0-ri+si+ti)
1943  hr(3)=+0.125*rp*sp*tm-0.125*sp*tm*(2.0-ri-si+ti)
1944  hr(4)=-0.125*rm*sp*tm+0.125*sp*tm*(2.0+ri-si+ti)
1945  hr(5)=-0.125*rm*sm*tp+0.125*sm*tp*(2.0+ri+si-ti)
1946  hr(6)=+0.125*rp*sm*tp-0.125*sm*tp*(2.0-ri+si-ti)
1947  hr(7)=+0.125*rp*sp*tp-0.125*sp*tp*(2.0-ri-si-ti)
1948  hr(8)=-0.125*rm*sp*tp+0.125*sp*tp*(2.0+ri-si-ti)
1949  hr(9 )=-0.50*ri*sm*tm
1950  hr(10)=+0.25*(1.0-si**2)*tm
1951  hr(11)=-0.50*ri*sp*tm
1952  hr(12)=-0.25*(1.0-si**2)*tm
1953  hr(13)=-0.50*ri*sm*tp
1954  hr(14)=+0.25*(1.0-si**2)*tp
1955  hr(15)=-0.50*ri*sp*tp
1956  hr(16)=-0.25*(1.0-si**2)*tp
1957  hr(17)=-0.25*sm*(1.0-ti**2)
1958  hr(18)=+0.25*sm*(1.0-ti**2)
1959  hr(19)=+0.25*sp*(1.0-ti**2)
1960  hr(20)=-0.25*sp*(1.0-ti**2)
1961  ! FOR S-COORDINATE
1962  hs(1)=-0.125*rm*sm*tm+0.125*rm*tm*(2.0+ri+si+ti)
1963  hs(2)=-0.125*rp*sm*tm+0.125*rp*tm*(2.0-ri+si+ti)
1964  hs(3)=+0.125*rp*sp*tm-0.125*rp*tm*(2.0-ri-si+ti)
1965  hs(4)=+0.125*rm*sp*tm-0.125*rm*tm*(2.0+ri-si+ti)
1966  hs(5)=-0.125*rm*sm*tp+0.125*rm*tp*(2.0+ri+si-ti)
1967  hs(6)=-0.125*rp*sm*tp+0.125*rp*tp*(2.0-ri+si-ti)
1968  hs(7)=+0.125*rp*sp*tp-0.125*rp*tp*(2.0-ri-si-ti)
1969  hs(8)=+0.125*rm*sp*tp-0.125*rm*tp*(2.0+ri-si-ti)
1970  hs(9)=-0.25*(1.0-ri**2)*tm
1971  hs(10)=-0.50*rp*si*tm
1972  hs(11)=+0.25*(1.0-ri**2)*tm
1973  hs(12)=-0.50*rm*si*tm
1974  hs(13)=-0.25*(1.0-ri**2)*tp
1975  hs(14)=-0.50*rp*si*tp
1976  hs(15)=+0.25*(1.0-ri**2)*tp
1977  hs(16)=-0.50*rm*si*tp
1978  hs(17)=-0.25*rm*(1.0-ti**2)
1979  hs(18)=-0.25*rp*(1.0-ti**2)
1980  hs(19)=+0.25*rp*(1.0-ti**2)
1981  hs(20)=+0.25*rm*(1.0-ti**2)
1982  ! FOR T-COORDINATE
1983  ht(1)=-0.125*rm*sm*tm+0.125*rm*sm*(2.0+ri+si+ti)
1984  ht(2)=-0.125*rp*sm*tm+0.125*rp*sm*(2.0-ri+si+ti)
1985  ht(3)=-0.125*rp*sp*tm+0.125*rp*sp*(2.0-ri-si+ti)
1986  ht(4)=-0.125*rm*sp*tm+0.125*rm*sp*(2.0+ri-si+ti)
1987  ht(5)=+0.125*rm*sm*tp-0.125*rm*sm*(2.0+ri+si-ti)
1988  ht(6)=+0.125*rp*sm*tp-0.125*rp*sm*(2.0-ri+si-ti)
1989  ht(7)=+0.125*rp*sp*tp-0.125*rp*sp*(2.0-ri-si-ti)
1990  ht(8)=+0.125*rm*sp*tp-0.125*rm*sp*(2.0+ri-si-ti)
1991  ht(9)=-0.25*(1.0-ri**2)*sm
1992  ht(10)=-0.25*rp*(1.0-si**2)
1993  ht(11)=-0.25*(1.0-ri**2)*sp
1994  ht(12)=-0.25*rm*(1.0-si**2)
1995  ht(13)=0.25*(1.0-ri**2)*sm
1996  ht(14)=0.25*rp*(1.0-si**2)
1997  ht(15)=0.25*(1.0-ri**2)*sp
1998  ht(16)=0.25*rm*(1.0-si**2)
1999  ht(17)=-0.5*rm*sm*ti
2000  ht(18)=-0.5*rp*sm*ti
2001  ht(19)=-0.5*rp*sp*ti
2002  ht(20)=-0.5*rm*sp*ti
2003  ! JACOBI MATRIX
2004  xj11=0.0
2005  xj21=0.0
2006  xj31=0.0
2007  xj12=0.0
2008  xj22=0.0
2009  xj32=0.0
2010  xj13=0.0
2011  xj23=0.0
2012  xj33=0.0
2013  do i=1,nn
2014  xj11=xj11+hr(i)*xx(i)
2015  xj21=xj21+hs(i)*xx(i)
2016  xj31=xj31+ht(i)*xx(i)
2017  xj12=xj12+hr(i)*yy(i)
2018  xj22=xj22+hs(i)*yy(i)
2019  xj32=xj32+ht(i)*yy(i)
2020  xj13=xj13+hr(i)*zz(i)
2021  xj23=xj23+hs(i)*zz(i)
2022  xj33=xj33+ht(i)*zz(i)
2023  enddo
2024  !DETERMINANT OF JACOBIAN
2025  det=xj11*xj22*xj33 &
2026  +xj12*xj23*xj31 &
2027  +xj13*xj21*xj32 &
2028  -xj13*xj22*xj31 &
2029  -xj12*xj21*xj33 &
2030  -xj11*xj23*xj32
2031  wg=wgt(lx)*wgt(ly)*wgt(lz)*det
2032  do i=1,nn
2033  vect(i)=vect(i)-h(i)*wg*val
2034  enddo
2035  enddo
2036  enddo
2037  enddo
2038  endif
2039  !*
2040  return
2041 
2042  end subroutine heat_dflux_362
2043  !----------------------------------------------------------------------*
2044  subroutine heat_dflux_731(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
2045  !----------------------------------------------------------------------*
2046  !**
2047  !** SET S3 DFLUX
2048  !**
2049  ! BF LTYPE=0 :BODY FLUX
2050  ! S1 LTYPE=1 :SURFACE FLUX
2051  use hecmw
2052  implicit none
2053  ! I/F VARIABLES
2054  integer(kind=kint) :: NN, LTYPE
2055  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
2056  ! LOCAL VARIABLES
2057  integer(kind=kint) :: I
2058  real(kind=kreal) :: v1x, v1y, v1z, v2x, v2y, v2z, v3x, v3y, v3z, aa
2059  !
2060  if( ltype.EQ.0 ) then
2061  thick = thick
2062  elseif( ltype.EQ.1 ) then
2063  thick = 1.0d0
2064  else
2065  thick = 0.0d0
2066  endif
2067 
2068  do i=1,nn
2069  vect(i)=0.0
2070  enddo
2071  !
2072  v1x=xx(2)-xx(1)
2073  v1y=yy(2)-yy(1)
2074  v1z=zz(2)-zz(1)
2075  v2x=xx(3)-xx(1)
2076  v2y=yy(3)-yy(1)
2077  v2z=zz(3)-zz(1)
2078  v3x= v1y*v2z-v1z*v2y
2079  v3y= v1z*v2x-v1x*v2z
2080  v3z= v1x*v2y-v1y*v2x
2081  aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
2082  vect(1)= -val*aa*thick/3.0
2083  vect(2)= -val*aa*thick/3.0
2084  vect(3)= -val*aa*thick/3.0
2085  !
2086  return
2087 
2088  end subroutine heat_dflux_731
2089  !----------------------------------------------------------------------*
2090  subroutine heat_dflux_741(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
2091  !----------------------------------------------------------------------*
2092  !**
2093  !** SET S4 DFLUX
2094  !**
2095  ! BF LTYPE=0 : BODY FLUX
2096  ! S1 LTYPE=1 : SURFACE FLUX
2097  use hecmw
2098  implicit none
2099  ! I/F VARIABLES
2100  integer(kind=kint) :: NN, LTYPE
2101  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
2102  ! LOCAL VARIABLES
2103  integer(kind=kint) :: I, IG1, IG2
2104  real(kind=kreal) :: ri, si, xsum
2105  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
2106  real(kind=kreal) :: h(4), hr(4), hs(4), ht(4), pl(4)
2107  real(kind=kreal) :: xg(2), wgt(2)
2108  !*************************
2109  ! GAUSS INTEGRATION POINT
2110  !*************************
2111  data xg / -0.5773502691896,0.5773502691896 /
2112  data wgt / 1.0, 1.0 /
2113  !
2114  ! SELECTION OF LOAD TYPE
2115  !
2116  if ( ltype.EQ.0 ) then
2117  thick = 1.0d0 * thick
2118  elseif( ltype.EQ.1 ) then
2119  thick = 1.0d0
2120  else
2121  thick = 0.0d0
2122  endif
2123  !
2124  do i=1,nn
2125  vect(i)=0.0
2126  enddo
2127  !
2128  do ig2=1,2
2129  si=xg(ig2)
2130  do ig1=1,2
2131  ri=xg(ig1)
2132  !
2133  h(1)=0.25*(1.0-ri)*(1.0-si)
2134  h(2)=0.25*(1.0+ri)*(1.0-si)
2135  h(3)=0.25*(1.0+ri)*(1.0+si)
2136  h(4)=0.25*(1.0-ri)*(1.0+si)
2137  hr(1)=-.25*(1.0-si)
2138  hr(2)= .25*(1.0-si)
2139  hr(3)= .25*(1.0+si)
2140  hr(4)=-.25*(1.0+si)
2141  hs(1)=-.25*(1.0-ri)
2142  hs(2)=-.25*(1.0+ri)
2143  hs(3)= .25*(1.0+ri)
2144  hs(4)= .25*(1.0-ri)
2145  !
2146  g1x=0.0
2147  g1y=0.0
2148  g1z=0.0
2149  g2x=0.0
2150  g2y=0.0
2151  g2z=0.0
2152  do i=1,nn
2153  g1x=g1x+hr(i)*xx(i)
2154  g1y=g1y+hr(i)*yy(i)
2155  g1z=g1z+hr(i)*zz(i)
2156  g2x=g2x+hs(i)*xx(i)
2157  g2y=g2y+hs(i)*yy(i)
2158  g2z=g2z+hs(i)*zz(i)
2159  enddo
2160  !
2161  g3x=g1y*g2z-g1z*g2y
2162  g3y=g1z*g2x-g1x*g2z
2163  g3z=g1x*g2y-g1y*g2x
2164  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
2165  do i=1,nn
2166  vect(i)=vect(i)-xsum*val*wgt(ig1)*wgt(ig2)*h(i)*thick
2167  enddo
2168  !
2169  enddo
2170  enddo
2171  return
2172 
2173  end subroutine heat_dflux_741
2174 end module m_heat_lib_dflux
m_heat_lib_dflux::heat_dflux_351
subroutine heat_dflux_351(NN, XX, YY, ZZ, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:962
m_heat_lib_dflux
This module provides subroutines for calculating distributed heat flux for various elements.
Definition: heat_LIB_DFLUX.f90:7
m_heat_lib_dflux::heat_dflux_242
subroutine heat_dflux_242(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:405
m_heat_lib_dflux::heat_dflux_241
subroutine heat_dflux_241(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:293
m_heat_lib_dflux::heat_dflux_741
subroutine heat_dflux_741(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:2091
m_heat_lib_dflux::heat_dflux_342
subroutine heat_dflux_342(NN, XX, YY, ZZ, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:683
m_heat_lib_dflux::heat_dflux_111
subroutine heat_dflux_111(NN, XX, YY, ZZ, ASECT, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:25
m_heat_lib_dflux::heat_dflux_341
subroutine heat_dflux_341(NN, XX, YY, ZZ, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:539
m_heat_lib_dflux::heat_dflux_362
subroutine heat_dflux_362(NN, XX, YY, ZZ, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:1732
m_heat_lib_dflux::heat_dflux_231
subroutine heat_dflux_231(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:64
m_heat_lib_dflux::heat_dflux_361
subroutine heat_dflux_361(NN, XX, YY, ZZ, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:1509
m_heat_lib_dflux::heat_dflux_731
subroutine heat_dflux_731(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:2045
hecmw
Definition: hecmw.f90:6
m_heat_lib_dflux::heat_dflux_352
subroutine heat_dflux_352(NN, XX, YY, ZZ, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:1175
m_heat_lib_dflux::heat_dflux_232
subroutine heat_dflux_232(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
Definition: heat_LIB_DFLUX.f90:155