*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