Initial commit from Polytechnique Montreal

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