path: root/Dragon/src/SAL_GEOMETRY_MOD.f90

!
!---------------------------------------------------------------------
!
!Purpose:
! Support module used to create a 2D geometry.
!
!Copyright:
! Copyright (C) 2001 Ecole Polytechnique de Montreal.
!
!Author(s):
! X. Warin
!
!---------------------------------------------------------------------
!
MODULE SAL_GEOMETRY_MOD

  USE SAL_GEOMETRY_TYPES
  USE PRECISION_AND_KINDS, ONLY : PDB, PI,TWOPI,HALFPI
  USE SAL_NUMERIC_MOD,    ONLY : SAL141
  USE SALGET_FUNS_MOD

CONTAINS

  SUBROUTINE SAL100(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! perform input data and first allocation for geometry OBJECT T_G_BASIC
    !
    !Parameters: input/output
    ! GG    geometry basic information.
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC,TYPGEO,NBFOLD,EPS,ANGGEO, &
         INDEX,KNDEX,PREC,ISPEC,LGSPEC
    USE SAL_TRACKING_TYPES, ONLY : PRTIND
    IMPLICIT NONE
    TYPE(T_G_BASIC) :: GG
    !*****
    ! local variable
    ! **************
    INTEGER, PARAMETER :: N_DATAIN=25, N_DATARE=20
    INTEGER, DIMENSION (N_DATAIN) :: DATAIN
    REAL,    DIMENSION (N_DATARE) :: DATARE
    INTEGER :: OK
    INTEGER, PARAMETER :: FOUT =6
    !*****
    !     TYPGEO  = type of geometry:
    !     NBFOLD  = n in angle definition of rotation or symmetry geometry
    !     NB_NODE = number of nodes
    !     NBELEM  = number of elements
    !
    CALL SALGET(DATAIN,6,F_GEO,FOUT0,'dimensions for geometry')
    TYPGEO=DATAIN(1)
    NBFOLD=DATAIN(2)
    GG%NB_NODE=DATAIN(3)
    GG%NB_ELEM=DATAIN(4)
    GG%NB_MACRO=DATAIN(5)
    GG%NB_FLUX=DATAIN(6)
    ANGGEO=0.0D0
    SELECT CASE(TYPGEO)
       CASE(1)
       IF(NBFOLD.EQ.3) THEN
         ANGGEO=TWOPI/8.0D0 ! half diagonal symmetry case
       ELSE
         ANGGEO=TWOPI/NBFOLD
       ENDIF
       CASE(2)
       ANGGEO=TWOPI/NBFOLD
       CASE(5:6)
       ANGGEO=HALFPI
       CASE(7)
       ANGGEO=PI*0.25
       CASE(8:11)
       ANGGEO=PI/3.
       CASE(12)
       ANGGEO=PI/6.
    END SELECT
    IF((PRTIND >= 1).AND.(TYPGEO.EQ.0)) THEN
      WRITE(FOUT,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD
    ELSE IF(PRTIND >= 1) THEN
      WRITE(FOUT,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD, &
      & ' ANGGEO=',ANGGEO,' radians'
    ENDIF
    LGSPEC=(TYPGEO/=0).AND.(NBFOLD==0)
    IF(LGSPEC) THEN
      IF(ISPEC==0) THEN
        WRITE(*,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD
        CALL XABORT('SAL100: TISO option is incompatible with the surfacic file')
      ENDIF
    ELSE
      IF(ISPEC==1) THEN
        WRITE(*,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD
        CALL XABORT('SAL100: TSPC option is incompatible with the surfacic file')
      ENDIF
    ENDIF
    !
    !*    read printing indexes for general domain data and topological deformations
    !     INDEX  =  to print general domain data
    !     KNDEX  =  to print motions of topological adjustment
    !     PREC   =  if 0 then read RPAR & BCDATA with e20.12)
    CALL SALGET(DATAIN,3,F_GEO,FOUT0,'index kndex prec')
    INDEX=DATAIN(1)  
    KNDEX=DATAIN(2) 
    PREC=DATAIN(3) 
    !
    !*    read epsilons for topological deformations
    !     EPS    = if the distance of two ends of elements < eps, 
    !              they will be united to one point
    CALL SALGET(DATARE,1,F_GEO,FOUT0,'eps')
    EPS=DATARE(1)
    !
    IF(PRTIND >= 1) THEN
        WRITE(FOUT,'(//,5X,''domain checkout:'',/, &
        & 5X,''elements are in contact for distance < '',1P,E12.4,//)') EPS
    ENDIF
    !
    ALLOCATE(GG%NUM_MERGE(GG%NB_NODE),STAT =OK)
    IF(OK /= 0) CALL XABORT('SAL100: failure to allocate NB_NODE')
    CALL SALGET(GG%NUM_MERGE,GG%NB_NODE,F_GEO,FOUT0,'FLUX INDEX PER NODE')
    IF(MAXVAL(GG%NUM_MERGE) /= GG%NB_FLUX) CALL XABORT('SAL100: inconsistent NBFLUX')
    ALLOCATE(GG%NAME_MACRO(GG%NB_MACRO),STAT =OK)
    IF(OK /= 0) CALL XABORT('SAL100: failure to allocate NB_MACRO')
    CALL SALGET(GG%NAME_MACRO,GG%NB_MACRO,F_GEO,FOUT0,'NAMES OF MACROS')
    ALLOCATE(GG%NUM_MACRO(GG%NB_FLUX))
    CALL SALGET(GG%NUM_MACRO,GG%NB_FLUX,F_GEO,FOUT0,'macro order number per flux region')
    !*    do the work (SAL100_2 is called here!):
    CALL SAL110(GG)
    !
  END SUBROUTINE SAL100
  !
  SUBROUTINE SAL110(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! constructs geometrical domain
    !
    !Parameters: input/output
    ! GG    geometry basic information.
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,NBFOLD,NIPAR,NRPAR,ALLSUR,NANIS,IC,ISPEC,ANGGEO
    USE SAL_TRACKING_TYPES, ONLY : PRTIND
    !****
    IMPLICIT NONE
    ! in variable
    ! ************
    TYPE(T_G_BASIC), INTENT(INOUT) :: GG

    ! local variable
    ! ***************
    INTEGER :: ELEM,OK,I,TYPE
    CHARACTER(LEN=4) :: HSYM
    REAL(PDB),PARAMETER :: CONV=PI/180._PDB
    INTEGER, PARAMETER :: FOUT =6
    !****
    ! allocate node arrays
    HSYM=' '
    IF((TYPGEO==0).OR.((TYPGEO==6).AND.(NBFOLD==0))) THEN
      CONTINUE
    ELSE IF((TYPGEO==1).AND.(NBFOLD==3)) THEN
      HSYM='QUAR'
    ELSE IF(((TYPGEO==1).AND.(NBFOLD==8)).OR.((TYPGEO==7).AND.(NBFOLD==0))) THEN
      HSYM='EIGH'
    ELSE IF(((TYPGEO==1).AND.(NBFOLD==6)).OR.((TYPGEO==8).AND.(NBFOLD==0))) THEN
      HSYM='SIXT'
    ELSE IF(((TYPGEO==1).AND.(NBFOLD==12)).OR.((TYPGEO==12).AND.(NBFOLD==0))) THEN
      HSYM='S30'
    ELSE IF(((TYPGEO==1).AND.(NBFOLD==4)).OR.((TYPGEO==3).AND.(NBFOLD==0))) THEN
      HSYM='SYME'
    ELSE IF((TYPGEO==1).AND.(NBFOLD==2)) THEN
      HSYM='SY1D'
    ELSE IF(((TYPGEO==2).AND.(NBFOLD==6)).OR.((TYPGEO==10).AND.(NBFOLD==0))) THEN
      HSYM='RA60'
    ELSE IF((TYPGEO==5).AND.(NBFOLD==0)) THEN
      HSYM='TRAN'
    ELSE IF((TYPGEO==9).AND.(NBFOLD==0)) THEN
      HSYM='TRAN'
    ELSE IF(((TYPGEO==2).AND.(NBFOLD==3)).OR.(TYPGEO==11).AND.(NBFOLD==0)) THEN
      HSYM='R120'
    ELSE
      WRITE(*,*) "TYPGEO=",TYPGEO," NBFOLD=",NBFOLD
      CALL XABORT('SAL110: non supported type of symmetry')
    ENDIF
    ALLOCATE (GG%IPAR(NIPAR,GG%NB_ELEM), GG%RPAR(NRPAR,GG%NB_ELEM), STAT=OK)
    IF(OK/=0) CALL XABORT('SAL110: not enough memory I,R')

    !*    read surfacic file
    CALL SALINP(GG)
    !
    !*    unite the ends of elements, redefine elements
    CALL SAL128(GG%RPAR,GG%IPAR,GG%NB_ELEM)
    !
    IF((ISPEC==0).AND.(IC.EQ.4)) THEN
      !*    unfold domain
      IF(HSYM=='QUAR') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: DIAG unfold"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        ANGGEO=2.0*ANGGEO
        CALL SALFOLD_1('DIAG',GG)
      ELSE IF(HSYM=='EIGH') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: DIAG + SYME unfold"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        ANGGEO=2.0*ANGGEO
        CALL SALFOLD_1('DIAG',GG)
        CALL SALFOLD_1('SYMX',GG)
        CALL SALFOLD_1('SYMY',GG)
      ELSE IF(HSYM=='SYME') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: SYME unfold"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        CALL SALFOLD_1('SYMX',GG)
        CALL SALFOLD_1('SYMY',GG)
      ELSE IF(HSYM=='SY1D') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: SY1D unfold"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        CALL SALFOLD_1('SYMX',GG)
      ELSE IF(HSYM=='SIXT') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: SA60 unfold"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        CALL SALFOLD_1('SA60',GG)
        CALL SALFOLD_1('SYMH',GG)
      ELSE IF(HSYM=='S30') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: S30 unfold"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        CALL SALFOLD_1('S30 ',GG)
        CALL SALFOLD_1('SB60',GG)
        CALL SALFOLD_1('SYMX',GG)
      ELSE IF(HSYM=='RA60') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: SR60 unfold with rotation"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        CALL SALFOLD_2('SR60',GG)
        CALL SALFOLD_2('R180',GG)
      ELSE IF(HSYM=='R120') THEN
        IF(PRTIND>0) WRITE(*,*) "SAL110: SR120 unfold with rotation"
        IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1')
        CALL SALFOLD_2('R120',GG)
      ELSE
        IF(PRTIND>0) WRITE(*,*) "SAL110: no unfold"
      ENDIF
      IF((TYPGEO/=0).AND.(NBFOLD==0)) THEN
        TYPGEO=6
      ELSE
        TYPGEO=0; NBFOLD=0
      ENDIF
    ENDIF
    IF(PRTIND>0) WRITE(FOUT,*) 'SAL110: after unfolding -- NB_ELEM=',GG%NB_ELEM, &
    & ' NB_PERIM=',GG%NBBCDA

    IF(PRTIND>5) THEN
      !*    print surfacic file
      WRITE(FOUT,'(5H--cut,75(1H-))')
      WRITE(FOUT,'(5HBEGIN)')
      WRITE(FOUT,'(42H* typgeo nbfold nbnode nbelem nbmacr nbreg)')
      WRITE(FOUT,'(6I7)') TYPGEO,NBFOLD,GG%NB_NODE,GG%NB_ELEM,GG%NB_MACRO,GG%NB_NODE
      WRITE(FOUT,'(20H* index  kndex  prec)')
      WRITE(FOUT,'(4I7)') 0,0,1
      WRITE(FOUT,'(18H* eps         eps0)')
      WRITE(FOUT,'(1P,2E18.9)') 1.0E-03,1.0E-05
      WRITE(FOUT,'(20H* num_of_region/mesh)')
      WRITE(FOUT,'(10I7)') (GG%NUM_MERGE(I),I=1,GG%NB_NODE)
      WRITE(FOUT,'(13H* macro names)')
      WRITE(FOUT,'(4(3x,a10,2x))') (GG%NAME_MACRO(I),I=1,GG%NB_MACRO)
      WRITE(FOUT,'(35H* macro_order_index_per_flux_region)')
      WRITE(FOUT,'(10I7)') (GG%NUM_MACRO(I),I=1,GG%NB_FLUX)
      DO ELEM=1,GG%NB_ELEM
        TYPE=GG%IPAR(1,ELEM)
        WRITE(FOUT,'(7h elem =,I6)') ELEM
        WRITE(FOUT,'(22H*type    node-   node+)')
        WRITE(FOUT,'(3I6)') (GG%IPAR(I,ELEM),I=1,3)
        WRITE(FOUT,'(63H*cx            cy            ex_or_R       ey_or_theta1  theta2)')
        IF(TYPE<=2) THEN
          WRITE(FOUT,'(1P,5E18.9)') (GG%RPAR(I,ELEM),I=1,5)
        ELSE IF(TYPE==3) THEN
          WRITE(FOUT,'(1P,5E18.9)') (GG%RPAR(I,ELEM),I=1,3),GG%RPAR(4,ELEM)/CONV, &
                                 (GG%RPAR(5,ELEM)-GG%RPAR(4,ELEM))/CONV
        ENDIF
      ENDDO
      WRITE(FOUT,'(40H*defaul  nbbcda  allsur  divsur  ndivsur)')
      WRITE(FOUT,'(1P,5I8)') GG%DEFAUL,GG%NBBCDA,ALLSUR,0,0
      WRITE(FOUT,'(17H*albedo  deltasur)')
      WRITE(FOUT,'(1P,2E18.9)') GG%ALBEDO,0.0
      DO ELEM=1,GG%NBBCDA
        WRITE(FOUT,'(37H particular boundary condition number,i12)') ELEM
        WRITE(FOUT,'(13H*type    nber)')
        WRITE(FOUT,'(1P,2I8)') GG%BCDATAREAD(ELEM)%SALTYPE,GG%BCDATAREAD(ELEM)%NBER
        WRITE(FOUT,'(14H*elems(1,nber))')
        WRITE(FOUT,'(1P,10I8)') (GG%BCDATAREAD(ELEM)%ELEMNB(I),I=1,GG%BCDATAREAD(ELEM)%NBER)
        IF(GG%BCDATAREAD(ELEM)%SALTYPE==0) THEN
          WRITE(FOUT,'(7H*albedo)')
          WRITE(FOUT,'(1P,E18.9)') GG%BCDATAREAD(ELEM)%BCDATA(6)
        ELSE
          WRITE(FOUT,'(22H*cx      cy      angle)')
          WRITE(FOUT,'(1P,3E18.9)') GG%BCDATAREAD(ELEM)%BCDATA(1:2),GG%BCDATAREAD(ELEM)%BCDATA(5)*180._PDB/PI
        ENDIF
      ENDDO
    ENDIF

    ! allocate media and element arrays 
    ALLOCATE (GG%VOL_NODE(GG%NB_NODE),GG%PPERIM_NODE(GG%NB_NODE+1),GG%IBC2_ELEM(GG%NB_ELEM), &
         GG%ISURF2_ELEM(GG%NB_ELEM),GG%MED(GG%NB_NODE), STAT=OK)
    IF(OK/=0) CALL XABORT('SAL110: not enough memory VOL')
    GG%ISURF2_ELEM(:GG%NB_ELEM)=0

    !*    2D boundary conditions and macro contacts:
    !     - defines NB_BC2, NBSUR2
    !     - defines surface strctures for each 2D macro
    !     - defines perimeter structure for each 2D macro
    !     - read 2D boundary conditions
    CALL SAL130(GG)
    !
    !*    topological check
    CALL SAL140(GG%NB_NODE,GG%RPAR,GG%IPAR,GG%PPERIM_NODE,GG%PERIM_NODE)
    !
    !*    volumes, surfaces, put local nbers in node, and read media:
    CALL SAL160(GG)
    IF(PRTIND>5) THEN
      WRITE(FOUT,'(12H* mil(nbreg))')
      WRITE(FOUT,'(10I7)') (GG%MED(I),I=1,GG%NB_NODE)
      WRITE(FOUT,'(3HEND)')
      WRITE(FOUT,'(5H--cut,75(1H-))')
    ENDIF
    !
    !*    printout basic domain
    CALL SAL170(GG)
    !
  END SUBROUTINE SAL110
  !
  SUBROUTINE SALINP(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! read surfacic file
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : ALLSUR,PREC
    !****
    IMPLICIT NONE
    ! in variable
    ! ************
    TYPE(T_G_BASIC), INTENT(INOUT) :: GG
    !****
    INTEGER, PARAMETER :: N_DATAIN=25
    INTEGER, DIMENSION (N_DATAIN) :: DATAIN
    INTEGER   :: ELEM,I,TYPE,NBER
    INTEGER, PARAMETER, DIMENSION(0:4) :: READ_BC_LEN=(/1,1,3,3,3/)
    INTEGER, PARAMETER :: FOUT =6
    REAL(PDB) :: ANGLE,BCDATA_TDT(3)
    !
    !*    read element data
    DO ELEM=1,GG%NB_ELEM
       CALL SAL126(GG%RPAR(:,ELEM),GG%IPAR(:,ELEM))
    ENDDO
    !*    read
    !     DEFAUL = default bc condition
    !               0 = all external surfaces are isotropically reflected
    !                   with the same ALBEDO = ALBED0
    !               1 = all external surfaces are reflected
    !     NBBCDA = number of groups of bc conditions
    !
    CALL SALGET(DATAIN,3,F_GEO,FOUT0,'general bc data')
    GG%DEFAUL=DATAIN(1)
    GG%NBBCDA=DATAIN(2)
    ALLSUR=DATAIN(3)
    !     we can define only two default bc's
    IF(GG%DEFAUL>4.OR.GG%DEFAUL<0) THEN
         WRITE(FOUT,'(8H defaul=,I5)') GG%DEFAUL
         CALL XABORT('SALINP: wrong default bc type')
    ENDIF
    !*    read albedo : defaul bcdata
    CALL SALGET(GG%ALBEDO,F_GEO,FOUT0,PREC,'GENERAL ALBEDO')
    !
    !*    read detailed bcdata if required (motions)
    LBCDIAG=.FALSE.
    IF(GG%NBBCDA>0)THEN
       ALLOCATE(GG%BCDATAREAD(GG%NBBCDA))
       DO I=1,GG%NBBCDA
          GG%BCDATAREAD(I)%BCDATA(:6)=0.0
          CALL SALGET(DATAIN,2,F_GEO,FOUT0,'SPECIFIC BC: TYPE NBER')
          TYPE=DATAIN(1)
          NBER=DATAIN(2)
          GG%BCDATAREAD(I)%SALTYPE=TYPE
          GG%BCDATAREAD(I)%NBER=NBER
          IF(TYPE>5.OR.TYPE<0) CALL XABORT('SALINP: wrong bc type')
          IF(NBER>GG%NB_ELEM) CALL XABORT('SALINP: bc def exceeds nber of elements')
          !
          !           read order numbers of elements affected
          ALLOCATE(GG%BCDATAREAD(I)%ELEMNB(NBER))
          CALL SALGET(GG%BCDATAREAD(I)%ELEMNB,NBER,F_GEO,FOUT0,'BC ELEMENTS')
          !           read bc motion
          CALL SALGET(BCDATA_TDT,READ_BC_LEN(TYPE),F_GEO,FOUT0,PREC,'data for specific bc condition')
          IF(READ_BC_LEN(TYPE).EQ.1) THEN
            GG%BCDATAREAD(I)%BCDATA(1:5)=0._PDB
            GG%BCDATAREAD(I)%BCDATA(6)=BCDATA_TDT(1)
          ELSE
            GG%BCDATAREAD(I)%BCDATA(1:2)=BCDATA_TDT(1:2)
            ANGLE=BCDATA_TDT(3)*PI/180._PDB
            GG%BCDATAREAD(I)%BCDATA(3)=COS(ANGLE)
            GG%BCDATAREAD(I)%BCDATA(4)=SIN(ANGLE)
            GG%BCDATAREAD(I)%BCDATA(5)=ANGLE
            GG%BCDATAREAD(I)%BCDATA(6)=1._PDB
            LBCDIAG=LBCDIAG.OR.((GG%BCDATAREAD(I)%BCDATA(1)==0._PDB).AND.(GG%BCDATAREAD(I)%BCDATA(2)==0._PDB) &
                           .AND.(GG%BCDATAREAD(I)%BCDATA(5)==PI/4._PDB))
          ENDIF
       ENDDO
    ENDIF
  END SUBROUTINE SALINP
  !
  SUBROUTINE SAL126(RPAR,IPAR)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! reads one element data
    !
    !Parameters: input/output
    ! RPAR  floating point geometry descriptors
    ! IPAR  integer geometry descriptors
    !
    !---------------------------------------------------------------------
    !
    !****
    USE SAL_GEOMETRY_TYPES,  ONLY : PREC,G_ELE_TYPE
    !****
    IMPLICIT NONE
    ! in variable
    INTEGER,   INTENT(OUT), DIMENSION(:) :: IPAR
    REAL(PDB), INTENT(INOUT), DIMENSION(:) :: RPAR
    !****
    ! local variable
    INTEGER   :: TYPE,NBER,IAUX
    REAL(PDB) :: PHI1,PHI2,ANGMAX,DELPHI
    !****
    !     IPAR(1)          = 1 (segment), 2 (circle), 3(arc of circle)
    !     IPAR(2)          = order nber of node in - side
    !     IPAR(3)          = order nber of node in + side
    !
    !*    read integer element data
    CALL SALGET(IPAR,3,F_GEO,FOUT0,'integer descriptors')
    TYPE=IPAR(1)
    !
    !     ** segment:
    !        R = C+T*F  with T in (0,1)
    !        RPAR(1~5)  = CX CY FX FY F
    !
    !     ** arc of circle:
    !        R = C+R*F(THETA) with
    !            THETA in (THETA1 < THETA2)
    !            with THETA1 in (0,360)
    !        RPAR(1~5) = CX CY R THETA1 THETA2 (in degrees)
    !
    !*    read real element data

    IF(TYPE == G_ELE_TYPE(1)) THEN
      NBER=4
    ELSE IF(TYPE == G_ELE_TYPE(2)) THEN
      NBER=3
    ELSE IF(TYPE == G_ELE_TYPE(3)) THEN
      NBER=5
      ANGMAX=360._PDB
    ELSE
      WRITE(FOUT0,'(1X,''==> SAL126: unknown type '',I3)') TYPE
      CALL XABORT('SAL126: unknown element type')
    ENDIF

    CALL SALGET(RPAR,NBER,F_GEO,FOUT0,PREC,'real descriptors')

    IF(TYPE == G_ELE_TYPE(1)) THEN
      !        segment: compute length
      RPAR(5)=SQRT(RPAR(3)*RPAR(3)+RPAR(4)*RPAR(4))
      RPAR(6)=0._PDB
    ELSE IF(TYPE == G_ELE_TYPE(2)) THEN
      !        full circle: set angles
      RPAR(4)=0._PDB
      RPAR(5)=TWOPI
      RPAR(6)=0._PDB
    ELSE IF((TYPE == G_ELE_TYPE(3)).OR.(TYPE == G_ELE_TYPE(4))) THEN
      !        check angles
      PHI1=RPAR(NBER-1)
      DELPHI=RPAR(NBER)
      !        order angles in increasing values:
      IF(DELPHI>0._PDB)THEN
         IF(DELPHI>ANGMAX)THEN
            WRITE(FOUT0,'(1X,''==> SAL126: DELPHI = '',1P,E12.4, &
                 &'' > '',1P,E12.4,''  FOR TYPE'',I3)')DELPHI,ANGMAX,TYPE
            CALL XABORT('SAL126: invalid value of delphi')
         ENDIF
         PHI2=PHI1+DELPHI
      ELSE
         IF(DELPHI<-ANGMAX)THEN
            WRITE(FOUT0,'(1X,''==> SAL126: DELPHI = '',1P,E12.4, &
                 &'' < '',1P,E12.4,''  FOR TYPE'',I3)')DELPHI,-ANGMAX,TYPE
            CALL XABORT('SAL126: invalid value of delphi')
         ENDIF
         PHI2=PHI1
         PHI1=PHI1+DELPHI
      ENDIF
      IF(TYPE==G_ELE_TYPE(3))THEN
         !           arc of circle: put phi1 within 0 and 360.
         IF(PHI1>360._PDB)THEN
            IAUX=INT(PHI1/360._PDB)
            DELPHI=360._PDB*IAUX
            PHI2=PHI2-DELPHI
            PHI1=PHI1-DELPHI
         ELSEIF(PHI1<0._PDB)THEN
            IAUX=INT((-PHI1+1.D-7)/360._PDB)+1
            DELPHI=360._PDB*IAUX
            PHI2=PHI2+DELPHI
            PHI1=PHI1+DELPHI
         ENDIF
         RPAR(6)=0._PDB
      ELSE
         CALL XABORT('SAL126: unsupported option')
      ENDIF
      !        convert to radians
      RPAR(NBER-1)=PHI1*(TWOPI/360._PDB)
      RPAR(NBER)=PHI2*(TWOPI/360._PDB)
    ENDIF
    !
  END SUBROUTINE SAL126
  !
  SUBROUTINE SAL128(RPAR,IPAR,NB_ELEM)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! create library of ends of elements and redefine element discriptors
    !
    !Parameters: input/output
    ! RPAR     floating point geometry descriptors
    ! IPAR     integer geometry descriptors
    ! NB_ELEM  number of surfacic elements
    !
    !---------------------------------------------------------------------
    !
    !****
    USE SAL_GEOMETRY_TYPES, ONLY : G_ELE_TYPE,EPS,KNDEX,TYPGEO
    ! in variable
    !************
    IMPLICIT NONE
    INTEGER,   INTENT(IN), DIMENSION(:,:)     :: IPAR
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:)  :: RPAR
    INTEGER , INTENT(IN) :: NB_ELEM
    ! local variable
    ! ***************
    REAL(PDB), DIMENSION (:,:), ALLOCATABLE :: POINT    ! coordinates of the ends of elements
    INTEGER, DIMENSION (:,:), ALLOCATABLE  :: ELEM_NEW  ! point nbers of two ends of elements
    !****
    !     NN               = total number of points
    !     NBP              = count of number of ends having the same coordinates
    INTEGER   :: TYPE,NN,ELEM,I,P,POS,I1,I2,OK
    REAL(PDB) :: X(2),Y(2),D,DX,DY
    REAL      :: EPS2
    LOGICAL   :: LGONE
    INTEGER, DIMENSION(NB_ELEM)     :: NBP
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     - POINT(2,NB_ELEM)    = coordinates of the ends of elements
    !     - ELEM_NEW(2,NB_ELEM) = point nbers of two ends of elements
    !
    ALLOCATE(POINT(2,NB_ELEM),ELEM_NEW(2,NB_ELEM),STAT=OK)
    IF(OK/=0) CALL XABORT('SAL100_3_1: not enough memory I,CH')

    EPS2=EPS*EPS
    IF(KNDEX/=0)WRITE(FOUT,'(//,&
         &" reunite the ends of elements of distance less than ",E13.4)') EPS
    !
    !*    nn counts the nber of points
    NBP(1:NB_ELEM)=0
    NN=0
    DO ELEM=1,NB_ELEM
       TYPE=IPAR(1,ELEM)
       IF(TYPE/=G_ELE_TYPE(2)) THEN
          DO I=1,2
             !              get coordinates of end
             CALL SAL141(TYPE,RPAR(:,ELEM),X(I),Y(I),I)
             !              find if there is a defined neighbouring point
             LGONE=.FALSE.
             DO P=1,NN
                DX=POINT(1,P)-X(I)
                DY=POINT(2,P)-Y(I)
                D=DX*DX+DY*DY
                IF(D<EPS2) THEN
                   LGONE=.TRUE.
                   POS=P
                   EXIT
                ENDIF
             ENDDO
             IF(LGONE) THEN
                !                 choose average value
                X(I)=(X(I)+POINT(1,POS)*NBP(POS))/REAL(NBP(POS)+1)
                Y(I)=(Y(I)+POINT(2,POS)*NBP(POS))/REAL(NBP(POS)+1)
             ELSE
                !                 add points
                NN=NN+1
                IF(NN > NB_ELEM) CALL XABORT('SAL128: point overflow')
                POS=NN
             ENDIF
             POINT(1,POS)=X(I)
             POINT(2,POS)=Y(I)
             NBP(POS)=NBP(POS)+1
             !              define the end of element
             ELEM_NEW(I,ELEM)=POS
          ENDDO
       ELSE
          ELEM_NEW(1:2,ELEM)=0
       ENDIF
    ENDDO
    !*    adjust points on the axes
    SELECT CASE(TYPGEO)
       CASE(5:6)
       CALL SAL128_3(NN,POINT,EPS,EPS2)
       CASE(7)
       CALL SAL128_4(NN,POINT,EPS,EPS2)
       CASE(8,12) 
       CALL SAL128_5(NN,POINT,EPS,EPS2)
       CASE(9) 
       CALL SAL128_6(NN,POINT,EPS,EPS2)
       CASE(10) 
       CALL SAL128_7(NN,POINT,EPS,EPS2)
       CASE(11) 
       CALL SAL128_8(NN,POINT,EPS,EPS2)
    END SELECT

    !*    redefine the elements in rpar
    IF(KNDEX/=0)WRITE(FOUT,'(/," redefine elements : ",/)')
    DO ELEM=1,NB_ELEM
       TYPE=IPAR(1,ELEM)
       IF(TYPE/=G_ELE_TYPE(2)) THEN
          I1=ELEM_NEW(1,ELEM)
          I2=ELEM_NEW(2,ELEM)
          CALL SAL129(POINT(1,I1),POINT(2,I1),POINT(1,I2),POINT(2,I2), &
               TYPE,RPAR(:,ELEM),ELEM,KNDEX,FOUT)
       ENDIF
    ENDDO
    !
    DEALLOCATE(POINT,ELEM_NEW)
    !
  END SUBROUTINE SAL128
  !
  SUBROUTINE SAL128_3(NN,POINT,EPS,EPS2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! process retangular geometry with translation or symmetry boundary condition
    ! (TYPGEO=5,6): adjusts points on the axes, computes the lengths of the
    ! retangle sides; in case of with translation,element lengths on the opposite
    ! axes will be adjusted to be the same
    !
    !Parameters: input
    ! NN        total number of points
    ! EPS       when distance of points to axes less than EPS, displace
    !           the points onto the axes
    ! EPS2      variable set to EPS*EPS
    !
    !Parameters: input/output
    ! POINT     coordinates of points
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS,    ONLY : PDB
    USE SAL_GEOMETRY_TYPES,     ONLY : TYPGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)     :: NN
    REAL,      INTENT(IN)     :: EPS,EPS2
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:)  :: POINT
    !     DIMENSION        POINT(2,NN)
    !****
    INTEGER   :: NP,NAXIS,I,J,K,M,P,IPOINT(4,NN),BP(4),IP(4)
    LOGICAL   :: LGADJUST
    REAL(PDB) :: D,AUX,EX(4),EY(4),DIS(4),AX(4),AY(4),D_AXIS(4,NN)
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     flag to adjust the element lengths on the opposite axes
    LGADJUST=TYPGEO==5
    !     compute sides of the rectangle
    LX=0.; LY=0.
    DO P=1,NN
       IF(ABS(POINT(1,P))<EPS) THEN
          IF(LY<POINT(2,P)) LY=POINT(2,P)
       ENDIF
       IF(ABS(POINT(2,P))<EPS) THEN
          IF(LX<POINT(1,P)) LX=POINT(1,P)
       ENDIF
    ENDDO
    !*    axis directions and their distance to (0,0)
    NAXIS=4
    EX(1)=1.; EY(1)=0.; DIS(1)=0.
    EX(2)=1.; EY(2)=0.; DIS(2)=-LY
    EX(3)=0.; EY(3)=1.; DIS(3)=0.
    EX(4)=0.; EY(4)=1.; DIS(4)=LX
    NP=4
    AX(1)=0.; AY(1)=0.
    AX(2)=0.; AY(2)=LY
    AX(3)=LX; AY(3)=0.
    AX(4)=LX; AY(4)=LY
    !     beginning point of the axes
    IF(LGADJUST) THEN
       !     axis numbering:
       !                M=2
       !           ************
       !           *          *
       !      M=3  *          *  M=4
       !           ************
       !                M=1
       !
       BP(1)=1; BP(2)=2; BP(3)=1; BP(4)=3
       IP=0
    ENDIF
    ITER0:DO P=1,NN
       DO I=1,NP
          !           if it is the 4 corners of the rectangle
          D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2
          IF(D.LT.EPS2) THEN
             !              move point to the corners
             POINT(1,P)=AX(I) ; POINT(2,P)=AY(I)
             IF(LGADJUST) THEN
                !                 put distance to the axis origin
                SELECT CASE(I)
                   CASE(2)
                   !                    end of axis 3,keep its distance and number
                   IP(3)=IP(3)+1
                   D_AXIS(3,IP(3))=AY(I)
                   IPOINT(3,IP(3))=P
                   CASE(3)
                   !                    end of axis 1,keep its distance and number
                   IP(1)=IP(1)+1
                   D_AXIS(1,IP(1))=AX(I)
                   IPOINT(1,IP(1))=P
                   CASE(4)
                   !                    end of axes 2&4,keep its distance and number
                   IP(2)=IP(2)+1
                   D_AXIS(2,IP(2))=AX(I)
                   IPOINT(2,IP(2))=P
                   IP(4)=IP(4)+1
                   D_AXIS(4,IP(4))=AY(I)
                   IPOINT(4,IP(4))=P
                END SELECT
             ENDIF
             CYCLE ITER0
          ENDIF
       ENDDO
       DO I=1,NAXIS
          !           distance to the axis
          D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I)
          !           when aux>0,the point is near the half axis
          AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I)
          IF(ABS(D).LT.EPS.AND.AUX>=0) THEN
             !              move point to the axis
             POINT(1,P)=POINT(1,P)-D*EY(I)
             POINT(2,P)=POINT(2,P)+D*EX(I)
             IF(LGADJUST) THEN
                !                 compute distance to the axis origin
                D=SQRT((POINT(1,P)-AX(BP(I)))**2+(POINT(2,P)-AY(BP(I)))**2)
                IP(I)=IP(I)+1
                D_AXIS(I,IP(I))=D
                IPOINT(I,IP(I))=P
             ENDIF
             CYCLE ITER0
          ENDIF
       ENDDO
    ENDDO ITER0
    !
    IF(LGADJUST) THEN
       DO I=1,NAXIS,2
          IF(IP(I)/=IP(I+1)) &
               CALL XABORT('SAL128_3: axial points nber not the same,axis')
          !           sort the 'd_axis' table
          DO J=I,I+1
             DO P=1,IP(J)
                DO K=P+1,IP(J)
                   IF(D_AXIS(J,P)>D_AXIS(J,K)) THEN
                      D=D_AXIS(J,P)
                      D_AXIS(J,P)=D_AXIS(J,K)
                      D_AXIS(J,K)=D
                      M=IPOINT(J,P)
                      IPOINT(J,P)=IPOINT(J,K)
                      IPOINT(J,K)=M
                   ENDIF
                ENDDO
             ENDDO
          ENDDO
          DO P=1,IP(I)
             IF(ABS(D_AXIS(I,P)-D_AXIS(I+1,P))>EPS) THEN
                IF(KNDEX/=0) &
                     WRITE(FOUT,'(" warning: too great axial length difference",&
                     & 2(/,2X,"axis = ",I3," point = ",I3," d = ",E13.6))')&
                     I,P,D_AXIS(I,P),I+1,P,D_AXIS(I+1,P)
             ENDIF
             D=(D_AXIS(I,P)+D_AXIS(I+1,P))*0.5
             DO J=I,I+1
                K=IPOINT(J,P)
                POINT(1,K)=D*EX(J)+AX(BP(J))
                POINT(2,K)=D*EY(J)+AY(BP(J))
             ENDDO
          ENDDO
       ENDDO
    ENDIF
    !
  END SUBROUTINE SAL128_3
  !
  SUBROUTINE SAL128_4(NN,POINT,EPS,EPS2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! process a 1/8 square geometry (typgeo=7): adjusts points on the axes,
    ! computes the square side length
    !
    !Parameters: input
    ! NN        total number of points
    ! EPS       when distance of points to axes less than EPS, displace
    !           the points onto the axes
    ! EPS2      variable set to EPS*EPS
    !
    !Parameters: input/output
    ! POINT     coordinates of points
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS,    ONLY : PDB
    USE SAL_GEOMETRY_TYPES,     ONLY : ANGGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)     :: NN
    REAL,      INTENT(IN)     :: EPS,EPS2
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT
    !****
    INTEGER   :: NP,NAXIS,I,P
    REAL(PDB) :: D,AUX,EX(3),EY(3),DIS(3),AX(3),AY(3)
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     compute the square sides
    LX=0.; LY=0.
    DO P=1,NN
       IF(LY<POINT(2,P)) LY=POINT(2,P)
       IF(LX<POINT(1,P)) LX=POINT(1,P)
    ENDDO
    IF(KNDEX>0.AND.ABS(LY-LX)>EPS) &
         WRITE(FOUT,'(" warning: the square sides are not the same!", &
         &  " lx = ",E12.4," ly = ",E12.4)')LX,LY
    LX=(LX+LY)*0.5
    LY=LX
    !*    adjust points on the axes
    NAXIS=3
    !     axis directions
    EX(1)=1.; EY(1)=0.; DIS(1)=0.
    EX(2)=COS(ANGGEO); EY(2)=SIN(ANGGEO); DIS(2)=0.
    EX(3)=0.; EY(3)=1.; DIS(3)=LX
    NP=3
    !     axis ends
    AX(1)=0.; AY(1)=0.
    AX(2)=LX; AY(2)=0.
    AX(3)=LX; AY(3)=LY
    ITER0:DO P=1,NN
       DO I=1,NP
          D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2
          IF(D.LT.EPS2) THEN
             !              move point to the axis origin
             POINT(1,P)=AX(I) ; POINT(2,P)=AY(I)
             CYCLE ITER0
          ENDIF
       ENDDO
       DO I=1,NAXIS
          !           distance to the axis
          D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I)
          !           when aux>0,the point is near the half axis
          AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I)
          IF(ABS(D).LT.EPS.AND.AUX>=0) THEN
             !        move point to the axis
             POINT(1,P)=POINT(1,P)-D*EY(I)
             POINT(2,P)=POINT(2,P)+D*EX(I)
             CYCLE ITER0
          ENDIF
       ENDDO
    ENDDO ITER0
    !
  END SUBROUTINE SAL128_4
  !
  SUBROUTINE SAL128_5(NN,POINT,EPS,EPS2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! process a trianglar geometry (TYPGEO=8 or 12) with specular reflection:
    ! adjusts points on the axes, computes the triangular side length
    !
    !Parameters: input
    ! NN        total number of points
    ! EPS       when distance of points to axes less than EPS, displace
    !           the points onto the axes
    ! EPS2      variable set to EPS*EPS
    !
    !Parameters: input/output
    ! POINT     coordinates of points
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS,    ONLY : PDB
    USE SAL_GEOMETRY_TYPES,     ONLY : ANGGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)     :: NN
    REAL,      INTENT(IN)     :: EPS,EPS2
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:)  :: POINT
    !****
    INTEGER   :: NP,NAXIS,I,P
    REAL(PDB) :: D,AUX,EX(3),EY(3),DIS(3),AX(3),AY(3)
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     compute the triangular side
    LX=0.; LY=0.
    DO P=1,NN
       IF(ABS(POINT(2,P))<EPS) THEN
          IF(LX<POINT(1,P)) LX=POINT(1,P)
       ENDIF
       IF(LY<POINT(2,P)) LY=POINT(2,P)
    ENDDO
    IF(KNDEX>0.AND.ABS(LY-LX*SIN(ANGGEO))>EPS) &
         WRITE(FOUT,'(" warning: h is different from a*sin(anggeo)!", &
         &  " a = ",E12.4," h = ",E12.4)')LX,LY
    !     adjust h according a
    LY=LX*SIN(ANGGEO)
    !*    adjust points on the axes
    NAXIS=3
    !     axis directions
    EX(1)=1.; EY(1)=0.; DIS(1)=0.
    EX(2)=COS(ANGGEO); EY(2)=SIN(ANGGEO); DIS(2)=0.
    EX(3)=COS(ANGGEO*2.); EY(3)=SIN(ANGGEO*2.); DIS(3)=LY
    NP=3
    !     axis ends
    AX(1)=0.; AY(1)=0.
    AX(2)=LX; AY(2)=0.
    AX(3)=LX*0.5; AY(3)=LY
    ITER0:DO P=1,NN
       DO I=1,NP
          D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2
          IF(D.LT.EPS2) THEN
             !              move point to the axis origin
             POINT(1,P)=AX(I) ; POINT(2,P)=AY(I)
             CYCLE ITER0
          ENDIF
       ENDDO
       DO I=1,NAXIS
          !           distance to the axis
          D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I)
          !           when aux>0,the point is near the half axis
          AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I)
          IF(ABS(D).LT.EPS.AND.AUX>=0) THEN
             !              move point to the axis
             POINT(1,P)=POINT(1,P)-D*EY(I)
             POINT(2,P)=POINT(2,P)+D*EX(I)
             CYCLE ITER0
          ENDIF
       ENDDO
    ENDDO ITER0
    !
  END SUBROUTINE SAL128_5
  !
  SUBROUTINE SAL128_6(NN,POINT,EPS,EPS2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! process a hexagonal geometry with translations on all sides (typgeo=9):
    ! adjusts points on the axes, computes the lengths of hexagon sides;
    ! element lengths on the opposite axes will be adjusted to be the same
    !
    !Parameters: input
    ! NN        total number of points
    ! EPS       when distance of points to axes less than EPS, displace
    !           the points onto the axes
    ! EPS2      variable set to EPS*EPS
    !
    !Parameters: input/output
    ! POINT     coordinates of points
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS,    ONLY : PDB
    USE SAL_GEOMETRY_TYPES,     ONLY : ANGGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)     :: NN
    REAL,      INTENT(IN)     :: EPS,EPS2
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:)  :: POINT
    !****
    INTEGER   :: NP,NAXIS,I,J,K,M,P,IPOINT(6,NN),BP(6),IP(6)
    REAL(PDB) :: D,AUX,EX(6),EY(6),DIS(6),AX(6),AY(6),D_AXIS(6,NN)
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     compute sides of the hexagon: LX=SIDE LENGTH,LY=LX*SIN(A)
    LX=0.; LY=0.
    DO P=1,NN
       IF(ABS(POINT(2,P))<EPS) THEN
          IF(LX<POINT(1,P)) LX=POINT(1,P)
       ENDIF
       IF(LY<POINT(2,P)) LY=POINT(2,P)
    ENDDO
    IF(KNDEX>0.AND.ABS(LY-LX*SIN(ANGGEO))>EPS) &
         WRITE(FOUT,'(" warning: h is different from A*SIN(ANGGEO)!", &
         & " a = ",E12.4," h = ",E12.4)')LX,LY
    !*    axis directions and their distance to (0,0)
    NAXIS=6
    EX(1)=1.; EY(1)=0.; DIS(1)=LY
    EX(2)=1.; EY(2)=0.; DIS(2)=-LY
    EX(3)=COS(ANGGEO); EY(3)=SIN(ANGGEO); DIS(3)=LY
    EX(4)=EX(3);       EY(4)=EY(3);       DIS(4)=-LY
    EX(5)=COS(2.*ANGGEO); EY(5)=SIN(2.*ANGGEO); DIS(5)=LY
    EX(6)=EX(5);          EY(6)=EY(5);          DIS(6)=-LY
    NP=6
    AX(1)=-LX*0.5; AY(1)=-LY
    AX(2)=LX*0.5;  AY(2)=-LY
    AX(3)=LX;      AY(3)=0.
    AX(4)=LX*0.5;  AY(4)=LY
    AX(5)=-LX*0.5; AY(5)=LY
    AX(6)=-LX;     AY(6)=0.
    !     axis numbering:
    !               M=2
    !              *****
    !         M=4 *     *   M=5
    !            *   *   *
    !         M=6 *     *   M=3
    !              *****
    !               M=1
    !
    BP(1)=1; BP(2)=5; BP(3)=2; BP(4)=6; BP(5)=3; BP(6)=1
    IP=0
    ITER0:DO P=1,NN
       DO I=1,NP
          !           if it is the 6 corners of the hexagon
          D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2
          IF(D.LT.EPS2) THEN
             !              move point to the corner
             POINT(1,P)=AX(I) ; POINT(2,P)=AY(I)
             !              put distance to the axis origin
             J=0
             SELECT CASE(I)
                CASE(1)
                CYCLE ITER0
                CASE(2)
                !                 END OF AXIS 1
                J=1
                CASE(3)
                !                 END OF AXIS 3
                J=3
                CASE(4)
                !                 END OF AXES 2&5
                J=2
                CASE(5)
                !                 END OF AXIS 4
                J=4
                CASE(6)
                !                 END OF AXIS 6
                J=6
             END SELECT
             !              keep its distance and number
10           IP(J)=IP(J)+1
             D_AXIS(J,IP(J))=LX
             IPOINT(J,IP(J))=P
             IF(J==2) THEN
                J=5
                GOTO 10
             ENDIF
             CYCLE ITER0
          ENDIF
       ENDDO
       DO I=1,NAXIS
          !           distance to the axis
          D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I)
          !           when aux>0,the point is near the half axis
          AUX=(POINT(1,P)-AX(BP(I)))*EX(I)+(POINT(2,P)-AY(BP(I)))*EY(I)
          IF(ABS(D).LT.EPS.AND.AUX>=0) THEN
             !              move point to the axis
             POINT(1,P)=POINT(1,P)-D*EY(I)
             POINT(2,P)=POINT(2,P)+D*EX(I)
             !              compute distance to the axis origin
             D=SQRT((POINT(1,P)-AX(BP(I)))**2+(POINT(2,P)-AY(BP(I)))**2)
             IP(I)=IP(I)+1
             D_AXIS(I,IP(I))=D
             IPOINT(I,IP(I))=P
             CYCLE ITER0
          ENDIF
       ENDDO
    ENDDO ITER0
    !
    DO I=1,NAXIS,2
       IF(IP(I)/=IP(I+1)) &
            CALL XABORT('SAL128_6: axial points nber not the same,axis')
       !        sort the 'd_axis' table
       DO J=I,I+1
          DO P=1,IP(J)
             DO K=P+1,IP(J)
                IF(D_AXIS(J,P)>D_AXIS(J,K)) THEN
                   D=D_AXIS(J,P)
                   D_AXIS(J,P)=D_AXIS(J,K)
                   D_AXIS(J,K)=D
                   M=IPOINT(J,P)
                   IPOINT(J,P)=IPOINT(J,K)
                   IPOINT(J,K)=M
                ENDIF
             ENDDO
          ENDDO
       ENDDO
       DO P=1,IP(I)
          IF(ABS(D_AXIS(I,P)-D_AXIS(I+1,P))>EPS) THEN
             IF(KNDEX/=0) &
                  WRITE(FOUT,'(" warning: too great axial length difference",&
                  & 2(/,2X,"axis = ",I3," point = ",I3," d = ",E13.6))')&
                  I,P,D_AXIS(I,P),I+1,P,D_AXIS(I+1,P)
          ENDIF
          D=(D_AXIS(I,P)+D_AXIS(I+1,P))*0.5
          DO J=I,I+1
             K=IPOINT(J,P)
             POINT(1,K)=D*EX(J)+AX(BP(J))
             POINT(2,K)=D*EY(J)+AY(BP(J))
          ENDDO
       ENDDO
    ENDDO
    !
  END SUBROUTINE SAL128_6
  !
  SUBROUTINE SAL128_7(NN,POINT,EPS,EPS2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! process an isocel trianglar geometry (TYPGEO=10) with RA60 rotation:
    ! adjusts points on the axes, computes the triangular side length
    !
    !Parameters: input
    ! NN        total number of points
    ! EPS       when distance of points to axes less than EPS, displace
    !           the points onto the axes
    ! EPS2      variable set to EPS*EPS
    !
    !Parameters: input/output
    ! POINT     coordinates of points
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS,    ONLY : PDB
    USE SAL_GEOMETRY_TYPES,     ONLY : ANGGEO,LX=>LENGTHX,LY=>LENGTHY
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)     :: NN
    REAL,      INTENT(IN)     :: EPS,EPS2
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:)  :: POINT
    !****
    INTEGER   :: NP,NAXIS,I,P
    REAL(PDB) :: D,AUX,EX(3),EY(3),DIS(3),AX(3),AY(3)
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     compute the triangular side
    LX=0.; LY=0.
    DO P=1,NN
       IF(ABS(POINT(2,P))<EPS) THEN
          IF(LX<POINT(1,P)) LX=POINT(1,P)
       ENDIF
       IF(LY<POINT(2,P)) LY=POINT(2,P)
    ENDDO
    IF(KNDEX>0.AND.ABS(LY-LX*SIN(ANGGEO))>EPS) &
         WRITE(FOUT,'(" warning: h is different from a*sin(anggeo)!", &
         &  " a = ",E12.4," h = ",E12.4)')LX,LY
    !     adjust h according a
    LY=LX*SIN(ANGGEO)
    !*    adjust points on the axes
    NAXIS=3
    !     axis directions
    EX(1)=1.; EY(1)=0.; DIS(1)=0.
    EX(2)=COS(ANGGEO); EY(2)=SIN(ANGGEO); DIS(2)=0.
    EX(3)=COS(ANGGEO*2.); EY(3)=SIN(ANGGEO*2.); DIS(3)=LY
    NP=3
    !     axis ends
    AX(1)=0.; AY(1)=0.
    AX(2)=LX; AY(2)=0.
    AX(3)=LX*0.5; AY(3)=LY
    ITER0:DO P=1,NN
       DO I=1,NP
          D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2
          IF(D.LT.EPS2) THEN
             !              move point to the axis origin
             POINT(1,P)=AX(I) ; POINT(2,P)=AY(I)
             CYCLE ITER0
          ENDIF
       ENDDO
       DO I=1,NAXIS
          !           distance to the axis
          D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I)
          !           when aux>0,the point is near the half axis
          AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I)
          IF(ABS(D).LT.EPS.AND.AUX>=0) THEN
             !              move point to the axis
             POINT(1,P)=POINT(1,P)-D*EY(I)
             POINT(2,P)=POINT(2,P)+D*EX(I)
             CYCLE ITER0
          ENDIF
       ENDDO
    ENDDO ITER0
    !
  END SUBROUTINE SAL128_7
  !
  SUBROUTINE SAL128_8(NN,POINT,EPS,EPS2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! process R120 lozenge geometry (TYPGEO=11): adjusts points on the axes,
    ! computes the lozenge side length
    !
    !Parameters: input
    ! NN        total number of points
    ! EPS       when distance of points to axes less than EPS, displace
    !           the points onto the axes
    ! EPS2      variable set to EPS*EPS
    !
    !Parameters: input/output
    ! POINT     coordinates of points
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS,    ONLY : PDB
    USE SAL_GEOMETRY_TYPES,     ONLY : ANGGEO,LX=>LENGTHX,LY=>LENGTHY
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)     :: NN
    REAL,      INTENT(IN)     :: EPS,EPS2
    REAL(PDB), INTENT(INOUT), DIMENSION(:,:)  :: POINT
    !****
    INTEGER   :: NP,NAXIS,I,P
    REAL(PDB) :: D,AUX,EX(4),EY(4),DIS(4),AX(4),AY(4)
    INTEGER, PARAMETER :: FOUT =6
    !****
    !     compute the triangular side
    LX=0.; LY=0.
    DO P=1,NN
       IF(ABS(POINT(2,P))<EPS) THEN
          IF(LX<POINT(1,P)) LX=POINT(1,P)
       ENDIF
       IF(LY<POINT(2,P)) LY=POINT(2,P)
    ENDDO
    IF(ABS(LY-LX*SIN(ANGGEO))>EPS) &
         WRITE(FOUT,'(" warning: h is different from a*sin(anggeo)!", &
         &  " a = ",E12.4," h = ",E12.4)') LX,LY
    !     adjust h according to base
    LY=LX*SIN(ANGGEO)
    !*    adjust points on the axes
    NAXIS=4
    !     axis directions
    EX(1)=1.; EY(1)=0.; DIS(1)=0.
    EX(2)=1.; EY(2)=0.; DIS(2)=-LY
    EX(3)=COS(ANGGEO); EY(3)=SIN(ANGGEO); DIS(3)=0.
    EX(4)=COS(ANGGEO); EY(4)=SIN(ANGGEO); DIS(4)=LY
    NP=4
    !     axis ends
    AX(1)=0.    ; AY(1)=0.
    AX(2)=LX*0.5; AY(2)=LY
    AX(3)=LX    ; AY(3)=0.
    AX(4)=LX*1.5; AY(4)=LY
    ITER0:DO P=1,NN
       DO I=1,NP
          D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2
          IF(D.LT.EPS2) THEN
             !              move point to the axis origin
             POINT(1,P)=AX(I) ; POINT(2,P)=AY(I)
             CYCLE ITER0
          ENDIF
       ENDDO
       DO I=1,NAXIS
          !           distance to the axis
          D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I)
          !           when aux>0,the point is near the half axis
          AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I)
          IF(ABS(D).LT.EPS.AND.AUX>=0) THEN
             !        move point to the axis
             POINT(1,P)=POINT(1,P)-D*EY(I)
             POINT(2,P)=POINT(2,P)+D*EX(I)
             CYCLE ITER0
          ENDIF
       ENDDO
    ENDDO ITER0
    !
  END SUBROUTINE SAL128_8
  !
  SUBROUTINE SAL129(X1,Y1,X2,Y2,TYPE,RPAR,II,KNDEX,FOUT)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! recompute element data by giving two ends of the element
    !
    !Parameters: input
    ! X1      new position in X for the end 1
    ! Y1      new position in Y for the end 1
    ! X2      new position in X for the end 2
    ! Y2      new position in Y for the end 2
    ! TYPE    type of element  1 (segment)  3 (arc of circle)
    ! II      element to be changed
    ! KNDEX   if not 0 print out data when motion>EPS
    ! FOUT    printout file
    !
    !Parameters: input/output
    ! RPAR    descriptors of the element
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES,     ONLY : NRPAR,EPS,G_ELE_TYPE,G_ELE_LEN
    USE SAL_NUMERIC_MOD, ONLY : SALACO
    !****
    IMPLICIT NONE
    ! in variable
    ! ***********
    INTEGER, INTENT(IN)                    :: TYPE,KNDEX,FOUT,II
    REAL(PDB), INTENT(IN)                  :: X1,Y1,X2,Y2
    REAL(PDB), INTENT(INOUT), DIMENSION(:) :: RPAR
    !****

    ! local variable
    ! **************
    REAL(PDB) :: CX,CY,DELX,DELY,R,THETA,SX,SY,VX,VY,V2,DIST,MOV1,MOV2
    INTEGER   :: I
    REAL(PDB) :: RPAR_OLD(NRPAR)
    !****
    !     keep old parameters:
    RPAR_OLD(1:NRPAR)=RPAR(1:NRPAR)
    !
    MOV1=0._PDB; MOV2=0._PDB;
    IF(TYPE == G_ELE_TYPE(1)) THEN
      !*       segment :
      !        compute the motions of two ends
      DELX=X1-RPAR(1)
      DELY=Y1-RPAR(2)
      MOV1=SQRT(DELX*DELX+DELY*DELY)
      DELX=X2-(RPAR(1)+RPAR(3))
      DELY=Y2-(RPAR(2)+RPAR(4))
      MOV2=SQRT(DELX*DELX+DELY*DELY)
      !        set new ends, compute direction
      RPAR(1)=X1
      RPAR(2)=Y1
      RPAR(3)=X2-X1
      RPAR(4)=Y2-Y1
      RPAR(5)=SQRT(RPAR(3)*RPAR(3)+RPAR(4)*RPAR(4))
    ELSE IF(TYPE == G_ELE_TYPE(3)) THEN
      !*       arc of circle :
      !        get old radius
      CX=RPAR(1); CY=RPAR(2)
      !        compute the motion of two ends
      R=RPAR(3)
      THETA=RPAR(4)
      DELX=CX+R*COS(THETA)-X1
      DELY=CY+R*SIN(THETA)-Y1
      MOV1=SQRT(DELX*DELX+DELY*DELY)
      THETA=RPAR(5)
      DELX=CX+R*COS(THETA)-X2
      DELY=CY+R*SIN(THETA)-Y2
      MOV2=SQRT(DELX*DELX+DELY*DELY)
      !        compute new radius, new center and new angles
      SX=X1-CX
      SY=Y1-CY
      !        get demi-vector from (x1,y1) to (x2,y2)
      VX=(X2-X1)/2.0
      VY=(Y2-Y1)/2.0
      V2=VX*VX+VY*VY
      !        get center closest to old center, change radius
      !        compute new r
      R=SQRT(V2+(SX*VY-SY*VX)**2/V2)
      RPAR(3)=R
      !        compute new center
      DIST=1.0_PDB+(VX*SX+VY*SY)/V2
      !        get new center and compute angles
      CX=CX+VX*DIST
      CY=CY+VY*DIST
      RPAR(1)=CX
      RPAR(2)=CY
      THETA=SALACO((X1-CX)/R,Y1-CY)
      RPAR(4)=THETA
      THETA=SALACO((X2-CX)/R,Y2-CY)
      RPAR(5)=THETA
      !        if phi1 > phi2 put phi2 equal to phi2+twopi
      IF(RPAR(4)>RPAR(5)) RPAR(5)=RPAR(5)+TWOPI
    ENDIF
    IF(KNDEX/=0) THEN
       IF(MOV1>EPS.OR.MOV2>EPS) THEN
          WRITE(FOUT,'(/," old element ",I4," : ",6(1P,E13.4))') &
               II,(RPAR_OLD(I),I=1,G_ELE_LEN(TYPE))
          WRITE(FOUT,'(" new element ",I4," : ",6(1P,E13.4))') &
               II,(RPAR(I),I=1,G_ELE_LEN(TYPE))
       ENDIF
    ENDIF
    !
  END SUBROUTINE SAL129
  !
  SUBROUTINE SAL130(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! perform domain definition and set boundary conditions
    !   - computes node perimeters per 2d macro
    !   - reads boundary conditions
    !   - computes external perimeters 
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES,      ONLY : G_BC_TYPE,TYPGEO
    !****
    IMPLICIT NONE
    ! in variable
    ! ************
    TYPE(T_G_BASIC), INTENT(INOUT):: GG
    ! local variable
    !****
    INTEGER :: NN,OK
    INTEGER, DIMENSION(GG%NB_ELEM*2) :: AUX_ARR
    !****
    !*    compute node perimeters for the macro
    CALL SAL130_2(GG%NB_ELEM,GG%NB_NODE,GG%IPAR,GG%PPERIM_NODE, &
         GG%PERIM_NODE,AUX_ARR)
    !
    !*    - compute number of bc's per 2D macro,
    !       NB_BC2 counts total nber of 2D bc's
    !     - compute IBC2_ELEM, keep relative 2D bc nber to elements
    !     - allocation : 2D bc structures
    !                    2D perimeter structure for a macro
    !     - get list of elements in 2d macro perimeter

    CALL SAL130_4(GG%NB_ELEM,NN,GG%IPAR,GG%IBC2_ELEM,AUX_ARR)
    
    GG%NB_BC2=NN

    ALLOCATE (GG%TYPE_BC2(NN),GG%IDATA_BC2(NN),GG%PERIM_MAC2(NN), &
         GG%PPERIM_MAC2(7),STAT=OK)
    IF(OK/=0) CALL XABORT('SAL130: not enough memory I,R')
    GG%PPERIM_MAC2(:7)=0
    GG%PERIM_MAC2(1:NN)=AUX_ARR(1:NN)
    GG%NPERIM_MAC2=NN
    !
    !*    read boundary conditions:
    CALL SAL131(GG)
    !
    !     - define IELEM_SURF2
    IF(GG%NB_SURF2>0) THEN
       ! allocate
       ALLOCATE(GG%SURF2(GG%NB_SURF2),STAT = OK)
       IF(OK /= 0) CALL XABORT('SAL130: not enough memory I,R')
       CALL SAL130_6(GG%NB_SURF2,GG%IBC2_SURF2,GG%PERIM_MAC2,GG%IELEM_SURF2)
    ENDIF
    !
    !*    construct perimeter structures for rotative or symmetrical geometry
    SELECT CASE(TYPGEO)
       CASE(1:2)
       CALL SAL130_8(GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%PPERIM_MAC2,GG%DIST_AXIS, &
            GG%IBC2_ELEM,GG%TYPE_BC2,GG%IDATA_BC2,GG%BCDATA,GG%IPAR,GG%RPAR)
       CASE(5:12)
       CALL SAL130_10(GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%PPERIM_MAC2,GG%DIST_AXIS, &
            GG%IBC2_ELEM,GG%TYPE_BC2,GG%IDATA_BC2,GG%BCDATA,GG%IPAR,GG%RPAR)
    END SELECT
    !
  END SUBROUTINE SAL130
  !
  SUBROUTINE SAL130_2(NB_ELEM,NB_NODE,IPAR,PPERIM,PERIM,LIST)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! compute node perimeters for one 2D macro
    !
    !Parameters: input
    ! NB_ELEM   number of elements
    ! NB_NODE   number of nodes
    ! IPAR      integer descriptors for elements
    !
    !Parameters: output
    ! PPERIM    array pointer to list of elements in perimeter per node
    ! PERIM     elements in perimeters per node
    ! LIST      temporary array
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN)                    :: NB_ELEM,NB_NODE
    INTEGER, INTENT(IN),   DIMENSION(:,:)  :: IPAR
    INTEGER, INTENT(OUT),  DIMENSION(:)    :: LIST
    INTEGER, INTENT(OUT),  DIMENSION(:)    :: PPERIM
    INTEGER, POINTER,      DIMENSION(:)    :: PERIM
    !****
    INTEGER :: NT,NN,NODE,ELEM,OK
    CHARACTER(LEN=131) :: HSMG
    !****
    NT=0
    PPERIM(1)=1
    DO NODE=1,NB_NODE
       DO ELEM=1, NB_ELEM
          IF(IPAR(2,ELEM)==NODE.OR.IPAR(3,ELEM)==NODE) THEN
             NT=NT+1
             LIST(NT)=ELEM
          ENDIF
       ENDDO
       NN=NT+1-PPERIM(NODE)
       IF(NN>0) THEN
          PPERIM(NODE+1)=NT+1
       ELSE
          WRITE(HSMG,'(15HSAL130_2: node=,i5,19H without perimeter.)') NODE
          CALL XABORT(HSMG)
       ENDIF
    ENDDO
    IF(NT>0) THEN
       ALLOCATE (PERIM(NT), STAT=OK)
       IF(OK/=0) CALL XABORT('SAL130_2: not enough memory I')
       PERIM(1:NT)=LIST(1:NT)
    ENDIF
    !
  END SUBROUTINE SAL130_2
  !
  SUBROUTINE SAL130_4(NB_ELEM,NB_BC,IPAR,IBC2_ELEM,LIST_BC)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! count number of bc's in a 2D macro
    !
    !Parameters: input
    ! NB_ELEM   number of elements
    !
    !Parameters: input/output
    ! IPAR      integer descriptors for elements
    !
    !Parameters: output
    ! NB_BC     nber of bc's
    ! IBC2_ELEM relative 2D bc index per element
    ! LIST_BC   list of elements in boundary
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN)                 :: NB_ELEM
    INTEGER, INTENT(OUT)                :: NB_BC
    INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR
    INTEGER, INTENT(OUT),  DIMENSION(:) :: IBC2_ELEM,LIST_BC
    !****
    INTEGER :: ELEM
    !****
    !     initiation
    IBC2_ELEM(:)=0
    NB_BC=0
    DO ELEM=1, NB_ELEM
       IF(IPAR(2,ELEM)<=0.AND.IPAR(3,ELEM)<=0) THEN
          CALL XABORT('SAL130_4: two boundaries for element')
       ELSE IF(IPAR(2,ELEM)<=0.OR.IPAR(3,ELEM)<=0) THEN
          NB_BC=NB_BC+1
          LIST_BC(NB_BC)=ELEM
          IF(IBC2_ELEM(ELEM)/=0) CALL XABORT('SAL130_4: two surfaces to element')
          IBC2_ELEM(ELEM)=NB_BC
          !           set mark for the macro connection surface :
       ENDIF
    ENDDO
    !
  END SUBROUTINE SAL130_4
  !
  SUBROUTINE SAL130_6(NSURF,IBC2_SURF2,IELEM_BC2,IELEM_SURF2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! get element order indices per surface
    !
    !Parameters: input
    ! NSURF       number of surfaces
    ! IBC2_SURF2  2D bc order index per surface
    ! IELEM_BC2   element order index per bc
    !
    !Parameters: output
    ! IELEM_SURF2 element order index per surface
    !
    !---------------------------------------------------------------------
    !
    !****
    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE
    !****
    IMPLICIT NONE
    INTEGER, INTENT(IN)                 :: NSURF
    INTEGER, INTENT(IN),   DIMENSION(:) :: IBC2_SURF2,IELEM_BC2
    INTEGER, INTENT(OUT),  DIMENSION(:) :: IELEM_SURF2
    !****
    INTEGER :: SURF,IBC,ELEM
    !****
    IF(NSURF > 0) THEN
       DO SURF=1, NSURF
          !           get relative bc order number
          IBC=IBC2_SURF2(SURF)
          !           get relative element
          ELEM=IELEM_BC2(IBC)
          !           define ielem_surf2
          IELEM_SURF2(SURF)=ELEM
       ENDDO
    ENDIF
    !
  END SUBROUTINE SAL130_6
  !
  SUBROUTINE SAL130_8(NPERIM,PERIM,PPERIM,DIST_AXIS,IBC2_ELEM,TYPE_BC2,IDATA_BC2, &
       BCDATA,IPAR,RPAR)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! calculate PERIM_MAC2,DIST_AXIS for rotative or symmetrical geometry
    ! (TYPGEO=1 & 2)
    !
    !Parameters: input
    ! NPERIM     number of elements in perimeter
    ! IBC2_ELEM  2D bc order number per element
    ! TYPE_BC2   type of boundary conditions per 2D bc
    ! IDATA_BC2  position in the 'bcdata' table for each 2D bc
    ! BCDATA     table of boundary conditions data
    ! IPAR       integer element descriptor table
    ! RPAR       floating point element descriptor table
    !
    !Parameters: input/output
    ! PERIM      list of elements in perimeter,in return it will be in following
    !            order: (elems on axis 1)+(elems on axis 2)+(other elems)
    !
    !Parameters: output
    ! PPERIM     pointers to the table 'perim':
    !            (1): first elem on axis 1 (2): first elem on axis 2
    !            (3): first of other elems (4): NPERIM + 1
    ! DIST_AXIS  distance of points on axis to the center (0,0),in order of:
    !            (distances of points on axis 1)+(distances of points on axis 2)
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES,     ONLY : ANGGEO,EPS,G_BC_TYPE
    USE SAL_NUMERIC_MOD, ONLY : SAL141
    !****
    IMPLICIT NONE
    INTEGER, INTENT(IN)                  :: NPERIM
    INTEGER, INTENT(IN),    DIMENSION(:) :: IBC2_ELEM,TYPE_BC2,IDATA_BC2
    INTEGER, INTENT(INOUT), DIMENSION(:) :: PERIM
    INTEGER, INTENT(OUT),   DIMENSION(:) :: PPERIM
    REAL(PDB), POINTER,     DIMENSION(:) :: DIST_AXIS
    INTEGER, INTENT(IN),  DIMENSION(:,:) :: IPAR
    REAL(PDB), INTENT(INOUT),DIMENSION(:,:) :: RPAR
    REAL(PDB), INTENT(IN),DIMENSION(:,:) :: BCDATA
    !****
    !     LIST_ELEMS = table of elements,elements on axis are in increasing order
    !                  of distance to (0,0):
    !                  1=elements on axis 1; 2=elements on axis 2; 3=other elements
    !     AUX_DIST   = max distance of element ends to the beginnings of the axes:
    !                  1=distances on axis 1; 2=distances on axis 2
    INTEGER,   DIMENSION(NPERIM,3) :: LIST_ELEMS
    REAL(PDB), DIMENSION(NPERIM,2) :: AUX_DIST
    INTEGER   :: I,J,K,M,ELEM,TYPBC,IBC,IDATA,NBE(3),OK,NAXES
    REAL(PDB) :: ANGLE,X,Y,D
    !****
    NAXES=2
    !*    calculate number of elements on axis 1 & 2
    !     and their distances to the beginning of the axes
    NBE=0
    DO I=1,NPERIM
       ELEM=PERIM(I)
       IBC=IBC2_ELEM(ELEM)
       TYPBC=TYPE_BC2(IBC)
       IDATA=IDATA_BC2(IBC)
       !        default is the elements not on axis
       M=NAXES+1
       IF(TYPBC==G_BC_TYPE(3)) THEN
          !*          rotation:
          ANGLE=BCDATA(5,IDATA)
          IF(ABS(ANGLE-ANGGEO)<EPS) THEN
             !              element is on the axis 1 (angle=anggeo)
             M=1
          ELSEIF(ABS(ANGLE+ANGGEO)<EPS) THEN
             !              element is on the axis 2 (angle=-anggeo)
             M=2
          ENDIF
       ELSEIF(TYPBC==G_BC_TYPE(4)) THEN
          !*          symmetry:
          ANGLE=BCDATA(5,IDATA)
          IF(ABS(ANGLE)<EPS) THEN
             !              element is on the axis 1 (angle=0)
             M=1
          ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN
             !              element is on the axis 2 (angle=anggeo)
             M=2
          ENDIF
       ENDIF
       NBE(M)=NBE(M)+1
       LIST_ELEMS(NBE(M),M)=ELEM
       IF(M/=NAXES+1) THEN
          !           sort the element list according their distance to (0,0)
          D=0.
          DO K=1,2
             CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
             D=MAX(D,SQRT(X*X+Y*Y))
          ENDDO
          AUX_DIST(NBE(M),M)=D
          IF(NBE(M)>1) THEN
             DO K=1,NBE(M)-1
                IF(AUX_DIST(NBE(M),M)<AUX_DIST(K,M)) THEN
                   !                    insert the last element to position k
                   D=AUX_DIST(NBE(M),M)
                   ELEM=LIST_ELEMS(NBE(M),M)
                   DO J=NBE(M),K+1,-1
                      AUX_DIST(J,M)=AUX_DIST(J-1,M)
                      LIST_ELEMS(J,M)=LIST_ELEMS(J-1,M)
                   ENDDO
                   AUX_DIST(K,M)=D
                   LIST_ELEMS(K,M)=ELEM
                   EXIT
                ENDIF
             ENDDO
          ENDIF
       ENDIF
    ENDDO
    !     nber of blocks in 'perim'
    M=NAXES+1
    !*    organize elements as:
    !       (elems on axis 1)+(elems on axis 2)+(other elems)
    PPERIM(:)=0
    PPERIM(1)=1
    DO I=1,M
       PPERIM(I+1)=PPERIM(I)+NBE(I)
    ENDDO
    DO I=1,M
       PERIM(PPERIM(I):PPERIM(I+1)-1)=LIST_ELEMS(1:NBE(I),I)
    ENDDO
    !*    allocate and set dist_axis
    IF(PPERIM(NAXES+1)-1>0) THEN
       ALLOCATE(DIST_AXIS(PPERIM(NAXES+1)-1),STAT=OK)
       IF(OK.NE.0) CALL XABORT('SAL130_8: not enough memory R')
       DO I=1,NAXES
          DIST_AXIS(PPERIM(I):PPERIM(I+1)-1)=AUX_DIST(1:NBE(I),I)
       ENDDO
    ELSE
       NULLIFY(DIST_AXIS)
    ENDIF
    !
  end subroutine sal130_8
  !
  SUBROUTINE SAL130_10(NPERIM,PERIM,PPERIM,DIST_AXIS,IBC2_ELEM,TYPE_BC2,IDATA_BC2, &
       BCDATA,IPAR,RPAR)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! calculate PERIM_MAC2 and DIST_AXIS for the cyclical geometries:
    ! type 5&6: retangle with translations or symmetries on all sides
    ! type 7  : 1/8 assembly with symmetries on all sides
    ! type 8  : equilateral triangle with symmetries on all sides
    ! type 9  : hexagon with translations on all sides
    ! type 10 : hexagon with RA60 rotation
    ! type 11 : hexagon with R120 rotation
    ! type 12 : rectangular S30 triangle geometry
    !
    !Parameters: input
    ! NPERIM     number of elements in perimeter
    ! IBC2_ELEM  2D bc order number per element
    ! TYPE_BC2   type of boundary conditions per 2D bc
    ! IDATA_BC2  position in the 'bcdata' table for each 2D bc
    ! BCDATA     table of boundary conditions data
    ! IPAR       integer element descriptor table
    ! RPAR       floating point element descriptor table
    !
    !Parameters: input/output
    ! PERIM      list of elements in perimeter,in return it will be in following
    !            order: (elems on axis 1)+(elems on axis 2)+(other elems)
    !
    !Parameters: output
    ! PPERIM     pointers to the table 'perim':
    ! DIST_AXIS  distance of points on axis to the axis origin
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES,     ONLY : TYPGEO,ANGGEO,EPS,LX=>LENGTHX,LY=>LENGTHY,G_BC_TYPE
    USE SAL_NUMERIC_MOD, ONLY : SAL141
    !****
    IMPLICIT NONE
    INTEGER, INTENT(IN)                  :: NPERIM
    INTEGER, INTENT(IN),    DIMENSION(:) :: IBC2_ELEM,TYPE_BC2,IDATA_BC2
    INTEGER, INTENT(INOUT), DIMENSION(:) :: PERIM
    INTEGER, INTENT(OUT),   DIMENSION(:) :: PPERIM
    REAL(PDB), POINTER,     DIMENSION(:) :: DIST_AXIS
    INTEGER, INTENT(IN),  DIMENSION(:,:) :: IPAR
    REAL(PDB), INTENT(INOUT),DIMENSION(:,:) :: RPAR
    REAL(PDB), INTENT(IN),DIMENSION(:,:) :: BCDATA
    !****
    !     LIST_ELEMS = table of elements,elements on axis are in increasing order
    !                  of distance to axis origin:
    !                  1=elements on axis 1; 2=elements on axis 2; ...
    !     AUX_DIST   = max distance of element ends to axis origin:
    !                  1=distances on axis 1; 2=distances on axis 2; ...
    INTEGER,   DIMENSION(NPERIM,6) :: LIST_ELEMS
    REAL(PDB), DIMENSION(NPERIM,6) :: AUX_DIST
    INTEGER   :: I,J,K,M,ELEM,TYPBC,IBC,IDATA,NBE(6),OK,NAXES
    REAL(PDB) :: ANGLE,X,Y,D
    !****
    NAXES=0
    IF((TYPGEO==5).OR.(TYPGEO==6).OR.(TYPGEO==11)) THEN
       NAXES=4
    ELSEIF((TYPGEO==7).OR.(TYPGEO==8).OR.(TYPGEO==10).OR.(TYPGEO==12)) THEN
       NAXES=3
    ELSEIF(TYPGEO==9) THEN
       NAXES=6
    ENDIF
    !*    calculate number of elements on axes,
    !     and their distances to the origin of the axis
    NBE=0

    DO I=1,NPERIM
       ELEM=PERIM(I)
       IBC=IBC2_ELEM(ELEM)
       TYPBC=TYPE_BC2(IBC)
       IDATA=IDATA_BC2(IBC)
    ENDDO
    DO I=1,NPERIM
       ELEM=PERIM(I)
       IBC=IBC2_ELEM(ELEM)
       TYPBC=TYPE_BC2(IBC)
       IDATA=IDATA_BC2(IBC)
       ANGLE=BCDATA(5,IDATA)
       !        default is the elements not on axis
       M=NAXES+1
       SELECT CASE(TYPGEO)
          CASE(5)
          !*          rectangle with translations:
          !           axes definition:
          !                       M=3
          !                  ************
          !                  *          *
          !              M=2 *          * M=4
          !                  ************
          !                       M=1
          IF(TYPBC==G_BC_TYPE(2)) THEN
             IF(BCDATA(2,IDATA)>0) THEN
                M=1
             ELSEIF(BCDATA(2,IDATA)<0) THEN
                M=3
             ELSEIF(BCDATA(1,IDATA)>0) THEN
                M=2
             ELSEIF(BCDATA(1,IDATA)<0) THEN
                M=4
             ENDIF
          ELSE
             CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=5.')
          ENDIF
          CASE(6)
          !*          rectangle with symmetries:
          !           axes definition:
          !                       M=3
          !                  ************
          !                  *          *
          !              M=2 *          * M=4
          !                  ************
          !                       M=1
          IF(TYPBC==G_BC_TYPE(4)) THEN
             IF(ABS(ANGLE)<EPS) THEN
                IF(BCDATA(2,IDATA)>0) THEN
                   !                    cy>0: element is on the axis 3 (angle=0)
                   M=3
                ELSE
                   !                    cy=0: element is on the axis 1 (angle=0)
                   M=1
                ENDIF
             ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN
                IF(BCDATA(1,IDATA)>0) THEN
                   !                    cx>0: element is on the axis 4 (angle=anggeo)
                   M=4
                ELSE
                   !                    cx=0: element is on the axis 2 (angle=anggeo)
                   M=2
                ENDIF
             ENDIF
          ELSE
             CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=6.')
          ENDIF
          CASE(7,8,12)
          !*          triangles with symmetries:
          !           typgeo=7 axes definition:
          !                        *
          !                 M=2  * * M=3
          !                    *   *
          !                  *******
          !                     M=1
          !           typgeo=8 axes definition:
          !                        *
          !                M=2   *   *    M=3
          !                    *       *
          !                  *************
          !                        M=1
          IF(TYPBC==G_BC_TYPE(4)) THEN
             IF(ABS(ANGLE)<EPS) THEN
                !                 element is on the axis 1 (angle=0)
                M=1
             ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN
                !                 element is on the axis 2 (angle=anggeo)
                M=2
             ELSEIF(ABS(ANGLE-HALFPI)<EPS) THEN
                !                 typgeo=7:element is on the axis 3 (angle=pi/2)
                M=3
             ELSEIF(ABS(ANGLE-2.*ANGGEO)<EPS) THEN
                !                 typgeo=8:element is on the axis 3 (angle=2*anggeo)
                M=3
             ELSEIF(ABS(ANGLE-2.*PI/3.)<EPS) THEN
                !                 typgeo=12:element is on the axis 3 (angle=2*pi/3)
                M=3
             ELSE
                CALL XABORT('SAL130_10: boundary condition data error in element for TYPGEO=7,8,12.')
             ENDIF
          ELSE
             CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=7,8.')
          ENDIF
          CASE(9)
          !*          hexagon with translations:
          !           axes definition:
          !                                   M=4 (0,-2ly)
          !                                  *****
          !                 M=3 (3/2lx,-ly) *     *  M=5 (-3/2lx,-ly)
          !                                *       *
          !                 M=2 (3/2lx, ly) *     *  M=6 (-3/2lx, ly)
          !                                  *****
          !                                   M=1 (0, 2ly)
          !           origins of axes:
          !           axis 1&2: (-lx/2,-ly)
          !           axis 3:   (-lx  ,  0)
          !           axis 4:   (-lx/2, ly)
          !           axis 5:   ( lx  ,  0)
          !           axis 6:   ( lx/2,-ly)
          IF(TYPBC==G_BC_TYPE(2)) THEN
             IF(ABS(BCDATA(1,IDATA))<EPS) THEN
                IF(BCDATA(2,IDATA)>0) THEN
                   M=1
                ELSE
                   M=4
                ENDIF
             ELSEIF(BCDATA(1,IDATA)>0.AND.BCDATA(2,IDATA)>0) THEN
                M=2
             ELSEIF(BCDATA(1,IDATA)>0.AND.BCDATA(2,IDATA)<0) THEN
                M=3
             ELSEIF(BCDATA(1,IDATA)<0.AND.BCDATA(2,IDATA)>0) THEN
                M=6
             ELSEIF(BCDATA(1,IDATA)<0.AND.BCDATA(2,IDATA)<0) THEN
                M=5
             ENDIF
          ELSE
             CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=9.')
          ENDIF
          CASE(10)
          !*          triangles with rotations and translations:
          !           typgeo=10 axes definition:
          !                        *
          !                M=2   *   *    M=3
          !                    *       *
          !                  *************
          !                        M=1
          IF((TYPBC==G_BC_TYPE(2)).OR.(TYPBC==G_BC_TYPE(3))) THEN
             IF(ABS(ANGLE)<EPS) THEN
                !                 element is on the axis 1 (angle=0)
                M=1
             ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN
                !                 element is on the axis 2 (angle=anggeo)
                M=2
             ELSEIF(ABS(ANGLE-2.*ANGGEO)<EPS) THEN
                !                 element is on the axis 3 (angle=2*anggeo)
                M=3
             ELSE
                CALL XABORT('SAL130_10: boundary condition data error in element for TYPGEO=10.')
             ENDIF
          ELSE
             CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=10.')
          ENDIF
          CASE(11)
          !*          lozenge with rotations and translations:
          !           axes definition:
          !                          M=3
          !                     ************
          !                    *          *
          !              M=2  *          *  M=4
          !                  ************
          !                       M=1
          IF((TYPBC==G_BC_TYPE(2)).OR.(TYPBC==G_BC_TYPE(3))) THEN
             IF(ABS(ANGLE)<EPS) THEN
                IF(BCDATA(2,IDATA)>0) THEN
                   !                    cy>0: element is on the axis 3 (angle=0)
                   M=3
                ELSE
                   !                    cy=0: element is on the axis 1 (angle=0)
                   M=1
                ENDIF
             ELSE
                IF(BCDATA(1,IDATA)>0) THEN
                   !                    cx>0: element is on the axis 4 (angle>0)
                   M=4
                ELSE
                   !                    cx=0: element is on the axis 2 (angle>0)
                   M=2
                ENDIF
             ENDIF
          ELSE
             CALL XABORT('SAL130_10: TYPGEO=11:wrong boundary condition type.')
          ENDIF
          CASE DEFAULT
             CALL XABORT('SAL130_10: boundary condition not implemented.')
       END SELECT
       IF(M>=NAXES+1) CALL XABORT('SAL130_10: element not on axes')
       NBE(M)=NBE(M)+1
       LIST_ELEMS(NBE(M),M)=ELEM
       !        sort the element list according their distance to
       !        the origins of the axes
       D=0.
       SELECT CASE(TYPGEO)
          CASE(5:6)
          SELECT CASE(M)
             CASE(1:2)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT(X*X+Y*Y))
             ENDDO
             CASE(3)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,X)
             ENDDO
             CASE(4)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,Y)
             ENDDO
          END SELECT
          CASE(7)
          SELECT CASE(M)
             CASE(1:2)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT(X*X+Y*Y))
             ENDDO
             CASE(3)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,Y)
             ENDDO
          END SELECT
          CASE(8,10,12)
          SELECT CASE(M)
             CASE(1:2)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT(X*X+Y*Y))
             ENDDO
             CASE(3)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                !                 axis origin is (lx,0)
                D=MAX(D,SQRT((X-LX)*(X-LX)+Y*Y))
             ENDDO
          END SELECT
          CASE(9)
          !           origins of axes:
          !           axis 1&2: (-lx/2,-ly)
          !           axis 3:   (-lx  ,  0)
          !           axis 4:   (-lx/2, ly)
          !           axis 5:   ( lx  ,  0)
          !           axis 6:   ( lx/2,-ly)
          SELECT CASE(M)
             CASE(1:2)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X+LX*0.5)*(X+LX*0.5)+(Y+LY)*(Y+LY)))
             ENDDO
             CASE(3)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X+LX)*(X+LX)+Y*Y))
             ENDDO
             CASE(4)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X+LX*0.5)*(X+LX*0.5)+(Y-LY)*(Y-LY)))
             ENDDO
             CASE(5)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X-LX)*(X-LX)+Y*Y))
             ENDDO
             CASE(6)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X-LX*0.5)*(X-LX*0.5)+(Y+LY)*(Y+LY)))
             ENDDO
          END SELECT
          CASE(11)
          SELECT CASE(M)
             CASE(1:2)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT(X*X+Y*Y))
             ENDDO
             CASE(3)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X-LX*0.5)*(X-LX*0.5)+(Y-LY)*(Y-LY)))
             ENDDO
             CASE(4)
             DO K=1,2
                CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K)
                D=MAX(D,SQRT((X-LX)*(X-LX)+Y*Y))
             ENDDO
          END SELECT
       END SELECT
       AUX_DIST(NBE(M),M)=D
       IF(NBE(M)>1) THEN
          DO K=1,NBE(M)-1
             IF(AUX_DIST(NBE(M),M)<AUX_DIST(K,M)) THEN
                !                 insert the last element to position k
                D=AUX_DIST(NBE(M),M)
                ELEM=LIST_ELEMS(NBE(M),M)
                DO J=NBE(M),K+1,-1
                   AUX_DIST(J,M)=AUX_DIST(J-1,M)
                   LIST_ELEMS(J,M)=LIST_ELEMS(J-1,M)
                ENDDO
                AUX_DIST(K,M)=D
                LIST_ELEMS(K,M)=ELEM
                EXIT
             ENDIF
          ENDDO
       ENDIF
    ENDDO
    !     nber of blocks in 'perim'
    M=NAXES
    !*    organize elements as: (elems on axis 1)+(elems on axis 2)+...
    PPERIM(1)=1
    DO I=1,M
       PPERIM(I+1)=PPERIM(I)+NBE(I)
    ENDDO
    DO I=1,M
       PERIM(PPERIM(I):PPERIM(I+1)-1)=LIST_ELEMS(1:NBE(I),I)
    ENDDO
    !*    allocate and set dist_axis
    IF(PPERIM(NAXES+1)-1>0) THEN
       ALLOCATE(DIST_AXIS(PPERIM(NAXES+1)-1),STAT=OK)
       IF(OK.NE.0) CALL XABORT('SAL130_10: not enough memory R')
       DO I=1,NAXES
          DIST_AXIS(PPERIM(I):PPERIM(I+1)-1)=AUX_DIST(1:NBE(I),I)
       ENDDO
    ELSE
       NULLIFY(DIST_AXIS)
    ENDIF
    !
  END SUBROUTINE SAL130_10
  !
  SUBROUTINE SAL131(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! modifies GG%IPAR and constructs bcdata (modifies GG%IPAR and
    ! constructs BCDATA)
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_MAX_LEN,G_BC_LEN,G_BC_TYPE,NIPAR,EPS, &
                                   ANGGEO,TYPGEO,LENGTHX,LENGTHY,ALLSUR
    !****
    IMPLICIT NONE
    ! in variable
    ! ************
    TYPE(T_G_BASIC), INTENT(INOUT) :: GG
    !****
    LOGICAL   :: LGBC,LGALLS
    INTEGER   :: TYPE,II,NN,NBER,I,J,ITBC,IELEM
    REAL(PDB) :: ANGLE
    INTEGER   :: ELEM,IBC,OK,IDATA(2)
    INTEGER, POINTER, DIMENSION(:) :: TYPE_BC2,ISURF2_ELEM,IELEM_BC2
    !*****
    !**   TMP_BCDATA        = description of motions at the boundary
    !**   IPAR              = pointer to geometry descriptors
    INTEGER, DIMENSION(GG%NB_ELEM) :: AUX_ARR
    REAL(PDB), DIMENSION(G_BC_MAX_LEN,GG%NB_ELEM+5) :: TMP_BCDATA
    !*****
    !*    initialization
    IF(TYPGEO==1.OR.TYPGEO==2) IDATA(:)=0
    !*    BCDATA for surfaces of type G_BC_TYPE(-1)
    !
    LGALLS=ALLSUR/=0
    ITBC=0      ! the first bc data
    !
    !*    treat approximate boundary condictions
    IF(LGALLS)THEN
       ITBC=1
       !        ALL BC'S PRODUIT SURFACES
       IF(GG%DEFAUL==1)THEN
          !           SPECULAR REFLEXION -> ISOTROPIC REFLEXION WITH ALBEDO=1
          GG%DEFAUL=0 ; TMP_BCDATA(6,ITBC)=1._PDB
       ENDIF
    ENDIF
    !
    !*    put default value in all bc elements:
    CALL SAL131_2(GG%NB_ELEM,GG%DEFAUL,GG%IPAR,GG%IBC2_ELEM,GG%TYPE_BC2,GG%IDATA_BC2)
    !
    IF(GG%NBBCDA>0)THEN
       ANGLE=0._PDB
       DO I=1,GG%NBBCDA
          ITBC=ITBC+1
          IF(ITBC.GT.GG%NB_ELEM+5) CALL XABORT('SAL131: BCDATA overflow')
          TMP_BCDATA(:,ITBC)=GG%BCDATAREAD(I)%BCDATA(:)
          TYPE=GG%BCDATAREAD(I)%SALTYPE
          SELECT CASE(TYPE)
             CASE(1)
             IF(LGALLS)THEN
                !                specular reflexion -> isotropic reflexion with albedo=1
                TYPE=0 ; TMP_BCDATA(6,ITBC)=1._PDB
             ENDIF
             CASE(2)
             ANGLE=TMP_BCDATA(5,ITBC)
             SELECT CASE(TYPGEO)
                CASE(5) 
                !                adjust translation data according to the sides of rectangle
                IF(TMP_BCDATA(1,ITBC)/=0) TMP_BCDATA(1,ITBC)=SIGN(LENGTHX,TMP_BCDATA(1,ITBC))
                IF(TMP_BCDATA(2,ITBC)/=0) TMP_BCDATA(2,ITBC)=SIGN(LENGTHY,TMP_BCDATA(2,ITBC))
                CASE(9)
                !                adjust translation data according to the hexagonal sides
                IF(ABS(TMP_BCDATA(1,ITBC))<EPS) THEN
                   !                    axes 1 & 4
                   TMP_BCDATA(2,ITBC)=SIGN(2.*LENGTHY,TMP_BCDATA(2,ITBC))
                ELSE
                   TMP_BCDATA(1,ITBC)=SIGN(LENGTHX*1.5,TMP_BCDATA(1,ITBC))
                   TMP_BCDATA(2,ITBC)=SIGN(LENGTHY,TMP_BCDATA(2,ITBC))
                ENDIF
                TMP_BCDATA(3,ITBC)=COS(ANGLE)
                TMP_BCDATA(4,ITBC)=SIN(ANGLE)
                TMP_BCDATA(5,ITBC)=ANGLE
                CASE(11)
                !                adjust translation data according to the hexagonal sides
                IF(ABS(ANGLE)<EPS) THEN
                   IF(TMP_BCDATA(2,ITBC)>0.) THEN
                     !                    cy>0 and angle=0
                     TMP_BCDATA(1,ITBC)=LENGTHX/2.0
                     TMP_BCDATA(2,ITBC)=LENGTHY
                   ELSE
                     !                    cy=0 and angle=0
                     TMP_BCDATA(1,ITBC)=0.0
                     TMP_BCDATA(2,ITBC)=0.0
                   ENDIF
                ELSE
                   IF(TMP_BCDATA(1,ITBC)>0.) THEN
                     !                    cx>0 and angle=pi/3 rad
                     TMP_BCDATA(1,ITBC)=LENGTHX
                     TMP_BCDATA(2,ITBC)=0.0
                   ELSE
                     !                    cx=0 and angle=pi/3 rad
                     TMP_BCDATA(1,ITBC)=0.0
                     TMP_BCDATA(2,ITBC)=0.0
                   ENDIF
                ENDIF
                TMP_BCDATA(3,ITBC)=COS(ANGLE)
                TMP_BCDATA(4,ITBC)=SIN(ANGLE)
                TMP_BCDATA(5,ITBC)=ANGLE
             END SELECT
             IF(LGALLS)THEN
                !                 translation -> isotropic translation with albedo=1
                TYPE=TYPE+4
                !                 perform array shift
                TMP_BCDATA(:,ITBC)=CSHIFT(TMP_BCDATA(:,ITBC),1)
                !                 albedo=1
                TMP_BCDATA(6,ITBC)=1._PDB
             ENDIF
             CASE(3) 
             !              cases of rotation:
             !              read angle, compute cos and sin
             ANGLE=TMP_BCDATA(5,ITBC)
             IF(TYPGEO==1.OR.TYPGEO==2) THEN
                IF(ABS(ANGLE-ANGGEO)<EPS) THEN
                   !                    for the axis 1:keep anggeo,cos(anggeo),sin(anggeo)
                   ANGLE=ANGGEO
                   IF(IDATA(1)==0) IDATA(1)=ITBC
                ELSEIF(ABS(ANGLE+ANGGEO)<EPS) THEN
                   !                    for the axis 2:keep -anggeo,cos(-anggeo),sin(-anggeo)
                   ANGLE=-ANGGEO
                   IF(IDATA(2)==0) IDATA(2)=ITBC
                ELSE
                   CALL XABORT('SAL131: error in angle of rotative axis')
                ENDIF
             ENDIF
             TMP_BCDATA(3,ITBC)=COS(ANGLE)
             TMP_BCDATA(4,ITBC)=SIN(ANGLE)
             TMP_BCDATA(5,ITBC)=ANGLE
             IF(LGALLS)THEN
                !                 rotation -> isotropic rotation with albedo=1
                !                 axial symmetry -> isotropic axial symmetry with albedo=1
                TYPE=TYPE+4
                !                 perform array shift
                TMP_BCDATA(:,ITBC)=CSHIFT(TMP_BCDATA(:,ITBC),1)
                !                 albedo=1
                TMP_BCDATA(6,ITBC)=1._PDB
             ENDIF
             CASE(4)
             !              cases of specular symmetry:
             !              read center and angle, compute cos and sin
             ANGLE=TMP_BCDATA(5,ITBC)
             SELECT CASE(TYPGEO)
                CASE(1,2)
                IF(ABS(ANGLE-0.)<EPS) THEN
                   IF(IDATA(1)==0) IDATA(1)=ITBC
                ELSE
                   IF(IDATA(2)==0) IDATA(2)=ITBC
                ENDIF
                CASE(6:8)
                !                 adjust translation data according to the rectangle/triangle sides
                IF(TMP_BCDATA(1,ITBC)/=0) TMP_BCDATA(1,ITBC)=SIGN(LENGTHX,TMP_BCDATA(1,ITBC))
                IF(TMP_BCDATA(2,ITBC)/=0) TMP_BCDATA(2,ITBC)=SIGN(LENGTHY,TMP_BCDATA(2,ITBC))
             END SELECT
             TMP_BCDATA(3,ITBC)=COS(ANGLE)
             TMP_BCDATA(4,ITBC)=SIN(ANGLE)
             TMP_BCDATA(5,ITBC)=ANGLE
             IF(LGALLS)THEN
                !                 rotation -> isotropic rotation with albedo=1
                !                 axial symmetry -> isotropic axial symmetry with albedo=1
                TYPE=TYPE+4
                !                 perform array shift
                TMP_BCDATA(:,ITBC)=CSHIFT(TMP_BCDATA(:,ITBC),1)
                !                 albedo=1
                TMP_BCDATA(6,ITBC)=1._PDB
             ENDIF
          END SELECT
          !
          !           modify notation for boundary conditions
          NBER=GG%BCDATAREAD(I)%NBER
          DO J=1,NBER
             ELEM=GG%BCDATAREAD(I)%ELEMNB(J)
             IF(ELEM>GG%NB_ELEM.OR.ELEM<=0) CALL XABORT('SAL131: unknown bc element')
             !        get local surface nber
             IBC=GG%IBC2_ELEM(ELEM)
             LGBC=GG%IPAR(2,ELEM)<=0
             II=0
             IF(LGBC)THEN
                II=2
             ELSE
                LGBC=GG%IPAR(3,ELEM)<=0
                IF(LGBC) II=3
             ENDIF
             IF(.NOT.LGBC) THEN
                WRITE(*,*) 'elem :',ELEM
                WRITE(*,*) 'GG%IPAR(:,ELEM) :',GG%IPAR(:,ELEM)
                CALL XABORT('SAL131: wrong bc element')
             ENDIF
             !              put bc type
             GG%IPAR(II,ELEM)=G_BC_TYPE(TYPE)
             GG%TYPE_BC2(IBC)=G_BC_TYPE(TYPE)
             !              put bc data position :
             GG%IDATA_BC2(IBC)=ITBC
          ENDDO
       ENDDO
    ENDIF
    !
    !*    - set BCDATA position for surfaces of type G_BC_TYPE(-1)
    !     - compute the nber of surfaces (type -1,0,-12,-13,-14,-15) : nbsur2
    !     - allocate structures for the surfaces
    !     - compute surf_mac2
    TYPE_BC2=>GG%TYPE_BC2
    ISURF2_ELEM=>GG%ISURF2_ELEM
    IELEM_BC2=>GG%PERIM_MAC2
    NN=0
    DO IBC=1,GG%NB_BC2
       IF(TYPE_BC2(IBC)==G_BC_TYPE(-1)) THEN
          !              macro contact surfaces : set bcdata position to 1
          GG%IDATA_BC2(IBC)=1
       ENDIF
       !           relative element nber
       IELEM=IELEM_BC2(IBC)
       !           count 2D surfaces number
       IF(TYPE_BC2(IBC)==G_BC_TYPE(-1) .OR. TYPE_BC2(IBC)==G_BC_TYPE(0) .OR. &
          TYPE_BC2(IBC)==G_BC_TYPE(1)) THEN
          NN=NN+1
          AUX_ARR(NN)=IBC
          ISURF2_ELEM(IELEM)=NN
       ELSE
          ISURF2_ELEM(IELEM)=0
       ENDIF
    ENDDO
    GG%NB_SURF2=NN
    IF(NN>0) THEN
       ALLOCATE (GG%IBC2_SURF2(NN),GG%IELEM_SURF2(NN),STAT=OK)
       IF(OK/=0) CALL XABORT('SAL131: NOT ENOUGH MEMORY I,R')
       GG%IBC2_SURF2(1:NN)=AUX_ARR(1:NN)
    ELSE
       NULLIFY(GG%IBC2_SURF2,GG%IELEM_SURF2)
    ENDIF
    !
    !*    allocate idata_axis
    IF(TYPGEO==1.OR.TYPGEO==2) THEN
       DO I=1,2
          IF(IDATA(I)==0) THEN
             !              there is no elements on this axis,add a bcdata for this axis
             ITBC=ITBC+1
             ANGLE=0._PDB
             IF(TYPGEO==1) THEN
                !                 symmetry
                IF(I==1) THEN
                   ANGLE=0._PDB
                ELSE
                   ANGLE=ANGGEO
                ENDIF
             ELSEIF(TYPGEO==2) THEN
                !                 rotation
                IF(I==1) THEN
                   ANGLE=ANGGEO
                ELSE
                   ANGLE=-ANGGEO
                ENDIF
             ENDIF
             TMP_BCDATA(1,ITBC)=0._PDB
             TMP_BCDATA(2,ITBC)=0._PDB
             TMP_BCDATA(3,ITBC)=COS(ANGLE)
             TMP_BCDATA(4,ITBC)=SIN(ANGLE)
             TMP_BCDATA(5,ITBC)=ANGLE
             IDATA(I)=ITBC
          ENDIF
       ENDDO
       ALLOCATE(GG%IDATA_AXIS(2), STAT=OK)
       IF(OK.NE.0) CALL XABORT('SAL131: not enough memory I')
       GG%IDATA_AXIS=IDATA
    ENDIF
    !*    allocate bcdata
    ALLOCATE (GG%BCDATA(G_BC_MAX_LEN,ITBC), STAT=OK)
    IF(OK/=0) CALL XABORT('SAL131: not enough memory R')
    GG%BCDATA(:,1:ITBC)=TMP_BCDATA(:,1:ITBC)
    GG%NALBG=ITBC
    !
  END SUBROUTINE SAL131
  !
  SUBROUTINE SAL131_2(NB_ELEM,DEFAUL,IPAR,IBC2_ELEM,TYPE_BC2,IDATA_BC2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! put default bc type and BCDATA to a 2D macro
    !
    !Parameters: input
    ! NB_ELEM    nber of elements
    ! DEFAUL     default bc type
    ! IBC2_ELEM  relative 2D bc number per element
    !
    !Parameters: input/output
    ! IPAR       integer descriptors for elements
    !
    !Parameters: input
    ! TYPE_BC2    2D boundary condiction type
    ! IDATA_BC2   position in bcdata par 2D bc
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE
    !****
    IMPLICIT NONE
    INTEGER, INTENT(IN)                    :: NB_ELEM,DEFAUL
    INTEGER, INTENT(INOUT), DIMENSION(:,:) :: IPAR
    INTEGER, INTENT(IN),   DIMENSION(:)    :: IBC2_ELEM
    INTEGER, INTENT(OUT),  DIMENSION(:)    :: TYPE_BC2
    INTEGER, INTENT(OUT),  DIMENSION(:)    :: IDATA_BC2
    !****
    INTEGER :: ELEM,II,IBC
    LOGICAL :: LGBC
    !****
    !     initiation
    TYPE_BC2=G_BC_TYPE(-1)
    !
    DO ELEM=1,NB_ELEM
       LGBC=IPAR(2,ELEM)==0
       II=0
       IF(LGBC) THEN
          II=2
          IF(IPAR(3,ELEM)<=0) CALL XABORT('SAL131_2: element with 2 bc''s')
       ELSE
          LGBC=IPAR(3,ELEM)==0
          IF(LGBC) II=3
       ENDIF
       IF(LGBC) THEN
          !           put bc type value to ipar
          IPAR(II,ELEM)=G_BC_TYPE(DEFAUL)
          !           put bc type value to bc type structure
          IBC=IBC2_ELEM(ELEM)
          IF(IBC==0) CALL XABORT('SAL131_2: surf-element relation error')
          IF(TYPE_BC2(IBC)/=G_BC_TYPE(-1)) CALL XABORT('SAL131_2: two elements to a surface')
          TYPE_BC2(IBC)=G_BC_TYPE(DEFAUL)
          !           put position of "defaul" in bcdata table (use default albedo)
          IDATA_BC2(IBC)=0
       ENDIF
    ENDDO
    !
  END SUBROUTINE SAL131_2
  !
  SUBROUTINE SAL140(NB_NODE,RPAR,IPAR,PPERIM,PERIM)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! checks domain topology in a 2D macro
    !
    !Parameters: input
    ! NB_NODE    number of nodes in macro
    !
    !Parameters: input/output
    ! RPAR       floating point geometry descriptors
    ! IPAR       integer descriptors for elements
    ! PPERIM     pointer to list of elements in perimeter
    ! PERIM      list of perimeters
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : NRPAR,NIPAR
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)   :: NB_NODE
    INTEGER,   INTENT(IN),    DIMENSION (:)   :: PPERIM
    INTEGER,   INTENT(INOUT), DIMENSION (:,:) :: IPAR
    INTEGER,   INTENT(INOUT), DIMENSION (:)   :: PERIM
    REAL(PDB), INTENT(INOUT), DIMENSION (:,:) :: RPAR
    !****
    INTEGER   :: ELEM,NODE,ID,L,FIRST,LAST,NEXT,NLOOP, &
         NEWEND,KEEP,IEND,I1,I2
    REAL(PDB) :: X,Y,XNEW,YNEW,DIST
    LOGICAL   :: LGOPEN,LGON,LGERR
    INTEGER, PARAMETER      :: MXKEEP = 20
    REAL, DIMENSION(MXKEEP) :: KEEPD
    INTEGER, PARAMETER :: FOUT =6
    REAL(PDB),PARAMETER :: EPS2=1.0E-7_PDB
    !*****
    NEXT=0
    !*    checks topology of each node
    DO NODE=1,NB_NODE
       I1=PPERIM(NODE)
       !        ID counts elements in the perimeter already treated
       ID=I1-1
       I2=PPERIM(NODE+1)-1
       IF(I2<I1) CALL XABORT('SAL140: node without perimeter')
       !        mark elements in perimeter of node as untreated
       !        (except full circle: ipar=2)
       DO L=I1,I2
          ELEM=PERIM(L)
          IF(IPAR(1,ELEM)==2)THEN
             ID=ID+1
             IF(ID<L)THEN
                ! move circles at beginning
                PERIM(L)=PERIM(ID)
                PERIM(ID)=ELEM
             ENDIF
          ENDIF
       ENDDO
       ID=ID+1
       !        nloop counts nber of loops for node
       NLOOP=ID-1
       !        treat all elements in perimeter until they have been all done:
       DO WHILE (ID<I2)
          !           an element in the perimeter that has not been done
          FIRST=PERIM(ID)
          NLOOP=NLOOP+1
          !           get first end for this element:
          IEND=1
          !           determine the loop defined by element last:
          LGOPEN=.TRUE.
          LAST=FIRST
          DO WHILE(LGOPEN)
             !              get coordinates x and y of end of element last
             CALL SAL141(IPAR(1,LAST),RPAR(:,LAST),X,Y,IEND)
             !              get next element in the perimeter in touch with
             !              last element at its end iend
             LGON=.TRUE.
             L=ID+1
             KEEP=0
             DO WHILE(LGON.AND.L<=I2)
                NEXT=PERIM(L)
                !                 check whether next is really 'next' to last
                CALL SAL142(X,Y,XNEW,YNEW,IPAR(1,NEXT),RPAR(:,NEXT),NEWEND,EPS2,DIST)
                IF(KEEP<MXKEEP)KEEP=KEEP+1
                KEEPD(KEEP)=REAL(DIST)
                LGON=NEWEND<0
                IF(LGON)L=L+1
             ENDDO
             LGOPEN=.NOT.LGON
             IF(LGOPEN)THEN
                !                 replace last element in order in the perimeter
                ID=ID+1
                IF(ID<L)THEN
                   PERIM(L)=PERIM(ID)
                   PERIM(ID)=NEXT
                ENDIF
             ELSE
                LGERR=.FALSE.
                !                 try to close with first element of loop:
                CALL SAL142(X,Y,XNEW,YNEW,IPAR(1,FIRST),RPAR(:,FIRST),NEWEND,EPS2,DIST)
                LGOPEN=NEWEND<0
                IF(LGOPEN)THEN
                   !                    fatal error: cannot close loop for this node
                   WRITE(FOUT,'(//,''==> cannot close node '',i5, &
                        &''   AT ELEMENT'',I5,/)')NODE,LAST
                   LGERR=.TRUE.
                ELSEIF(FIRST==LAST)THEN
                   WRITE(FOUT,'(//,''==> node '',i5,''    with '', &
                        &''isolated element'',i5,/)')NODE,LAST
                   LGERR=.TRUE.
                ENDIF
                IF(LGERR) CALL XABORT('SAL140: node not closed')
                NEXT=FIRST
             ENDIF
             !              define last = next element and proceed
             LAST=NEXT
             IEND=3-NEWEND
          ENDDO
          ID=ID+1
       ENDDO
    ENDDO
    !
  END SUBROUTINE SAL140
  !
  SUBROUTINE SAL142(X,Y,XNEW,YNEW,TYPE,RPAR,IEND,EPS2,DIST)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! checks whether element is very close to point (X,Y)
    !
    !Parameters: input
    ! X        abscissa coordinate
    ! Y        ordinate coordinate
    ! TYPE     type of element  1 (segment), 3 (arc of circle)
    ! EPS2     criterium for closeness
    !
    !Parameters: input/output
    ! XNEW     abscissa coordinate of end of element close to point
    ! YNEW     ordinate coordinate of end of element close to point
    ! RPAR     floating point geometry descriptors
    ! IEND     =  1 (beginning is close to point dist < EPS2)
    !             2 (end is close to point DIST < EPS2)
    !            -1 (beginning is close to point dist > EPS2)
    !            -2 (end is close to point DIST > EPS2)
    ! DIST     distance from end of element to point (X,Y)
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : G_ELE_TYPE
    !****
    IMPLICIT NONE
    INTEGER,   INTENT(IN)               :: TYPE
    INTEGER,   INTENT(OUT)              :: IEND
    REAL(PDB), INTENT(IN)               :: X,Y,EPS2
    REAL(PDB), INTENT(OUT)              :: XNEW,YNEW,DIST
    REAL(PDB), INTENT(IN), DIMENSION(:) :: RPAR
    !     DIMENSION        RPAR(NRPAR)
    !****
    REAL(PDB) :: CX,CY,THETA,R,DIST2
    !****
    !*    function giving distance between two points:
    REAL(PDB) :: FUNC,X1,Y1,X2,Y2
    FUNC(X1,Y1,X2,Y2)=(X1-X2)**2+(Y1-Y2)**2
    !****
    CX=RPAR(1)
    CY=RPAR(2)
    DIST2=0._PDB
    IF(TYPE==G_ELE_TYPE(1))THEN
       !*       segment
       XNEW=CX
       YNEW=CY
       DIST=FUNC(XNEW,YNEW,X,Y)
       IF(DIST<=EPS2)THEN
          IEND=1
          RETURN
       ELSE
          XNEW=CX+RPAR(3)
          YNEW=CY+RPAR(4)
          DIST2=FUNC(XNEW,YNEW,X,Y)
       ENDIF
    ELSEIF(TYPE<=G_ELE_TYPE(3))THEN
       !*       arc of circle
       R=RPAR(3)
       THETA=RPAR(4)
       XNEW=CX+R*COS(THETA)
       YNEW=CY+R*SIN(THETA)
       DIST=FUNC(XNEW,YNEW,X,Y)
       IF(DIST<=EPS2)THEN
          IEND=1
          RETURN
       ELSE
          THETA=RPAR(5)
          XNEW=CX+R*COS(THETA)
          YNEW=CY+R*SIN(THETA)
          DIST2=FUNC(XNEW,YNEW,X,Y)
       ENDIF
    ELSE
       CALL XABORT('SAL142: not implemented')
    ENDIF
    !
    IF(DIST2<=EPS2)THEN
       DIST=DIST2
       IEND=2
       RETURN
    ELSE
       IF(DIST<=DIST2)THEN
          IEND=-1
       ELSE
          DIST=DIST2
          IEND=-2
       ENDIF
    ENDIF
    !
  END SUBROUTINE SAL142
  !
  SUBROUTINE SAL160(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! analyse domain definition: 2D volumes, surfaces
    !  - compute node volumes
    !  - compute areas of 2d surfaces
    !  - read medium data
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE,NBMED
    !****
    IMPLICIT NONE
    TYPE(T_G_BASIC), INTENT(INOUT)  :: GG
    !****
    INTEGER :: INB,IMED
    INTEGER, DIMENSION(GG%NB_NODE) :: DATAIN
    !****
    !     SUBROUTINE SAL160_2(NB_ELEM,IPAR,RPAR,VOL2,ISURF2_ELEM,NB_SURF2,SURF2)
    IF(GG%NB_SURF2==0) THEN
      CALL SAL160_2(GG%NB_ELEM,GG%IPAR,GG%RPAR,GG%VOL_NODE,GG%ISURF2_ELEM, &
         GG%NB_SURF2)
    ELSE
      CALL SAL160_2(GG%NB_ELEM,GG%IPAR,GG%RPAR,GG%VOL_NODE,GG%ISURF2_ELEM, &
         GG%NB_SURF2,GG%SURF2)
    ENDIF
    !
    !*    read medium per region
    CALL SALGET(DATAIN,GG%NB_NODE,F_GEO,FOUT0,'media per node')
  
    ! number of media fixed to maximum of datain
    NBMED = MAXVAL(DATAIN(1:GG%NB_NODE))
    !
    !*    check and define med for code regions
    DO INB=1,GG%NB_NODE
       IMED=DATAIN(INB)
       IF(IMED>NBMED.OR.IMED<0)THEN
          WRITE(*,*) 'medium : ',IMED
          WRITE(*,*) 'inb, nbmed, nbnode : ',INB,NBMED,GG%NB_NODE
          CALL XABORT('SAL160: wrong medium in a region')
       ENDIF
       GG%MED(INB)=IMED
    ENDDO
    !
  END SUBROUTINE SAL160
  !
  SUBROUTINE SAL160_2(NB_ELEM,IPAR,RPAR,VOL2,ISURF2_ELEM,NB_SURF2,SURF2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! compute 2D volumes and surfaces
    !  - compute node volumes
    !  - compute 2D surface areas
    !
    !Parameters: input
    ! NB_ELEM      number of elements
    ! NB_NODE      number of nodes
    ! IPAR         integer descriptors for elements
    ! RPAR         floating point descriptors for elements
    ! ISURF2_ELEM  2D surface nber per elem
    ! NB_SURF2     number of 2D surface
    !
    !Parameters: output
    ! VOL2         2D volumes of node 
    ! SURF2        2D areas of node
    !
    !---------------------------------------------------------------------
    !
    USE PRECISION_AND_KINDS, ONLY : PDB
    USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE
    !****
    IMPLICIT NONE
    INTEGER, INTENT(IN)                    :: NB_ELEM
    INTEGER, INTENT(IN),   DIMENSION(:,:)  :: IPAR
    REAL(PDB), INTENT(IN), DIMENSION(:,:)  :: RPAR
    REAL(PDB), INTENT(OUT),  DIMENSION(:)  :: VOL2
    INTEGER, INTENT(IN),   DIMENSION(:)    :: ISURF2_ELEM
    INTEGER, INTENT(IN)                    :: NB_SURF2
    REAL(PDB), INTENT(OUT), OPTIONAL,  DIMENSION(:)  :: SURF2
    !****
    INTEGER :: ELEM,NODEBC,NODE,ISURF
    LOGICAL :: LGBC
    REAL(PDB) :: AUX
    !****
    !     initiation
    VOL2=0._PDB
    IF(NB_SURF2 > 0) SURF2=0._PDB
    DO ELEM=1,NB_ELEM
       !*       compute volume of node and add to volume of region (local)
       CALL SAL161(IPAR(1,ELEM),RPAR(:,ELEM),AUX)
       NODEBC=IPAR(2,ELEM)
       LGBC=NODEBC<=0
       IF(.NOT.LGBC) VOL2(NODEBC)=VOL2(NODEBC)+AUX
       NODE=IPAR(3,ELEM)
       IF(NODE>0)THEN
          VOL2(NODE)=VOL2(NODE)-AUX
       ELSE
          IF(LGBC) CALL XABORT('SAL160_2: isolated element')
          LGBC=.TRUE.
          NODEBC=NODE
       ENDIF
       IF(LGBC) THEN
          IF(NODEBC==G_BC_TYPE(-1).OR.NODEBC==G_BC_TYPE(0).OR.NODEBC==G_BC_TYPE(1))THEN
             !  compute external surface (loaded in total order number)
             ISURF=ISURF2_ELEM(ELEM)
             IF(ISURF==0) CALL XABORT('SAL160_2: wrong relation of element - surf ')
             CALL SAL162(IPAR(1,ELEM),RPAR(:,ELEM),SURF2(ISURF))
          ENDIF
       ENDIF
    ENDDO
    !
  END SUBROUTINE SAL160_2
  !
  SUBROUTINE SAL161(TYPE,RPAR,VOL2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! computes 'volume' between an element and the x axis
    !
    !Parameters: input
    ! TYPE         type of element
    ! RPAR         floating point descriptors for elements
    !
    !Parameters: output
    ! VOL2         '2D area' between the element and the horizontal axis
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    INTEGER,   INTENT(IN)                 :: TYPE
    REAL(PDB), INTENT(IN),  DIMENSION (:) :: RPAR
    REAL(PDB), INTENT(OUT)                :: VOL2
    !     DIMENSION        RPAR(*)
    !****
    REAL(PDB) :: YC,EX,EY,R,PHI1,PHI2,COS1,COS2
    !****
    !     volume is added to node- and substracted from node+
    YC=RPAR(2)
    IF(TYPE==1)THEN
       !        segment:
       EX=RPAR(3)
       EY=RPAR(4)
       VOL2=EX*(YC+EY/2.)
    ELSEIF(TYPE==2)THEN
       !        whole circle
       VOL2=PI*RPAR(3)*RPAR(3)
    ELSEIF(TYPE==3)THEN
       !        arc of circle:
       R=RPAR(3)
       PHI1=RPAR(4)
       PHI2=RPAR(5)
       COS1=COS(PHI1)
       COS2=COS(PHI2)
       VOL2=R*(YC*(COS1-COS2)+(R/2.0)*(PHI2-PHI1+ &
            & (SIN(PHI1)*COS1-SIN(PHI2)*COS2)))
    ELSE
       CALL XABORT('SAL161: not implemented')
    ENDIF
    !
  END SUBROUTINE SAL161
  !
  SUBROUTINE SAL162(TYPE,RPAR,SURF2)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! computes surface of a boundary element
    !
    !Parameters: input
    ! TYPE         type of element
    ! RPAR         floating point descriptors for elements
    !
    !Parameters: output
    ! SURF2        '2D length' of element
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    INTEGER,   INTENT(IN)                :: TYPE
    REAL(PDB), INTENT(IN), DIMENSION (:) :: RPAR
    REAL(PDB), INTENT(OUT)               :: SURF2
    !****
    IF(TYPE==1)THEN
       !        segment
       SURF2=RPAR(5)
    ELSEIF(TYPE<=3)THEN
       !        arc of circle
       SURF2=(RPAR(5)-RPAR(4))*RPAR(3)
    ELSE
       CALL XABORT('SAL162: not implemented')
    ENDIF
    !
  END SUBROUTINE SAL162
  !
  SUBROUTINE SAL170(GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! prints out volume, surface and medium info for 2D macros
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES,      ONLY : G_BC_TYPE
    USE SAL_TRACKING_TYPES,      ONLY : PRTIND
    !****
    IMPLICIT NONE
    TYPE(T_G_BASIC), INTENT(INOUT)  :: GG
    INTEGER, PARAMETER :: FOUT =6

    !****
     IF(PRTIND == 1) THEN
         WRITE(FOUT,'(//,10x,''2D geometry'',/,10X,11(''=''),//, &
         &I8,''  regions'',/, &
         &I8,''  external surfaces'',//)') GG%NB_NODE,GG%NB_SURF2
     ENDIF
    !
  END SUBROUTINE SAL170
  !
  SUBROUTINE SALSYM(X1,Y1,X2,Y2,X0,Y0,X4,Y4)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! set the symmetric point (x4,y4) of (x0,y0) relative to symmetry axis
    ! (x1,y1)->(x2,y2)
    !
    !Parameters: input
    ! X1   abscissa coordinate of a point on the symmetry axis
    ! Y1   ordinate coordinate of a point on the symmetry axis
    ! X2   abscissa coordinate of a point on the symmetry axis
    ! Y2   ordinate coordinate of a point on the symmetry axis
    ! X0   abscissa coordinate of the point to mirror
    ! Y0   ordinate coordinate of the point to mirror
    !
    !Parameters: output
    ! X4   abscissa coordinate of the symmetric point
    ! Y4   ordinate coordinate of the symmetric point
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL(PDB),INTENT(IN)  :: X1,Y1,X2,Y2,X0,Y0
    REAL(PDB),INTENT(OUT) :: X4,Y4

    REAL(PDB) :: A,B,C,DEN,X3,Y3
    A=Y2-Y1; B=X1-X2; C=X2*Y1-X1*Y2;
    DEN=A*A+B*B;
    IF(DEN==0._PDB) CALL XABORT('SALSYM: division by zero')
    X3=(B*(B*X0-A*Y0)-A*C)/DEN;
    Y3=(A*(-B*X0+A*Y0)-B*C)/DEN;
    X4=X0+2._PDB*(X3-X0); Y4=Y0+2._PDB*(Y3-Y0);
  END SUBROUTINE SALSYM
  !
  SUBROUTINE SALROT(X0,Y0,THETA,X3,Y3,X4,Y4)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! set the symmetric point (x4,y4) of (x3,y3) relative to the rotation
    ! center (x1,y1) with angle theta.
    !
    !Parameters: input
    ! X0    abscissa coordinate of the rotation center
    ! Y0    ordinate coordinate of the rotation center
    ! THETA rotation angle
    ! X3    abscissa coordinate of the point to mirror
    ! Y3    ordinate coordinate of the point to mirror
    !
    !Parameters: output
    ! X4    abscissa coordinate of the symmetric point
    ! Y4    ordinate coordinate of the symmetric point
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL(PDB),INTENT(IN)  :: X0,Y0,THETA,X3,Y3
    REAL(PDB),INTENT(OUT) :: X4,Y4

    REAL(PDB) :: A,B
    IF(THETA==0.0_PDB) THEN
      X4=X3; Y4=Y3
    ELSE
      A=X3-X0; B=Y3-Y0
      X4=X0+A*COS(THETA)-B*SIN(THETA); Y4=Y0+A*SIN(THETA)+B*COS(THETA)
    ENDIF
  END SUBROUTINE SALROT
  !
  SUBROUTINE SALFOLD_0(GG,IPASS,IB,NBBCDA,ALIGN,LFOLD,IFOLD)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! unfold the domain with reflection relative to axis AXIS_XY
    !
    !Parameters: input
    ! IB     actual unfolding axis
    ! NBBCDA number of perimeters before unfolding
    ! ALIGN  unfolding axes
    ! LFOLD  identification flag to all unfolding axes
    !
    !Parameters: input/output
    ! GG     geometry descriptor
    !
    !Parameters: output
    ! IFOLD folded element indices corresponding to unfolded ones
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,LENGTHX,LENGTHY
    USE SAL_NUMERIC_MOD,  ONLY : FINDLC,DET_ROSETTA
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: IPASS,IB,NBBCDA
    LOGICAL, DIMENSION(NBBCDA), INTENT(IN) :: LFOLD
    REAL(PDB), DIMENSION(3,3,NBBCDA), INTENT(INOUT) :: ALIGN
    TYPE(T_G_BASIC), INTENT(INOUT) :: GG
    INTEGER, DIMENSION(:), INTENT(OUT) :: IFOLD
    !
    ! AXIS_XY values of AXIS_X1, AXIS_Y1, AXIS_X2 and AXIS_Y2 for the
    !         reflecting axis
    INTEGER :: ELEM,TYPE,OK,TMP_NB_ELEM,TMP_NBBCDA,I,J,IBC,INDBC,IAUX
    INTEGER, DIMENSION(3) :: IPAR_TMP
    REAL(PDB), DIMENSION(4) :: AXIS_XY
    REAL(PDB), DIMENSION(6) :: RPAR_TMP
    REAL(PDB),PARAMETER :: EPS=1.0E-5_PDB
    REAL(PDB) :: X1,X2,X4,Y1,Y2,Y4,DX4,DY4,RAD,THETA1,THETA2,X1B,Y1B,X4B, &
                 Y4B,XMIN,YMIN,XMAX,YMAX,PHI1,PHI2,DELPHI,DET1,DET2
    !
    ! allocatable arrays
    INTEGER, ALLOCATABLE, DIMENSION(:) :: PERIM_ELEM
    LOGICAL, ALLOCATABLE, DIMENSION(:) :: ISPERIM
    INTEGER, POINTER, DIMENSION(:,:) :: TMP_IPAR
    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: I2
    REAL(PDB), ALLOCATABLE, DIMENSION(:) :: ANGLE,ALBEDO
    REAL(PDB), POINTER, DIMENSION(:,:) :: TMP_RPAR
    REAL(PDB), ALLOCATABLE, DIMENSION(:,:,:) :: ALIGN2
    TYPE(T_SALBCDATA), POINTER, DIMENSION(:) :: TMP_BCDATAREAD
    !
    ! compute size of the unfold geometry
    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
    DO ELEM=1,GG%NB_ELEM
      TYPE=GG%IPAR(1,ELEM)
      IF(TYPE==1) THEN
        X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
        XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
        X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
        XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2);
      ENDIF
    ENDDO
    LENGTHX=XMAX-XMIN; LENGTHY=YMAX-YMIN;
    !
    ! allocate new surfacic element containers
    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
    ALLOCATE(TMP_IPAR(NIPAR,3*GG%NB_ELEM), TMP_RPAR(NRPAR,3*GG%NB_ELEM), &
             I2(2,GG%NB_ELEM), STAT=OK)
    IF(OK/=0) CALL XABORT('SALFOLD_0: not enough memory')
    TMP_IPAR(:,:)=0; TMP_RPAR(:,:)=0._PDB;
    !
    ! loop over old elements
    TMP_NB_ELEM=0
    THETA1=0._PDB; THETA2=0._PDB;
    I2(:2,:GG%NB_ELEM)=0
    AXIS_XY(1)=ALIGN(1,1,IB) ; AXIS_XY(2)=ALIGN(1,2,IB)
    AXIS_XY(3)=ALIGN(2,1,IB) ; AXIS_XY(4)=ALIGN(2,2,IB)
    OUT1: DO ELEM=1,GG%NB_ELEM
      TYPE=GG%IPAR(1,ELEM)
      X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); RAD=GG%RPAR(3,ELEM)
      IF(TYPE==1) THEN
        XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
        X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM)
        ! Cycle if this element is sitting on an unfolding axe
        DO IBC=1,NBBCDA
          IF(.NOT.LFOLD(IBC)) CYCLE
          IF((IPASS.EQ.1).AND.(IBC.NE.IB)) CYCLE
          ALIGN(3,1,IBC)=X1; ALIGN(3,2,IBC)=Y1
          DET1 = DET_ROSETTA(ALIGN(1,1,IBC),3)
          ALIGN(3,1,IBC)=X2; ALIGN(3,2,IBC)=Y2;
          DET2 = DET_ROSETTA(ALIGN(1,1,IBC),3)
          IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) CYCLE OUT1
        ENDDO
        !
        CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1,Y1,X4,Y4)
        CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X2,Y2,DX4,DY4)
      ELSE IF(TYPE==2) THEN
        CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1,Y1,X4,Y4)
        THETA1=0._PDB; THETA2=0._PDB;
      ELSE IF(TYPE==3) THEN
        CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1,Y1,X4,Y4)
        X1B=X1+RAD*COS(GG%RPAR(4,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(4,ELEM));
        CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1B,Y1B,X4B,Y4B)
        IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN
          THETA1=PI/2._PDB
        ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN
          THETA1=3._PDB*PI/2._PDB
        ELSE IF(X4B-X4 > 0._PDB) THEN
          THETA1=ATAN((Y4B-Y4)/(X4B-X4))
        ELSE
          THETA1=ATAN((Y4B-Y4)/(X4B-X4))+PI
        ENDIF
        X1B=X1+RAD*COS(GG%RPAR(5,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(5,ELEM));
        CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1B,Y1B,X4B,Y4B)
        IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN
          THETA2=PI/2._PDB
        ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN
          THETA2=3._PDB*PI/2._PDB
        ELSE IF(X4B-X4 > 0._PDB) THEN
          THETA2=ATAN((Y4B-Y4)/(X4B-X4))
        ELSE
          THETA2=ATAN((Y4B-Y4)/(X4B-X4))+PI
        ENDIF
      ELSE
        WRITE(*,*) " elem=",ELEM," type=",TYPE
        CALL XABORT('SALFOLD_0: invalid type of surfacic element')
      ENDIF
      IPAR_TMP(:3)=0
      RPAR_TMP(:6)=0_PDB
      RPAR_TMP(1)=X4; RPAR_TMP(2)=Y4;
      IPAR_TMP(1)=TYPE;
      IF(TYPE==1) THEN
        RPAR_TMP(3)=DX4-X4; RPAR_TMP(4)=DY4-Y4;
        RPAR_TMP(5)=SQRT(RPAR_TMP(3)**2+RPAR_TMP(4)**2)
        XMIN=MIN(XMIN,X4); YMIN=MIN(YMIN,Y4); XMAX=MAX(XMAX,X4); YMAX=MAX(YMAX,Y4);
        XMIN=MIN(XMIN,DX4); YMIN=MIN(YMIN,DY4); XMAX=MAX(XMAX,DX4); YMAX=MAX(YMAX,DY4);
        IPAR_TMP(2)=GG%IPAR(3,ELEM); IPAR_TMP(3)=GG%IPAR(2,ELEM);
      ELSE IF((TYPE==2).OR.(TYPE==3)) THEN
        RPAR_TMP(3)=GG%RPAR(3,ELEM) ! RADIUS
        IF(THETA2>THETA1) THETA1=THETA1+2._PDB*PI
        PHI1=THETA2; DELPHI=THETA1-THETA2;
        IF(DELPHI>0._PDB)THEN
          PHI2=PHI1+DELPHI
        ELSE
          PHI2=PHI1
          PHI1=PHI1+DELPHI
        ENDIF
        IF(TYPE==3)THEN
          ! arc of circle: put phi1 within 0 and 2*pi
          IF(PHI1>2._PDB*PI)THEN
            IAUX=INT(PHI1/(2._PDB*PI))
            DELPHI=(2._PDB*PI)*IAUX
            PHI1=PHI1-DELPHI ; PHI2=PHI2-DELPHI
          ELSEIF(PHI1<0._PDB)THEN
            IAUX=INT((-PHI1+1.D-7)/(2._PDB*PI))+1
            DELPHI=(2._PDB*PI)*IAUX
            PHI1=PHI1+DELPHI ; PHI2=PHI2+DELPHI
          ENDIF
        ENDIF
        RPAR_TMP(4)=PHI1; RPAR_TMP(5)=PHI2; ! ANGLES
        IPAR_TMP(2)=GG%IPAR(2,ELEM); IPAR_TMP(3)=GG%IPAR(3,ELEM)
      ENDIF
      RPAR_TMP(6)=0._PDB
      IF(IPASS==2) THEN
        ! remove identical elements at pass 2
        DO I=1,TMP_NB_ELEM
          IF((ABS(TMP_RPAR(1,I)-RPAR_TMP(1))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(2,I)-RPAR_TMP(2))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(3,I)-RPAR_TMP(3))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(4,I)-RPAR_TMP(4))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(5,I)-RPAR_TMP(5))<=10.0*EPS)) THEN
             CYCLE OUT1
          ENDIF
        ENDDO
      ENDIF
      TMP_NB_ELEM=TMP_NB_ELEM+1
      IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_0: tmp_nb_elem overflow(1)')
      TMP_IPAR(:3,TMP_NB_ELEM)=IPAR_TMP(:3)
      TMP_RPAR(:6,TMP_NB_ELEM)=RPAR_TMP(:6)
      I2(2,ELEM)=TMP_NB_ELEM
      IF(IPASS==2) THEN
        ! remove identical elements at pass 2
        DO I=1,TMP_NB_ELEM
          IF((ABS(TMP_RPAR(1,I)-GG%RPAR(1,ELEM))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(2,I)-GG%RPAR(2,ELEM))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(3,I)-GG%RPAR(3,ELEM))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(4,I)-GG%RPAR(4,ELEM))<=10.0*EPS).AND. &
             (ABS(TMP_RPAR(5,I)-GG%RPAR(5,ELEM))<=10.0*EPS)) THEN
             CYCLE OUT1
          ENDIF
        ENDDO
      ENDIF
      TMP_NB_ELEM=TMP_NB_ELEM+1
      IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_0: tmp_nb_elem overflow(2)')
      TMP_IPAR(:,TMP_NB_ELEM)=GG%IPAR(:,ELEM)
      TMP_RPAR(:,TMP_NB_ELEM)=GG%RPAR(:,ELEM)
      I2(1,ELEM)=TMP_NB_ELEM
    ENDDO OUT1
    DEALLOCATE(GG%IPAR,GG%RPAR)
    DO ELEM=1,GG%NB_ELEM
      IF(I2(1,ELEM).GT.2*GG%NB_ELEM) CALL XABORT('SALFOLD_0: IFOLD overflow')
      IF(I2(1,ELEM).NE.0) IFOLD(I2(1,ELEM))=ELEM
      IF(I2(2,ELEM).NE.0) IFOLD(I2(2,ELEM))=ELEM
    ENDDO
    GG%IPAR=>TMP_IPAR; GG%RPAR=>TMP_RPAR;
    GG%NB_ELEM=TMP_NB_ELEM
    !
    ! loop over boundary conditions
    ALLOCATE(ISPERIM(GG%NB_ELEM),ALIGN2(3,3,GG%NB_ELEM),ANGLE(GG%NB_ELEM), &
    & ALBEDO(GG%NB_ELEM),PERIM_ELEM(GG%NB_ELEM))
    ALIGN2(:3,3,:GG%NB_ELEM)=1.0_PDB
    PERIM_ELEM(:GG%NB_ELEM)=0
    ISPERIM(:GG%NB_ELEM)=.FALSE.
    TMP_NBBCDA=0
    DO IBC=1,GG%NBBCDA
      DO I=1,GG%BCDATAREAD(IBC)%NBER
        INDBC=GG%BCDATAREAD(IBC)%ELEMNB(I)
        IF(INDBC==0) CYCLE
        IF(I2(1,INDBC)/=0) ISPERIM(I2(1,INDBC))=.TRUE.
        IF(I2(2,INDBC)/=0) ISPERIM(I2(2,INDBC))=.TRUE.
      ENDDO
    ENDDO
    ITER0: DO ELEM=1,GG%NB_ELEM
      IF(.NOT.ISPERIM(ELEM)) CYCLE
      X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
      X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
      DO J=1,TMP_NBBCDA
        ALIGN2(3,1,J)=X1; ALIGN2(3,2,J)=Y1;
        DET1 = DET_ROSETTA(ALIGN2(1,1,J),3)
        ALIGN2(3,1,J)=X2; ALIGN2(3,2,J)=Y2;
        DET2 = DET_ROSETTA(ALIGN2(1,1,J),3)
        IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) THEN
          PERIM_ELEM(ELEM) = J
          CYCLE ITER0
        ENDIF
      ENDDO
      TMP_NBBCDA=TMP_NBBCDA+1
      PERIM_ELEM(ELEM) = TMP_NBBCDA
      ANGLE(TMP_NBBCDA)=ATAN((Y2-Y1)/(X2-X1))
      IF(ABS(ANGLE(TMP_NBBCDA)).LE.1.0E-5) ANGLE(TMP_NBBCDA)=0.0
      ALIGN2(1,1,TMP_NBBCDA)=X1; ALIGN2(1,2,TMP_NBBCDA)=Y1
      ALIGN2(2,1,TMP_NBBCDA)=X2; ALIGN2(2,2,TMP_NBBCDA)=Y2
      ! Recover albedo from folded geometry
      ALBEDO(TMP_NBBCDA)=1.0
      DO IBC=1,GG%NBBCDA
        J = FINDLC(GG%BCDATAREAD(IBC)%ELEMNB,ELEM)
        IF(J.EQ.1) THEN
          ALBEDO(TMP_NBBCDA)=GG%BCDATAREAD(IBC)%BCDATA(6)
          EXIT
        ENDIF
      ENDDO
    ENDDO ITER0
    ALLOCATE(TMP_BCDATAREAD(TMP_NBBCDA))
    DO IBC=1,TMP_NBBCDA
      TMP_BCDATAREAD(IBC)%NBER = COUNT(PERIM_ELEM(:GG%NB_ELEM) == IBC)
      ALLOCATE(TMP_BCDATAREAD(IBC)%ELEMNB(TMP_BCDATAREAD(IBC)%NBER))
      TMP_BCDATAREAD(IBC)%SALTYPE = 0
      J=0
      DO I=1,GG%NB_ELEM
        IF(PERIM_ELEM(I) == IBC) THEN
          J=J+1
          TMP_BCDATAREAD(IBC)%ELEMNB(J) = I
        ENDIF
      ENDDO
      TMP_BCDATAREAD(IBC)%BCDATA(1) = ALIGN2(1,1,IBC)
      TMP_BCDATAREAD(IBC)%BCDATA(2) = ALIGN2(1,2,IBC)
      TMP_BCDATAREAD(IBC)%BCDATA(3) = COS(ANGLE(IBC))
      TMP_BCDATAREAD(IBC)%BCDATA(4) = SIN(ANGLE(IBC))
      TMP_BCDATAREAD(IBC)%BCDATA(5) = ANGLE(IBC)
      TMP_BCDATAREAD(IBC)%BCDATA(6) = ALBEDO(IBC)
    ENDDO
    DEALLOCATE(I2,PERIM_ELEM,ALBEDO,ANGLE,ALIGN2,ISPERIM)
    DEALLOCATE(GG%BCDATAREAD)
    GG%BCDATAREAD=>TMP_BCDATAREAD
    GG%NBBCDA=TMP_NBBCDA
    GG%ALBEDO=1.D0
  END SUBROUTINE SALFOLD_0
  !
  SUBROUTINE SALFOLD_1(HSYM,GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! unfold the domain with reflection
    !
    !Parameters: input
    ! HSYM: type of symmetry: 'DIAG': diagonal symmetry; 'SYMX': symmetry
    !       relative to X axis; 'SYMY': symmetry relative to Y axis;
    !       'SA60': SA60 symmetry; 'S30': S30 symmetry
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,LENGTHX,LENGTHY
    USE SAL_NUMERIC_MOD,  ONLY : FINDLC,DET_ROSETTA
    IMPLICIT NONE
    CHARACTER(LEN=4),INTENT(IN)  :: HSYM
    TYPE(T_G_BASIC), INTENT(INOUT) :: GG
    !
    INTEGER :: ELEM,TYPE,OK,TMP_NB_ELEM,TMP_NBBCDA,I,J,IBC,INDBC,IAUX,ISYM,NSYM
    REAL(PDB),PARAMETER :: EPS=1.0E-5_PDB
    REAL(PDB) :: AXIS_X1(2),AXIS_X2(2),AXIS_Y1(2),AXIS_Y2(2)
    REAL(PDB) :: X1,X2,X4,Y1,Y2,Y4,DX4,DY4,RAD,THETA1,THETA2,X1B,Y1B,X4B, &
                 Y4B,XMIN,YMIN,XMAX,YMAX,PHI1,PHI2,DELPHI,DET1,DET2,ALIGN3(3,3)
    LOGICAL :: NOCOPY(2)
    !
    ! allocatable arrays
    INTEGER, ALLOCATABLE, DIMENSION(:) :: PERIM_ELEM
    LOGICAL, ALLOCATABLE, DIMENSION(:) :: ISPERIM
    INTEGER, POINTER, DIMENSION(:,:) :: TMP_IPAR
    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: I2
    REAL(PDB), ALLOCATABLE, DIMENSION(:) :: ANGLE,ALBEDO
    REAL(PDB), POINTER, DIMENSION(:,:) :: TMP_RPAR
    REAL(PDB), ALLOCATABLE, DIMENSION(:,:,:) :: ALIGN
    TYPE(T_SALBCDATA), POINTER, DIMENSION(:) :: TMP_BCDATAREAD
    !
    ! compute size of the unfold geometry
    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
    DO ELEM=1,GG%NB_ELEM
      TYPE=GG%IPAR(1,ELEM)
      IF(TYPE==1) THEN
        X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
        XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
        X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
        XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2);
      ENDIF
    ENDDO
    LENGTHX=XMAX-XMIN; LENGTHY=YMAX-YMIN;
    !
    NSYM=1
    IF(HSYM=='SYMX') THEN ! define symmetry axis
      AXIS_X1(1)=0._PDB; AXIS_X2(1)=100._PDB; AXIS_Y1(1)=0._PDB; AXIS_Y2(1)=0._PDB;
    ELSE IF(HSYM=='SYMY') THEN
      AXIS_X1(1)=0._PDB; AXIS_X2(1)=0._PDB; AXIS_Y1(1)=0._PDB; AXIS_Y2(1)=100._PDB;
    ELSE IF(HSYM=='DIAG') THEN
      AXIS_X1(1)=0._PDB; AXIS_X2(1)=100._PDB; AXIS_Y1(1)=0._PDB; AXIS_Y2(1)=100._PDB;
    ELSE IF(HSYM=='SA60') THEN
      ! the hexagon side is on south
      NSYM=2
      AXIS_X1(1)=XMIN; AXIS_Y1(1)=YMIN; AXIS_X2(1)=XMIN+0.5_PDB*LENGTHX
      AXIS_Y2(1)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX
      AXIS_X1(2)=XMIN+LENGTHX; AXIS_Y1(2)=YMIN; AXIS_X2(2)=XMIN+0.5_PDB*LENGTHX
      AXIS_Y2(2)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX
    ELSE IF(HSYM=='SB60') THEN
      ! the hexagon side is on north-east
      NSYM=2
      AXIS_X1(1)=XMIN; AXIS_Y1(1)=YMIN; AXIS_X2(1)=XMIN+0.5_PDB*LENGTHX
      AXIS_Y2(1)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX
      AXIS_X1(2)=XMIN; AXIS_Y1(2)=YMIN; AXIS_X2(2)=XMIN
      AXIS_Y2(2)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX
    ELSE IF(HSYM=='S30') THEN
      AXIS_X1(1)=XMIN; AXIS_Y1(1)=YMIN; AXIS_X2(1)=XMIN+0.75_PDB*LENGTHX
      AXIS_Y2(1)=YMIN+0.25_PDB*SQRT(3._PDB)*LENGTHX
    ELSE IF(HSYM=='SYMH') THEN
      AXIS_X1(1)=XMIN ; AXIS_Y1(1)=YMIN+0.25_PDB*SQRT(3._PDB)*LENGTHX
      AXIS_X2(1)=XMIN+LENGTHX ; AXIS_Y2(1)=AXIS_Y1(1)
    ELSE
      CALL XABORT('SALFOLD_1: invalid type of symmetry axis')
    ENDIF
    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
    !
    ! allocate new surfacic element containers
    ALLOCATE(TMP_IPAR(NIPAR,3*GG%NB_ELEM), TMP_RPAR(NRPAR,3*GG%NB_ELEM), &
             I2(3,GG%NB_ELEM), STAT=OK)
    IF(OK/=0) CALL XABORT('SALFOLD_1: not enough memory')
    TMP_IPAR(:,:)=0; TMP_RPAR(:,:)=0._PDB;
    !
    ! loop over old elements
    TMP_NB_ELEM=0
    THETA1=0._PDB; THETA2=0._PDB;
    DO ELEM=1,GG%NB_ELEM
      I2(:,ELEM)=0
      TYPE=GG%IPAR(1,ELEM)
      NOCOPY(:2)=.FALSE.
      DO ISYM=1,NSYM
        X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); RAD=GG%RPAR(3,ELEM)
        IF(TYPE==1) THEN
          XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
          X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
          CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1,Y1,X4,Y4)
          CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X2,Y2,DX4,DY4)
          NOCOPY(ISYM)=ABS(X1-X4)<10.0*EPS .AND. ABS(Y1-Y4)<10.0*EPS .AND. &
            ABS(X2-DX4)<10.0*EPS .AND. ABS(Y2-DY4)<10.0*EPS
          IF(NOCOPY(ISYM)) CYCLE
          IF(HSYM=='SB60') THEN
            IF(ISYM==1) THEN
              IF((ABS(Y1)<10.0*EPS).AND.(ABS(Y2)<10.0*EPS)) CYCLE
            ELSE
              IF((ABS(2._PDB*ABS(Y1-YMIN)-ABS(X4-X1)*SQRT(3._PDB))<10.0*EPS .AND. &
                & ABS(Y1-Y4)<10.0*EPS .AND. ABS(2._PDB*ABS(Y2-YMIN)-ABS(DX4-X2)*SQRT(3._PDB))<10.0*EPS &
                & .AND. ABS(Y2-DY4)<10.0*EPS)) CYCLE
              ALIGN3(1,2)=YMIN ; ALIGN3(2,2)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX
              ALIGN3(1,1)=-XMIN
              ALIGN3(2,1)=-XMIN-0.5_PDB*LENGTHX
              ALIGN3(:3,3)=1.0_PDB
              ALIGN3(3,1)=X4; ALIGN3(3,2)=Y4;
              DET1 = DET_ROSETTA(ALIGN3(1,1),3)
              ALIGN3(3,1)=DX4; ALIGN3(3,2)=DY4;
              DET2 = DET_ROSETTA(ALIGN3(1,1),3)
              IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) CYCLE
            ENDIF
          ENDIF
        ELSE IF(TYPE==2) THEN
          CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1,Y1,X4,Y4)
          THETA1=0._PDB; THETA2=0._PDB;
        ELSE IF(TYPE==3) THEN
          CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1,Y1,X4,Y4)
          X1B=X1+RAD*COS(GG%RPAR(4,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(4,ELEM));
          CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1B,Y1B,X4B,Y4B)
          IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN
            THETA1=PI/2._PDB
          ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN
            THETA1=3._PDB*PI/2._PDB
          ELSE IF(X4B-X4 > 0._PDB) THEN
            THETA1=ATAN((Y4B-Y4)/(X4B-X4))
          ELSE
            THETA1=ATAN((Y4B-Y4)/(X4B-X4))+PI
          ENDIF
          X1B=X1+RAD*COS(GG%RPAR(5,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(5,ELEM));
          CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1B,Y1B,X4B,Y4B)
          IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN
            THETA2=PI/2._PDB
          ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN
            THETA2=3._PDB*PI/2._PDB
          ELSE IF(X4B-X4 > 0._PDB) THEN
            THETA2=ATAN((Y4B-Y4)/(X4B-X4))
          ELSE
            THETA2=ATAN((Y4B-Y4)/(X4B-X4))+PI
          ENDIF
        ELSE
          WRITE(*,*) " elem=",ELEM," type=",TYPE," isym=",ISYM
          CALL XABORT('SALFOLD_1: invalid type of surfacic element')
        ENDIF
        TMP_NB_ELEM=TMP_NB_ELEM+1
        IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_1: tmp_nb_elem overflow(1)')
        I2(ISYM+1,ELEM)=TMP_NB_ELEM
        TMP_RPAR(1,TMP_NB_ELEM)=X4; TMP_RPAR(2,TMP_NB_ELEM)=Y4;
        TMP_IPAR(1,TMP_NB_ELEM)=TYPE;
        IF(TYPE==1) THEN
          TMP_RPAR(3,TMP_NB_ELEM)=DX4-X4; TMP_RPAR(4,TMP_NB_ELEM)=DY4-Y4;
          TMP_RPAR(5,TMP_NB_ELEM)=SQRT(TMP_RPAR(3,TMP_NB_ELEM)**2+TMP_RPAR(4,TMP_NB_ELEM)**2)
          XMIN=MIN(XMIN,X4); YMIN=MIN(YMIN,Y4); XMAX=MAX(XMAX,X4); YMAX=MAX(YMAX,Y4);
          XMIN=MIN(XMIN,DX4); YMIN=MIN(YMIN,DY4); XMAX=MAX(XMAX,DX4); YMAX=MAX(YMAX,DY4);
          TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(3,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(2,ELEM);
        ELSE IF((TYPE==2).OR.(TYPE==3)) THEN
          TMP_RPAR(3,TMP_NB_ELEM)=GG%RPAR(3,ELEM) ! RADIUS
          IF(THETA2>THETA1) THETA1=THETA1+2._PDB*PI
          PHI1=THETA2; DELPHI=THETA1-THETA2;
          IF(DELPHI>0._PDB)THEN
            PHI2=PHI1+DELPHI
          ELSE
            PHI2=PHI1
            PHI1=PHI1+DELPHI
          ENDIF
          IF(TYPE==3)THEN
            ! arc of circle: put phi1 within 0 and 2*pi
            IF(PHI1>2._PDB*PI)THEN
              IAUX=INT(PHI1/(2._PDB*PI))
              DELPHI=(2._PDB*PI)*IAUX
              PHI1=PHI1-DELPHI ; PHI2=PHI2-DELPHI
            ELSEIF(PHI1<0._PDB)THEN
              IAUX=INT((-PHI1+1.D-7)/(2._PDB*PI))+1
              DELPHI=(2._PDB*PI)*IAUX
              PHI1=PHI1+DELPHI ; PHI2=PHI2+DELPHI
            ENDIF
          ENDIF
          TMP_RPAR(4,TMP_NB_ELEM)=PHI1; TMP_RPAR(5,TMP_NB_ELEM)=PHI2; ! ANGLES
          TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(2,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(3,ELEM)
        ENDIF
        TMP_RPAR(6,TMP_NB_ELEM)=0._PDB
      ENDDO
      IF((.NOT.NOCOPY(1)).AND.(.NOT.NOCOPY(2))) THEN
        TMP_NB_ELEM=TMP_NB_ELEM+1
        IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_1: tmp_nb_elem overflow(2)')
        TMP_IPAR(:,TMP_NB_ELEM)=GG%IPAR(:,ELEM)
        TMP_RPAR(:,TMP_NB_ELEM)=GG%RPAR(:,ELEM)
        I2(1,ELEM)=TMP_NB_ELEM
      ENDIF
    ENDDO
    DEALLOCATE(GG%IPAR,GG%RPAR)
    GG%IPAR=>TMP_IPAR; GG%RPAR=>TMP_RPAR;
    GG%NB_ELEM=TMP_NB_ELEM
    !
    ! translate the domain
    DO ELEM=1,GG%NB_ELEM
      GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-XMIN
      GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-YMIN
    ENDDO
    LENGTHX=XMAX-XMIN ; LENGTHY=YMAX-YMIN ;
    !
    ! loop over boundary conditions
    TMP_NBBCDA=0
    ALLOCATE(ISPERIM(GG%NB_ELEM),ALIGN(3,3,GG%NB_ELEM),ANGLE(GG%NB_ELEM), &
    & ALBEDO(GG%NB_ELEM),PERIM_ELEM(GG%NB_ELEM))
    ALIGN(:3,3,:GG%NB_ELEM)=1.0_PDB
    PERIM_ELEM(:GG%NB_ELEM)=0
    ISPERIM(:GG%NB_ELEM)=.FALSE.
    DO IBC=1,GG%NBBCDA
      DO I=1,GG%BCDATAREAD(IBC)%NBER
        INDBC=GG%BCDATAREAD(IBC)%ELEMNB(I)
        IF(INDBC==0) CYCLE
        DO ISYM=1,NSYM+1
          IF(I2(ISYM,INDBC)/=0) ISPERIM(I2(ISYM,INDBC))=.TRUE.
        ENDDO
      ENDDO
    ENDDO
    ITER0: DO ELEM=1,GG%NB_ELEM
      IF(.NOT.ISPERIM(ELEM)) CYCLE
      X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
      X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
      DO J=1,TMP_NBBCDA
        ALIGN(3,1,J)=X1; ALIGN(3,2,J)=Y1;
        DET1 = DET_ROSETTA(ALIGN(1,1,J),3)
        ALIGN(3,1,J)=X2; ALIGN(3,2,J)=Y2;
        DET2 = DET_ROSETTA(ALIGN(1,1,J),3)
        IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) THEN
          PERIM_ELEM(ELEM) = J
          CYCLE ITER0
        ENDIF
      ENDDO
      TMP_NBBCDA=TMP_NBBCDA+1
      PERIM_ELEM(ELEM) = TMP_NBBCDA
      ANGLE(TMP_NBBCDA)=ATAN((Y2-Y1)/(X2-X1))
      IF(ABS(ANGLE(TMP_NBBCDA)).LE.1.0E-5) ANGLE(TMP_NBBCDA)=0.0
      ALIGN(1,1,TMP_NBBCDA)=X1; ALIGN(1,2,TMP_NBBCDA)=Y1
      ALIGN(2,1,TMP_NBBCDA)=X2; ALIGN(2,2,TMP_NBBCDA)=Y2
      ! Recover albedo from folded geometry
      ALBEDO(TMP_NBBCDA)=1.0
      DO IBC=1,GG%NBBCDA
        J = FINDLC(GG%BCDATAREAD(IBC)%ELEMNB,ELEM)
        IF(J.EQ.1) THEN
          ALBEDO(TMP_NBBCDA)=GG%BCDATAREAD(IBC)%BCDATA(6)
          EXIT
        ENDIF
      ENDDO
    ENDDO ITER0
    ALLOCATE(TMP_BCDATAREAD(TMP_NBBCDA))
    DO IBC=1,TMP_NBBCDA
      TMP_BCDATAREAD(IBC)%NBER = COUNT(PERIM_ELEM(:GG%NB_ELEM) == IBC)
      ALLOCATE(TMP_BCDATAREAD(IBC)%ELEMNB(TMP_BCDATAREAD(IBC)%NBER))
      TMP_BCDATAREAD(IBC)%SALTYPE = 0
      J=0
      DO I=1,GG%NB_ELEM
        IF(PERIM_ELEM(I) == IBC) THEN
          J=J+1
          TMP_BCDATAREAD(IBC)%ELEMNB(J) = I
        ENDIF
      ENDDO
      TMP_BCDATAREAD(IBC)%BCDATA(1) = ALIGN(1,1,IBC)
      TMP_BCDATAREAD(IBC)%BCDATA(2) = ALIGN(1,2,IBC)
      TMP_BCDATAREAD(IBC)%BCDATA(3) = COS(ANGLE(IBC))
      TMP_BCDATAREAD(IBC)%BCDATA(4) = SIN(ANGLE(IBC))
      TMP_BCDATAREAD(IBC)%BCDATA(5) = ANGLE(IBC)
      TMP_BCDATAREAD(IBC)%BCDATA(6) = ALBEDO(IBC)
    ENDDO
    DEALLOCATE(I2,PERIM_ELEM,ALBEDO,ANGLE,ALIGN,ISPERIM)
    DEALLOCATE(GG%BCDATAREAD)
    GG%BCDATAREAD=>TMP_BCDATAREAD
    GG%NBBCDA=TMP_NBBCDA
    GG%ALBEDO=1.D0
  END SUBROUTINE SALFOLD_1
  !
  SUBROUTINE SALFOLD_2(HSYM,GG)
    !
    !---------------------------------------------------------------------
    !
    !Purpose:
    ! unfold the domain with rotation
    !
    !Parameters: input
    ! HSYM: type of symmetry: SR60': dual 60-degree rotation; R180: 180-degree
    !       rotation; R120: dual 120-degree rotation
    !
    !Parameters: input/output
    ! GG    geometry descriptor
    !
    !---------------------------------------------------------------------
    !
    USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,LENGTHX,LENGTHY
    IMPLICIT NONE
    CHARACTER(LEN=4),INTENT(IN)  :: HSYM
    TYPE(T_G_BASIC), INTENT(INOUT) :: GG
    !
    INTEGER :: ELEM,ELEM2,TYPE,OK,TMP_NB_ELEM,ISYM,NSYM,IAUX
    REAL(PDB) :: THROT(3),X1,X2,X4,Y1,Y2,Y4,DX4,DY4,RAD,THETA1,THETA2,X1B,Y1B,X4B, &
                 Y4B,XMIN,YMIN,XMAX,YMAX,CENTER_X,CENTER_Y,DELPHI
    REAL(PDB),PARAMETER :: EPS=1.0E-5_PDB
    !
    ! allocatable arrays
    LOGICAL, ALLOCATABLE, DIMENSION(:) :: ELEM_DUP
    INTEGER, POINTER, DIMENSION(:,:) :: TMP_IPAR
    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: I2
    REAL(PDB), POINTER, DIMENSION(:,:) :: TMP_RPAR
    !
    ! compute size of the unfold geometry
    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
    DO ELEM=1,GG%NB_ELEM
      TYPE=GG%IPAR(1,ELEM)
      IF(TYPE==1) THEN
        X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
        XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
        X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
        XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2);
      ENDIF
    ENDDO
    LENGTHX=XMAX-XMIN; LENGTHY=YMAX-YMIN;
    !
    NSYM=0
    IF(HSYM=='SR60') THEN ! define rotation center
      NSYM=3
      CENTER_X=XMIN+0.5_PDB*LENGTHX; CENTER_Y=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX
      THROT(1)=0.0_PDB; THROT(2)=-PI/3.0_PDB; THROT(3)=PI/3.0_PDB
    ELSE IF(HSYM=='R180') THEN
      NSYM=2
      CENTER_X=XMIN+0.5_PDB*LENGTHX; CENTER_Y=YMIN+0.25_PDB*SQRT(3._PDB)*LENGTHX
      THROT(1)=0.0_PDB; THROT(2)=PI; 
    ELSE IF(HSYM=='R120') THEN
      NSYM=3
      CENTER_X=XMIN+0.5_PDB*LENGTHX/1.5_PDB; CENTER_Y=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX/1.5_PDB
      THROT(1)=0.0_PDB; THROT(2)=-2.0_PDB*PI/3.0_PDB; THROT(3)=2.0_PDB*PI/3.0_PDB
    ELSE
      CALL XABORT('SALFOLD_2: invalid type of symmetry axis')
    ENDIF
    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
    !
    ! allocate new surfacic element containers
    ALLOCATE(I2(NSYM,GG%NB_ELEM))
    ALLOCATE(TMP_IPAR(NIPAR,3*GG%NB_ELEM), TMP_RPAR(NRPAR,3*GG%NB_ELEM), STAT=OK)
    IF(OK/=0) CALL XABORT('SALFOLD_2: not enough memory')
    I2(:NSYM,:GG%NB_ELEM)=0
    TMP_IPAR(:,:)=0; TMP_RPAR(:,:)=0._PDB;
    !
    ! loop over old elements
    TMP_NB_ELEM=0
    THETA1=0._PDB; THETA2=0._PDB;
    DO ELEM=1,GG%NB_ELEM
      TYPE=GG%IPAR(1,ELEM)
      DO ISYM=1,NSYM
        X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); RAD=GG%RPAR(3,ELEM)
        IF(TYPE==1) THEN
          XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
          X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
          CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1,Y1,X4,Y4)
          CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X2,Y2,DX4,DY4)
        ELSE IF(TYPE==2) THEN
          CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1,Y1,X4,Y4)
          THETA1=0._PDB; THETA2=0._PDB;
        ELSE IF(TYPE==3) THEN
          DELPHI=GG%RPAR(5,ELEM)-GG%RPAR(4,ELEM)
          X1B=X1+RAD*COS(GG%RPAR(4,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(4,ELEM));
          CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1,Y1,X4,Y4)
          CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1B,Y1B,X4B,Y4B)
          IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN
            THETA1=PI/2._PDB
          ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN
            THETA1=3._PDB*PI/2._PDB
          ELSE IF(X4B-X4 > 0._PDB) THEN
            THETA1=ATAN((Y4B-Y4)/(X4B-X4))
          ELSE
            THETA1=ATAN((Y4B-Y4)/(X4B-X4))+PI
          ENDIF
          THETA2=THETA1+DELPHI
          ! put THETA1 within 0 and 2*pi
          IF(THETA1>2._PDB*PI)THEN
            IAUX=INT(THETA1/(2._PDB*PI))
            DELPHI=(2._PDB*PI)*IAUX
            THETA1=THETA1-DELPHI ; THETA2=THETA2-DELPHI
          ELSEIF(THETA1<0._PDB)THEN
            IAUX=INT((-THETA1+1.D-7)/(2._PDB*PI))+1
            DELPHI=(2._PDB*PI)*IAUX
            THETA1=THETA1+DELPHI ; THETA2=THETA2+DELPHI
          ENDIF
        ELSE
          WRITE(*,*) " elem=",ELEM," type=",TYPE," isym=",ISYM
          CALL XABORT('SALFOLD_2: invalid type of surfacic element')
        ENDIF
        TMP_NB_ELEM=TMP_NB_ELEM+1
        IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_2: TMP_NB_ELEM overflow')
        I2(ISYM,ELEM)=TMP_NB_ELEM
        TMP_RPAR(1,TMP_NB_ELEM)=X4; TMP_RPAR(2,TMP_NB_ELEM)=Y4;
        TMP_IPAR(1,TMP_NB_ELEM)=TYPE;
        IF(TYPE==1) THEN
          TMP_RPAR(3,TMP_NB_ELEM)=DX4-X4; TMP_RPAR(4,TMP_NB_ELEM)=DY4-Y4;
          TMP_RPAR(5,TMP_NB_ELEM)=SQRT(TMP_RPAR(3,TMP_NB_ELEM)**2+TMP_RPAR(4,TMP_NB_ELEM)**2)
          XMIN=MIN(XMIN,X4); YMIN=MIN(YMIN,Y4); XMAX=MAX(XMAX,X4); YMAX=MAX(YMAX,Y4);
          XMIN=MIN(XMIN,DX4); YMIN=MIN(YMIN,DY4); XMAX=MAX(XMAX,DX4); YMAX=MAX(YMAX,DY4);
          TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(2,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(3,ELEM)
        ELSE IF((TYPE==2).OR.(TYPE==3)) THEN
          TMP_RPAR(3,TMP_NB_ELEM)=GG%RPAR(3,ELEM) ! RADIUS
          TMP_RPAR(4,TMP_NB_ELEM)=THETA1; TMP_RPAR(5,TMP_NB_ELEM)=THETA2; ! ANGLES
          TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(2,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(3,ELEM)
        ENDIF
        TMP_RPAR(6,TMP_NB_ELEM)=0._PDB
      ENDDO
    ENDDO
    DEALLOCATE(GG%IPAR,GG%RPAR)
    !
    ! eliminate duplicate elements
    ALLOCATE(ELEM_DUP(TMP_NB_ELEM))
    ELEM_DUP(:TMP_NB_ELEM)=.FALSE.
    DO ELEM=1,TMP_NB_ELEM
      TYPE=TMP_IPAR(1,ELEM)
      IF((TYPE==1).AND.((TMP_IPAR(2,ELEM)==0).OR.(TMP_IPAR(3,ELEM)==0))) THEN
        X4=TMP_RPAR(1,ELEM); Y4=TMP_RPAR(2,ELEM)
        DX4=X4+TMP_RPAR(3,ELEM); DY4=Y4+TMP_RPAR(4,ELEM)
        DO ELEM2=ELEM+1,TMP_NB_ELEM
          IF(TMP_IPAR(1,ELEM2)/=1) CYCLE
          X1=TMP_RPAR(1,ELEM2); Y1=TMP_RPAR(2,ELEM2)
          X2=X1+TMP_RPAR(3,ELEM2); Y2=Y1+TMP_RPAR(4,ELEM2)
          IF(((ABS(X1-X4)<10.0*EPS  .AND. ABS(Y1-Y4)<10.0*EPS  .AND. &
               ABS(X2-DX4)<10.0*EPS .AND. ABS(Y2-DY4)<10.0*EPS)).OR. &
             ((ABS(X1-DX4)<10.0*EPS .AND. ABS(Y1-DY4)<10.0*EPS .AND. &
               ABS(X2-X4)<10.0*EPS  .AND. ABS(Y2-Y4)<10.0*EPS))) THEN
            ELEM_DUP(ELEM2)=.TRUE.
            IF(TMP_IPAR(2,ELEM)==0) TMP_IPAR(2,ELEM)=MAX(TMP_IPAR(2,ELEM2),TMP_IPAR(3,ELEM2))
            IF(TMP_IPAR(3,ELEM)==0) TMP_IPAR(3,ELEM)=MAX(TMP_IPAR(2,ELEM2),TMP_IPAR(3,ELEM2))
            IF(TMP_IPAR(2,ELEM)==TMP_IPAR(3,ELEM)) ELEM_DUP(ELEM)=.TRUE.
            EXIT
          ENDIF
        ENDDO
      ENDIF
    ENDDO
    !
    ELEM=1
    DO WHILE(ELEM<=TMP_NB_ELEM)
      IF(ELEM_DUP(ELEM)) THEN
        TMP_NB_ELEM=TMP_NB_ELEM-1
        DO ELEM2=ELEM,TMP_NB_ELEM
          TMP_IPAR(:,ELEM2)=TMP_IPAR(:,ELEM2+1); TMP_RPAR(:,ELEM2)=TMP_RPAR(:,ELEM2+1)
          ELEM_DUP(ELEM2)=ELEM_DUP(ELEM2+1)
        ENDDO
      ELSE
        ELEM=ELEM+1
      ENDIF
    ENDDO
    GG%IPAR=>TMP_IPAR; GG%RPAR=>TMP_RPAR;
    GG%NB_ELEM=TMP_NB_ELEM
    !
    ! translate the domain
    DO ELEM=1,GG%NB_ELEM
      GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-XMIN
      GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-YMIN
    ENDDO
    LENGTHX=XMAX-XMIN ; LENGTHY=YMAX-YMIN ;
    GG%NBBCDA=0; GG%ALBEDO=1.D0
    !
    DEALLOCATE(ELEM_DUP,I2)
  END SUBROUTINE SALFOLD_2
END MODULE SAL_GEOMETRY_MOD