Initial commit from Polytechnique Montreal

1 files changed, 570 insertions, 0 deletions

diff --git a/Trivac/src/FLDTRS.f b/Trivac/src/FLDTRS.f
new file mode 100755
index 0000000..3cd6367
--- /dev/null
+++ b/Trivac/src/FLDTRS.f
@@ -0,0 +1,570 @@
+*DECK FLDTRS
+      SUBROUTINE FLDTRS(NAMP,IPTRK,IPSYS,LL4,S1,F1,NADI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform NADI inner iterations with the ADI preconditionning.  Special
+* version for Thomas-Raviart or Raviart-Thomas-Schneider basis.
+*
+*Copyright:
+* Copyright (C) 2006 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
+*
+*Reference:
+* A. Hebert, "A Raviart-Thomas-Schneider implementation of the
+* simplified Pn method in 3-D hexagonal geometry,"  PHYSOR 2010 -
+* Int. Conf. on Advances in Reactor Physics to Power the Nuclear
+* Renaissance, May 9-14, Pittsburgh, Pennsylvania, 2010.
+*
+*Parameters: input
+* NAMP    name of the ADI-splitted matrix.
+* IPTRK   L_TRACK pointer to the tracking information.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* LL4     order of the matrix.
+* S1      source term of the linear system.
+* F1      initial solution of the linear system.
+* NADI    number of inner ADI iterations.
+*
+*Parameters: output
+* F1      solution of the linear system after NADI iterations.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPSYS
+      CHARACTER NAMP*12
+      INTEGER LL4,NADI
+      REAL S1(LL4),F1(LL4)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (NSTATE=40)
+      CHARACTER NAMT*12
+      INTEGER ITP(NSTATE)
+      LOGICAL LMUX,DIAG
+      REAL, DIMENSION(:), ALLOCATABLE :: FL,FW,FX,FY,FZ,T,GAR
+      INTEGER C11W_LEN,C11X_LEN,C11Y_LEN,C11Z_LEN
+      INTEGER, DIMENSION(:), POINTER :: KN,IPERT,IPBW,MUW,IPVW,NBLW,
+     1 LBLW,IPBX,MUX,IPVX,NBLX,LBLX,IPBY,MUY,IPVY,NBLY,LBLY,IPBZ,MUZ,
+     2 IPVZ,NBLZ,LBLZ
+      REAL, DIMENSION(:), POINTER :: TF,DIFF,BW,C11W,BX,C11X,BY,C11Y,
+     1 BZ,C11Z
+      DOUBLE PRECISION, DIMENSION(:), POINTER :: CTRAN
+      TYPE(C_PTR) KN_PTR,IPERT_PTR,DIFF_PTR,TF_PTR,CTRAN_PTR,
+     1 BW_PTR,C11W_PTR,IPBW_PTR,MUW_PTR,IPVW_PTR,NBLW_PTR,LBLW_PTR,
+     2 BX_PTR,C11X_PTR,IPBX_PTR,MUX_PTR,IPVX_PTR,NBLX_PTR,LBLX_PTR,
+     3 BY_PTR,C11Y_PTR,IPBY_PTR,MUY_PTR,IPVY_PTR,NBLY_PTR,LBLY_PTR,
+     4 BZ_PTR,C11Z_PTR,IPBZ_PTR,MUZ_PTR,IPVZ_PTR,NBLZ_PTR,LBLZ_PTR
+*----
+*  INITIALIZATION
+*----
+      NAMT=NAMP
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ITP)
+      IELEM=ITP(9)
+      ISEG=ITP(17)
+      LTSW=ITP(19)
+      LL4F=ITP(25)
+      LL4W=ITP(26)
+      LL4X=ITP(27)
+      LL4Y=ITP(28)
+      LL4Z=ITP(29)
+      NLF=ITP(30)
+      IOFW=LL4F
+      IOFX=LL4F+LL4W
+      IOFY=LL4F+LL4W+LL4X
+      IOFZ=LL4F+LL4W+LL4X+LL4Y
+      CALL LCMLEN(IPTRK,'MUX',IDUM,ITYLCM)
+      LMUX=(IDUM.NE.0).AND.(ITYLCM.EQ.1)
+      DIAG=(LL4Y.GT.0).AND.(.NOT.LMUX)
+      CALL LCMGPD(IPSYS,'TF'//NAMT,TF_PTR)
+      CALL C_F_POINTER(TF_PTR,TF,(/ LL4F*NLF/2 /))
+*
+      NULLIFY(IPBW)
+      NULLIFY(IPVW)
+      NULLIFY(BW)
+      IF(LL4W.GT.0) THEN
+         ISPLH=ITP(13)
+         LX=ITP(14)
+         LZ=ITP(16)
+         NBLOS=LX*LZ/3
+         CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
+         CALL LCMGPD(IPTRK,'CTRAN',CTRAN_PTR)
+         CALL LCMGPD(IPTRK,'KN',KN_PTR)
+         CALL LCMGPD(IPTRK,'IPERT',IPERT_PTR)
+         CALL LCMGPD(IPSYS,'DIFF'//NAMT,DIFF_PTR)
+         CALL C_F_POINTER(CTRAN_PTR,CTRAN,(/ ((IELEM+1)*IELEM)**2 /))
+         CALL C_F_POINTER(KN_PTR,KN,(/ MAXKN /))
+         CALL C_F_POINTER(IPERT_PTR,IPERT,(/ NBLOS /))
+         CALL C_F_POINTER(DIFF_PTR,DIFF,(/ NBLOS /))
+*
+         ALLOCATE(FW(LL4W))
+         CALL LCMGPD(IPTRK,'IPBBW',IPBW_PTR)
+         CALL LCMLEN(IPSYS,'WB',LENWB,ITYL)
+         IF(LENWB.EQ.0)THEN
+           CALL LCMGPD(IPTRK,'WB',BW_PTR)
+         ELSE
+           CALL LCMGPD(IPSYS,'WB',BW_PTR)
+         ENDIF
+         CALL C_F_POINTER(IPBW_PTR,IPBW,(/ 2*IELEM*LL4W /))
+         CALL C_F_POINTER(BW_PTR,BW,(/ 2*IELEM*LL4W /))
+         CALL LCMLEN(IPSYS,'WI'//NAMT,C11W_LEN,ITYLCM)
+         CALL LCMGPD(IPSYS,'WI'//NAMT,C11W_PTR)
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A W-ORIENTED LINEAR SYSTEM.
+            CALL LCMGPD(IPTRK,'MUW',MUW_PTR)
+            CALL C_F_POINTER(MUW_PTR,MUW,(/ LL4W /))
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A W-ORIENTED LINEAR SYSTEM.
+            CALL LCMGET(IPTRK,'LL4VW',LL4VW)
+            CALL LCMGPD(IPTRK,'MUVW',MUW_PTR)
+            CALL LCMGPD(IPTRK,'IPVW',IPVW_PTR)
+            CALL LCMLEN(IPTRK,'NBLW',LONW,ITYLCM)
+            CALL LCMGPD(IPTRK,'NBLW',NBLW_PTR)
+            CALL LCMGPD(IPTRK,'LBLW',LBLW_PTR)
+            CALL C_F_POINTER(MUW_PTR,MUW,(/ LL4VW/ISEG /))
+            CALL C_F_POINTER(IPVW_PTR,IPVW,(/ LL4W /))
+            CALL C_F_POINTER(NBLW_PTR,NBLW,(/ LONW /))
+            CALL C_F_POINTER(LBLW_PTR,LBLW,(/ LONW /))
+         ENDIF
+         CALL C_F_POINTER(C11W_PTR,C11W,(/ C11W_LEN /))
+      ENDIF
+      ALLOCATE(FX(LL4X))
+      DO 10 I0=1,LL4X
+      FX(I0)=F1(IOFX+I0)
+   10 CONTINUE
+      CALL LCMGPD(IPTRK,'IPBBX',IPBX_PTR)
+      CALL LCMLEN(IPSYS,'XB',LENXB,ITYL)
+      IF(LENXB.EQ.0) THEN
+        CALL LCMGPD(IPTRK,'XB',BX_PTR)
+      ELSE
+        CALL LCMGPD(IPSYS,'XB',BX_PTR)
+      ENDIF
+      CALL C_F_POINTER(IPBX_PTR,IPBX,(/ 2*IELEM*LL4X /))
+      CALL C_F_POINTER(BX_PTR,BX,(/ 2*IELEM*LL4X /))
+      NULLIFY(IPVX)
+      IF(DIAG) THEN
+         CALL LCMLEN(IPSYS,'YI'//NAMT,C11X_LEN,ITYLCM)
+         CALL LCMGPD(IPSYS,'YI'//NAMT,C11X_PTR)
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A X-ORIENTED LINEAR SYSTEM.
+            CALL LCMGPD(IPTRK,'MUY',MUX_PTR)
+            CALL C_F_POINTER(MUX_PTR,MUX,(/ LL4X /))
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A X-ORIENTED LINEAR SYSTEM.
+            CALL LCMGET(IPTRK,'LL4VY',LL4VX)
+            CALL LCMGPD(IPTRK,'MUVY',MUX_PTR)
+            CALL LCMGPD(IPTRK,'IPVY',IPVX_PTR)
+            CALL LCMLEN(IPTRK,'NBLY',LONX,ITYLCM)
+            CALL LCMGPD(IPTRK,'NBLY',NBLX_PTR)
+            CALL LCMGPD(IPTRK,'LBLY',LBLX_PTR)
+            CALL C_F_POINTER(MUX_PTR,MUX,(/ LL4VX/ISEG /))
+            CALL C_F_POINTER(IPVX_PTR,IPVX,(/ LL4X /))
+            CALL C_F_POINTER(NBLX_PTR,NBLX,(/ LONX /))
+            CALL C_F_POINTER(LBLX_PTR,LBLX,(/ LONX /))
+         ENDIF
+      ELSE
+         CALL LCMLEN(IPSYS,'XI'//NAMT,C11X_LEN,ITYLCM)
+         CALL LCMGPD(IPSYS,'XI'//NAMT,C11X_PTR)
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A X-ORIENTED LINEAR SYSTEM.
+            CALL LCMGPD(IPTRK,'MUX',MUX_PTR)
+            CALL C_F_POINTER(MUX_PTR,MUX,(/ LL4X /))
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A X-ORIENTED LINEAR SYSTEM.
+            CALL LCMGET(IPTRK,'LL4VX',LL4VX)
+            CALL LCMGPD(IPTRK,'MUVX',MUX_PTR)
+            CALL LCMGPD(IPTRK,'IPVX',IPVX_PTR)
+            CALL LCMLEN(IPTRK,'NBLX',LONX,ITYLCM)
+            CALL LCMGPD(IPTRK,'NBLX',NBLX_PTR)
+            CALL LCMGPD(IPTRK,'LBLX',LBLX_PTR)
+            CALL C_F_POINTER(MUX_PTR,MUX,(/ LL4VX/ISEG /))
+            CALL C_F_POINTER(IPVX_PTR,IPVX,(/ LL4X /))
+            CALL C_F_POINTER(NBLX_PTR,NBLX,(/ LONX /))
+            CALL C_F_POINTER(LBLX_PTR,LBLX,(/ LONX /))
+         ENDIF
+      ENDIF
+      CALL C_F_POINTER(C11X_PTR,C11X,(/ C11X_LEN /))
+      NULLIFY(IPBY)
+      NULLIFY(IPVY)
+      NULLIFY(BY)
+      IF(LL4Y.GT.0) THEN
+         ALLOCATE(FY(LL4Y))
+         DO 20 I0=1,LL4Y
+         FY(I0)=F1(IOFY+I0)
+   20    CONTINUE
+         CALL LCMGPD(IPTRK,'IPBBY',IPBY_PTR)
+         CALL LCMLEN(IPSYS,'YB',LENYB,ITYL)
+         IF(LENYB.EQ.0) THEN
+            CALL LCMGPD(IPTRK,'YB',BY_PTR)
+         ELSE
+            CALL LCMGPD(IPSYS,'YB',BY_PTR)
+         ENDIF
+         CALL C_F_POINTER(IPBY_PTR,IPBY,(/ 2*IELEM*LL4Y /))
+         CALL C_F_POINTER(BY_PTR,BY,(/ 2*IELEM*LL4Y /))
+         CALL LCMLEN(IPSYS,'YI'//NAMT,C11Y_LEN,ITYLCM)
+         CALL LCMGPD(IPSYS,'YI'//NAMT,C11Y_PTR)
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A Y-ORIENTED LINEAR SYSTEM.
+            CALL LCMGPD(IPTRK,'MUY',MUY_PTR)
+            CALL C_F_POINTER(MUY_PTR,MUY,(/ LL4Y /))
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A Y-ORIENTED LINEAR SYSTEM.
+            CALL LCMGET(IPTRK,'LL4VY',LL4VY)
+            CALL LCMGPD(IPTRK,'MUVY',MUY_PTR)
+            CALL LCMGPD(IPTRK,'IPVY',IPVY_PTR)
+            CALL LCMLEN(IPTRK,'NBLY',LONY,ITYLCM)
+            CALL LCMGPD(IPTRK,'NBLY',NBLY_PTR)
+            CALL LCMGPD(IPTRK,'LBLY',LBLY_PTR)
+            CALL C_F_POINTER(MUY_PTR,MUY,(/ LL4VY/ISEG /))
+            CALL C_F_POINTER(IPVY_PTR,IPVY,(/ LL4Y /))
+            CALL C_F_POINTER(NBLY_PTR,NBLY,(/ LONY /))
+            CALL C_F_POINTER(LBLY_PTR,LBLY,(/ LONY /))
+         ENDIF
+         CALL C_F_POINTER(C11Y_PTR,C11Y,(/ C11Y_LEN /))
+      ENDIF
+      NULLIFY(IPBZ)
+      NULLIFY(IPVZ)
+      NULLIFY(BZ)
+      IF(LL4Z.GT.0) THEN
+         ALLOCATE(FZ(LL4Z))
+         DO 30 I0=1,LL4Z
+         FZ(I0)=F1(IOFZ+I0)
+   30    CONTINUE
+         CALL LCMGPD(IPTRK,'IPBBZ',IPBZ_PTR)
+         CALL LCMLEN(IPSYS,'ZB',LENZB,ITYL)
+         IF(LENZB.EQ.0) THEN
+            CALL LCMGPD(IPTRK,'ZB',BZ_PTR)
+         ELSE
+            CALL LCMGPD(IPSYS,'ZB',BZ_PTR)
+         ENDIF
+         CALL C_F_POINTER(IPBZ_PTR,IPBZ,(/ 2*IELEM*LL4Z /))
+         CALL C_F_POINTER(BZ_PTR,BZ,(/ 2*IELEM*LL4Z /))
+         CALL LCMLEN(IPSYS,'ZI'//NAMT,C11Z_LEN,ITYLCM)
+         CALL LCMGPD(IPSYS,'ZI'//NAMT,C11Z_PTR)
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A Z-ORIENTED LINEAR SYSTEM.
+            CALL LCMGPD(IPTRK,'MUZ',MUZ_PTR)
+            CALL C_F_POINTER(MUZ_PTR,MUZ,(/ LL4Z /))
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A Z-ORIENTED LINEAR SYSTEM.
+            CALL LCMGET(IPTRK,'LL4VZ',LL4VZ)
+            CALL LCMGPD(IPTRK,'MUVZ',MUZ_PTR)
+            CALL LCMGPD(IPTRK,'IPVZ',IPVZ_PTR)
+            CALL LCMLEN(IPTRK,'NBLZ',LONZ,ITYLCM)
+            CALL LCMGPD(IPTRK,'NBLZ',NBLZ_PTR)
+            CALL LCMGPD(IPTRK,'LBLZ',LBLZ_PTR)
+            CALL C_F_POINTER(MUZ_PTR,MUZ,(/ LL4VZ/ISEG /))
+            CALL C_F_POINTER(IPVZ_PTR,IPVZ,(/ LL4Z /))
+            CALL C_F_POINTER(NBLZ_PTR,NBLZ,(/ LONZ /))
+            CALL C_F_POINTER(LBLZ_PTR,LBLZ,(/ LONZ /))
+         ENDIF
+         CALL C_F_POINTER(C11Z_PTR,C11Z,(/ C11Z_LEN /))
+      ENDIF
+      ALLOCATE(FL(LL4F))
+*----
+*  W DIRECTION
+*----
+      IF(ISEG.GT.0) ALLOCATE(T(ISEG))
+      DO 520 IADI=1,NADI
+      IF(LL4W.GT.0) THEN
+         DO 40 I0=1,LL4F
+         FL(I0)=S1(I0)
+   40    CONTINUE
+         DO 60 I0=1,LL4X
+         DO 50 J0=1,2*IELEM
+         JJ=IPBX((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 60
+         FL(JJ)=FL(JJ)-BX((I0-1)*2*IELEM+J0)*FX(I0)
+   50    CONTINUE
+   60    CONTINUE
+         DO 80 I0=1,LL4Y
+         DO 70 J0=1,2*IELEM
+         JJ=IPBY((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 80
+         FL(JJ)=FL(JJ)-BY((I0-1)*2*IELEM+J0)*FY(I0)
+   70    CONTINUE
+   80    CONTINUE
+         DO 100 I0=1,LL4Z
+         DO 90 J0=1,2*IELEM
+         JJ=IPBZ((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 100
+         FL(JJ)=FL(JJ)-BZ((I0-1)*2*IELEM+J0)*FZ(I0)
+   90    CONTINUE
+  100    CONTINUE
+         DO 130 I0=1,LL4W
+         GGW=-S1(IOFW+I0)
+         DO 110 J0=1,2*IELEM
+         JJ=IPBW((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 120
+         GGW=GGW+BW((I0-1)*2*IELEM+J0)*FL(JJ)/TF(JJ)
+  110    CONTINUE
+  120    FW(I0)=GGW
+  130    CONTINUE
+*
+*        PIOLAT TRANSFORM TERM.
+         CALL FLDPWY(LL4W,LL4X,LL4Y,NBLOS,IELEM,CTRAN,IPERT,KN,DIFF,
+     1   FY,FW)
+         CALL FLDPWX(LL4W,LL4X,NBLOS,IELEM,CTRAN,IPERT,KN,DIFF,FX,FW)
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A W-ORIENTED LINEAR SYSTEM.
+            CALL ALLDLS(LL4W,MUW,C11W,FW)
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A W-ORIENTED LINEAR SYSTEM.
+            ALLOCATE(GAR(LL4VW))
+            GAR(:LL4VW)=0.0
+            DO 140 I=1,LL4W
+            GAR(IPVW(I))=FW(I)
+  140       CONTINUE
+            CALL ALVDLS(LTSW,MUW,C11W,GAR,ISEG,LONW,NBLW,LBLW,T)
+            DO 150 I=1,LL4W
+            FW(I)=GAR(IPVW(I))
+  150       CONTINUE
+            DEALLOCATE(GAR)
+         ENDIF
+      ENDIF
+*----
+*  X DIRECTION
+*----
+      DO 160 I0=1,LL4F
+      FL(I0)=S1(I0)
+  160 CONTINUE
+      DO 180 I0=1,LL4W
+      DO 170 J0=1,2*IELEM
+      JJ=IPBW((I0-1)*2*IELEM+J0)
+      IF(JJ.EQ.0) GO TO 180
+      FL(JJ)=FL(JJ)-BW((I0-1)*2*IELEM+J0)*FW(I0)
+  170 CONTINUE
+  180 CONTINUE
+      DO 200 I0=1,LL4Y
+      DO 190 J0=1,2*IELEM
+      JJ=IPBY((I0-1)*2*IELEM+J0)
+      IF(JJ.EQ.0) GO TO 200
+      FL(JJ)=FL(JJ)-BY((I0-1)*2*IELEM+J0)*FY(I0)
+  190 CONTINUE
+  200 CONTINUE
+      DO 220 I0=1,LL4Z
+      DO 210 J0=1,2*IELEM
+      JJ=IPBZ((I0-1)*2*IELEM+J0)
+      IF(JJ.EQ.0) GO TO 220
+      FL(JJ)=FL(JJ)-BZ((I0-1)*2*IELEM+J0)*FZ(I0)
+  210 CONTINUE
+  220 CONTINUE
+      DO 250 I0=1,LL4X
+      GGX=-S1(IOFX+I0)
+      DO 230 J0=1,2*IELEM
+      JJ=IPBX((I0-1)*2*IELEM+J0)
+      IF(JJ.EQ.0) GO TO 240
+      GGX=GGX+BX((I0-1)*2*IELEM+J0)*FL(JJ)/TF(JJ)
+  230 CONTINUE
+  240 FX(I0)=GGX
+  250 CONTINUE
+      IF(LL4W.GT.0) THEN
+*        PIOLAT TRANSFORM TERM.
+         CALL FLDPXW(LL4W,LL4X,NBLOS,IELEM,CTRAN,IPERT,KN,DIFF,FW,FX)
+         CALL FLDPXY(LL4W,LL4X,LL4Y,NBLOS,IELEM,CTRAN,IPERT,KN,DIFF,
+     1   FY,FX)
+      ENDIF
+      IF(ISEG.EQ.0) THEN
+*        SCALAR SOLUTION FOR A X-ORIENTED LINEAR SYSTEM.
+         CALL ALLDLS(LL4X,MUX,C11X,FX)
+      ELSE IF(ISEG.GT.0) THEN
+*        SUPERVECTORIAL SOLUTION FOR A X-ORIENTED LINEAR SYSTEM.
+         ALLOCATE(GAR(LL4VX))
+         GAR(:LL4VX)=0.0
+         DO 260 I=1,LL4X
+         GAR(IPVX(I))=FX(I)
+  260    CONTINUE
+         CALL ALVDLS(LTSW,MUX,C11X,GAR,ISEG,LONX,NBLX,LBLX,T)
+         DO 270 I=1,LL4X
+         FX(I)=GAR(IPVX(I))
+  270    CONTINUE
+         DEALLOCATE(GAR)
+      ENDIF
+*----
+*  Y DIRECTION
+*----
+      IF(LL4Y.GT.0) THEN
+         DO 280 I0=1,LL4F
+         FL(I0)=S1(I0)
+  280    CONTINUE
+         DO 300 I0=1,LL4W
+         DO 290 J0=1,2*IELEM
+         JJ=IPBW((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 300
+         FL(JJ)=FL(JJ)-BW((I0-1)*2*IELEM+J0)*FW(I0)
+  290    CONTINUE
+  300    CONTINUE
+         DO 320 I0=1,LL4X
+         DO 310 J0=1,2*IELEM
+         JJ=IPBX((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 320
+         FL(JJ)=FL(JJ)-BX((I0-1)*2*IELEM+J0)*FX(I0)
+  310    CONTINUE
+  320    CONTINUE
+         DO 340 I0=1,LL4Z
+         DO 330 J0=1,2*IELEM
+         JJ=IPBZ((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 340
+         FL(JJ)=FL(JJ)-BZ((I0-1)*2*IELEM+J0)*FZ(I0)
+  330    CONTINUE
+  340    CONTINUE
+         DO 370 I0=1,LL4Y
+         GGY=-S1(IOFY+I0)
+         DO 350 J0=1,2*IELEM
+         JJ=IPBY((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 360
+         GGY=GGY+BY((I0-1)*2*IELEM+J0)*FL(JJ)/TF(JJ)
+  350    CONTINUE
+  360    FY(I0)=GGY
+  370    CONTINUE
+         IF(LL4W.GT.0) THEN
+*           PIOLAT TRANSFORM TERM.
+            CALL FLDPYX(LL4W,LL4X,LL4Y,NBLOS,IELEM,CTRAN,IPERT,KN,
+     1      DIFF,FX,FY)
+            CALL FLDPYW(LL4W,LL4X,LL4Y,NBLOS,IELEM,CTRAN,IPERT,KN,
+     1      DIFF,FW,FY)
+         ENDIF
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A Y-ORIENTED LINEAR SYSTEM.
+            CALL ALLDLS(LL4Y,MUY,C11Y,FY)
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A Y-ORIENTED LINEAR SYSTEM.
+            ALLOCATE(GAR(LL4VY))
+            GAR(:LL4VY)=0.0
+            DO 380 I=1,LL4Y
+            GAR(IPVY(I))=FY(I)
+  380       CONTINUE
+            CALL ALVDLS(LTSW,MUY,C11Y,GAR,ISEG,LONY,NBLY,LBLY,T)
+            DO 390 I=1,LL4Y
+            FY(I)=GAR(IPVY(I))
+  390       CONTINUE
+            DEALLOCATE(GAR)
+         ENDIF
+      ENDIF
+*----
+*  Z DIRECTION
+*----
+      IF(LL4Z.GT.0) THEN
+         DO 400 I0=1,LL4F
+         FL(I0)=S1(I0)
+  400    CONTINUE
+         DO 420 I0=1,LL4W
+         DO 410 J0=1,2*IELEM
+         JJ=IPBW((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 420
+         FL(JJ)=FL(JJ)-BW((I0-1)*2*IELEM+J0)*FW(I0)
+  410    CONTINUE
+  420    CONTINUE
+         DO 440 I0=1,LL4X
+         DO 430 J0=1,2*IELEM
+         JJ=IPBX((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 440
+         FL(JJ)=FL(JJ)-BX((I0-1)*2*IELEM+J0)*FX(I0)
+  430    CONTINUE
+  440    CONTINUE
+         DO 460 I0=1,LL4Y
+         DO 450 J0=1,2*IELEM
+         JJ=IPBY((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 460
+         FL(JJ)=FL(JJ)-BY((I0-1)*2*IELEM+J0)*FY(I0)
+  450    CONTINUE
+  460    CONTINUE
+         DO 490 I0=1,LL4Z
+         GGZ=-S1(IOFZ+I0)
+         DO 470 J0=1,2*IELEM
+         JJ=IPBZ((I0-1)*2*IELEM+J0)
+         IF(JJ.EQ.0) GO TO 480
+         GGZ=GGZ+BZ((I0-1)*2*IELEM+J0)*FL(JJ)/TF(JJ)
+  470    CONTINUE
+  480    FZ(I0)=GGZ
+  490    CONTINUE
+         IF(ISEG.EQ.0) THEN
+*           SCALAR SOLUTION FOR A Z-ORIENTED LINEAR SYSTEM.
+            CALL ALLDLS(LL4Z,MUZ,C11Z,FZ)
+         ELSE IF(ISEG.GT.0) THEN
+*           SUPERVECTORIAL SOLUTION FOR A Z-ORIENTED LINEAR SYSTEM.
+            ALLOCATE(GAR(LL4VZ))
+            GAR(:LL4VZ)=0.0
+            DO 500 I=1,LL4Z
+            GAR(IPVZ(I))=FZ(I)
+  500       CONTINUE
+            CALL ALVDLS(LTSW,MUZ,C11Z,GAR,ISEG,LONZ,NBLZ,LBLZ,T)
+            DO 510 I=1,LL4Z
+            FZ(I)=GAR(IPVZ(I))
+  510       CONTINUE
+            DEALLOCATE(GAR)
+         ENDIF
+      ENDIF
+  520 CONTINUE
+      IF(ISEG.GT.0) DEALLOCATE(T)
+      DEALLOCATE(FL)
+*----
+*  COMPUTE FLUX AND RECOVER CURRENTS
+*----
+      DO 530 I0=1,LL4F
+      F1(I0)=S1(I0)
+  530 CONTINUE
+      DO 550 J0=1,LL4W
+      DO 540 I0=1,2*IELEM
+      II=IPBW((J0-1)*2*IELEM+I0)
+      IF(II.EQ.0) GO TO 550
+      F1(II)=F1(II)-BW((J0-1)*2*IELEM+I0)*FW(J0)
+  540 CONTINUE
+  550 CONTINUE
+      DO 570 J0=1,LL4X
+      DO 560 I0=1,2*IELEM
+      II=IPBX((J0-1)*2*IELEM+I0)
+      IF(II.EQ.0) GO TO 570
+      F1(II)=F1(II)-BX((J0-1)*2*IELEM+I0)*FX(J0)
+  560 CONTINUE
+  570 CONTINUE
+      DO 590 J0=1,LL4Y
+      DO 580 I0=1,2*IELEM
+      II=IPBY((J0-1)*2*IELEM+I0)
+      IF(II.EQ.0) GO TO 590
+      F1(II)=F1(II)-BY((J0-1)*2*IELEM+I0)*FY(J0)
+  580 CONTINUE
+  590 CONTINUE
+      DO 610 J0=1,LL4Z
+      DO 600 I0=1,2*IELEM
+      II=IPBZ((J0-1)*2*IELEM+I0)
+      IF(II.EQ.0) GO TO 610
+      F1(II)=F1(II)-BZ((J0-1)*2*IELEM+I0)*FZ(J0)
+  600 CONTINUE
+  610 CONTINUE
+      DO 620 I0=1,LL4F
+      F1(I0)=F1(I0)/TF(I0)
+  620 CONTINUE
+      IF(LL4W.GT.0) THEN
+         DO 630 I0=1,LL4W
+         F1(IOFW+I0)=FW(I0)
+  630    CONTINUE
+         DEALLOCATE(FW)
+      ENDIF
+      DO 640 I0=1,LL4X
+      F1(IOFX+I0)=FX(I0)
+  640 CONTINUE
+      DEALLOCATE(FX)
+      IF(LL4Y.GT.0) THEN
+         DO 650 I0=1,LL4Y
+         F1(IOFY+I0)=FY(I0)
+  650    CONTINUE
+         DEALLOCATE(FY)
+      ENDIF
+      IF(LL4Z.GT.0) THEN
+         DO 660 I0=1,LL4Z
+         F1(IOFZ+I0)=FZ(I0)
+  660    CONTINUE
+         DEALLOCATE(FZ)
+      ENDIF
+      RETURN
+      END