Initial commit from Polytechnique Montreal

1 files changed, 286 insertions, 0 deletions

diff --git a/Dragon/src/SYBILT.f b/Dragon/src/SYBILT.f
new file mode 100644
index 0000000..d1db3f7
--- /dev/null
+++ b/Dragon/src/SYBILT.f
@@ -0,0 +1,286 @@
+*DECK SYBILT
+      SUBROUTINE SYBILT(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Sybil-type (2D) tracking operator.
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+*Author(s): A. Hebert
+*
+*Parameters: input/output
+* NENTRY  number of LCM objects or files used by the operator.
+* HENTRY  name of each LCM object or file:
+*         HENTRY(1) creation or modification type(L_TRACK);
+*         HENTRY(2) read-only type(L_GEOM).
+* IENTRY  type of each LCM object or file:
+*         =1 LCM memory object; =2 XSM file; =3 sequential binary file;
+*         =4 sequential ascii file.
+* JENTRY  access of each LCM object or file:
+*         =0 the LCM object or file is created;
+*         =1 the LCM object or file is open for modifications;
+*         =2 the LCM object or file is open in read-only mode.
+* KENTRY  LCM object address or file unit number.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER      NENTRY,IENTRY(NENTRY),JENTRY(NENTRY)
+      TYPE(C_PTR)  KENTRY(NENTRY)
+      CHARACTER    HENTRY(NENTRY)*12
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (NSTATE=40,IOUT=6)
+      CHARACTER TEXT4*4,TEXT12*12,TITLE*72,HSIGN*12
+      DOUBLE PRECISION DFLOTT
+      INTEGER ITITL(18),IGP(NSTATE),ISTATE(NSTATE),JQUA2(2),IPARAM(15)
+*----
+*  PARAMETER VALIDATION
+*----
+      IF(NENTRY.LE.1) CALL XABORT('SYBILT: TWO PARAMETERS EXPECTED.')
+      IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2)) CALL XABORT('SYBILT: L'
+     1 //'CM OBJECT EXPECTED AT LHS.')
+      IF((JENTRY(1).NE.0).AND.(JENTRY(1).NE.1)) CALL XABORT('SYBILT: E'
+     1 //'NTRY IN CREATE OR MODIFICATION MODE EXPECTED.')
+      IF((JENTRY(2).NE.2).OR.((IENTRY(2).NE.1).AND.(IENTRY(2).NE.2)))
+     1 CALL XABORT('SYBILT: LCM OBJECT IN READ-ONLY MODE EXPECTED AT R'
+     1 //'HS.')
+      CALL LCMGTC(KENTRY(2),'SIGNATURE',12,HSIGN)
+      IF(HSIGN.NE.'L_GEOM') THEN
+         TEXT12=HENTRY(2)
+         CALL XABORT('SYBILT: SIGNATURE OF '//TEXT12//' IS '//HSIGN//
+     1   '. L_GEOM EXPECTED.')
+      ENDIF
+      HSIGN='L_TRACK'
+      CALL LCMPTC(KENTRY(1),'SIGNATURE',12,HSIGN)
+      HSIGN='SYBIL'
+      CALL LCMPTC(KENTRY(1),'TRACK-TYPE',12,HSIGN)
+      CALL LCMGET(KENTRY(2),'STATE-VECTOR',ISTATE)
+*
+      IMPX=1
+      TITLE=' '
+      IF(JENTRY(1).EQ.0) THEN
+         MAXPTS=ISTATE(6)
+         MAXZ=10000
+         MAXJ=MAX(18,4*MAXPTS)
+         MULTC=3
+         IWIGN=3
+         INORM=0
+         IQW=0
+         JQUA1=5
+         JQUA2(1)=3
+         JQUA2(2)=3
+         IQUA10=5
+         IBIHET=2
+      ELSE IF(JENTRY(1).EQ.1) THEN
+         CALL LCMGTC(KENTRY(1),'SIGNATURE',12,HSIGN)
+         IF(HSIGN.NE.'L_TRACK') THEN
+            TEXT12=HENTRY(1)
+            CALL XABORT('SYBILT: SIGNATURE OF '//TEXT12//' IS '//HSIGN//
+     1      '. L_TRACK EXPECTED.')
+         ENDIF
+         CALL LCMGTC(KENTRY(1),'TRACK-TYPE',12,HSIGN)
+         IF(HSIGN.NE.'SYBIL') THEN
+            TEXT12=HENTRY(1)
+            CALL XABORT('SYBILT: TRACK-TYPE OF '//TEXT12//' IS '//HSIGN
+     1      //'. SYBIL EXPECTED.')
+         ENDIF
+         CALL LCMGET(KENTRY(1),'STATE-VECTOR',IGP)
+         MAXPTS=IGP(1)
+         ITG=IGP(6)
+         MAXZ=IGP(7)
+         MAXJ=IGP(8)
+         IF(ITG.EQ.2) THEN
+            CALL LCMSIX(KENTRY(1),'PURE-GEOM',1)
+            CALL LCMGET(KENTRY(1),'PARAM',IPARAM)
+            CALL LCMSIX(KENTRY(1),' ',2)
+            JQUA1=IPARAM(3)
+         ELSE IF(ITG.EQ.3) THEN
+            CALL LCMSIX(KENTRY(1),'DOITYOURSELF',1)
+            CALL LCMGET(KENTRY(1),'PARAM',IPARAM)
+            CALL LCMSIX(KENTRY(1),' ',2)
+            JQUA1=IPARAM(2)
+         ELSE IF(ITG.EQ.4) THEN
+            CALL LCMSIX(KENTRY(1),'EURYDICE',1)
+            CALL LCMGET(KENTRY(1),'PARAM',IPARAM)
+            CALL LCMSIX(KENTRY(1),' ',2)
+            MULTC=IPARAM(2)
+            IWIGN=IPARAM(3)
+            INORM=IPARAM(12)
+            IQW=IPARAM(13)
+            JQUA1=IPARAM(10)
+            JQUA2(1)=IPARAM(8)
+            JQUA2(2)=IPARAM(9)
+         ENDIF
+         IF(IGP(40).NE.0) THEN
+            CALL LCMSIX(KENTRY(1),'BIHET',1)
+            CALL LCMGET(KENTRY(1),'PARAM',IPARAM)
+            CALL LCMSIX(KENTRY(1),' ',2)
+            IBIHET=IPARAM(6)
+            IQUA10=IPARAM(8)
+         ELSE
+            IBIHET=0
+            IQUA10=0
+         ENDIF
+         CALL LCMLEN(KENTRY(1),'TITLE',LENGT,ITYLCM)
+         IF(LENGT.GT.0) CALL LCMGTC(KENTRY(1),'TITLE',72,TITLE)
+      ENDIF
+      IRECT=0
+      IHALT=0
+      ILIGN=0
+      FRTM=0.05
+      EPSJ=0.5E-5
+   20 CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT)
+   21 CONTINUE
+      IF(INDIC.EQ.10) GO TO 30
+      IF(INDIC.NE.3) CALL XABORT('SYBILT: CHARACTER DATA EXPECTED.')
+      IF(TEXT4.EQ.'EDIT') THEN
+         CALL REDGET(INDIC,IMPX,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED(1).')
+      ELSE IF(TEXT4.EQ.'MAXR') THEN
+         CALL REDGET(INDIC,MAXPTS,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED(2).')
+         MAXJ=MAX(MAXJ,4*MAXPTS)
+      ELSE IF(TEXT4.EQ.'MAXJ') THEN
+*        ALLOCATED STORAGE FOR STORING THE INTERFACE CURRENTS.
+         CALL REDGET(INDIC,MAXJ,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED.')
+      ELSE IF(TEXT4.EQ.'MAXZ') THEN
+*        ALLOCATED STORAGE FOR STORING THE TRACKING INFORMATION.
+         CALL REDGET(INDIC,MAXZ,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED.')
+      ELSE IF(TEXT4.EQ.'TITL') THEN
+         CALL REDGET(INDIC,NITMA,FLOTT,TITLE,DFLOTT)
+         IF(INDIC.NE.3) CALL XABORT('SYBILT: TITLE EXPECTED.')
+      ELSE IF(TEXT4.EQ.'ROTH') THEN
+         MULTC=1
+      ELSE IF(TEXT4.EQ.'ROT+') THEN
+         MULTC=2
+      ELSE IF(TEXT4.EQ.'DP00') THEN
+         MULTC=3
+      ELSE IF(TEXT4.EQ.'DP01') THEN
+         MULTC=4
+      ELSE IF(TEXT4.EQ.'ASKE') THEN
+*        USE ASKEW CYLINDERIZATION IN EURYDICE.
+         IWIGN=1
+      ELSE IF(TEXT4.EQ.'WIGN') THEN
+*        USE ASKEW CYLINDERIZATION IN EURYDICE.
+         IWIGN=2
+      ELSE IF(TEXT4.EQ.'SANC') THEN
+*        USE SANCHEZ CYLINDERIZATION IN EURYDICE.
+         IWIGN=3
+      ELSE IF(TEXT4.EQ.'HALT') THEN
+*        STOP AT THE END OF TRACKING.
+         IHALT=1
+      ELSE IF(TEXT4.EQ.'LIGN') THEN
+*        PRINT THE TRACKING INFORMATION.
+         ILIGN=1
+      ELSE IF(TEXT4.EQ.'RENO') THEN
+*        NORMALIZE TRACKS.
+         INORM=0
+      ELSE IF(TEXT4.EQ.'NORE') THEN
+*        DO NOT NORMALIZE TRACKS.
+         INORM=1
+      ELSE IF(TEXT4.EQ.'RECT') THEN
+*        DO NOT CONSIDER THE SYMMETRIES OF SQUARE CELLS IN EURYDICE.
+         IRECT=1
+      ELSE IF(TEXT4.EQ.'GAUS') THEN
+*        USE GAUSS QUADRATURES IN ANGLE AND SPACE.
+         IQW=0
+      ELSE IF(TEXT4.EQ.'EQW') THEN
+*        USE EQUAL WEIGHT QUADRATURES IN ANGLE AND SPACE.
+         IQW=1
+      ELSE IF(TEXT4.EQ.'QUA1') THEN
+*        1-D QUADRATURE PARAMETER.
+         CALL REDGET(INDIC,JQUA1,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED.')
+      ELSE IF(TEXT4.EQ.'QUA2') THEN
+*        2-D QUADRATURE PARAMETERS (ANGLE AND SPACE).
+         CALL REDGET(INDIC,JQUA2(1),FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED.')
+         CALL REDGET(INDIC,JQUA2(2),FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED.')
+      ELSE IF(TEXT4.EQ.'EPSJ') THEN
+         CALL REDGET(INDIC,NITMA,EPSJ,TEXT4,DFLOTT)
+         IF (INDIC.NE.2) CALL XABORT('SYBILT: REAL DATA EXPECTED.')
+      ELSE IF(TEXT4.EQ.'QUAB') THEN
+         CALL REDGET(INDIC,IQUA10,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('SYBILT: INTEGER DATA EXPECTED.')
+      ELSE IF(TEXT4.EQ.'SAPO') THEN
+         IBIHET=1
+      ELSE IF(TEXT4.EQ.'HEBE') THEN
+         IBIHET=2
+      ELSE IF(TEXT4.EQ.'SLSI') THEN
+         IBIHET=3
+         CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT)
+         IF (INDIC.NE.2) GOTO 21
+         FRTM=FLOTT
+      ELSE IF(TEXT4.EQ.'SLSS') THEN
+         IBIHET=4
+         CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT)
+         IF (INDIC.NE.2) GOTO 21
+         FRTM=FLOTT
+      ELSE IF(TEXT4.EQ.';') THEN
+         GO TO 30
+      ELSE
+         CALL XABORT('SYBILT: '//TEXT4//' IS AN INVALID KEY WORD.')
+      ENDIF
+      GO TO 20
+*
+   30 IF(TITLE.NE.' ') THEN
+         READ(TITLE,'(18A4)') (ITITL(I),I=1,18)
+         CALL LCMPUT(KENTRY(1),'TITLE',18,3,ITITL)
+      ENDIF
+      TEXT12=HENTRY(2)
+      READ(TEXT12,'(3A4)') (ITITL(I),I=1,3)
+      CALL LCMPUT(KENTRY(1),'LINK.GEOM',3,3,ITITL)
+      IF(IMPX.GT.1) WRITE(IOUT,100) TITLE
+*
+      IF(MAXPTS.EQ.0) CALL XABORT('SYBILT: MAXPTS NOT DEFINED.')
+      CALL LCMPUT(KENTRY(1),'EPSJ',1,2,EPSJ)
+      CALL SYBTRK (KENTRY(1),KENTRY(2),IMPX,MAXPTS,MAXJ,MAXZ,MULTC,
+     1 IWIGN,IHALT,ILIGN,INORM,IRECT,IQW,JQUA1,JQUA2,IQUA10,IBIHET,FRTM)
+*
+      IF(IMPX.GT.0) THEN
+         CALL LCMGET(KENTRY(1),'STATE-VECTOR',IGP)
+         WRITE(IOUT,110) (IGP(I),I=1,9),IGP(40)
+         WRITE(IOUT,120) EPSJ
+      ENDIF
+      RETURN
+*
+  100 FORMAT(1H1,
+     1 36H SSSSS  YY  YY  BBBBBB  IIIIII LL   ,10X,2H1 ,83(1H*)/
+     2 38H SSSSSSS YY  YY  BBBBBBB IIIIII LL    ,8X,3H11 ,45(1H*),
+     3 38H MULTIGROUP VERSION.  A. HEBERT (1987)/
+     4 48H SS   SS YY  YY  BB   BB   II   LL           111/
+     5 48H  SSS     Y  Y   BBBBBB    II   LL      ===   11/
+     6 48H    SSS    YY    BBBBBB    II   LL      ===   11/
+     7 48H SS   SS   YY    BB   BB   II   LL            11/
+     8 50H SSSSSSS   YY    BBBBBBB IIIIII LLLLLLL     111111/
+     9 50H  SSSSS    YY    BBBBBB  IIIIII LLLLLLL     111111//
+     1 1X,A72//)
+  110 FORMAT(/14H STATE VECTOR:/
+     1 7H NREG  ,I7,22H   (NUMBER OF REGIONS)/
+     2 7H NUN   ,I7,28H   (NUMBER OF FLUX UNKNOWNS)/
+     3 7H ILK   ,I7,39H   (0=LEAKAGE PRESENT/1=LEAKAGE ABSENT)/
+     4 7H NBMIX ,I7,36H   (MAXIMUM NUMBER OF MIXTURES USED)/
+     5 7H NSURF ,I7,29H   (NUMBER OF OUTER SURFACES)/
+     6 7H ITG   ,I7,21H   (TYPE OF GEOMETRY)/
+     7 5H MAXZ  ,I9,20H   (TRACKING LENGTH)/
+     8 5H MAXJ  ,I9,41H   (INTERNAL STORAGE LENGTH FOR CURRENTS)/
+     9 7H NUNCUR,I7,47H   (NUMBER OF ADDITIONAL INTERFACE CURRENT COMP,
+     1 7HONENTS)/
+     2 7H IBIHET,I7,46H   (0/1=DOUBLE HETEROGENEITY IS NOT/IS ACTIVE))
+  120 FORMAT(5H EPSJ,1P,E11.2,33H (FLUX-CURRENT ITERATION EPSILON))
+      END