103 lines
3.3 KiB
Fortran
103 lines
3.3 KiB
Fortran
SUBROUTINE WGTJH1
|
|
C =================
|
|
C
|
|
C Angle quadrature weights
|
|
C from Hummer, Kunasz, & Kunasz, 1973, Comp. Phys. Comm. 6, 38
|
|
C
|
|
C The present version of this routine assumes that there are
|
|
C impact rays tangent to every depth layers (i.e. NREXT=ND)
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'WINCOM.FOR'
|
|
PARAMETER (UN=1., TWO=2., HALF=0.5)
|
|
PARAMETER (SIX=6.)
|
|
PARAMETER (C03=UN/3.,D03=2./3.,C04=UN/4.,C06=UN/6.)
|
|
PARAMETER (C24=UN/24.,C45=UN/45.,D45=2./45.,C72=UN/72.)
|
|
DIMENSION WAJ(MKU),WBJ(MKU),AHH(MKU,4)
|
|
DIMENSION BMUH(MKU),BMUHP(MKU),WAH(MKU),WBH(MKU)
|
|
DIMENSION WSD(MKU),WSU(MKU),WSL(MKU),WUU(MKU)
|
|
DIMENSION WTD(MKU),WTU(MKU),WTL(MKU)
|
|
C
|
|
DO 100 ID=1,ND
|
|
DO IU=ID+1,KMU
|
|
AHH(IU,1)=BMU(IU,ID)-BMU(IU-1,ID)
|
|
AHH(IU,2)=AHH(IU,1)*AHH(IU,1)
|
|
AHH(IU,3)=AHH(IU,2)*AHH(IU,1)
|
|
AHH(IU,4)=AHH(IU,3)*AHH(IU,1)
|
|
BMUH(IU)=BMU(IU,ID)*AHH(IU,1)
|
|
BMUHP(IU)=BMU(IU-1,ID)*AHH(IU,1)
|
|
END DO
|
|
C
|
|
C Weights for J
|
|
C
|
|
WAJ(ID)=HALF*AHH(ID+1,1)
|
|
WAJ(KMU)=HALF*AHH(KMU,1)
|
|
WBJ(ID)=-C24*AHH(ID+1,3)
|
|
WBJ(KMU)=-C24*AHH(KMU,3)
|
|
WSL(ID+1)=C06*AHH(ID+1,1)
|
|
WSU(KMU-1)=0.
|
|
WSD(ID)=C03*AHH(ID+1,1)
|
|
WSD(KMU)=UN
|
|
WTL(ID+1)=UN/AHH(ID+1,1)
|
|
WTU(KMU-1)=0.
|
|
WTD(ID)=-WTL(ID+1)
|
|
WTD(KMU)=0.
|
|
DO IU=ID+1,KMU-1
|
|
WAJ(IU)=HALF*(AHH(IU,1)+AHH(IU+1,1))
|
|
WBJ(IU)=-C24*(AHH(IU+1,3)+AHH(IU,3))
|
|
AH1=SIX/(AHH(IU,1)+AHH(IU+1,1))
|
|
WSL(IU+1)=C06*AH1*AHH(IU+1,1)
|
|
WSU(IU-1)=UN-WSL(IU+1)
|
|
WSD(IU)=TWO
|
|
WTL(IU+1)=AH1/AHH(IU+1,1)
|
|
WTU(IU-1)=AH1/AHH(IU,1)
|
|
WTD(IU)=-SIX/AHH(IU,1)/AHH(IU+1,1)
|
|
END DO
|
|
NMUD=KMU-ID+1
|
|
CALL TRIDAG(WSL,WSD,WSU,WBJ,WUU,NMUD)
|
|
WMUJ(ID,ID)=WAJ(ID)+WTD(ID)*WUU(ID)+WTU(ID)*WUU(ID+1)
|
|
WMUJ(KMU,ID)=WAJ(KMU)+WTL(KMU)*WUU(KMU-1)+WTD(KMU)*WUU(KMU)
|
|
DO IU=ID+1,KMU-1
|
|
WMUJ(IU,ID)=WAJ(IU)+WTL(IU)*WUU(IU-1)+
|
|
* WTD(IU)*WUU(IU)+WTU(IU)*WUU(IU+1)
|
|
END DO
|
|
C
|
|
C Weights for emergent flux H
|
|
C
|
|
IF(ID.GT.1) GO TO 100
|
|
WAH(ID)=HALF*BMUH(ID+1)-C03*AHH(ID+1,2)
|
|
WAH(KMU)=HALF*BMUHP(KMU)+C03*AHH(KMU,2)
|
|
WBH(ID)=AHH(ID+1,3)*(C45*AHH(ID+1,1)-C24*BMU(ID+1,ID))
|
|
WBH(KMU)=-AHH(KMU,3)*(C45*AHH(KMU,1)+C24*BMU(KMU-1,ID))
|
|
WSL(ID+1)=0.
|
|
WSD(ID)=UN
|
|
WTL(ID+1)=0.
|
|
WTD(ID)=0.
|
|
DO IU=ID+1,KMU-1
|
|
WAH(IU)=HALF*(BMUH(IU+1)+BMUHP(IU))-
|
|
* C03*(AHH(IU+1,2)-AHH(IU,2))
|
|
WBH(IU)=-C24*(BMUH(IU+1)*AHH(IU+1,2)+BMUHP(IU)*AHH(IU,2))+
|
|
* C45*(AHH(IU+1,4)-AHH(IU,4))
|
|
END DO
|
|
CALL TRIDAG(WSL,WSD,WSU,WBH,WUU,NMUD)
|
|
WMUH(ID)=WAH(ID)+WTD(ID)*WUU(ID)+WTU(ID)*WUU(ID+1)
|
|
WMUH(KMU)=WAH(KMU)+WTL(KMU)*WUU(KMU-1)+WTD(KMU)*WUU(KMU)
|
|
DO IU=ID+1,KMU-1
|
|
WMUH(IU)=WAH(IU)+WTL(IU)*WUU(IU-1)+
|
|
* WTD(IU)*WUU(IU)+WTU(IU)*WUU(IU+1)
|
|
END DO
|
|
C
|
|
100 CONTINUE
|
|
C
|
|
C Weights for H are overwritten by trapezoidal weigths
|
|
C
|
|
id=1
|
|
wmuh(1)=bmu(1,id)*(bmu(2,id)-bmu(1,id))*half
|
|
wmuh(kmu)=bmu(kmu,id)*(bmu(kmu,id)-bmu(kmu-1,id))*half
|
|
do iu=2,kmu-1
|
|
wmuh(iu)=bmu(iu,id)*(bmu(iu+1,id)-bmu(iu-1,id))*half
|
|
end do
|
|
c
|
|
RETURN
|
|
END
|