FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
heat_LIB_FILM.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 !-------------------------------------------------------------------------------
8 contains
9  !***********************************************************************
10  ! FILM_231(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
11  ! FILM_232(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
12  ! FILM_241(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
13  ! FILM_242(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
14  ! FILM_341(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
15  ! FILM_342(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
16  ! FILM_351(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
17  ! FILM_352(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
18  ! FILM_361(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
19  ! FILM_362(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
20  ! FILM_731(NN,XX,YY,ZZ,LTYPE,HH,SINK,TERM1,TERM2)
21  ! FILM_741(NN,XX,YY,ZZ,LTYPE,HH,SINK,TERM1,TERM2)
22  !----------------------------------------------------------------------*
23  subroutine heat_film_231(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
24  !----------------------------------------------------------------------*
25  !**
26  !** SET 231 FILM
27  !**
28  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
29  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
30  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
31  use hecmw
32  implicit real(kind=kreal)(a-h,o-z)
33  ! I/F VARIABLES
34  integer(kind=kint) NOD(MM)
35  real(kind=kreal) xx(nn),yy(nn),zz(nn),term1(mm*mm),term2(mm)
36  real(kind=kreal) xg(2),wgt(2),h(2),hr(2)
37  data wgt/1.0,1.0/
38  data xg/-0.5773502691896, 0.5773502691896/
39  !
40  if(ltype.EQ.1) then
41  nod(1)=1
42  nod(2)=2
43  else if(ltype.EQ.2) then
44  nod(1)=2
45  nod(2)=3
46  else if(ltype.EQ.3) then
47  nod(1)=3
48  nod(2)=1
49  endif
50  !
51  ic = 0
52  do ip = 1, 2
53  term2(ip) = 0.0
54  do jp = 1, 2
55  ic = ic + 1
56  term1(ic) = 0.0
57  !** INTEGRATION OVER SURFACE
58  do lx=1,2
59  ri=xg(lx)
60  h(1)=0.5*(1.0-ri)
61  h(2)=0.5*(1.0+ri)
62  hr(1)=-0.5
63  hr(2)= 0.5
64  gx=0.0
65  gy=0.0
66  do i=1,2
67  gx=gx+hr(i)*xx(nod(i))
68  gy=gy+hr(i)*yy(nod(i))
69  enddo
70  xsum = dsqrt( gx*gx+gy*gy )*thick
71  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
72  if( ip.EQ.jp ) then
73  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
74  endif
75  enddo
76  enddo
77  enddo
78  !
79  return
80 
81  end subroutine heat_film_231
82  !----------------------------------------------------------------------*
83  subroutine heat_film_232(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
84  !----------------------------------------------------------------------*
85  !**
86  !** SET 232 FILM
87  !**
88  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
89  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
90  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
91  use hecmw
92  implicit real(kind=kreal) (a - h, o - z)
93  ! I/F VARIABLES
94  integer(kind=kint) :: NOD(MM)
95  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
96  real(kind=kreal) :: xg(3), wgt(3), h(3), hr(3)
97  !*************************
98  ! GAUSS INTEGRATION POINT
99  !*************************
100  xg(1) = -0.7745966692
101  xg(2) = 0.0
102  xg(3) = 0.7745966692
103  wgt(1) = 0.5555555555
104  wgt(2) = 0.8888888888
105  wgt(3) = 0.5555555555
106 
107  if(ltype.EQ.1) then
108  nod(1)=1
109  nod(2)=2
110  nod(3)=4
111  else if(ltype.EQ.2) then
112  nod(1)=2
113  nod(2)=3
114  nod(3)=5
115  else if(ltype.EQ.3) then
116  nod(1)=3
117  nod(2)=1
118  nod(3)=6
119  endif
120 
121  ic = 0
122  do ip = 1, 3
123  term2(ip) = 0.0
124  do jp = 1, 3
125  ic = ic + 1
126  term1(ic) = 0.0
127  !** INTEGRATION OVER SURFACE
128  do lx=1,3
129  ri=xg(lx)
130  h(1)=-0.5*(1.0-ri)*ri
131  h(2)= 0.5*(1.0+ri)*ri
132  h(3)=1.0-ri*ri
133  hr(1)= ri-0.5
134  hr(2)= ri+0.5
135  hr(3)=-2.0*ri
136  gx=0.0
137  gy=0.0
138  do i=1,3
139  gx=gx+hr(i)*xx(nod(i))
140  gy=gy+hr(i)*yy(nod(i))
141  enddo
142  xsum = dsqrt( gx*gx+gy*gy )*thick
143  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
144  if( ip.EQ.jp ) then
145  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
146  endif
147  enddo
148  enddo
149  enddo
150 
151  return
152 
153  end subroutine heat_film_232
154  !----------------------------------------------------------------------*
155  subroutine heat_film_241(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
156  !----------------------------------------------------------------------*
157  !**
158  !** SET 241 FILM
159  !**
160  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
161  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
162  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
163  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
164  use hecmw
165  implicit real(kind=kreal) (a - h, o - z)
166  ! I/F VARIABLES
167  integer(kind=kint) :: NOD(MM)
168  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
169  real(kind=kreal) :: xg(2), wgt(2), h(2), hr(2)
170  data wgt/1.0, 1.0/
171  data xg/-0.5773502691896, 0.5773502691896/
172 
173  if(ltype.EQ.1) then
174  nod(1)=1
175  nod(2)=2
176  else if(ltype.EQ.2) then
177  nod(1)=2
178  nod(2)=3
179  else if(ltype.EQ.3) then
180  nod(1)=3
181  nod(2)=4
182  else if(ltype.EQ.4) then
183  nod(1)=4
184  nod(2)=1
185  endif
186 
187  ic = 0
188  do ip = 1, 2
189  term2(ip) = 0.0
190  do jp = 1, 2
191  ic = ic + 1
192  term1(ic) = 0.0
193  !** INTEGRATION OVER SURFACE
194  do lx=1,2
195  ri=xg(lx)
196  h(1)=0.5*(1.0-ri)
197  h(2)=0.5*(1.0+ri)
198  hr(1)=-0.5
199  hr(2)= 0.5
200  gx=0.0
201  gy=0.0
202  do i=1,2
203  gx=gx+hr(i)*xx(nod(i))
204  gy=gy+hr(i)*yy(nod(i))
205  enddo
206  xsum = dsqrt( gx*gx+gy*gy )*thick
207  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
208  if( ip.EQ.jp ) then
209  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
210  endif
211  enddo
212  enddo
213  enddo
214 
215  return
216 
217  end subroutine heat_film_241
218  !----------------------------------------------------------------------*
219  subroutine heat_film_242(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
220  !----------------------------------------------------------------------*
221  !**
222  !** SET 242 FILM
223  !**
224  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
225  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
226  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
227  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
228  use hecmw
229  implicit real(kind=kreal) (a - h, o - z)
230  ! I/F VARIABLES
231  integer(kind=kint) :: NOD(MM)
232  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
233  real(kind=kreal) :: xg(3), wgt(3), h(3), hr(3)
234  !*************************
235  ! GAUSS INTEGRATION POINT
236  !*************************
237  xg(1) = -0.7745966692
238  xg(2) = 0.0
239  xg(3) = 0.7745966692
240  wgt(1) = 0.5555555555
241  wgt(2) = 0.8888888888
242  wgt(3) = 0.5555555555
243 
244  if(ltype.EQ.1) then
245  nod(1)=1
246  nod(2)=2
247  nod(3)=5
248  else if(ltype.EQ.2) then
249  nod(1)=2
250  nod(2)=3
251  nod(3)=6
252  else if(ltype.EQ.3) then
253  nod(1)=3
254  nod(2)=4
255  nod(3)=7
256  else if(ltype.EQ.4) then
257  nod(1)=4
258  nod(2)=1
259  nod(3)=8
260  endif
261 
262  ic = 0
263  do ip = 1, 3
264  term2(ip) = 0.0
265  do jp = 1, 3
266  ic = ic + 1
267  term1(ic) = 0.0
268  !** INTEGRATION OVER SURFACE
269  do lx=1,3
270  ri=xg(lx)
271  h(1)=-0.5*(1.0-ri)*ri
272  h(2)= 0.5*(1.0+ri)*ri
273  h(3)=1.0-ri*ri
274  hr(1)= ri-0.5
275  hr(2)= ri+0.5
276  hr(3)=-2.0*ri
277  gx=0.0
278  gy=0.0
279  do i=1,3
280  gx=gx+hr(i)*xx(nod(i))
281  gy=gy+hr(i)*yy(nod(i))
282  enddo
283  xsum = dsqrt( gx*gx+gy*gy )*thick
284  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
285  if( ip.EQ.jp ) then
286  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
287  endif
288  enddo
289  enddo
290  enddo
291 
292  return
293 
294  end subroutine heat_film_242
295  !----------------------------------------------------------------------*
296  subroutine heat_film_341(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
297  !----------------------------------------------------------------------*
298  !**
299  !** SET 341 FILM
300  !**
301  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
302  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
303  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
304  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
305  use hecmw
306  implicit real(kind=kreal) (a - h, o - z)
307  ! I/F VARIABLES
308  integer(kind=kint) :: NOD(MM)
309  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
310 
311  !
312  if ( ltype.EQ.1 ) then
313  nod(1) = 1
314  nod(2) = 2
315  nod(3) = 3
316  else if( ltype.EQ.2 ) then
317  nod(1) = 4
318  nod(2) = 2
319  nod(3) = 1
320  else if( ltype.EQ.3 ) then
321  nod(1) = 4
322  nod(2) = 3
323  nod(3) = 2
324  else if( ltype.EQ.4 ) then
325  nod(1) = 4
326  nod(2) = 1
327  nod(3) = 3
328  endif
329  !
330  ax = xx( nod(2) ) - xx( nod(1) )
331  ay = yy( nod(2) ) - yy( nod(1) )
332  az = zz( nod(2) ) - zz( nod(1) )
333  bx = xx( nod(3) ) - xx( nod(1) )
334  by = yy( nod(3) ) - yy( nod(1) )
335  bz = zz( nod(3) ) - zz( nod(1) )
336  !
337  aa = dsqrt( (ay*bz-az*by)**2+(az*bx-ax*bz)**2+(ax*by-ay*bx)**2 )/6.0
338  !
339  ic = 0
340  do ip = 1, 3
341  term2(ip) = - aa*hh*sink
342  do jp = 1, 3
343  ic = ic + 1
344  term1(ic) = - aa*hh/3.0
345  enddo
346  enddo
347  !
348  return
349 
350  end subroutine heat_film_341
351  !----------------------------------------------------------------------*
352  subroutine heat_film_342(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
353  !----------------------------------------------------------------------*
354  !**
355  !** SET 342 FILM
356  !**
357  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
358  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
359  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
360  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
361  use hecmw
362  implicit real(kind=kreal) (a - h, o - z)
363  ! I/F VARIABLES
364  integer(kind=kint) :: NOD(MM)
365  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
366  real(kind=kreal) :: xg(3), wgt(3), h(6), hl1(6), hl2(6), hl3(6)
367  !*************************
368  ! GAUSS INTEGRATION POINT
369  !*************************
370  xg(1) = -0.7745966692
371  xg(2) = 0.0
372  xg(3) = 0.7745966692
373  wgt(1) = 0.5555555555
374  wgt(2) = 0.8888888888
375  wgt(3) = 0.5555555555
376  !
377  if(ltype.EQ.1) then
378  nod(1)=1
379  nod(2)=2
380  nod(3)=3
381  nod(4)=5
382  nod(5)=6
383  nod(6)=7
384  else if(ltype.EQ.2) then
385  nod(1)=4
386  nod(2)=2
387  nod(3)=1
388  nod(4)=9
389  nod(5)=5
390  nod(6)=8
391  else if(ltype.EQ.3) then
392  nod(1)=4
393  nod(2)=3
394  nod(3)=2
395  nod(4)=10
396  nod(5)=6
397  nod(6)=9
398  else if(ltype.EQ.4) then
399  nod(1)=4
400  nod(2)=1
401  nod(3)=3
402  nod(4)=8
403  nod(5)=7
404  nod(6)=10
405  endif
406  !
407  ic = 0
408  do ip = 1, 6
409  term2(ip) = 0.0
410  do jp = 1, 6
411  ic = ic + 1
412  term1(ic) = 0.0
413  ! INTEGRATION OVER SURFACE
414  do l2=1,3
415  xl2=xg(l2)
416  x2 =(xl2+1.0)*0.5
417  do l1=1,3
418  xl1=xg(l1)
419  x1=0.5*(1.0-x2)*(xl1+1.0)
420  ! INTERPOLATION FUNCTION
421  x3=1.0-x1-x2
422  h(1)= x1*(2.0*x1-1.)
423  h(2)= x2*(2.0*x2-1.)
424  h(3)= x3*(2.0*x3-1.)
425  h(4)= 4.0*x1*x2
426  h(5)= 4.0*x2*x3
427  h(6)= 4.0*x1*x3
428  ! DERIVATIVE OF INTERPOLATION FUNCTION
429  ! FOR L1-COORDINATE
430  hl1(1)=4.0*x1-1.0
431  hl1(2)= 0.0
432  hl1(3)= 0.0
433  hl1(4)= 4.0*x2
434  hl1(5)= 0.0
435  hl1(6)= 4.0*x3
436  ! FOR L2-COORDINATE
437  hl2(1)= 0.0
438  hl2(2)= 4.0*x2-1.0
439  hl2(3)= 0.0
440  hl2(4)= 4.0*x1
441  hl2(5)= 4.0*x3
442  hl2(6)= 0.0
443  ! FOR L3-COORDINATE
444  hl3(1)= 0.0
445  hl3(2)= 0.0
446  hl3(3)= 4.0*x3-1.0
447  hl3(4)= 0.0
448  hl3(5)= 4.0*x2
449  hl3(6)= 4.0*x1
450  !
451  g1x=0.0
452  g1y=0.0
453  g1z=0.0
454  g2x=0.0
455  g2y=0.0
456  g2z=0.0
457  do i=1,6
458  g1x=g1x+(hl1(i)-hl3(i))*xx(nod(i))
459  g1y=g1y+(hl1(i)-hl3(i))*yy(nod(i))
460  g1z=g1z+(hl1(i)-hl3(i))*zz(nod(i))
461  g2x=g2x+(hl2(i)-hl3(i))*xx(nod(i))
462  g2y=g2y+(hl2(i)-hl3(i))*yy(nod(i))
463  g2z=g2z+(hl2(i)-hl3(i))*zz(nod(i))
464  enddo
465  g3x=g1y*g2z-g1z*g2y
466  g3y=g1z*g2x-g1x*g2z
467  g3z=g1x*g2y-g1y*g2x
468  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
469  g3x=g3x/xsum
470  g3y=g3y/xsum
471  g3z=g3z/xsum
472  !JACOBI MATRIX
473  xj11=g1x
474  xj12=g1y
475  xj13=g1z
476  xj21=g2x
477  xj22=g2y
478  xj23=g2z
479  xj31=g3x
480  xj32=g3y
481  xj33=g3z
482  !DETERMINANT OF JACOBIAN
483  det=xj11*xj22*xj33 &
484  +xj12*xj23*xj31 &
485  +xj13*xj21*xj32 &
486  -xj13*xj22*xj31 &
487  -xj12*xj21*xj33 &
488  -xj11*xj23*xj32
489  !
490  wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25*hh
491  term1(ic) = term1(ic)-wg*h(ip)*h(jp)
492  if( ip.EQ.jp ) term2(ip) = term2(ip)-wg*h(jp)*sink
493  !
494  enddo
495  enddo
496  enddo
497  enddo
498  !
499  return
500 
501  end subroutine heat_film_342
502  !----------------------------------------------------------------------*
503  subroutine heat_film_351(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
504  !----------------------------------------------------------------------*
505  !**
506  !** SET 351 FILM
507  !**
508  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
509  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
510  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
511  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
512  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
513  use hecmw
514  implicit real(kind=kreal) (a - h, o - z)
515  ! I/F VARIABLES
516  integer(kind=kint) :: NOD(MM)
517  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
518  ! LOCAL VARIABLES
519  real(kind=kreal) :: h(6), hr(6), hs(6), ht(6)
520  real(kind=kreal) :: xg(2), wgt(2)
521  !*************************
522  ! GAUSS INTEGRATION POINT
523  !*************************
524  data xg/-0.5773502691896,0.5773502691896/
525  data wgt/1.0,1.0/
526 
527  isuf = 0
528  !
529  if( ltype.EQ.1 ) then
530  nod(1) = 1
531  nod(2) = 2
532  nod(3) = 3
533  isuf = 1
534  else if( ltype.EQ.2 ) then
535  nod(1) = 6
536  nod(2) = 5
537  nod(3) = 4
538  isuf = 1
539  else if( ltype.EQ.3 ) then
540  nod(1) = 4
541  nod(2) = 5
542  nod(3) = 2
543  nod(4) = 1
544  isuf = 2
545  else if( ltype.EQ.4 ) then
546  nod(1) = 5
547  nod(2) = 6
548  nod(3) = 3
549  nod(4) = 2
550  isuf = 2
551  else if( ltype.EQ.5 ) then
552  nod(1) = 6
553  nod(2) = 4
554  nod(3) = 1
555  nod(4) = 3
556  isuf = 2
557  endif
558  !
559  if( isuf.EQ.2 ) then
560 
561  ic = 0
562  do ip = 1, 4
563  term2(ip) = 0.0
564  do jp = 1, 4
565  ic = ic + 1
566  term1(ic) = 0.0
567  !
568  do ig2=1,2
569  si=xg(ig2)
570  do ig1=1,2
571  ri=xg(ig1)
572  !
573  h(1)=0.25*(1.0-ri)*(1.0-si)
574  h(2)=0.25*(1.0+ri)*(1.0-si)
575  h(3)=0.25*(1.0+ri)*(1.0+si)
576  h(4)=0.25*(1.0-ri)*(1.0+si)
577  hr(1)=-.25*(1.0-si)
578  hr(2)= .25*(1.0-si)
579  hr(3)= .25*(1.0+si)
580  hr(4)=-.25*(1.0+si)
581  hs(1)=-.25*(1.0-ri)
582  hs(2)=-.25*(1.0+ri)
583  hs(3)= .25*(1.0+ri)
584  hs(4)= .25*(1.0-ri)
585  !
586  g1x=0.0
587  g1y=0.0
588  g1z=0.0
589  g2x=0.0
590  g2y=0.0
591  g2z=0.0
592  do i = 1, 4
593  g1x=g1x+hr(i)*xx(nod(i))
594  g1y=g1y+hr(i)*yy(nod(i))
595  g1z=g1z+hr(i)*zz(nod(i))
596  g2x=g2x+hs(i)*xx(nod(i))
597  g2y=g2y+hs(i)*yy(nod(i))
598  g2z=g2z+hs(i)*zz(nod(i))
599  enddo
600  g3x=g1y*g2z-g1z*g2y
601  g3y=g1z*g2x-g1x*g2z
602  g3z=g1x*g2y-g1y*g2x
603  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
604  !
605  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
606  if( ip.EQ.jp ) then
607  term2(ip)=term2(ip)-xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
608  endif
609  enddo
610  enddo
611  enddo
612  enddo
613  elseif( isuf.EQ.1 ) then
614  ax = xx( nod(2) ) - xx( nod(1) )
615  ay = yy( nod(2) ) - yy( nod(1) )
616  az = zz( nod(2) ) - zz( nod(1) )
617  bx = xx( nod(3) ) - xx( nod(1) )
618  by = yy( nod(3) ) - yy( nod(1) )
619  bz = zz( nod(3) ) - zz( nod(1) )
620  aa = dsqrt( (ay*bz-az*by)**2+(az*bx-ax*bz)**2+(ax*by-ay*bx)**2 )/6.0
621  !
622  ic = 0
623  do ip = 1, 3
624  term2(ip) = - aa*hh*sink
625  do jp = 1, 3
626  ic = ic + 1
627  term1(ic) = - aa*hh/3.0
628  enddo
629  enddo
630 
631  endif
632 
633  return
634 
635  end subroutine heat_film_351
636  !----------------------------------------------------------------------*
637  subroutine heat_film_352(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
638  !----------------------------------------------------------------------*
639  !**
640  !** SET 352 FILM
641  !**
642  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
643  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
644  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
645  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
646  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
647  ! F6 LTYPE=6 : FILM IN NORMAL-DIRECTION FOR FACE-6
648  use hecmw
649  implicit real(kind=kreal) (a - h, o - z)
650  ! I/F VARIABLES
651  integer(kind=kint) :: NOD(MM)
652  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
653  ! LOCAL VARIABLES
654  real(kind=kreal) :: h(8), hr(8), hs(8), ht(8), hl1(6), hl2(6), hl3(6)
655  real(kind=kreal) :: xg(3), wgt(3)
656  !*************************
657  ! GAUSS INTEGRATION POINT
658  !*************************
659  xg(1) = -0.7745966692
660  xg(2) = 0.0
661  xg(3) = 0.7745966692
662  wgt(1) = 0.5555555555
663  wgt(2) = 0.8888888888
664  wgt(3) = 0.5555555555
665  isuf = 0
666  !
667  if( ltype.EQ.1 ) then
668  nod(1) = 1
669  nod(2) = 2
670  nod(3) = 3
671  nod(4) = 7
672  nod(5) = 8
673  nod(6) = 9
674  isuf = 1
675  else if( ltype.EQ.2 ) then
676  nod(1) = 6
677  nod(2) = 5
678  nod(3) = 4
679  nod(4) = 11
680  nod(5) = 10
681  nod(6) = 12
682  isuf = 1
683  else if( ltype.EQ.3 ) then
684  nod(1) = 4
685  nod(2) = 5
686  nod(3) = 2
687  nod(4) = 1
688  nod(5) = 10
689  nod(6) = 14
690  nod(7) = 7
691  nod(8) = 13
692  isuf = 2
693  else if( ltype.EQ.4 ) then
694  nod(1) = 5
695  nod(2) = 6
696  nod(3) = 3
697  nod(4) = 2
698  nod(5) = 11
699  nod(6) = 15
700  nod(7) = 8
701  nod(8) = 14
702  isuf = 2
703  else if( ltype.EQ.5 ) then
704  nod(1) = 6
705  nod(2) = 4
706  nod(3) = 1
707  nod(4) = 3
708  nod(5) = 12
709  nod(6) = 13
710  nod(7) = 9
711  nod(8) = 15
712  isuf = 2
713  endif
714  !
715  if( isuf.EQ.2 ) then
716  ic = 0
717  do ip = 1,8
718  term2(ip) = 0.0
719  do jp = 1,8
720  ic = ic + 1
721  term1(ic) = 0.0
722  do ig2=1,3
723  si=xg(ig2)
724  do ig1=1,3
725  ri=xg(ig1)
726  rp=1.0+ri
727  sp=1.0+si
728  rm=1.0-ri
729  sm=1.0-si
730  h(1)=0.25*rm*sm*(-1.0-ri-si)
731  h(2)=0.25*rp*sm*(-1.0+ri-si)
732  h(3)=0.25*rp*sp*(-1.0+ri+si)
733  h(4)=0.25*rm*sp*(-1.0-ri+si)
734  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
735  h(6)=0.5*(1.0-si*si)*(1.0+ri)
736  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
737  h(8)=0.5*(1.0-si*si)*(1.0-ri)
738  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
739  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
740  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
741  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
742  hr(5)=-ri*(1.0-si)
743  hr(6)= .5*(1.0-si*si)
744  hr(7)=-ri*(1.0+si)
745  hr(8)=-.5*(1.0-si*si)
746  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
747  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
748  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
749  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
750  hs(5)=-.5*(1.0-ri*ri)
751  hs(6)=-si*(1.0+ri)
752  hs(7)= .5*(1.0-ri*ri)
753  hs(8)=-si*(1.0-ri)
754  g1x=0.0
755  g1y=0.0
756  g1z=0.0
757  g2x=0.0
758  g2y=0.0
759  g2z=0.0
760  do i = 1,8
761  g1x=g1x+hr(i)*xx(nod(i))
762  g1y=g1y+hr(i)*yy(nod(i))
763  g1z=g1z+hr(i)*zz(nod(i))
764  g2x=g2x+hs(i)*xx(nod(i))
765  g2y=g2y+hs(i)*yy(nod(i))
766  g2z=g2z+hs(i)*zz(nod(i))
767  enddo
768  g3x=g1y*g2z-g1z*g2y
769  g3y=g1z*g2x-g1x*g2z
770  g3z=g1x*g2y-g1y*g2x
771  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
772  wg=xsum*wgt(ig1)*wgt(ig2)*hh
773  term1(ic)=term1(ic)-wg*h(ip)*h(jp)
774  if( ip.EQ.jp ) term2(ip)=term2(ip)-wg*h(jp)*sink
775  enddo
776  enddo
777  enddo
778  enddo
779  !
780  elseif( isuf.EQ.1 ) then
781  ic = 0
782  do ip = 1, 6
783  term2(ip) = 0.0
784  do jp = 1, 6
785  ic = ic + 1
786  term1(ic) = 0.0
787  ! INTEGRATION OVER SURFACE
788  do l2=1,3
789  xl2=xg(l2)
790  x2 =(xl2+1.0)*0.5
791  do l1=1,3
792  xl1=xg(l1)
793  x1=0.5*(1.0-x2)*(xl1+1.0)
794  ! INTERPOLATION FUNCTION
795  x3=1.0-x1-x2
796  h(1)= x1*(2.0*x1-1.)
797  h(2)= x2*(2.0*x2-1.)
798  h(3)= x3*(2.0*x3-1.)
799  h(4)= 4.0*x1*x2
800  h(5)= 4.0*x2*x3
801  h(6)= 4.0*x1*x3
802  ! DERIVATIVE OF INTERPOLATION FUNCTION
803  ! FOR L1-COORDINATE
804  hl1(1)=4.0*x1-1.0
805  hl1(2)= 0.0
806  hl1(3)= 0.0
807  hl1(4)= 4.0*x2
808  hl1(5)= 0.0
809  hl1(6)= 4.0*x3
810  ! FOR L2-COORDINATE
811  hl2(1)= 0.0
812  hl2(2)= 4.0*x2-1.0
813  hl2(3)= 0.0
814  hl2(4)= 4.0*x1
815  hl2(5)= 4.0*x3
816  hl2(6)= 0.0
817  ! FOR L3-COORDINATE
818  hl3(1)= 0.0
819  hl3(2)= 0.0
820  hl3(3)= 4.0*x3-1.0
821  hl3(4)= 0.0
822  hl3(5)= 4.0*x2
823  hl3(6)= 4.0*x1
824  !
825  g1x=0.0
826  g1y=0.0
827  g1z=0.0
828  g2x=0.0
829  g2y=0.0
830  g2z=0.0
831  do i=1,6
832  g1x=g1x+(hl1(i)-hl3(i))*xx(nod(i))
833  g1y=g1y+(hl1(i)-hl3(i))*yy(nod(i))
834  g1z=g1z+(hl1(i)-hl3(i))*zz(nod(i))
835  g2x=g2x+(hl2(i)-hl3(i))*xx(nod(i))
836  g2y=g2y+(hl2(i)-hl3(i))*yy(nod(i))
837  g2z=g2z+(hl2(i)-hl3(i))*zz(nod(i))
838  enddo
839  g3x=g1y*g2z-g1z*g2y
840  g3y=g1z*g2x-g1x*g2z
841  g3z=g1x*g2y-g1y*g2x
842  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
843  g3x=g3x/xsum
844  g3y=g3y/xsum
845  g3z=g3z/xsum
846  !JACOBI MATRIX
847  xj11=g1x
848  xj12=g1y
849  xj13=g1z
850  xj21=g2x
851  xj22=g2y
852  xj23=g2z
853  xj31=g3x
854  xj32=g3y
855  xj33=g3z
856  !DETERMINANT OF JACOBIAN
857  det=xj11*xj22*xj33 &
858  +xj12*xj23*xj31 &
859  +xj13*xj21*xj32 &
860  -xj13*xj22*xj31 &
861  -xj12*xj21*xj33 &
862  -xj11*xj23*xj32
863  wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25*hh
864  term1(ic) = term1(ic)-wg*h(ip)*h(jp)
865  if( ip.EQ.jp ) term2(ip) = term2(ip)-wg*h(jp)*sink
866  !
867  enddo
868  enddo
869  enddo
870  enddo
871  endif
872  !
873  return
874 
875  end subroutine heat_film_352
876  !----------------------------------------------------------------------*
877  subroutine heat_film_361(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
878  !----------------------------------------------------------------------*
879  !**
880  !** SET 361 FILM
881  !**
882  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
883  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
884  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
885  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
886  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
887  ! F6 LTYPE=6 : FILM IN NORMAL-DIRECTION FOR FACE-6
888  use hecmw
889  implicit real(kind=kreal) (a - h, o - z)
890  ! I/F VARIABLES
891  integer(kind=kint) :: NOD(MM)
892  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
893  ! LOCAL VARIABLES
894  real(kind=kreal) :: h(4), hr(4), hs(4), ht(4)
895  real(kind=kreal) :: xg(2), wgt(2)
896  !*************************
897  ! GAUSS INTEGRATION POINT
898  !*************************
899  data xg/-0.5773502691896,0.5773502691896/
900  data wgt/1.0,1.0/
901  !
902  if( ltype.EQ.1 ) then
903  nod(1) = 1
904  nod(2) = 2
905  nod(3) = 3
906  nod(4) = 4
907  else if( ltype.EQ.2 ) then
908  nod(1) = 8
909  nod(2) = 7
910  nod(3) = 6
911  nod(4) = 5
912  else if( ltype.EQ.3 ) then
913  nod(1) = 5
914  nod(2) = 6
915  nod(3) = 2
916  nod(4) = 1
917  else if( ltype.EQ.4 ) then
918  nod(1) = 6
919  nod(2) = 7
920  nod(3) = 3
921  nod(4) = 2
922  else if( ltype.EQ.5 ) then
923  nod(1) = 7
924  nod(2) = 8
925  nod(3) = 4
926  nod(4) = 3
927  else if( ltype.EQ.6 ) then
928  nod(1) = 8
929  nod(2) = 5
930  nod(3) = 1
931  nod(4) = 4
932  endif
933  !
934  ic = 0
935  do ip = 1, 4
936  term2(ip) = 0.0
937  do jp = 1, 4
938  ic = ic + 1
939  term1(ic) = 0.0
940  do ig2=1,2
941  si=xg(ig2)
942  do ig1=1,2
943  ri=xg(ig1)
944  h(1)=0.25*(1.0-ri)*(1.0-si)
945  h(2)=0.25*(1.0+ri)*(1.0-si)
946  h(3)=0.25*(1.0+ri)*(1.0+si)
947  h(4)=0.25*(1.0-ri)*(1.0+si)
948  hr(1)=-.25*(1.0-si)
949  hr(2)= .25*(1.0-si)
950  hr(3)= .25*(1.0+si)
951  hr(4)=-.25*(1.0+si)
952  hs(1)=-.25*(1.0-ri)
953  hs(2)=-.25*(1.0+ri)
954  hs(3)= .25*(1.0+ri)
955  hs(4)= .25*(1.0-ri)
956  g1x=0.0
957  g1y=0.0
958  g1z=0.0
959  g2x=0.0
960  g2y=0.0
961  g2z=0.0
962  do i = 1, 4
963  g1x=g1x+hr(i)*xx(nod(i))
964  g1y=g1y+hr(i)*yy(nod(i))
965  g1z=g1z+hr(i)*zz(nod(i))
966  g2x=g2x+hs(i)*xx(nod(i))
967  g2y=g2y+hs(i)*yy(nod(i))
968  g2z=g2z+hs(i)*zz(nod(i))
969  enddo
970  g3x=g1y*g2z-g1z*g2y
971  g3y=g1z*g2x-g1x*g2z
972  g3z=g1x*g2y-g1y*g2x
973  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
974  !
975  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
976  if( ip.EQ.jp ) then
977  term2(ip)=term2(ip)-xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
978  endif
979  !
980  enddo
981  enddo
982  enddo
983  enddo
984 
985  return
986 
987  end subroutine heat_film_361
988  !----------------------------------------------------------------------*
989  subroutine heat_film_362(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
990  !----------------------------------------------------------------------*
991  !**
992  !** SET 362 FILM
993  !**
994  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
995  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
996  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
997  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
998  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
999  ! F6 LTYPE=6 : FILM IN NORMAL-DIRECTION FOR FACE-6
1000  use hecmw
1001  implicit real(kind=kreal) (a - h, o - z)
1002  ! I/F VARIABLES
1003  integer(kind=kint) :: NOD(MM)
1004  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
1005  ! LOCAL VARIABLES
1006  real(kind=kreal) :: h(8), hr(8), hs(8), ht(8)
1007  real(kind=kreal) :: xg(3), wgt(3)
1008  !*************************
1009  ! GAUSS INTEGRATION POINT
1010  !*************************
1011  xg(1) = -0.7745966692
1012  xg(2) = 0.0
1013  xg(3) = 0.7745966692
1014  wgt(1) = 0.5555555555
1015  wgt(2) = 0.8888888888
1016  wgt(3) = 0.5555555555
1017  !
1018  if( ltype.EQ.1 ) then
1019  nod(1) = 1
1020  nod(2) = 2
1021  nod(3) = 3
1022  nod(4) = 4
1023  nod(5) = 9
1024  nod(6) = 10
1025  nod(7) = 11
1026  nod(8) = 12
1027  else if( ltype.EQ.2 ) then
1028  nod(1) = 8
1029  nod(2) = 7
1030  nod(3) = 6
1031  nod(4) = 5
1032  nod(5) = 15
1033  nod(6) = 14
1034  nod(7) = 13
1035  nod(8) = 16
1036  else if( ltype.EQ.3 ) then
1037  nod(1) = 5
1038  nod(2) = 6
1039  nod(3) = 2
1040  nod(4) = 1
1041  nod(5) = 13
1042  nod(6) = 18
1043  nod(7) = 9
1044  nod(8) = 17
1045  else if( ltype.EQ.4 ) then
1046  nod(1) = 6
1047  nod(2) = 7
1048  nod(3) = 3
1049  nod(4) = 2
1050  nod(5) = 14
1051  nod(6) = 19
1052  nod(7) = 10
1053  nod(8) = 18
1054  else if( ltype.EQ.5 ) then
1055  nod(1) = 7
1056  nod(2) = 8
1057  nod(3) = 4
1058  nod(4) = 3
1059  nod(5) = 15
1060  nod(6) = 20
1061  nod(7) = 11
1062  nod(8) = 19
1063  else if( ltype.EQ.6 ) then
1064  nod(1) = 8
1065  nod(2) = 5
1066  nod(3) = 1
1067  nod(4) = 4
1068  nod(5) = 16
1069  nod(6) = 17
1070  nod(7) = 12
1071  nod(8) = 20
1072  endif
1073 
1074  ic = 0
1075  do ip = 1,8
1076  term2(ip) = 0.0
1077  do jp = 1,8
1078  ic = ic + 1
1079  term1(ic) = 0.0
1080  !
1081  do ig2=1,3
1082  si=xg(ig2)
1083  do ig1=1,3
1084  ri=xg(ig1)
1085 
1086  rp=1.0+ri
1087  sp=1.0+si
1088  rm=1.0-ri
1089  sm=1.0-si
1090  h(1)=0.25*rm*sm*(-1.0-ri-si)
1091  h(2)=0.25*rp*sm*(-1.0+ri-si)
1092  h(3)=0.25*rp*sp*(-1.0+ri+si)
1093  h(4)=0.25*rm*sp*(-1.0-ri+si)
1094  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
1095  h(6)=0.5*(1.0-si*si)*(1.0+ri)
1096  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
1097  h(8)=0.5*(1.0-si*si)*(1.0-ri)
1098  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1099  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1100  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
1101  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
1102  hr(5)=-ri*(1.0-si)
1103  hr(6)= .5*(1.0-si*si)
1104  hr(7)=-ri*(1.0+si)
1105  hr(8)=-.5*(1.0-si*si)
1106  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1107  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1108  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
1109  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
1110  hs(5)=-.5*(1.0-ri*ri)
1111  hs(6)=-si*(1.0+ri)
1112  hs(7)= .5*(1.0-ri*ri)
1113  hs(8)=-si*(1.0-ri)
1114 
1115  g1x=0.0
1116  g1y=0.0
1117  g1z=0.0
1118  g2x=0.0
1119  g2y=0.0
1120  g2z=0.0
1121 
1122  do i = 1,8
1123  g1x=g1x+hr(i)*xx(nod(i))
1124  g1y=g1y+hr(i)*yy(nod(i))
1125  g1z=g1z+hr(i)*zz(nod(i))
1126  g2x=g2x+hs(i)*xx(nod(i))
1127  g2y=g2y+hs(i)*yy(nod(i))
1128  g2z=g2z+hs(i)*zz(nod(i))
1129  enddo
1130 
1131  g3x=g1y*g2z-g1z*g2y
1132  g3y=g1z*g2x-g1x*g2z
1133  g3z=g1x*g2y-g1y*g2x
1134  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1135  !
1136  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
1137  if( ip.EQ.jp ) then
1138  term2(ip) = term2(ip)- xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
1139  endif
1140  !
1141  enddo
1142  enddo
1143  enddo
1144  enddo
1145 
1146  return
1147 
1148  end subroutine heat_film_362
1149  !----------------------------------------------------------------------*
1150  subroutine heat_film_731( NN,XX,YY,ZZ,HH,SINK,TERM1,TERM2 )
1151  !----------------------------------------------------------------------*
1152  !**
1153  !** SET 731 FILM
1154  !**
1155  ! F1 LTYPE=1 : SURFACE FILM
1156  use hecmw
1157  implicit real(kind=kreal) (a - h, o - z)
1158  ! I/F VARIABLES
1159  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), hh, sink, term1(nn*nn), term2(nn)
1160  !
1161  ax = xx(2) - xx(1)
1162  ay = yy(2) - yy(1)
1163  az = zz(2) - zz(1)
1164  bx = xx(3) - xx(1)
1165  by = yy(3) - yy(1)
1166  bz = zz(3) - zz(1)
1167  aa = dsqrt( (ay*bz-az*by)**2+(az*bx-ax*bz)**2+(ax*by-ay*bx)**2 )/6.0
1168  !
1169  ic = 0
1170  do ip = 1, nn
1171  term2(ip) = - aa*hh*sink
1172  do jp = 1, nn
1173  ic = ic + 1
1174  term1(ic) = - aa*hh/3.0
1175  enddo
1176  enddo
1177 
1178  return
1179 
1180  end subroutine heat_film_731
1181  !----------------------------------------------------------------------*
1182  subroutine heat_film_741( NN,XX,YY,ZZ,HH,SINK,TERM1,TERM2 )
1183  !----------------------------------------------------------------------*
1184  !**
1185  !** SET 741 FILM
1186  !**
1187  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
1188  use hecmw
1189  implicit real(kind=kreal) (a - h, o - z)
1190  ! I/F VARIABLES
1191  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(nn*nn), term2(nn)
1192  ! LOCAL VARIABLES
1193  real(kind=kreal) :: h(4), hr(4), hs(4), ht(4)
1194  real(kind=kreal) :: xg(2), wgt(2)
1195  !*************************
1196  ! GAUSS INTEGRATION POINT
1197  !*************************
1198  data xg/-0.5773502691896,0.5773502691896/
1199  data wgt/1.0,1.0/
1200  !
1201  ic = 0
1202  do ip = 1, nn
1203  term2(ip) = 0.0
1204  do jp = 1, nn
1205  ic = ic + 1
1206  term1(ic) = 0.0
1207  !
1208  do ig2=1,2
1209  si=xg(ig2)
1210  do ig1=1,2
1211  ri=xg(ig1)
1212  !
1213  h(1)=0.25*(1.0-ri)*(1.0-si)
1214  h(2)=0.25*(1.0+ri)*(1.0-si)
1215  h(3)=0.25*(1.0+ri)*(1.0+si)
1216  h(4)=0.25*(1.0-ri)*(1.0+si)
1217  hr(1)=-.25*(1.0-si)
1218  hr(2)= .25*(1.0-si)
1219  hr(3)= .25*(1.0+si)
1220  hr(4)=-.25*(1.0+si)
1221  hs(1)=-.25*(1.0-ri)
1222  hs(2)=-.25*(1.0+ri)
1223  hs(3)= .25*(1.0+ri)
1224  hs(4)= .25*(1.0-ri)
1225  !
1226  g1x=0.0
1227  g1y=0.0
1228  g1z=0.0
1229  g2x=0.0
1230  g2y=0.0
1231  g2z=0.0
1232  do i = 1, nn
1233  g1x=g1x+hr(i)*xx(i)
1234  g1y=g1y+hr(i)*yy(i)
1235  g1z=g1z+hr(i)*zz(i)
1236  g2x=g2x+hs(i)*xx(i)
1237  g2y=g2y+hs(i)*yy(i)
1238  g2z=g2z+hs(i)*zz(i)
1239  enddo
1240  g3x=g1y*g2z-g1z*g2y
1241  g3y=g1z*g2x-g1x*g2z
1242  g3z=g1x*g2y-g1y*g2x
1243  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1244  !
1245  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
1246  if( ip.EQ.jp ) then
1247  term2(ip)=term2(ip)-xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
1248  endif
1249  !
1250  enddo
1251  enddo
1252  enddo
1253  enddo
1254 
1255  return
1256 
1257  end subroutine heat_film_741
1258 end module m_heat_lib_film
Definition: hecmw.f90:6
This module provides subroutines to generate heat transfer boundary.
subroutine heat_film_242(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_232(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_342(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_352(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_341(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_731(NN, XX, YY, ZZ, HH, SINK, TERM1, TERM2)
subroutine heat_film_241(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_231(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_351(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_741(NN, XX, YY, ZZ, HH, SINK, TERM1, TERM2)
subroutine heat_film_362(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
subroutine heat_film_361(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)