mt04pp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt04pp.f
00002 !
00096                      SUBROUTINE MT04PP
00097 !                    *****************
00098 !
00099      &( T,XM,XMUL,SU,SV,SW,U,V,W,X,Y,Z,SURFAC,IKLE,NELEM,NELMAX,FORMUL)
00100 !
00101 !***********************************************************************
00102 ! BIEF   V6P3                                   21/08/2010
00103 !***********************************************************************
00104 !
00105 !
00106 !
00107 !
00108 !
00109 !
00110 !
00111 !
00112 !
00113 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00114 !| FORMUL         |-->| FORMULA DESCRIBING THE RESULTING MATRIX
00115 !| IKLE           |-->| CONNECTIVITY TABLE.
00116 !| NELEM          |-->| NUMBER OF ELEMENTS
00117 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00118 !| SU             |-->| STRUCTURE OF FUNCTIONS U
00119 !| SV             |-->| STRUCTURE OF FUNCTIONS V
00120 !| SW             |-->| STRUCTURE OF FUNCTIONS W
00121 !| SURFAC         |-->| AREA OF 2D ELEMENTS
00122 !| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
00123 !| U              |-->| FUNCTION USED IN THE FORMULA
00124 !| V              |-->| FUNCTION USED IN THE FORMULA
00125 !| W              |-->| FUNCTION USED IN THE FORMULA
00126 !| X              |-->| ABSCISSAE OF POINTS, PER ELEMENT
00127 !| Y              |-->| ORDINATES OF POINTS, PER ELEMENT
00128 !| Z              |-->| ELEVATIONS OF POINTS, STILL PER POINTS !!!
00129 !| XM             |<->| OFF-DIAGONAL TERMS
00130 !| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
00131 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00132 !
00133       USE BIEF, EX_MT04PP => MT04PP
00134 !
00135       IMPLICIT NONE
00136       INTEGER LNG,LU
00137       COMMON/INFO/LNG,LU
00138 !
00139 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00140 !
00141       INTEGER, INTENT(IN) :: NELEM,NELMAX
00142       INTEGER, INTENT(IN) :: IKLE(NELMAX,6)
00143 !
00144       DOUBLE PRECISION, INTENT(INOUT) :: T(NELMAX,6),XM(NELMAX,30)
00145       DOUBLE PRECISION, INTENT(IN)    :: SURFAC(NELMAX)
00146 !
00147       DOUBLE PRECISION, INTENT(IN) :: XMUL
00148       DOUBLE PRECISION, INTENT(IN) :: U(*),V(*),W(*)
00149 !
00150 !     STRUCTURES OF      U, V, W
00151       TYPE(BIEF_OBJ), INTENT(IN) :: SU,SV,SW
00152 !
00153       DOUBLE PRECISION, INTENT(IN) :: X(NELMAX,6),Y(NELMAX,6),Z(*)
00154 !
00155       CHARACTER(LEN=16), INTENT(IN) :: FORMUL
00156 !
00157 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00158 !
00159 !     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
00160 !
00161       DOUBLE PRECISION X2,X3,Y2,Y3
00162       DOUBLE PRECISION U0,V0,F0,G0,HH,C,DX,SURNORMU
00163 !
00164       INTEGER I1,I2,I3,I4,I5,I6,IELEM
00165 !
00166       INTRINSIC SQRT
00167 !
00168 !**********************************************************************
00169 !
00170       IF(SU%ELM.NE.41) THEN
00171         IF (LNG.EQ.1) WRITE(LU,1000) SU%ELM
00172         IF (LNG.EQ.2) WRITE(LU,1001) SU%ELM
00173 1000    FORMAT(1X,'MT04PP (BIEF) : TYPE DE U NON PREVU : ',I6)
00174 1001    FORMAT(1X,'MT04PP (BIEF) : TYPE OF U NOT IMPLEMENTED: ',I6)
00175         CALL PLANTE(1)
00176         STOP
00177       ENDIF
00178       IF(SV%ELM.NE.41) THEN
00179         IF (LNG.EQ.1) WRITE(LU,2000) SV%ELM
00180         IF (LNG.EQ.2) WRITE(LU,2001) SV%ELM
00181 2000    FORMAT(1X,'MT04PP (BIEF) : TYPE DE V NON PREVU : ',I6)
00182 2001    FORMAT(1X,'MT04PP (BIEF) : TYPE OF V NOT IMPLEMENTED: ',I6)
00183         CALL PLANTE(1)
00184         STOP
00185       ENDIF
00186 !     HERE WSCONV WHICH IS IN FACT DEFINED PER LAYER
00187       IF(SW%ELM.NE.41) THEN
00188         IF (LNG.EQ.1) WRITE(LU,3000) SW%ELM
00189         IF (LNG.EQ.2) WRITE(LU,3001) SW%ELM
00190 3000    FORMAT(1X,'MT04PP (BIEF) : TYPE DE W NON PREVU : ',I6)
00191 3001    FORMAT(1X,'MT04PP (BIEF) : TYPE OF W NOT IMPLEMENTED: ',I6)
00192         CALL PLANTE(1)
00193         STOP
00194       ENDIF
00195 !
00196       IF(FORMUL(1:7).EQ.'MAUGUG2') THEN
00197 !
00198 !     LOOP ON THE 3D ELEMENTS
00199 !
00200       DO IELEM = 1 , NELEM
00201 !
00202       I1 = IKLE(IELEM,1)
00203       I2 = IKLE(IELEM,2)
00204       I3 = IKLE(IELEM,3)
00205       I4 = IKLE(IELEM,4)
00206       I5 = IKLE(IELEM,5)
00207       I6 = IKLE(IELEM,6)
00208 !
00209 !     X2  =  X(I2) - X(I1)
00210 !     X3  =  X(I3) - X(I1)
00211 !     Y2  =  Y(I2) - Y(I1)
00212 !     Y3  =  Y(I3) - Y(I1)
00213       X2  =  X(IELEM,2)
00214       X3  =  X(IELEM,3)
00215       Y2  =  Y(IELEM,2)
00216       Y3  =  Y(IELEM,3)
00217 !
00218       U0 = (U(I1)+U(I2)+U(I3)+U(I4)+U(I5)+U(I6))/6.D0
00219       V0 = (V(I1)+V(I2)+V(I3)+V(I4)+V(I5)+V(I6))/6.D0
00220 !     AVERAGE HEIGHT OF PRISM
00221       HH = (Z(I4)-Z(I1)+Z(I5)-Z(I2)+Z(I6)-Z(I3))/3.D0
00222       C  = XMUL*HH/24.D0/SURFAC(IELEM)
00223 !
00224       T(IELEM,1)=2*C*(U0*Y3-U0*Y2-V0*X3+V0*X2)**2
00225       T(IELEM,2)=2*C*(-U0*Y3+V0*X3)**2
00226       T(IELEM,3)=2*C*(-U0*Y2+V0*X2)**2
00227       T(IELEM,4)=T(IELEM,1)
00228       T(IELEM,5)=T(IELEM,2)
00229       T(IELEM,6)=T(IELEM,3)
00230 !
00231       XM(IELEM,01)= 2*(-U0*Y3+V0*X3)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00232       XM(IELEM,02)=-2*(-U0*Y2+V0*X2)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00233       XM(IELEM,03)=(U0*Y3-U0*Y2-V0*X3+V0*X2)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00234       XM(IELEM,04)=(-U0*Y3+V0*X3)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00235       XM(IELEM,05)=-(-U0*Y2+V0*X2)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00236       XM(IELEM,06)=-2*(-U0*Y2+V0*X2)*(-U0*Y3+V0*X3)*C
00237       XM(IELEM,07)=(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-U0*Y3+V0*X3)*C
00238       XM(IELEM,08)=(-U0*Y3+V0*X3)*(-U0*Y3+V0*X3)*C
00239       XM(IELEM,09)=-(-U0*Y2+V0*X2)*(-U0*Y3+V0*X3)*C
00240       XM(IELEM,10)=-(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-U0*Y2+V0*X2)*C
00241       XM(IELEM,11)=-(-U0*Y3+V0*X3)*(-U0*Y2+V0*X2)*C
00242       XM(IELEM,12)=(-U0*Y2+V0*X2)*(-U0*Y2+V0*X2)*C
00243       XM(IELEM,13)= 2*(-U0*Y3+V0*X3)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00244       XM(IELEM,14)=-2*(-U0*Y2+V0*X2)*(U0*Y3-U0*Y2-V0*X3+V0*X2)*C
00245       XM(IELEM,15)=-2*(-U0*Y2+V0*X2)*(-U0*Y3+V0*X3)*C
00246 !
00247       ENDDO
00248 !
00249       ELSEIF(FORMUL(1:7).EQ.'MAUGUG1') THEN
00250 !
00251       DO IELEM = 1 , NELEM
00252 !
00253       I1 = IKLE(IELEM,1)
00254       I2 = IKLE(IELEM,2)
00255       I3 = IKLE(IELEM,3)
00256       I4 = IKLE(IELEM,4)
00257       I5 = IKLE(IELEM,5)
00258       I6 = IKLE(IELEM,6)
00259 !
00260 !     X2  =  X(I2) - X(I1)
00261 !     X3  =  X(I3) - X(I1)
00262 !     Y2  =  Y(I2) - Y(I1)
00263 !     Y3  =  Y(I3) - Y(I1)
00264       X2  =  X(IELEM,2)
00265       X3  =  X(IELEM,3)
00266       Y2  =  Y(IELEM,2)
00267       Y3  =  Y(IELEM,3)
00268 !
00269 !     VELOCITIES CONSIDERED CONSTANT
00270       U0 = (U(I1)+U(I2)+U(I3)+U(I4)+U(I5)+U(I6))/6.D0
00271       V0 = (V(I1)+V(I2)+V(I3)+V(I4)+V(I5)+V(I6))/6.D0
00272 !
00273       SURNORMU=1.D0/MAX(SQRT(U0**2+V0**2),1.D-8)
00274       DX=SQRT(2.D0*SURFAC(IELEM))
00275 !     HERE F0 AND G0 MUST BE
00276 !     DX*U/(2*NORME(U)) AND DY*V/(2*NORME(U)) BUT CONSIDERS DY=DX
00277 !     APPROXIMATED BY SQUARE ROOT OF (2*TRIANGLE AREA)
00278       F0 = 0.5D0*DX*U0*SURNORMU
00279       G0 = 0.5D0*DX*V0*SURNORMU
00280 !     AVERAGE HEIGHT OF PRISM
00281       HH = (Z(I4)-Z(I1)+Z(I5)-Z(I2)+Z(I6)-Z(I3))/3.D0
00282       C  = XMUL*HH/24.D0/SURFAC(IELEM)
00283 !
00284       T(IELEM,1)=2*(U0*Y3-U0*Y2-V0*X3+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00285       T(IELEM,2)=2*(-U0*Y3+V0*X3)*(-F0*Y3+G0*X3)*C
00286       T(IELEM,3)=2*(-U0*Y2+V0*X2)*(-F0*Y2+G0*X2)*C
00287       T(IELEM,4)=T(IELEM,1)
00288       T(IELEM,5)=T(IELEM,2)
00289       T(IELEM,6)=T(IELEM,3)
00290 !
00291       XM(IELEM,01)=2*(-U0*Y3+V0*X3)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00292       XM(IELEM,02)=-2*(-U0*Y2+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00293       XM(IELEM,03)=(U0*Y3-U0*Y2-V0*X3+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00294       XM(IELEM,04)=(-U0*Y3+V0*X3)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00295       XM(IELEM,05)=-(-U0*Y2+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00296       XM(IELEM,06)=-2*(-U0*Y2+V0*X2)*(-F0*Y3+G0*X3)*C
00297       XM(IELEM,07)=(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y3+G0*X3)*C
00298       XM(IELEM,08)=(-U0*Y3+V0*X3)*(-F0*Y3+G0*X3)*C
00299       XM(IELEM,09)=-(-U0*Y2+V0*X2)*(-F0*Y3+G0*X3)*C
00300       XM(IELEM,10)=-(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y2+G0*X2)*C
00301       XM(IELEM,11)=-(-U0*Y3+V0*X3)*(-F0*Y2+G0*X2)*C
00302       XM(IELEM,12)=(-U0*Y2+V0*X2)*(-F0*Y2+G0*X2)*C
00303       XM(IELEM,13)=2*(-U0*Y3+V0*X3)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00304       XM(IELEM,14)=-2*(-U0*Y2+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00305       XM(IELEM,15)=-2*(-U0*Y2+V0*X2)*(-F0*Y3+G0*X3)*C
00306 !
00307       XM(IELEM,16)= 2*(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y3+G0*X3)*C
00308       XM(IELEM,17)= -2*(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y2+G0*X2)*C
00309       XM(IELEM,21)= -2*(-U0*Y3+V0*X3)*(-F0*Y2+G0*X2)*C
00310       XM(IELEM,18)=(U0*Y3-U0*Y2-V0*X3+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00311       XM(IELEM,22)= (-U0*Y3+V0*X3)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00312       XM(IELEM,25)= -(-U0*Y2+V0*X2)*(F0*Y3-F0*Y2-G0*X3+G0*X2)*C
00313       XM(IELEM,19)= (U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y3+G0*X3)*C
00314       XM(IELEM,23)= (-U0*Y3+V0*X3)*(-F0*Y3+G0*X3)*C
00315       XM(IELEM,26)= -(-U0*Y2+V0*X2)*(-F0*Y3+G0*X3)*C
00316       XM(IELEM,28)= 2*(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y3+G0*X3)*C
00317       XM(IELEM,20)= -(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y2+G0*X2)*C
00318       XM(IELEM,24)= -(-U0*Y3+V0*X3)*(-F0*Y2+G0*X2)*C
00319       XM(IELEM,27)= (-U0*Y2+V0*X2)*(-F0*Y2+G0*X2)*C
00320       XM(IELEM,29)= -2*(U0*Y3-U0*Y2-V0*X3+V0*X2)*(-F0*Y2+G0*X2)*C
00321       XM(IELEM,30)= -2*(-U0*Y3+V0*X3)*(-F0*Y2+G0*X2)*C
00322 !
00323       ENDDO
00324 !
00325       ELSE
00326         IF (LNG.EQ.1) WRITE(LU,4000) FORMUL
00327         IF (LNG.EQ.2) WRITE(LU,4001) FORMUL
00328 4000    FORMAT(1X,'MT04PP (BIEF) : FORMULE NON PREVUE : ',A16)
00329 4001    FORMAT(1X,'MT04PP (BIEF) : UNEXPECTED FORMULA: ',A16)
00330         CALL PLANTE(1)
00331         STOP
00332       ENDIF
00333 !
00334 !-----------------------------------------------------------------------
00335 !
00336       RETURN
00337       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0