*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