diff options
Diffstat (limited to 'Trivac/src/FLDTRS.f')
| -rwxr-xr-x | Trivac/src/FLDTRS.f | 570 |
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 |