summaryrefslogtreecommitdiff
path: root/Dragon/src/SAL_AUX_MOD.f90
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SAL_AUX_MOD.f90')
-rw-r--r--Dragon/src/SAL_AUX_MOD.f90496
1 files changed, 496 insertions, 0 deletions
diff --git a/Dragon/src/SAL_AUX_MOD.f90 b/Dragon/src/SAL_AUX_MOD.f90
new file mode 100644
index 0000000..e64416a
--- /dev/null
+++ b/Dragon/src/SAL_AUX_MOD.f90
@@ -0,0 +1,496 @@
+!
+!---------------------------------------------------------------------
+!
+!Purpose:
+! Support subroutines for isotropic and specular boundary conditions cases
+!
+!Copyright:
+! Copyright (C) 2001 Ecole Polytechnique de Montreal.
+!
+!Author(s):
+! X. Warin
+!
+!---------------------------------------------------------------------
+!
+MODULE SAL_AUX_MOD
+
+ USE PRECISION_AND_KINDS, ONLY : PDB,SMALL,PI,TWOPI,HALFPI,INFINITY
+ USE SAL_NUMERIC_MOD, ONLY : SAL141
+
+CONTAINS
+ SUBROUTINE SAL231(RTRACK,ITRACK,DELX,EX0,EY0,ANGLE)
+ !
+ !---------------------------------------------------------------------
+ !
+ !Purpose:
+ ! print out trajectory information
+ !
+ !Parameters: input
+ ! RTRACK floating point vectors to store trajectory information
+ ! ITRACK integer vectors to store trajectory information
+ ! DELX initial point of trajectory (D=0)
+ ! EX0 first direction cosine
+ ! EY0 second direction cosine
+ ! ANGLE track angle
+ !
+ !---------------------------------------------------------------------
+ !
+ USE PRECISION_AND_KINDS, ONLY : PDB
+ USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2
+ !**
+ IMPLICIT NONE
+ INTEGER, INTENT(IN), DIMENSION(:) :: ITRACK
+ REAL(PDB), INTENT(IN), DIMENSION(:) :: RTRACK
+ REAL(PDB), INTENT(IN) :: DELX,EX0,EY0
+ REAL(PDB), INTENT(IN) :: ANGLE
+ !**
+ INTEGER :: I,KM,K,II,JSURF,JPHI,JPSI,NTRACK,NBTOT,I0
+ REAL(PDB) :: ANG0
+ REAL(PDB), PARAMETER :: SMALLT=1.E-10
+ INTEGER, PARAMETER :: FOUT =6
+ !**
+ ANG0=ANGLE
+ NTRACK=ITRACK(1)
+ NBTOT=ITRACK(2)
+ WRITE(FOUT,'(//,3X,"TRAJECTORY",/,3X,"==========", &
+ & //3X,"DELX = ",1P,E12.4,6X,"EX EY = ",1P,2E12.4, &
+ & /,3X,"WITH SMALL = ",1P,E12.4,5X,"(",1P,E12.4," DEGREES )", &
+ & /,3X,"ANGLE # AND WEIGHT ",I6,1P,E12.4, &
+ & /,3X,"NBER OF SUB-TRAJ = ",I6, &
+ & /,3X,"NBER OF ELEM IN TRAJ = ",I6,/)') &
+ DELX,EX0,EY0,SMALLT,ANG0,ITRACK(7),RTRACK(7),NBTOT,NTRACK
+ WRITE(FOUT,'(/,20X,"SURF",3X,"PHI",3X,"PSI",5X,"SINPHI",7X, &
+ & "COSPHI",/,20X,4("-"),2(3X,3("-")),5X,6("-"),7X,6("-"))')
+ IF(ITRACK(5)/=0)THEN
+ JSURF=1
+ JPHI=ITRACK(5)
+ JPSI=JPHI
+ WRITE(FOUT,'(3X,A14,3I6,2X,1P,2E13.4)')'LEFT SURFACE',JSURF,JPHI,JPSI, &
+ RTRACK(3),RTRACK(5)
+ ENDIF
+ IF(ITRACK(6)/=0)THEN
+ JSURF=1
+ JPHI=ITRACK(6)
+ JPSI=JPHI
+ WRITE(FOUT,'(3X,A14,3I6,2X,1P,2E13.4)')'RIGHT SURFACE',JSURF,JPHI,JPSI, &
+ RTRACK(4),RTRACK(6)
+ ENDIF
+ IF(NTRACK/=0)THEN
+ I0=NTRACK+NNN
+ ! print sub-trajectories information
+ WRITE(FOUT,'(/," SUB-TRAJECTORIES:")')
+ WRITE(FOUT,'(/," NBER",2X,"NBER OF ELEM",2X," ANGLE", &
+ &/," ----",2X,"------------",2X," -----")')
+ DO I=1,NBTOT
+ WRITE(FOUT,'(I6,"*",4X,I6,5X,I6)')I,ITRACK(I0+2*I-1),ITRACK(I0+2*I)
+ ENDDO
+ ! print trajectory
+ WRITE(FOUT,'(/," TRAJECTORY:",/," NBER",3(5X,"REG",2X," LENGTH",3X,"ELEM"),/, &
+ & 3X,"----",3(5X,"---",3X,"------",3X,"----"),1X,/)')
+ DO I=1,NTRACK,3
+ II=I+NNN
+ KM=MIN(I+2,NTRACK)+NNN
+ WRITE(FOUT,'(1P,I6,"*",3(I7,E10.2,I7))')I, &
+ (ITRACK(K),RTRACK(K),ITRACK(K+NMAX2),K=II,KM)
+ ENDDO
+ ELSE
+ WRITE(FOUT,'(1X,"==> Track without intersections")')
+ ENDIF
+ WRITE(FOUT,'(/)')
+ !
+ END SUBROUTINE SAL231
+ !
+ SUBROUTINE SAL232(ITRACK,RTRACK,FACNRM,GG,SURFN,CURRN)
+ !
+ !---------------------------------------------------------------------
+ !
+ !Purpose:
+ ! computes numerical volumes: uses local macro arrays
+ !
+ !Parameters: input
+ ! ITRACK integer vectors to store trajectory information
+ ! RTRACK floating point vectors to store trajectory information
+ !
+ !Parameters: input/output
+ ! FACNRM numerical volumes per direction
+ ! SURFN numerical areas
+ ! CURRN numerical currents
+ ! GG geometry basic information.
+ !
+ !---------------------------------------------------------------------
+ !
+ USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC
+ USE SAL_TRACKING_TYPES, ONLY : NNN
+ !***
+ IMPLICIT NONE
+ INTEGER, INTENT(IN), DIMENSION(:) :: ITRACK
+ REAL(PDB), INTENT(IN), DIMENSION(:) :: RTRACK
+ REAL(PDB), INTENT(INOUT), DIMENSION(:), OPTIONAL :: SURFN,CURRN
+ REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: FACNRM
+ TYPE(T_G_BASIC) :: GG
+ ! DIMENSION ITRACK(*),RTRACK(*),FACNRM(NBREG,NPHI),
+ ! SURFN(NCURR,2),CURRN(NCURR,2)
+ !***
+ INTEGER :: LASTI,II,I,ICURR,NPHI,IPHI,C,P
+ REAL :: WT,WR
+ !***
+ ! total weight and space weight
+ NPHI=SIZE(FACNRM,2)
+ WT=REAL(RTRACK(7)) ; WR=REAL(RTRACK(8))
+ LASTI=ITRACK(1)
+ C=0; P=LASTI+NNN; IPHI=0
+ DO II=1+NNN,LASTI+NNN
+ IF((II-NNN)>C) THEN
+ C=C+ITRACK(P+1)
+ IPHI=ITRACK(P+2)
+ IF(IPHI>NPHI) IPHI=IPHI-NPHI
+ P=P+2
+ ENDIF
+ I=ITRACK(II)
+ IF(IPHI==0) CALL XABORT('SAL232: invalid IPHI')
+ FACNRM(I,IPHI)=FACNRM(I,IPHI)+RTRACK(II)*WR
+ ENDDO
+ IF(GG%NB_SURF2/=0)THEN
+ ICURR=ITRACK(5)
+ IF(ICURR>0) THEN
+ ! left surface: convert into 2d horizontal currents:
+ SURFN(ICURR)=SURFN(ICURR)+WT/RTRACK(5)
+ CURRN(ICURR)=CURRN(ICURR)+WT
+ ENDIF
+ ICURR=ITRACK(6)
+ IF(ICURR>0) THEN
+ ! right surface: convert into 2d horizontal currents:
+ SURFN(ICURR)=SURFN(ICURR)+WT/RTRACK(6)
+ CURRN(ICURR)=CURRN(ICURR)+WT
+ ENDIF
+ ENDIF
+ !
+ END SUBROUTINE SAL232
+ !
+ SUBROUTINE SAL235(NPIECE,THETA0,EX0,EY0,IPAR,RPAR,PEREXT,NPERIM)
+ !
+ !---------------------------------------------------------------------
+ !
+ !Purpose:
+ ! computes total perimeter projection on line orthogonal to
+ ! trajectory and with origin the center of coordinates
+ !
+ !Parameters: input
+ ! THETA0 THETA- and THETA+ for this trajectory
+ ! EX0 first direction cosine
+ ! EY0 second direction cosine
+ ! IPAR integer element data
+ ! RPAR floating point element data
+ ! PEREXT macro perimeter
+ ! NPERIM number of elements in perimeter
+ !
+ !Parameters: output
+ ! NPIECE number of pieces
+ !
+ !---------------------------------------------------------------------
+ !
+ USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,G_ELE_TYPE
+ USE SAL_TRACKING_TYPES, ONLY : DPIECE
+ !***
+ IMPLICIT NONE
+ REAL(PDB), INTENT(IN) :: EX0,EY0
+ REAL(PDB), INTENT(IN), DIMENSION(:) :: THETA0
+ INTEGER, INTENT(OUT) :: NPIECE
+ INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR
+ REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR
+ INTEGER, INTENT(IN), DIMENSION(:) :: PEREXT
+ INTEGER, INTENT(IN) :: NPERIM
+ !***
+ REAL(PDB) :: X,Y,RAD,DCENT,THETA1,THETA2,THETAM,THETA,DAUX,DMIN,DMAX
+ INTEGER :: L,ELEM,TYPE,IEND,ISIDE
+ LOGICAL :: LGONE
+ REAL(PDB), PARAMETER, DIMENSION(2) :: SIGNV = (/-1._PDB,1._PDB/)
+ INTEGER, PARAMETER :: FOUT =6
+ !***
+ LGONE=.TRUE.
+ DMIN=0._PDB; DMAX=0._PDB;
+ DO L=1,NPERIM
+ ELEM=PEREXT(L)
+ ! treat element
+ TYPE=IPAR(1,ELEM)
+ IF(TYPE==G_ELE_TYPE(2))THEN
+ ! circle:
+ RAD=RPAR(3,ELEM)
+ DCENT=RPAR(1,ELEM)*EY0-RPAR(2,ELEM)*EX0
+ ! project tangents to circle
+ IF(LGONE)THEN
+ DMIN=DCENT-RAD
+ DMAX=DCENT+RAD
+ LGONE=.FALSE.
+ ELSE
+ DMIN=MIN(DMIN,DCENT-RAD)
+ DMAX=MAX(DMAX,DCENT+RAD)
+ ENDIF
+ ELSE
+ DO IEND=1,2
+ CALL SAL141(TYPE,RPAR(:,ELEM),X,Y,IEND)
+ ! project end of element
+ DAUX=X*EY0-Y*EX0
+ IF(LGONE)THEN
+ DMIN=DAUX
+ DMAX=DAUX
+ LGONE=.FALSE.
+ ELSEIF(DAUX<DMIN)THEN
+ DMIN=DAUX
+ ELSEIF(DAUX>DMAX)THEN
+ DMAX=DAUX
+ ENDIF
+ ENDDO
+ IF(TYPE==G_ELE_TYPE(3))THEN
+ ! treat tangent to arc of circles
+ RAD=RPAR(3,ELEM)
+ DCENT=RPAR(1,ELEM)*EY0-RPAR(2,ELEM)*EX0
+ THETA1=RPAR(4,ELEM)
+ THETA2=RPAR(5,ELEM)
+ THETAM=THETA2-TWOPI
+ DO ISIDE=1,2
+ THETA=THETA0(ISIDE)
+ IF((THETA>THETA1.AND.THETA<THETA2).OR.THETA<THETAM) THEN
+ ! check projection of tangent
+ DAUX=DCENT+SIGNV(ISIDE)*RAD
+ IF(DAUX<DMIN)THEN
+ DMIN=DAUX
+ ELSEIF(DAUX>DMAX)THEN
+ DMAX=DAUX
+ ENDIF
+ ENDIF
+ ENDDO
+ ELSE
+ IF(TYPE/=G_ELE_TYPE(1)) CALL XABORT('SAL235: not implemented')
+ ENDIF
+ ENDIF
+ ENDDO
+ NPIECE=2
+ DPIECE(1)=DMIN
+ DPIECE(2)=DMAX
+ !
+ END SUBROUTINE SAL235
+ !
+ SUBROUTINE SAL237(EX0,EY0,MQ,NQ,PROJTAB,AXIS)
+ !
+ !---------------------------------------------------------------------
+ !
+ !Purpose:
+ ! computes geometry outline projections on the two symmetrical axis
+ !
+ !Parameters: input
+ ! EX0,EY0 horizontal tracking angle cosines
+ ! MQ,NQ cyclic tracking: for a rectangular geometry,
+ ! the length of track is SQRT((MQ*A)**2+(NQ*B)**2)
+ ! where A and B are rectangular sides.
+ !
+ !---------------------------------------------------------------------
+ !
+ USE PRECISION_AND_KINDS, ONLY : PDB
+ USE SAL_GEOMETRY_TYPES, ONLY : LENGTHX,LENGTHY,TYPGEO
+ USE SAL_TRACKING_TYPES, ONLY : DELR,LENGTH_INV_CYCL
+ !***
+ IMPLICIT NONE
+ REAL(PDB), INTENT(IN) :: EX0,EY0
+ INTEGER, INTENT(IN) :: MQ,NQ
+ !***
+ INTEGER :: IAXIS,NPIECE,AXIS(2),IMQ
+ REAL(PDB) :: X1,X2,DX,DR,R1,R2,PROJTAB(6),NNQ
+ REAL, PARAMETER :: EPS3 = 1.0E-3
+ !***
+ ! minimum radial interval to contain one trajectory
+ IF(TYPGEO.LE.7) THEN
+ ! Cartesian geometry
+ IF(NQ>0) THEN ! angle different from 0
+ X2=LENGTHX/NQ
+ IF(MQ>0) THEN
+ X1=0.
+ ELSE
+ X1=LENGTHX-X2; X2=LENGTHX
+ ENDIF
+ R1=X1*EY0; R2=X2*EY0
+ NPIECE=INT((R2-R1)/DELR)
+ IF(NPIECE==0) NPIECE=1
+ DR=(R2-R1)/NPIECE
+ DX=DR/EY0
+ IAXIS=1
+ ELSE ! angle equal to 0
+ X1=LENGTHY-LENGTHY/ABS(MQ)
+ X2=LENGTHY
+ R1=X1*EX0; R2=X2*EX0
+ NPIECE=INT((R2-R1)/DELR)
+ IF(NPIECE==0) NPIECE=1
+ DR=(R2-R1)/NPIECE
+ DX=DR/EX0
+ IF(TYPGEO.EQ.7) DX=DX*SQRT(2.0)
+ IAXIS=2
+ ENDIF
+ IF(TYPGEO.EQ.7) DR=DR/2.0
+ PROJTAB(:)=(/EX0,EY0,X1,X2,DX,DR/)
+ ! cyclical track length
+ IF(TYPGEO.EQ.5) THEN
+ ! translation
+ X1=LENGTHX*ABS(MQ); X2=LENGTHY*ABS(NQ)
+ ELSE
+ ! specular reflexion
+ X1=2.*LENGTHX*ABS(MQ); X2=2.*LENGTHY*ABS(NQ)
+ ENDIF
+ ELSE IF(TYPGEO.GE.8) THEN
+ ! hexagonal geometry
+ NNQ=(NQ-ABS(MQ))/2
+ X2=3.*LENGTHX/REAL(ABS(MQ)+2*NNQ)
+ IF(X2<=LENGTHX+EPS3) THEN
+ IAXIS=1
+ IF(MQ>0) THEN
+ X1=0.
+ ELSE
+ X1=LENGTHX-X2; X2=LENGTHX
+ ENDIF
+ ELSE
+ X1=0.
+ IF((TYPGEO==8).OR.(TYPGEO==10).OR.(TYPGEO==12)) THEN
+ ! MQ must be positive
+ IAXIS=2
+ X2=3.*LENGTHX/(MQ-NNQ)
+ ELSE
+ IF(MQ>0) THEN
+ IAXIS=2
+ X2=3.*LENGTHX/(2*MQ+NNQ)
+ ELSE
+ IAXIS=6
+ X2=3.*LENGTHX/(2*ABS(MQ)+NNQ)
+ ENDIF
+ ENDIF
+ ENDIF
+ R1=X1*EY0; R2=X2*EY0
+ NPIECE=INT((R2-R1)/DELR)
+ IF(NPIECE==0) NPIECE=1
+ DR=(R2-R1)/NPIECE
+ DX=DR/EY0
+ ! empirical correction of track weight in hexagonal cases (don't ask why)
+ IF(TYPGEO==8) THEN
+ IMQ=ABS(MQ)
+ IF(((IMQ==1).AND.(NQ==15)).OR.((IMQ==7).AND.(NQ==9)).OR.((IMQ==8).AND.(NQ==6))) THEN
+ DR=DR/3.0
+ ELSE IF(((IMQ==1).AND.(NQ==9)).OR.((IMQ==4).AND.(NQ==6)).OR.((IMQ==5).AND.(NQ==3))) THEN
+ DR=DR/3.0
+ ELSE IF(((IMQ==1).AND.(NQ==7)).OR.((IMQ==3).AND.(NQ==5)).OR.((IMQ==4).AND.(NQ==2))) THEN
+ DR=DR*5.0/12.0
+ ELSE IF(((IMQ==1).AND.(NQ==5)).OR.((IMQ==2).AND.(NQ==4)).OR.((IMQ==3).AND.(NQ==1))) THEN
+ DR=DR*4.0/9.0
+ ELSE IF(((IMQ==2).AND.(NQ==8)).OR.((IMQ==3).AND.(NQ==7)).OR.((IMQ==5).AND.(NQ==1))) THEN
+ DR=DR*7.0/15.0
+ ELSE IF(((IMQ==4).AND.(NQ==14)).OR.((IMQ==5).AND.(NQ==13)).OR.((IMQ==9).AND.(NQ==1))) THEN
+ DR=13.0*DR/27.0
+ ENDIF
+ ELSE IF(TYPGEO==9) THEN
+ IF(ABS(MQ)/NQ > 1) DR=(0.5+1.5*ABS(MQ)/NQ)*DR
+ ELSE IF(TYPGEO==10) THEN
+ IMQ=ABS(MQ)
+ IF(((IMQ==1).AND.(NQ==15)).OR.((IMQ==8).AND.(NQ==6))) THEN
+ DR=0.25*DR
+ ELSE IF(((IMQ==1).AND.(NQ==9)).OR.((IMQ==5).AND.(NQ==3))) THEN
+ DR=0.25*DR
+ ELSE IF(((IMQ==1).AND.(NQ==7)).OR.((IMQ==4).AND.(NQ==2))) THEN
+ DR=DR*5.0/14.0
+ ELSE IF(((IMQ==1).AND.(NQ==5)).OR.((IMQ==3).AND.(NQ==1))) THEN
+ DR=DR*2.0/5.0
+ ELSE IF(((IMQ==2).AND.(NQ==8)).OR.((IMQ==5).AND.(NQ==1))) THEN
+ DR=DR*7.0/16.0
+ ELSE IF(((IMQ==4).AND.(NQ==14)).OR.((IMQ==9).AND.(NQ==1))) THEN
+ DR=13.0*DR/28.0
+ ELSE
+ DR=0.5*DR
+ ENDIF
+ ELSE IF(TYPGEO==11) THEN
+ IMQ=ABS(MQ)
+ IF(((IMQ==1).AND.(NQ==15)).OR.((IMQ==8).AND.(NQ==6))) THEN
+ DR=0.5*DR
+ ELSE IF(((IMQ==1).AND.(NQ==9)).OR.((IMQ==5).AND.(NQ==3))) THEN
+ DR=0.5*DR
+ ELSE IF(((IMQ==1).AND.(NQ==7)).OR.((IMQ==4).AND.(NQ==2))) THEN
+ DR=0.7742663247*DR
+ ELSE IF(((IMQ==1).AND.(NQ==5)).OR.((IMQ==3).AND.(NQ==1))) THEN
+ DR=0.8257638060*DR
+ ELSE IF(((IMQ==2).AND.(NQ==8)).OR.((IMQ==5).AND.(NQ==1))) THEN
+ DR=0.8863607851*DR
+ ELSE IF(((IMQ==4).AND.(NQ==14)).OR.((IMQ==9).AND.(NQ==1))) THEN
+ DR=0.9327596923*DR
+ ENDIF
+ ELSE IF(TYPGEO==12) THEN
+ IMQ=ABS(MQ)
+ IF(((IMQ==1).AND.(NQ==15)).OR.((IMQ==7).AND.(NQ==9)).OR.((IMQ==8).AND.(NQ==6))) THEN
+ DR=DR/6.0
+ ELSE IF(((IMQ==1).AND.(NQ==9)).OR.((IMQ==4).AND.(NQ==6)).OR.((IMQ==5).AND.(NQ==3))) THEN
+ DR=DR/6.0
+ ELSE IF(((IMQ==1).AND.(NQ==7)).OR.((IMQ==3).AND.(NQ==5)).OR.((IMQ==4).AND.(NQ==2))) THEN
+ DR=DR*5.0/24.0
+ ELSE IF(((IMQ==1).AND.(NQ==5)).OR.((IMQ==2).AND.(NQ==4)).OR.((IMQ==3).AND.(NQ==1))) THEN
+ DR=DR*2.0/9.0
+ ELSE IF(((IMQ==2).AND.(NQ==8)).OR.((IMQ==3).AND.(NQ==7)).OR.((IMQ==5).AND.(NQ==1))) THEN
+ DR=0.232645602188*DR
+ ELSE IF(((IMQ==4).AND.(NQ==14)).OR.((IMQ==5).AND.(NQ==13)).OR.((IMQ==9).AND.(NQ==1))) THEN
+ DR=13.0*DR/54.0
+ ENDIF
+ ENDIF
+ PROJTAB(:)=(/EX0,EY0,X1,X2,DX,DR/)
+ ! cyclical track length
+ X1=LENGTHX*ABS(MQ)*1.5D0; X2=LENGTHX*(2*NNQ+ABS(MQ))*SQRT(3.D0)/2.D0
+ ENDIF
+ AXIS(:)=(/NPIECE,IAXIS/)
+ LENGTH_INV_CYCL=1./SQRT(X1*X1+X2*X2)
+ END SUBROUTINE SAL237
+ !
+ SUBROUTINE SAL220_1(ANGLE)
+ !
+ !---------------------------------------------------------------------
+ !
+ !Purpose:
+ ! computes unit vectors for the two rotative axis
+ !
+ !Parameters: input
+ ! ANGLE angle between the two rotative axis
+ !
+ !---------------------------------------------------------------------
+ !
+ USE PRECISION_AND_KINDS, ONLY : PDB
+ USE SAL_GEOMETRY_TYPES, ONLY : LENGTHX,LENGTHY,TYPGEO
+ USE SAL_TRACKING_TYPES, ONLY : HX,HY,BX,BY
+ !****
+ IMPLICIT NONE
+ REAL(PDB), INTENT(IN) :: ANGLE
+ !****
+ ! unit vector for the axis 1:
+ HX(1)=1.; HY(1)=0.
+ ! unit vector for the axis 2:
+ HX(2)=COS(ANGLE); HY(2)=SIN(ANGLE)
+ BX(:)=0.; BY(:)=0.
+ IF((TYPGEO.EQ.5).OR.(TYPGEO.EQ.6)) THEN
+ HX(3:4)=HX(1:2); HY(3:4)=HY(1:2)
+ BX(4)=LENGTHX
+ BY(3)=LENGTHY
+ ELSE IF(TYPGEO.EQ.7) THEN
+ HX(3)=0.; HY(3)=1.
+ BX(3)=LENGTHX
+ ELSE IF((TYPGEO.EQ.8).OR.(TYPGEO.EQ.10)) THEN
+ HX(3)=COS(ANGLE*2.); HY(3)=SIN(ANGLE*2.)
+ BX(3)=LENGTHX
+ ELSE IF(TYPGEO.EQ.9) THEN
+ HX(1)=1.; HY(1)=0.; BX(1)=-LENGTHX*0.5; BY(1)=-LENGTHY
+ HX(4)=1.; HY(4)=0.; BX(4)=-LENGTHX*0.5; BY(4)=LENGTHY
+ HX(2)=COS(ANGLE*2.); HY(2)=SIN(ANGLE*2.); BX(2)=BX(1); BY(2)=BY(1)
+ HX(5)=HX(2); HY(5)=HY(2); BX(5)=LENGTHX; BY(5)=0.
+ HX(3)=COS(ANGLE); HY(3)=SIN(ANGLE); BX(3)=-LENGTHX; BY(3)=0.
+ HX(6)=HX(3); HY(6)=HY(3); BX(6)=LENGTHX*0.5; BY(6)=-LENGTHY
+ ELSE IF(TYPGEO.EQ.11) THEN
+ BX(3)=LENGTHX*0.5; BY(3)=LENGTHY
+ HX(3)=1.; HY(3)=0.
+ BX(4)=LENGTHX
+ HX(4)=COS(ANGLE); HY(4)=SIN(ANGLE)
+ ELSE IF(TYPGEO.EQ.12) THEN
+ HX(3)=COS(PI/2.0+ANGLE); HY(3)=SIN(PI/2.0+ANGLE)
+ BX(3)=LENGTHX
+ ENDIF
+ !
+ END SUBROUTINE SAL220_1
+END MODULE SAL_AUX_MOD
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:31:10 +0000