Initial commit from Polytechnique Montreal

1 files changed, 1141 insertions, 0 deletions

diff --git a/Dragon/src/SAL_TRAJECTORY_MOD.f90 b/Dragon/src/SAL_TRAJECTORY_MOD.f90
new file mode 100644
index 0000000..e108fcb
--- /dev/null
+++ b/Dragon/src/SAL_TRAJECTORY_MOD.f90
@@ -0,0 +1,1141 @@
+!
+!---------------------------------------------------------------------
+!
+!Purpose:
+! Support module to compute a single track.
+!
+!Copyright:
+! Copyright (C) 2001 Ecole Polytechnique de Montreal.
+!
+!Author(s):
+! X. Warin
+!
+!---------------------------------------------------------------------
+!
+MODULE SAL_TRAJECTORY_MOD
+
+  USE PRECISION_AND_KINDS, ONLY : PDB,TWOPI,SMALL,PI
+  USE SAL_TRACKING_TYPES,  ONLY : NBER,LGOK,COSINE,AX,AY,EX,EY,ALPHA,F0,AT,BT, &
+       LGTYPE,R,D0,TORIG
+  USE SAL_NUMERIC_MOD,     ONLY : SALACO
+
+CONTAINS
+
+  SUBROUTINE SALTRA(DANGLT,NPERIM_MAC2,PERIM_MAC2,ISURF2_ELEM,IPAR,RPAR,PPERIM, &
+  PERIM,IBC2_ELEM,IDATA_BC2,BCDATA,PPERIM_MAC2,DIST_AXIS)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! computes intersection of trajectory (T): R=A+D*E with a mesh composed
+    ! of nodes and elements
+    !
+    !                   begin
+    !                     |
+    !    first entrance (LGMORE=f), get entry point
+    !                     |
+    !    get boundary condition at entry point and treat it
+    !                     |
+    !    if there is a left domain,do left tracking,inverse the trajectory
+    !                     |
+    !    if LGGEO1=t,test if there is a left re-entry point,keep it
+    !                     |
+    !            do basic tracking
+    !                     |
+    !    if LGGEO1=t,test if there is a right re-entry point,keep it
+    !                     |
+    !                    end
+    !
+    !Parameters: input
+    ! DANGLT       angle cosines
+    ! NPERIM_MAC2  number of elements composing perimeter of domain
+    ! PERIM_MAC2   elements composing perimeter of domain
+    ! ISURF2_ELEM  relative 2D surf number per elem
+    ! IPAR         integer geometry descriptors
+    ! RPAR         REAL(PDB) geometry descriptors
+    ! PPERIM       array pointer to elements in the perimeter of nodes
+    ! PERIM        array of elements in perimeter of nodes
+    !
+    !Parameters: input (optional data for cyclic tracking)
+    ! IBC2_ELEM    relative 2D boundary condition indices per element
+    ! IDATA_BC2    position of data per 2D boundary condition
+    ! BCDATA       table of boundary condition descriptor
+    ! PPERIM_MAC2  pointer to 'perim' and 'dist_axis':
+    !              PPERIM_MAC2(1): beginning of elements on axis 1
+    !              PPERIM_MAC2(2): beginning of elements on axis 2
+    !              PPERIM_MAC2(3): beginning of elements not on axis
+    ! DIST_AXIS    distance of points on this axis to the center (0,0)
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_GEOMETRY_TYPES,  ONLY : ISPEC
+    USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,ANGTAB,ELMTAB,CNT,CNT0,NB_TOT,DNEW,DINIT, &
+                                   NNEW,LNEW,IERR,LGMORE,DD0,NTRACK,EPS1,EX0,EY0,LGOK,IPART,DELX, &
+                                   N_AXIS
+    IMPLICIT NONE
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: DANGLT
+    INTEGER,   INTENT(IN)                  :: NPERIM_MAC2
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PERIM_MAC2,PPERIM,PERIM,ISURF2_ELEM
+    INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: RPAR
+    INTEGER,   INTENT(IN), DIMENSION(:), OPTIONAL :: IBC2_ELEM,IDATA_BC2,PPERIM_MAC2
+    REAL(PDB), INTENT(IN), DIMENSION(:), OPTIONAL :: DIST_AXIS
+    REAL(PDB), INTENT(IN), DIMENSION(:,:), OPTIONAL  :: BCDATA
+    !***
+    INTEGER   :: IOUT,LEN,P1,P2,IPHI
+    LOGICAL   :: LGLEFT 
+    REAL(PDB) :: RADIA
+    !***
+    !     initiate nber of sub-trajectories
+    NB_TOT=0
+    IF(ISPEC == 0) THEN
+      !*    compute entry point for a basic trajectory
+      CALL SAL240_3(PERIM_MAC2,NPERIM_MAC2,IPAR,RPAR)
+    ELSE
+      EX=EX0; EY=EY0
+      LGOK=.FALSE.
+      !     initiate all elements status to untreated
+      IPART(1,:)=-1
+      !     radia of the axis
+      IF(PPERIM_MAC2(N_AXIS+1)==PPERIM_MAC2(N_AXIS)) THEN
+        RADIA=0.
+      ELSE
+        RADIA=DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1)
+      ENDIF
+      IF(DELX<RADIA) THEN
+        !        point 'DELX' is on one of the elements on the axis
+        P1=PPERIM_MAC2(N_AXIS); P2=PPERIM_MAC2(N_AXIS+1)-1
+        CALL SAL241_2(P2-P1+1,PERIM_MAC2(P1:P2),DIST_AXIS(P1:P2),IPAR)
+        LGOK=.TRUE.
+      ELSE
+        WRITE(*,*) 'PPERIM_MAC2(N_AXIS+1),PPERIM_MAC2(N_AXIS) :',PPERIM_MAC2(N_AXIS+1),PPERIM_MAC2(N_AXIS)
+        WRITE(*,*) 'DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1) :',DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1)
+        WRITE(*,*) 'DELX :',DELX
+        CALL XABORT('SALTRA: Cant find entry point')
+      ENDIF
+    ENDIF
+    DINIT=DNEW
+    !*    treat boundary condition on entry point
+    CALL SAL245(ISURF2_ELEM,IPAR,RPAR,IOUT,LNEW,NNEW,DNEW,LGLEFT)
+    !*    abort if there is a left domain
+    IF(LGLEFT) CALL XABORT('SALTRA: LGLEFT True')
+    !*    track the basic domain
+    IF(ISPEC == 0) THEN
+      IPHI=ANGLE_TO_NUMBER(EX0,EY0,DANGLT)
+      CALL SAL240_4(PERIM,PPERIM,IPAR,RPAR,IPHI,ISURF2_ELEM,IOUT)
+    ELSE
+      CALL SAL240_4_2(DANGLT,PPERIM,PERIM,IPAR,RPAR,IDATA_BC2, &
+         BCDATA,PERIM_MAC2,PPERIM_MAC2,DIST_AXIS,IBC2_ELEM)
+    ENDIF
+    !     trajectory entered into the element joint
+    IF(IERR==0) RETURN
+    !     change DD0
+    IF(LGMORE) DD0=DNEW+1.001*EPS1
+    !*    put NB_TOT, ANGTAB to ITRAC2
+    !     total length of the trajectory
+    NTRACK=CNT-(CNT0+NNN)
+    ITRAC2(CNT0+1)=NTRACK
+    !     total nber of the sub-trajectories
+    ITRAC2(CNT0+2)=NB_TOT
+    !     put ANGTAB
+    LEN=2*NB_TOT
+    IF((CNT+LEN)>=NMAX2) THEN
+       CALL XABORT('SALTRA: Buffer overflow')
+    ELSE
+       ITRAC2(CNT+1:CNT+LEN)=ANGTAB(1:LEN)
+       ITRAC2(CNT+NMAX2+1:CNT+NMAX2+LEN)=ELMTAB(1:LEN)
+       CNT=CNT+LEN
+    ENDIF
+    !
+  END SUBROUTINE SALTRA
+
+  SUBROUTINE SAL241_2(NPERIM,PERIM,DIST_AXIS,IPAR)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! compute an entry point on the axial element
+    !
+    !Parameters: input
+    ! NPERIM          = number of elements on this axis
+    ! PERIM           = elements on this axis in perimeter
+    ! DIST_AXIS       = distance of points on this axis to the center (0,0)
+    ! IPAR            = integer descriptor of elements
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : IPART,N_AXIS,DNEW,DELX,NNEW,LNEW,COSINE,AX, &
+                                   AY,HX,HY,BX,BY,EX,EY
+    INTEGER,   INTENT(IN)                  :: NPERIM
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PERIM
+    REAL(PDB), INTENT(IN), DIMENSION(:)    :: DIST_AXIS
+    INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
+    INTEGER    :: I,J
+    !***
+    LNEW=0
+    !*    compute crossed element
+    DO I=1,NPERIM
+       IF(DELX<=DIST_AXIS(I)) THEN
+          LNEW=PERIM(I)
+          EXIT
+       ENDIF
+    ENDDO
+    IF(LNEW==0) CALL XABORT('SAL241_2: Error of distances on the axis')
+    !*    get entered node
+    NNEW=IPAR(2,LNEW)
+    IF(NNEW<0) NNEW=IPAR(3,LNEW)
+    IF(NNEW<0) CALL XABORT('SAL241_2: Error of element data')
+    !*    compute DNEW at entry point
+    DNEW=DELX*(EX*HX(N_AXIS)+EY*HY(N_AXIS))
+    IF(N_AXIS>2) DNEW=DNEW+BX(N_AXIS)*EX+BY(N_AXIS)*EY
+    !*    compute COSINE
+    COSINE=ABS(HX(N_AXIS)*EY-HY(N_AXIS)*EX)
+    !*    set all elements in this axis to be 'treated (0)'
+    !     others are 'untreated'
+    IPART(1,:)=-1
+    DO I=1,NPERIM
+       J=PERIM(I)
+       IPART(1,J)=0
+    ENDDO
+    !*    initial point
+    AX=BX(N_AXIS)+DELX*HX(N_AXIS)-DNEW*EX
+    AY=BY(N_AXIS)+DELX*HY(N_AXIS)-DNEW*EY
+    !
+  END SUBROUTINE SAL241_2
+  !
+  SUBROUTINE SAL240_3(PERIM_MAC2,NPERIM_MAC2,IPAR,RPAR)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! compute entry point for a basic trajectory
+    !
+    !Parameters: input
+    ! PERIM_MAC2   elements composing perimeter of domain
+    ! NPERIM_MAC2  number of elements composing perimeter of domain
+    ! IPAR         integer geometry descriptors
+    ! RPAR         floating point geometry descriptors
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : IPART,EX0,EY0,EX,EY,DD0,LGOK,LGMORE,NBER, &
+         DNEW,NNEW,LNEW,DINIT
+    IMPLICIT NONE
+    INTEGER,   INTENT(IN)                  :: NPERIM_MAC2
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PERIM_MAC2
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: RPAR
+    INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
+    !***
+    EX=EX0; EY=EY0
+    !     enter trajectory with vectors a and e and initial distance dd0
+    !     initiate all elements status to untreated
+    IPART(1,:)=-1
+    CALL SAL241(PERIM_MAC2,NPERIM_MAC2,IPAR,RPAR,DD0,-100,DNEW,NNEW,LNEW)
+    DINIT=DNEW
+    !     if we have sevaral entry points
+    LGMORE=NBER>3
+    !*    if not succeed
+    IF(.NOT.LGOK) CALL XABORT('SAL240_3: Could not enter domain')
+    !
+  END SUBROUTINE SAL240_3
+  !
+  SUBROUTINE SAL245(ISURF2_ELEM,IPAR,RPAR,IOUT,LOLD,NOLD,DOLD,LGLEFT)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! boundary condition treatment at entry point
+    !
+    !Parameters: input
+    ! ISURF2_ELEM  relative 2D surf nber per elem
+    ! IPAR         integer geometry descriptors
+    ! RPAR         floating point geometry descriptors
+    !
+    !Parameters: input/output
+    ! LOLD         entry point
+    ! NOLD         node index
+    ! DOLD         distance
+    ! LGLEFT       flag set to TRUE if there is a left domain
+    !
+    !Parameters: output
+    ! IOUT         (to load outside surface) = 5 (left trajectory)
+    !                                          6 (right trajectory)
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : ITRAC2,RTRAC2,CNT0,COSINE,EX,EY
+    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE,ISPEC
+    IMPLICIT NONE
+    INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: ISURF2_ELEM
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: RPAR
+    INTEGER,   INTENT(INOUT)   :: LOLD,NOLD
+    REAL(PDB), INTENT(INOUT)   :: DOLD
+    LOGICAL, INTENT(INOUT)     :: LGLEFT
+    INTEGER, INTENT(OUT)       :: IOUT
+    !***
+    INTEGER   :: BCIN,SURF
+    !***
+    !     initiate surface info
+    ITRAC2(CNT0+5)=0
+    ITRAC2(CNT0+6)=0
+    !*    get boundary condition at entry surface
+    BCIN=IPAR(2,LOLD)
+    IF(NOLD==BCIN)BCIN=IPAR(3,LOLD)
+    !     add entry cosine anyhow (for characteristics)
+    RTRAC2(CNT0+2)=COSINE
+    LGLEFT=.FALSE.
+    IF(ISPEC == 1) RETURN
+    IF(BCIN>=G_BC_TYPE(0).OR.BCIN==G_BC_TYPE(-1).OR.BCIN==G_BC_TYPE(1))THEN
+       !*       trajectory enters through vacuum:
+       !         - compute angle, trajectory-normal and store cos and sin
+       !           (inverse vector E in order to pass outgoing direction)
+       CALL SAL247_1(RPAR(:,LOLD),IPAR(:,LOLD),DOLD,RTRAC2(CNT0+3),RTRAC2(CNT0+5),-EX,-EY)
+       !*       get index of the entry surface
+       SURF=ISURF2_ELEM(LOLD)
+       IF(SURF/=0) ITRAC2(CNT0+5)=SURF ! store 2D surface index
+       IOUT=6
+    ELSE
+       CALL XABORT('SAL245: Reversed direction')
+    ENDIF
+    !
+  END SUBROUTINE SAL245
+  !
+  SUBROUTINE SAL240_4(PERIM,PPERIM,IPAR,RPAR,IPHI,ISURF2_ELEM,IOUT)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! track a trajectory until leaving the domain
+    !         begin: having (nnew,dnew,lnew)
+    !               |
+    !     (1)write horizontal angle information
+    !               |
+    !     (2)compute successive crossed nodes
+    !               |
+    !     (3)boundary condition:
+    !        re-entry: get re-entry point,go to (1)
+    !        go out: compute surface number,end
+    !
+    !Parameters: input
+    ! PERIM        array of elements in perimeter of nodes
+    ! PPERIM       array pointer to elements in the perimeter of nodes
+    ! IPAR         integer geometry descriptors
+    ! RPAR         floating point geometry descriptors
+    ! IPHI         angular index of the track
+    ! ISURF2_ELEM  relative 2D surf number per elem
+    ! IOUT         (to load outside surface) = 5 (left trajectory)
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,RTRAC2,ANGTAB,ELMTAB,PRTIND,CNT, &
+         CNT0,EX,EY,LGOK,NB_TOT,NB_MAX,DNEW,NNEW,LNEW,IERR
+    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE
+    IMPLICIT NONE
+    ! IN VARIABLE
+    INTEGER,   INTENT(IN)                  :: IPHI
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PPERIM,PERIM
+    INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: RPAR
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: ISURF2_ELEM
+    INTEGER,   INTENT(IN)                  :: IOUT
+    !***
+    INTEGER    :: NOLD,LOLD,P1,P2,CNT1,SURF
+    REAL(PDB)  :: DOLD,LENGTH
+    LOGICAL    :: LGON
+    INTEGER, PARAMETER :: FOUT =6
+    !***
+    LGON = .TRUE.
+    !     initiate counter
+    CNT1=CNT
+    EXTERIOR : DO WHILE(LGON)
+       NB_TOT=NB_TOT+1
+       IF(NB_TOT > NB_MAX) CALL XABORT('SAL240_4: NB_TOT overflow')
+       !        keep horizontal phi nber in angtab
+       IF(NB_TOT>1) CALL XABORT('SAL240_4: Angtab overflow')
+       ANGTAB(2*NB_TOT)=IPHI
+       ELMTAB(2*NB_TOT-1)=LNEW ; ELMTAB(2*NB_TOT)=0 ;
+       !
+       !*       track a sub-trajectory         
+       INTERIOR: DO WHILE(NNEW>0)
+          !           UPDATE DATA TO COMPUTE NEXT NODE:
+          DOLD=DNEW
+          LOLD=LNEW
+          NOLD=NNEW
+          !           crossing NODE NOLD
+          !           input: trajectory (T):R=A+D*E => A = (AX,AY), E = (EX,EY)
+          !           DOLD = D at last intersection
+          !           NOLD = NODE just entered
+          P1=PPERIM(NOLD); P2=PPERIM(NOLD+1)-1
+          CALL SAL241(PERIM(P1:P2),P2-P1+1,IPAR,RPAR,DOLD,NOLD,DNEW,NNEW,LNEW)
+          !           at return from SAL241:
+          !           DNEW        = D at point exiting node
+          !           COSINE      = cosine of trajectory with exiting normal
+          !           NNEW        = new node entered
+          !           LNEW        = element crossed when exiting node
+          !           NBER        = number of intersections with perimeter
+          !           LGOK        = .TRUE. if trajectory exits the node
+          IF(.NOT.LGOK) THEN
+             IERR=0
+             IF(PRTIND>0)WRITE(FOUT,'("SAL240_4 ==> couldnt exit node ",I5)') NOLD
+             RETURN
+          ENDIF
+          !           store data
+          LENGTH=DNEW-DOLD
+          !           store new length in track arrays
+          CNT=CNT+1
+          IF(CNT>=NMAX2) CALL XABORT('SAL240_4: NMAX2 overflow')
+          RTRAC2(CNT)=LENGTH
+          ITRAC2(CNT+NMAX2)=LNEW
+          ITRAC2(CNT)=NOLD
+       END DO INTERIOR
+       !
+       !*       exiting motif and analyzing bc condition
+       LGON =(NNEW<G_BC_TYPE(1)).AND.(NNEW>=G_BC_TYPE(5))
+       !        STORE NBER OF REGIONS TO ANGTAB
+       ANGTAB(2*NB_TOT-1)=CNT-CNT1
+       CNT1=CNT
+       IF(LGON) THEN
+          CALL XABORT('SAL240_4: Lgon is true')
+       ELSE
+          !           compute leaving surface
+          IF(NNEW<=0)THEN
+             !        we got a surface: end of trajectory
+             !        vacuum bd condition. put a marker (surface number) to allow
+             !        psi and pss computation
+             !        store exiting surface, cosphi and sinphi
+             !        compute angle trajectory-normal
+             CALL SAL247_1(RPAR(:,LNEW),IPAR(:,LNEW),DNEW, &
+                  RTRAC2(CNT0+IOUT-2),RTRAC2(CNT0+IOUT),EX,EY)
+             !        store 2D surface-cone nber
+             SURF=ISURF2_ELEM(LNEW)
+             IF(SURF/=0) ITRAC2(CNT0+IOUT)=SURF
+          ENDIF
+       ENDIF
+       !
+    ENDDO EXTERIOR
+    !     set success flag
+    IERR=1
+    !
+  END SUBROUTINE SAL240_4
+  !
+  SUBROUTINE SAL240_4_2(DANGLT,PPERIM,PERIM,IPAR,RPAR,IDATA_BC2,BCDATA, &
+       PERIM_MAC2,PPERIM_MAC2,DIST_AXIS,IBC2_ELEM)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! track a cyclic trajectory until the track length equal to the predefined
+    ! length of this trajectory
+    !      begin: having (nnew,dnew,lnew)
+    !               |
+    !     (1)write horizontal angle information
+    !               |
+    !     (2)compute successive crossed nodes
+    !               |
+    !     (3)if abs(total_length*length_inv_cycl-1)>eps1 :
+    !        yes:continue, get re-entry point,go to (1)
+    !        no: end of tracking
+    !
+    !Parameters: input
+    ! DANGLT       angle cosines
+    ! PERIM        array of elements in perimeter of nodes
+    ! PPERIM       array pointer to elements in the perimeter of nodes
+    ! IPAR         integer geometry descriptors
+    ! RPAR         real geometry descriptors
+    ! IDATA_BC2    position of bc data per 2D boundary conditions
+    ! PERIM_MAC2   elements composing perimeter of domain
+    ! PPERIM_MAC2  pointer to 'PERIM' and 'DIST_AXIS':
+    !                  PPERIM_MAC2(1):beginning of elements on axis 1
+    !                  PPERIM_MAC2(2):beginning of elements on axis 2
+    !                  PPERIM_MAC2(3):beginning of elements not on axis
+    ! DIST_AXIS    distance of points on this axis to the center (0,0)
+    ! BCDATA       table of bc descriptor
+    ! IBC2_ELEM    relative 2D bc nber per elem
+    !
+    !Parameters: output
+    ! NB_TOT       total number of sub-trajectories
+    ! ANGTAB       table of {N_K,ANGLE_K} (1<K<NB_TOT)
+    !              N_K: nber of regions in kth sub-trajectory
+    !              ANGLE_K: angle nber of kth sub-trajectory
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,RTRAC2,ANGTAB,ELMTAB,N_AXIS,CNT, &
+         EX,EY,AX,AY,LGOK,LENGTH_INV_CYCL,NB_TOT,NB_MAX,DNEW,NNEW,LNEW,IERR,EPS1,TORIG, &
+         N_AXIS_KEEP,IMPX
+    USE SAL_GEOMETRY_TYPES,      ONLY : G_BC_TYPE
+    IMPLICIT NONE
+    ! in variable
+    !************
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: DANGLT
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PPERIM,PERIM
+    INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: RPAR
+    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: BCDATA
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PPERIM_MAC2,IBC2_ELEM
+    INTEGER,   INTENT(IN), DIMENSION(:)    :: PERIM_MAC2
+    REAL(PDB), INTENT(IN), DIMENSION(:),OPTIONAL :: DIST_AXIS
+    INTEGER,   INTENT(IN), DIMENSION(:),OPTIONAL :: IDATA_BC2
+    ! local variable
+    INTEGER    :: NOLD,LOLD,P1,P2,CNT1,IDATA,OK
+    REAL(PDB)  :: DOLD,LENGTH,LENGTH_TOT
+    LOGICAL    :: LGON
+    INTEGER,   POINTER, DIMENSION(:) :: ITRAC3
+    REAL(PDB), POINTER, DIMENSION(:) :: RTRAC3
+    INTEGER, PARAMETER :: FOUT =6
+    !***
+    lgon=.true.
+    !     initiate counter
+    CNT1=CNT
+    LENGTH_TOT=0.
+    EXTERIOR : DO WHILE(LGON)
+       NB_TOT=NB_TOT+1
+       IF(NB_TOT > NB_MAX) CALL XABORT('SAL240_4_2: NB_TOT overflow')
+       N_AXIS_KEEP(NB_TOT)=N_AXIS
+       !        set horizontal angle index in angtab
+       ANGTAB(2*NB_TOT)=ANGLE_TO_NUMBER(EX,EY,DANGLT)
+       ELMTAB(2*NB_TOT-1)=LNEW ; ELMTAB(2*NB_TOT)=0 ;
+       TORIG(1,NB_TOT)=AX+DNEW*EX ; TORIG(2,NB_TOT)=AY+DNEW*EY ;
+       IF(IMPX > 4) WRITE(6,*) 'SAL240_4_2: beginning of track=',TORIG(:2,NB_TOT)
+       !
+       !*       track a sub-trajectory
+       INTERIOR: DO WHILE(NNEW.GT.0)
+          !           update data to compute next node:
+          DOLD=DNEW
+          LOLD=LNEW
+          NOLD=NNEW
+          !           crossing node NOLD
+          !           input: trajectory (t):r=a+d*e => a = (ax,ay), e = (ex,ey)
+          !           DOLD = d at last intersection
+          !           NOLD = node just entered
+          P1=PPERIM(NOLD); P2=PPERIM(NOLD+1)-1
+          CALL SAL241(PERIM(P1:P2),P2-P1+1,IPAR,RPAR,DOLD,NOLD,DNEW,NNEW,LNEW)
+          !           at return from SAL241:
+          !           DNEW        = d at point exiting node
+          !           COSINE      = cosine of trajectory with exiting normal
+          !           NNEW        = new node entered
+          !           LNEW        = element crossed when exiting node
+          !           NBER        = nber of intersections with perimeter
+          !           LGOK        = .true. if trajectory exits the node
+          IF(.NOT.LGOK) THEN
+             IERR=0
+             WRITE(FOUT,'(" SAL240_4_2 ==> couldnt exit node ",I5)') NOLD
+             RETURN
+          ENDIF
+          !           store data
+          LENGTH=DNEW-DOLD
+          !           add to total length
+          LENGTH_TOT=LENGTH_TOT+LENGTH
+          CNT=CNT+1
+          IF(CNT>=NMAX2) THEN
+             ALLOCATE(ITRAC3(4*NMAX2),RTRAC3(2*NMAX2),STAT=OK)
+             IF(OK/=0) CALL XABORT('SAL240_4_2: NMAX2 overflow.')
+             RTRAC3(:NMAX2)=RTRAC2(:NMAX2)
+             ITRAC3(:2*NMAX2)=ITRAC2(:2*NMAX2)
+             DEALLOCATE(RTRAC2,ITRAC2)
+             RTRAC2=>RTRAC3
+             ITRAC2=>ITRAC3
+             NMAX2=2*NMAX2
+          ENDIF
+          RTRAC2(CNT)=LENGTH
+          ITRAC2(CNT+NMAX2)=LNEW
+          ITRAC2(CNT)=NOLD
+       ENDDO INTERIOR
+       !
+       !*       exiting motif and analyzing bc condition
+       LGON=NNEW<=G_BC_TYPE(1).AND.NNEW>=G_BC_TYPE(5).AND.(ABS(LENGTH_TOT*LENGTH_INV_CYCL-1.)>EPS1)
+       !        store nber of regions to angtab
+       ANGTAB(2*NB_TOT-1)=CNT-CNT1
+       CNT1=CNT
+       IF(LGON)THEN
+           !     treat boundary condition and get new entry point
+           IF(PRESENT(IDATA_BC2).AND.PRESENT(DIST_AXIS)) THEN
+               IDATA=IDATA_BC2(IBC2_ELEM(LNEW))
+               !     treat bondary condition
+               CALL SAL247_3(BCDATA(:,IDATA))
+               !     at return from SAL247_3:
+               !        AX AY        = entry point at boundary
+               !        EX EY        = new direction
+               !        compute entry point
+               IF(NNEW==G_BC_TYPE(3).OR.NNEW==G_BC_TYPE(4).OR.NNEW==G_BC_TYPE(2)) THEN
+                 P1=PPERIM_MAC2(N_AXIS); P2=PPERIM_MAC2(N_AXIS+1)-1
+                 !        re-entry point is on the axis
+                 CALL SAL241_2(P2-P1+1,PERIM_MAC2(P1:P2),DIST_AXIS(P1:P2),IPAR)
+               ENDIF
+            ELSE
+               CALL XABORT('SAL240_4_2: missing IDATA_BC2 or DIST_AXIS argument')
+            ENDIF
+       ENDIF
+       !
+    ENDDO EXTERIOR
+    !     set success flag
+    IERR=1
+    !
+  END SUBROUTINE SAL240_4_2
+  !
+  SUBROUTINE SAL247_3(BCDATA)
+    !
+    !----------------------------------------------------------------------
+    !
+    !Purpose:
+    ! Treatment of boundary conditions for the cyclic cases
+    !                                                        side 4
+    !          side 3                                        ------
+    !        ----------                  /\          side 3 /      \ side 5
+    !        |        |                 /  \               /        \
+    ! side 2 |        | side 4  side 2 /    \ side 3       \        /
+    !        |        |               /      \       side 2 \      / side 6
+    !        ----------               --------               ------
+    !          side 1                  side 1                side 1
+    !
+    !Parameters: input
+    ! BCDATA  boundary condition descriptor
+    !
+    !----------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : NNEW,DNEW,AX,AY,BX,BY,EX,EY,DELX,N_AXIS
+    USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,EPS,LX=>LENGTHX
+    IMPLICIT NONE
+    REAL(PDB), INTENT(IN), DIMENSION(:) :: BCDATA
+    !***
+    REAL(PDB) :: AUX,ACX,ACY,COSTHE,SINTHE,CX,CY,EX0
+    !***
+    !     boundary point:
+    AX=AX+DNEW*EX
+    AY=AY+DNEW*EY
+    CX=BCDATA(1) ; CY=BCDATA(2)
+    COSTHE=BCDATA(3)
+    SINTHE=BCDATA(4)
+    SELECT CASE(NNEW)
+       CASE(-2,-3)
+       !           translation/rotation for a segment (get displacement data)
+       !           a = displacement vector
+       ! new axis
+       SELECT CASE(TYPGEO)
+       CASE(5)
+          IF (CY<0) THEN ! end axe 2
+             N_AXIS=1
+             DELX=AX
+          ELSEIF (CY>0) THEN ! end axe 1
+             N_AXIS=3
+             DELX=AX
+          ELSEIF (CX<0) THEN ! end axe 4
+             N_AXIS=2
+             DELX=AY
+          ELSEIF (CX>0) THEN ! end axe 3
+             N_AXIS=4
+             DELX=AY
+          ENDIF
+       CASE(9)
+          IF (ABS(BCDATA(1))<EPS) THEN ! axes 1&4
+             IF (CY<0) THEN
+                N_AXIS=1
+             ELSE
+                N_AXIS=4
+             ENDIF
+          ELSEIF (CX>0) THEN
+             IF (CY>0) THEN
+                N_AXIS=5
+             ELSE
+                N_AXIS=6
+             ENDIF
+          ELSEIF (CX<0) THEN
+             IF (CY>0) THEN
+                N_AXIS=3
+             ELSE
+                N_AXIS=2
+             ENDIF
+          ENDIF
+          ACX=AX+CX ; ACY=AY+CY
+          DELX=SQRT((ACX-BX(N_AXIS))**2+(ACY-BY(N_AXIS))**2)
+       CASE(10)
+          ACX=LX-AX
+          EX0=EX
+          IF(CX>0) THEN
+             !               cx>0: axis 3
+             EX=0.5*EX+0.5*SQRT(3.)*EY
+             EY=0.5*EY-0.5*SQRT(3.)*EX0
+             N_AXIS=2
+             DELX=SQRT((AX-CX)*(AX-CX)+AY*AY)
+          ELSEIF(BCDATA(5)>0.) THEN
+             !               cx=0, angle>0: axis 2 (axis of angle pi/3)
+             EX=0.5*EX-0.5*SQRT(3.)*EY
+             EY=0.5*EY+0.5*SQRT(3.)*EX0
+             N_AXIS=3
+             DELX=SQRT(AX*AX+AY*AY)
+          ELSE
+             !               cx=0, angle=0: axis 1 (axis X)
+             EX=-EX
+             EY=-EY
+             N_AXIS=1
+             DELX=ACX
+          ENDIF
+       CASE(11)
+          ACX=AX-CX ; ACY=AY-CY
+          AUX=2._PDB*(ACX*COSTHE+ACY*SINTHE)
+          AX=AUX*COSTHE-ACX+CX
+          AY=AUX*SINTHE-ACY+CY
+          DELX=LX-SQRT((AX-CX)*(AX-CX)+(AY-CY)*(AY-CY))
+          EX0=EX
+          IF(ABS(BCDATA(5))<EPS) THEN ! exit side is on axes 1 or 3
+             EX=-0.5*EX+0.5*SQRT(3.)*EY
+             EY=-0.5*EY-0.5*SQRT(3.)*EX0
+             AX=AX+1.5*DELX
+             AY=AY+0.5*SQRT(3.)*DELX
+             IF(CY>0.) THEN
+                N_AXIS=2 ! cy>0: exit element is on axis 3 (angle=0)
+             ELSE
+                N_AXIS=4 ! cy=0: exit element is on axis 1 (angle=0)
+             ENDIF
+          ELSE ! exit side is on axes 2 or 4
+             EX=-0.5*EX-0.5*SQRT(3.)*EY
+             EY=-0.5*EY+0.5*SQRT(3.)*EX0
+             AX=AX-1.5*(LX-DELX)
+             AY=AY-0.5*SQRT(3.)*(LX-DELX)
+             IF(CX>0.) THEN
+                N_AXIS=1 ! cx>0: exit element is on axis 4 (angle>0)
+             ELSE
+                N_AXIS=3 ! cx=0: exit element is on axis 2 (angle>0)
+             ENDIF
+          ENDIF
+       CASE DEFAULT
+          CALL XABORT('SAL247_3: option not available(1)')
+       END SELECT
+       CASE(-4)
+       !           symmetry with respect an axis (get axis data)
+       !           axis: r=c+t*f (f unit vector of angle theta)
+       ACX=AX-CX ; ACY=AY-CY
+       AUX=2._PDB*(ACX*COSTHE+ACY*SINTHE)
+       AX=AUX*COSTHE-ACX+CX
+       AY=AUX*SINTHE-ACY+CY
+       AUX=2._PDB*(EX*COSTHE+EY*SINTHE)
+       EX=AUX*COSTHE-EX
+       EY=AUX*SINTHE-EY
+       SELECT CASE(TYPGEO)
+       CASE(6)
+          IF(BCDATA(5)>0.) THEN
+             !                 vertical axes
+             IF(CX>0.) THEN
+                N_AXIS=4
+             ELSE
+                N_AXIS=2
+             ENDIF
+             DELX=AY
+          ELSE
+             !                 horizontal axes
+             IF(CY>0.) THEN
+                N_AXIS=3
+             ELSE
+                N_AXIS=1
+             ENDIF
+             DELX=AX
+          ENDIF
+       CASE(7) 
+          IF(CX>0) THEN
+             !                 cx>0: axis 3 (vertical axis)
+             N_AXIS=3
+             DELX=AY
+          ELSEIF(BCDATA(5)>0.) THEN
+             !                 cx=0, angle>0: axis 2 (axis of angle pi/4)
+             N_AXIS=2
+             DELX=SQRT(AX*AX+AY*AY)
+          ELSE
+             !                 cx=0, angle=0: axis 1 (axis x)
+             N_AXIS=1
+             DELX=SQRT(AX*AX+AY*AY)
+          ENDIF
+       CASE(8,12)
+          IF(CX>0) THEN
+             !               cx>0: axis 3
+             N_AXIS=3
+             DELX=SQRT((AX-CX)*(AX-CX)+AY*AY)
+          ELSEIF(BCDATA(5)>0.) THEN
+             !               cx=0, angle>0: axis 2 (axis of angle pi/3 or 2*pi/3)
+             N_AXIS=2
+             DELX=SQRT(AX*AX+AY*AY)
+          ELSE
+             !               cx=0, angle=0: axis 1 (axis X)
+             N_AXIS=1
+             DELX=AX
+          ENDIF
+       CASE DEFAULT
+          CALL XABORT('SAL247_3: option not available(2)')
+       END SELECT
+       CASE DEFAULT
+       CALL XABORT('SAL247_3: option not available(3)')
+    END SELECT
+    !
+  END SUBROUTINE SAL247_3
+  !
+  INTEGER FUNCTION ANGLE_TO_NUMBER(EX,EY,DANGLT)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! search order number of a horizontal angle in the angular quadrature
+    ! formula set
+    !
+    !Parameters: input
+    ! EX       angle cosine
+    ! EY       angle sine
+    ! DANGLT   angle cosines table
+    !
+    !---------------------------------------------------------------------
+    !
+    IMPLICIT NONE
+    REAL(PDB), INTENT(IN)               :: EX,EY
+    REAL(PDB), INTENT(IN), DIMENSION(:,:) :: DANGLT
+    !***
+    INTEGER   :: I,NPHI
+    REAL(PDB) :: EXREF,EYREF
+    CHARACTER(LEN=131) :: HSMG
+    !***
+    NPHI=SIZE(DANGLT,2)
+    ANGLE_TO_NUMBER=0
+    DO I=1,NPHI
+       EXREF=DANGLT(1,I) ; EYREF=DANGLT(2,I) ;
+       IF((ABS(EX-EXREF)<1.E-3).AND.(ABS(EY-EYREF)<1.E-3)) THEN
+         ANGLE_TO_NUMBER=I
+         GO TO 10
+       ELSE IF((ABS(EX-EXREF)<1.E-3).AND.(ABS(EY+EYREF)<1.E-3)) THEN
+         ANGLE_TO_NUMBER=2*NPHI-I+1
+         GO TO 10
+       ENDIF
+    ENDDO
+    WRITE(6,'(/29H ANGLE_TO_NUMBER: QUADRATURE:)')
+    DO I=1,2*NPHI
+       WRITE(6,'(1X,I5,1P,2E12.4)') I,DANGLT(:2,I)
+    ENDDO
+    WRITE(HSMG,'(47HANGLE_TO_NUMBER: UNABLE TO FIND ANGULAR COSINES,1P,2E12.4, &
+    & 26H INTO QUADRATURE SELECTED.)') EX,EY
+    CALL XABORT(HSMG)
+    10 RETURN
+    !
+  END FUNCTION ANGLE_TO_NUMBER
+  !
+  SUBROUTINE SAL247_1(RPAR,IPAR,D,SINPHI,COSPHI,EX,EY)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! computes COSPHI and SINPHI for the intersection of the trajectory
+    ! with element LNEW of descriptors RPAR and IPAR. SINPHI is
+    ! computed with respect to directions outgoing with the trajectory
+    !
+    !Parameters: input
+    ! RPAR         floating point geometry descriptors
+    ! IPAR         integer geometry descriptors
+    ! D            distance measured along the trajectory
+    ! EX           first exiting unit vector in along trajectory
+    ! EY           second exiting unit vector in along trajectory
+    !
+    !Parameters: output
+    ! SINPHI       cosine at intersection
+    ! COSPHI       sine at intersection
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_TRACKING_TYPES, ONLY : AX,AY
+    IMPLICIT NONE
+    INTEGER,   INTENT(IN), DIMENSION(:) :: IPAR
+    REAL(PDB), INTENT(IN), DIMENSION(:) :: RPAR
+    REAL(PDB), INTENT(IN)               :: D,EX,EY
+    REAL(PDB), INTENT(OUT)              :: SINPHI,COSPHI
+    !>    AX AY            = components of origin of trajectory
+    !***
+    INTEGER   :: TYPE
+    REAL(PDB) :: NX,NY
+    !***
+    NX=0._PDB
+    NY=0._PDB
+    TYPE=IPAR(1)
+    SELECT CASE(TYPE)
+       CASE(1)
+       !        TYPE=1=>  segment (s): R=C+T*F with T in (0,1)
+       !        RPAR(1),RPAR(2) = C = (CX,CY)
+       !        RPAR(3),RPAR(4) = F = (FX,FY)
+       !        components of normal
+       NX=RPAR(4)/RPAR(5)
+       NY=-RPAR(3)/RPAR(5)
+       CASE(2,3)
+       !        TYPE=2,3=> arc of circle
+       NX=(AX+D*EX-RPAR(1))/RPAR(3)
+       NY=(AY+D*EY-RPAR(2))/RPAR(3)
+    END SELECT
+    SINPHI=NX*EY-NY*EX
+    IF((EX*NX+EY*NY)<0._PDB) SINPHI=-SINPHI
+    COSPHI=SQRT(1._PDB-SINPHI*SINPHI)
+    !
+  END SUBROUTINE SAL247_1
+  !
+  SUBROUTINE SAL241(PERIM,NPERIM,IPAR,RPAR,DOLD,NOLD,DNEW,NNEW,LNEW)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! computes intersection of trajectory (t): R=A+D*E with a perimeter
+    ! composed of the elements given in array perim
+    !
+    !Parameters: input
+    ! PERIM        elements in the perimeter of a node
+    ! NPERIM       number of elements in the perimeter
+    ! RPAR         floating point geometry descriptors
+    ! IPAR         integer geometry descriptors
+    ! DOLD         value of D at current position
+    ! NOLD         current entered node
+    !
+    !Parameters: output
+    ! DNEW         value of D at the intersection
+    ! NNEW         node index that is enter after intersection
+    ! LNEW         element index that is intersected
+    !
+    !---------------------------------------------------------------------
+    !
+    USE SAL_GEOMETRY_TYPES,  ONLY : NRPAR,NIPAR
+    USE SAL_TRACKING_TYPES,  ONLY : NRPART,NIPART,RPART,IPART,EPS1
+    !***
+    IMPLICIT NONE
+    INTEGER,   INTENT(IN), DIMENSION(:)   :: PERIM
+    INTEGER,   INTENT(IN), DIMENSION(:,:) :: IPAR
+    REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR
+    INTEGER,   INTENT(IN)                 :: NPERIM,NOLD
+    INTEGER,   INTENT(OUT)                :: NNEW,LNEW
+    REAL(PDB), INTENT(IN)                 :: DOLD
+    REAL(PDB), INTENT(OUT)                :: DNEW
+    !***
+    REAL(PDB) :: D
+    INTEGER   :: I,L,INDEX,N,TYPE,NBINTE
+    !     INFTY is used to initialize search for minimum distance
+    REAL(PDB), PARAMETER :: INFTY=1.E+10
+    INTEGER, PARAMETER :: FOUT =6
+    !****
+    !     initialize distance for intersection
+    DNEW=INFTY
+    NBER=0
+    !     intersection:
+    DO I=1,NPERIM
+       L=PERIM(I)
+       !        get order nber of element in the perimeter
+       NBINTE=IPART(1,L)
+       IF(NBINTE<0)THEN
+          ! compute and store intersections
+          TYPE=IPAR(1,L)
+          IF(TYPE==1)THEN
+             !              segment
+             CALL SAL242(RPAR(:,L),IPAR(:,L),RPART(:,L),IPART(2:,L),NBINTE)
+          ELSEIF(TYPE<=3)THEN
+             !              arc of circle or circle
+             CALL SAL243(RPAR(:,L),IPAR(:,L),RPART(:,L),IPART(2:,L),NBINTE)
+          ELSE
+             CALL XABORT('SAL241: Not implemented')
+          ENDIF
+          IPART(1,L)=NBINTE
+       ENDIF
+       !
+       IF(NBINTE/=0)THEN
+          DO INDEX=1,NBINTE
+             D=RPART(INDEX,L)
+             N=IPART(1+INDEX,L)
+             ! analyzes feasability of intersection to eliminate conflicts due
+             ! to concavities and crossing of mesh points
+             IF(N/=NOLD.AND.D>(DOLD-EPS1))THEN
+                NBER=NBER+1
+                IF(D<(DNEW-EPS1))THEN
+                   DNEW=D
+                   COSINE=RPART(INDEX+2,L)
+                   NNEW=N
+                   LNEW=L
+                ELSEIF(D<(DNEW+EPS1))THEN
+                   ! case of two close intersections : liquidate smaller
+                   NBER=NBER-1
+                   IF(D>DNEW)THEN
+                      DNEW=D
+                      COSINE=RPART(INDEX+2,L)
+                      NNEW=N
+                      LNEW=L
+                   ENDIF
+                ENDIF
+             ENDIF
+          ENDDO
+       ENDIF
+    ENDDO
+    !
+    LGOK=DNEW/=INFTY
+    IF(LGOK)THEN
+       ! eliminate intersection
+       IPART(1,LNEW)=IPART(1,LNEW)-1
+       IF(IPART(1,LNEW)==1.AND.RPART(1,LNEW)==DNEW)THEN
+          ! FOR AN ELEMENT WITH TWO INTERSECTIONS, MOVE 2ND
+          ! INTERSECTION INTO FIRST IF FIRST HAS BEEN TAKEN
+          RPART(1,LNEW)=RPART(2,LNEW)
+          RPART(3,LNEW)=RPART(4,LNEW)
+          IPART(2,LNEW)=IPART(3,LNEW)
+       ENDIF
+    ELSE
+       ! print out problem
+       IF(NOLD/=-100) THEN
+          WRITE(FOUT,*)'Problem in SAL241: NOLD, DOLD, NBINTE ',NOLD,DOLD,NBINTE
+       ENDIF
+    ENDIF
+  END SUBROUTINE SAL241
+  !
+  SUBROUTINE SAL242(RPAR,IPAR,D,NODE,NBINTE)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! analysis of the intersection of the trajectory (T):R=A+D*E and
+    ! a segment
+    !
+    !Parameters: input
+    ! RPAR         floating point geometry descriptors
+    ! IPAR         integer geometry descriptors
+    !
+    !Parameters: output
+    ! D            value of D at intersection
+    ! NODE         order nber of node entered after intersection
+    ! NBINTE       number of intersections
+    !
+    !---------------------------------------------------------------------
+    !
+    IMPLICIT NONE
+    INTEGER,   INTENT(IN),  DIMENSION(:) :: IPAR
+    INTEGER,   INTENT(OUT), DIMENSION(:) :: NODE
+    INTEGER,   INTENT(OUT) :: NBINTE
+    REAL(PDB), INTENT(IN),  DIMENSION(:) :: RPAR
+    REAL(PDB), INTENT(OUT), DIMENSION(:) :: D
+    !     DIMENSION        RPAR(*),IPAR(*),NODE(*),D(*)
+    !***
+    REAL(PDB) :: CAX,CAY,FX,FY,A,DELTAM,DELTA
+    REAL(PDB), PARAMETER :: EPS2=0.
+    !***
+    !     TYPE=1=>  segment (S): R=C+T*F with T in (0,1)
+    !     RPAR(1),RPAR(2) = C = (CX,CY)
+    !     RPAR(3),RPAR(4) = F = (FX,FY)
+    !     components of vector F
+    FX=RPAR(3)
+    FY=RPAR(4)
+    !     DELTA=F X E
+    DELTA=FX*EY-FY*EX
+    IF(DELTA/=0._PDB)THEN
+       !        components of vector CA=C-A
+       CAX=RPAR(1)-AX
+       CAY=RPAR(2)-AY
+       !        A=AC X E
+       A=CAY*EX-CAX*EY
+       DELTAM=DELTA-A
+       IF(DELTAM>(-EPS2).AND.A>(-EPS2))THEN
+          !           crossing into + halfspace
+          NODE(1)=IPAR(3)
+       ELSEIF(DELTAM<=EPS2.AND.A<=EPS2)THEN
+          !           crossing into - halfspace
+          NODE(1)=IPAR(2)
+       ELSE
+          !           out-of-range crossing
+          NBINTE=0
+          RETURN
+       ENDIF
+       !        compute distance to intersection along trajectory
+       !        D = (AC X F)/DELTA
+       D(1)=(CAY*FX-CAX*FY)/DELTA
+       D(3)=ABS(DELTA)
+       NBINTE=1
+    ELSE
+       !        A/=0 => out-of-range crossing (infinity)
+       !        A==0 => trajectory coincides with segment
+       !        in any case neglect intersection
+       NBINTE=0
+    ENDIF
+  END SUBROUTINE SAL242
+  !
+  SUBROUTINE SAL243(RPAR,IPAR,D,NODE,NBINTE)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! analysis of the intersection of the trajectory (T):R=A+D*E and
+    ! a circle or arc of circle
+    !
+    !Parameters: input
+    ! RPAR         floating point geometry descriptors
+    ! IPAR         integer geometry descriptors
+    !
+    !Parameters: output
+    ! D            value of D at intersection
+    ! NODE         order nber of node entered after intersection
+    ! NBINTE       number of intersections
+    !
+    !---------------------------------------------------------------------
+    !
+    IMPLICIT NONE
+    INTEGER,   INTENT(IN),  DIMENSION(:) :: IPAR
+    INTEGER,   INTENT(OUT), DIMENSION(:) :: NODE
+    INTEGER,   INTENT(OUT) :: NBINTE
+    REAL(PDB), INTENT(IN),  DIMENSION(:) :: RPAR
+    REAL(PDB), INTENT(OUT), DIMENSION(:) :: D
+    !     DIMENSION        RPAR(*),IPAR(*),NODE(*),D(*)
+    !***
+    INTEGER   :: TYPE,I
+    REAL(PDB) :: CAX,CAY,RAD,RAD2,RHOMI2,DMIN,THETA,THETA1,THETA2, &
+         COSTHE,DELTA
+    REAL(PDB), PARAMETER :: EPS2=0._PDB
+    !***
+    !     TYPE=2,3=> arc of circle (C): R=C+R*F(THETA),
+    !                THETA in (THETA1,THETA2)
+    !     RPAR(1),RPAR(2) = C = (CX,CY)
+    !     RPAR(3)         = R = RADIUS
+    !     RPAR(4),RPAR(5) = (THETA1,THETA2) in (0,2PI) with THETA1<THETA2
+    !                        and THETA1<0 if arc crosses THETA=0
+    !     LGTYPE = .TRUE.  => THETA1 > 0
+    !            = .FALSE. => THETA1 < 0
+    !     components of vector CA=C-A
+    CAX=RPAR(1)-AX
+    CAY=RPAR(2)-AY
+    TYPE=IPAR(1)
+    !     value of R2
+    RAD=RPAR(3)
+    RAD2=RAD**2
+    !     RHOMI2=(CA X E)**2
+    RHOMI2=(CAX*EY-CAY*EX)**2
+    IF(RAD2>=RHOMI2)THEN
+       DELTA=SQRT(RAD2-RHOMI2)
+       !        tangent point = two very close points (to avoid infinite loop)
+       if(delta==0._pdb) delta=small
+       !        DMIN=CA*E
+       DMIN=CAX*EX+CAY*EY
+       !        D(1) D(2) = min and max distances for the two intersections
+       D(1)=DMIN-DELTA
+       D(2)=DMIN+DELTA
+       D(3)=DELTA/RAD
+       D(4)=D(3)
+       IF(TYPE==2)THEN
+          !           full circle. both intersections are possible
+          NODE(1)=IPAR(2)
+          NODE(2)=IPAR(3)
+          NBINTE=2
+       ELSE
+          !           analysis for arc of circle:
+          THETA1=RPAR(4)
+          THETA2=RPAR(5)
+          NBINTE=0
+          !           compute angles for closest and farthest intersections
+          !           and check feasability
+          DO I=1,2
+             COSTHE=(D(I)*EX-CAX)/RAD
+             THETA=SALACO(COSTHE,D(I)*EY-CAY)
+             IF( (((THETA1-THETA)<EPS2).AND.((THETA-THETA2)<EPS2)).OR. &
+                 ((THETA-THETA2)<(EPS2-TWOPI)) )THEN
+                NBINTE=NBINTE+1
+                IF(NBINTE/=I)D(NBINTE)=D(I)
+                NODE(NBINTE)=IPAR(I+1)
+             END IF
+          ENDDO
+       ENDIF
+    ELSE
+       NBINTE=0
+    ENDIF
+    !
+  END SUBROUTINE SAL243
+  !
+END MODULE SAL_TRAJECTORY_MOD