From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Trivac/src/NSSANM2.f90 | 603 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 603 insertions(+) create mode 100755 Trivac/src/NSSANM2.f90 (limited to 'Trivac/src/NSSANM2.f90') diff --git a/Trivac/src/NSSANM2.f90 b/Trivac/src/NSSANM2.f90 new file mode 100755 index 0000000..da08e30 --- /dev/null +++ b/Trivac/src/NSSANM2.f90 @@ -0,0 +1,603 @@ +subroutine NSSANM2(nunkn,nx,ny,ll4f,ll4x,ng,bndtl,npass,nmix,idl,kn,iqfr, & +& qfr,mat,xxx,yyy,keff,diff,sigr,chi,sigf,scat,fd,savg) +! +!----------------------------------------------------------------------- +! +!Purpose: +! Compute the ANM boundary fluxes and currents using a solution of +! one- and two-node relations in Cartesian 2D geometry. +! +!Copyright: +! Copyright (C) 2022 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 +! nunkn number of unknowns per energy group. +! nx number of x-nodes in the nodal calculation. +! ny number of y-nodes in the nodal calculation. +! ll4f number of averaged flux unknowns. +! ll4x number of X-directed net currents. +! ng number of energy groups. +! bndtl set to 'flat' or 'quadratic'. +! npass number of transverse current iterations. +! nmix number of material mixtures in the nodal calculation. +! idl position of averaged fluxes in unknown vector. +! kn element-ordered interface net current unknown list. +! iqfr node-ordered physical albedo indices. +! qfr albedo function information. +! mat material mixture index in eacn node. +! xxx Cartesian coordinates along the X axis. +! yyy Cartesian coordinates along the Y axis. +! keff effective multiplication facctor. +! diff diffusion coefficients +! sigr removal cross sections. +! chi fission spectra. +! sigf nu times fission cross section. +! scat scattering cross section. +! fd discontinuity factors. +! savg nodal fluxes and net currents. +! +!Parameters: output +! savg nodal fluxes, boundary fluxes and net currents. +! +!----------------------------------------------------------------------- + ! + !---- + ! subroutine arguments + !---- + integer,intent(in) :: nunkn,nx,ny,ll4f,ll4x,ng,npass,nmix,idl(nx,ny),kn(6,nx,ny), & + & iqfr(6,nx,ny),mat(nx,ny) + real,intent(in) :: qfr(6,nx,ny,ng),xxx(nx+1),yyy(ny+1),keff,diff(nmix,ng), & + & sigr(nmix,ng),chi(nmix,ng),sigf(nmix,ng),scat(nmix,ng,ng),fd(nmix,4,ng,ng) + real, dimension(nunkn,ng),intent(inout) :: savg + character(len=12), intent(in) :: bndtl + !---- + ! local and allocatable arrays + !---- + real :: xyz(4) + real, allocatable, dimension(:) :: work1,work2,work4,work5 + real, allocatable, dimension(:,:) :: A,Lambda,work3 + real(kind=8), allocatable, dimension(:,:) :: LLR + real(kind=8), allocatable, dimension(:,:,:,:) :: Lx,Rx,Ly,Ry + !---- + ! scratch storage allocation + !---- + allocate(A(ng,ng+1),Lambda(ng,ng),LLR(ng,8*ng)) + allocate(work1(ng),work2(ng),work3(ng,ng),work4(ng),work5(ng)) + allocate(Lx(ng,8*ng,nx,ny),Rx(ng,8*ng,nx,ny)) + allocate(Ly(ng,8*ng,nx,ny),Ry(ng,8*ng,nx,ny)) + !---- + ! compute 2D ANM coupling matrices for each single node + !---- + do j=1,ny + dely=yyy(j+1)-yyy(j) + do i=1,nx + delx=xxx(i+1)-xxx(i) + ibm=mat(i,j) + if(ibm == 0) cycle + work1(:ng)=diff(ibm,:ng) + work2(:ng)=sigr(ibm,:ng) + work3(:ng,:ng)=scat(ibm,:ng,:ng) + work4(:ng)=chi(ibm,:ng) + work5(:ng)=sigf(ibm,:ng) + ! + kk1=iqfr(1,i,j) + kk2=iqfr(2,i,j) + xyz(2:3)=xxx(i:i+1)-xxx(i) + if(kk1 == -2) then + ! reflection boundary condition + xyz(1)=2.0*xyz(2)-xyz(3) + else if(kk1 < 0) then + ! zero/void/albedo boundary condition + xyz(1)=-99999. + else + ! left neighbour + xyz(1)=xxx(i-1)-xxx(i) + endif + if(kk2 == -2) then + ! reflection boundary condition + xyz(4)=2.0*xyz(3)-xyz(2) + else if(kk2 < 0) then + ! zero/void/albedo boundary condition + xyz(4)=-99999. + else + ! right neighbour + xyz(4)=xxx(i+2)-xxx(i) + endif + call NSSLR2(keff,ng,bndtl,xyz,dely,work1,work2,work3,work4,work5,Lx(1,1,i,j),Rx(1,1,i,j)) + ! + kk3=iqfr(3,i,j) + kk4=iqfr(4,i,j) + xyz(2:3)=yyy(j:j+1)-yyy(j) + if(kk3 == -2) then + ! reflection boundary condition + xyz(1)=2.0*xyz(2)-xyz(3) + else if(kk3 < 0) then + ! zero/void/albedo boundary condition + xyz(1)=-99999. + else + ! left neighbour + xyz(1)=yyy(j-1)-yyy(j) + endif + if(kk4 == -2) then + ! reflection boundary condition + xyz(4)=2.0*xyz(3)-xyz(2) + else if(kk4 < 0) then + ! zero/void/albedo boundary condition + xyz(4)=-99999.0 + else + ! right neighbour + xyz(4)=yyy(j+2)-yyy(j) + endif + call NSSLR2(keff,ng,bndtl,xyz,delx,work1,work2,work3,work4,work5, & + & Ly(1,1,i,j),Ry(1,1,i,j)) + enddo + enddo + !---- + ! perform transverse current iterations + !---- + do ipass=1,npass + !---- + ! one- and two-node relations along X axis + !---- + do j=1,ny + nxmin=1 + do i=1,nx + if(mat(i,j) > 0) exit + nxmin=i+1 + enddo + if(nxmin > nx) cycle + nxmax=nx + do i=nx,1,-1 + if(mat(i,j) > 0) exit + nxmax=i-1 + enddo + ! one-node relation at left + ind1=idl(nxmin,j) + if(ind1 == 0) call XABORT('NSSANM2: invalid idl index.(1)') + iqf1=iqfr(1,nxmin,j) + jxm=kn(1,nxmin,j) ; jxp=kn(2,nxmin,j) ; jym=kn(3,nxmin,j) ; jyp=kn(4,nxmin,j) + jym_p=0 ; jyp_p=0 + if(nxmin < nx) then + jym_p=kn(3,nxmin+1,j) ; jyp_p=kn(4,nxmin+1,j) + endif + A(:ng,:ng+1)=0.0 + LLR(:ng,:8*ng)=0.0 + if((iqf1 > 0).or.(iqf1 == -1)) then + ! physical albedo + Lambda(:ng,:ng)=0.0 + do ig=1,ng + Lambda(ig,ig)=qfr(1,nxmin,j,ig) + enddo + LLR(:ng,:8*ng)=matmul(Lambda(:ng,:ng),Lx(:ng,:8*ng,nxmin,j)) + A(:ng,:ng)=-real(LLR(:ng,ng+1:2*ng),4) + do ig=1,ng + A(ig,ig)=-1.0+A(ig,ig) + enddo + else if(iqf1 == -2) then + ! zero net current + do ig=1,ng + A(ig,ig)=1.0 + enddo + else if(iqf1 == -3) then + ! zero flux + LLR(:ng,:8*ng)=Lx(:ng,:8*ng,nxmin,j) + A(:ng,:ng)=real(-LLR(:ng,ng+1:2*ng),4) + else if(iqf1 == -4) then + call XABORT('NSSANM2: SYME boundary condition is not supported.(1)') + else + call XABORT('NSSANM2: illegal left X-boundary condition.') + endif + if(iqf1 /= -2) then + A(:ng,ng+1)=real(matmul(LLR(:ng,:ng),savg(ind1,:ng)),4) + do ig=1,ng + do jg=1,ng + if(jym /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,3*ng+jg)*savg(5*ll4f+ll4x+jym,jg),4) + if(jym_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,4*ng+jg)*savg(5*ll4f+ll4x+jym_p,jg),4) + if(jyp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,6*ng+jg)*savg(5*ll4f+ll4x+jyp,jg),4) + if(jyp_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,7*ng+jg)*savg(5*ll4f+ll4x+jyp_p,jg),4) + enddo + enddo + endif + call ALSB(ng,1,A,ier,ng) + if(ier /= 0) call XABORT('NSSANM2: singular matrix.(1)') + savg(5*ll4f+jxm,:ng)=A(:ng,ng+1) + ! + ! two-node relations + do i=nxmin,nxmax-1 + ind1=idl(i,j) + if(ind1 == 0) call XABORT('NSSANM2: invalid idl index.(2)') + ind2=idl(i+1,j) + if(ind2 == 0) call XABORT('NSSANM2: invalid idl index.(3)') + if(kn(1,i+1,j) /= kn(2,i,j)) call XABORT('NSSANM2: invalid kn index.(1)') + if(iqfr(2,i,j) /= 0) call XABORT('NSSANM2: invalid iqfr index.(1)') + if(iqfr(1,i+1,j) /= 0) call XABORT('NSSANM2: invalid iqfr index.(2)') + jxm=kn(1,i,j) ; jxp=kn(2,i,j) ; jym=kn(3,i,j) ; jyp=kn(4,i,j) + jym_m=0 ; jyp_m=0 ; jym_pp=0 ; jyp_pp=0 + if((i == 1).and.(iqfr(1,1,j) == -2)) then + jym_m=kn(3,1,j) ; jyp_m=kn(4,1,j) + else if(i > 1) then + jym_m=kn(3,i-1,j) ; jyp_m=kn(4,i-1,j) + endif + jym_p=kn(3,i+1,j) ; jyp_p=kn(4,i+1,j) + if((i == nx-1).and.(iqfr(2,nx,j) == -2)) then + jym_pp=kn(3,nx,j) ; jyp_pp=kn(4,nx,j) + else if(i < nx-1) then + jym_pp=kn(3,i+2,j) ; jyp_pp=kn(4,i+2,j) + endif + ! + A(:ng,:ng+1)=0.0 + ! node i + LLR(:ng,:8*ng)=matmul(fd(mat(i,j),2,:ng,:ng),Rx(:ng,:8*ng,i,j)) + do ig=1,ng + A(:ng,ig)=A(:ng,ig)+real(LLR(:ng,ng+ig),4) + do jg=1,ng + A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,jg)*savg(ind1,jg),4) + if(jym_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,2*ng+jg)*savg(5*ll4f+ll4x+jym_m,jg),4) + if(jym /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,3*ng+jg)*savg(5*ll4f+ll4x+jym,jg),4) + if(jym_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,4*ng+jg)*savg(5*ll4f+ll4x+jym_p,jg),4) + if(jyp_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,5*ng+jg)*savg(5*ll4f+ll4x+jyp_m,jg),4) + if(jyp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,6*ng+jg)*savg(5*ll4f+ll4x+jyp,jg),4) + if(jyp_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,7*ng+jg)*savg(5*ll4f+ll4x+jyp_p,jg),4) + enddo + enddo + ! node i+1 + LLR(:ng,:8*ng)=matmul(fd(mat(i+1,j),1,:ng,:ng),Lx(:ng,:8*ng,i+1,j)) + do ig=1,ng + A(:ng,ig)=A(:ng,ig)+real(-LLR(:ng,ng+ig),4) + do jg=1,ng + A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,jg)*savg(ind2,jg),4) + if(jym /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,2*ng+jg)*savg(5*ll4f+ll4x+jym,jg),4) + if(jym_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,3*ng+jg)*savg(5*ll4f+ll4x+jym_p,jg),4) + if(jym_pp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,4*ng+jg)*savg(5*ll4f+ll4x+jym_pp,jg),4) + if(jyp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,5*ng+jg)*savg(5*ll4f+ll4x+jyp,jg),4) + if(jyp_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,6*ng+jg)*savg(5*ll4f+ll4x+jyp_p,jg),4) + if(jyp_pp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,7*ng+jg)*savg(5*ll4f+ll4x+jyp_pp,jg),4) + enddo + enddo + call ALSB(ng,1,A,ier,ng) + if(ier /= 0) call XABORT('NSSANM2: singular matrix.(2)') + if(jxp /= 0) savg(5*ll4f+jxp,:ng)=A(:ng,ng+1) + enddo + ! + ! one-node relation at right + ind1=idl(nxmax,j) + if(ind1 == 0) call XABORT('NSSANM2: invalid idl index.(4)') + iqf2=iqfr(2,nxmax,j) + jxm=kn(1,nxmax,j) ; jxp=kn(2,nxmax,j) ; jym=kn(3,nxmax,j) ; jyp=kn(4,nxmax,j) + jym_m=0 ; jyp_m=0 + if(nxmax > 1) then + jym_m=kn(3,nxmax-1,j) ; jyp_m=kn(4,nxmax-1,j) + endif + A(:ng,:ng+1)=0.0 + LLR(:ng,:8*ng)=0.0 + if((iqf2 > 0).or.(iqf2 == -1)) then + ! physical albedo + Lambda(:ng,:ng)=0.0 + do ig=1,ng + Lambda(ig,ig)=qfr(2,nxmax,j,ig) + enddo + LLR(:ng,:8*ng)=matmul(Lambda(:ng,:ng),Rx(:ng,:8*ng,nxmax,j)) + A(:ng,:ng)=real(LLR(:ng,ng+1:2*ng),4) + do ig=1,ng + A(ig,ig)=-1.0+A(ig,ig) + enddo + else if(iqf2 == -2) then + ! zero net current + do ig=1,ng + A(ig,ig)=1.0 + enddo + else if(iqf2 == -3) then + ! zero flux + LLR(:ng,:8*ng)=Rx(:ng,:8*ng,nxmax,j) + A(:ng,:ng)=real(LLR(:ng,ng+1:2*ng),4) + else if(iqf2 == -4) then + call XABORT('NSSANM2: SYME boundary condition is not supported.(2)') + else + call XABORT('NSSANM2: illegal right X-boundary condition.') + endif + if(iqf2 /= -2) then + A(:ng,ng+1)=real(matmul(-LLR(:ng,:ng),savg(ind1,:ng)),4) + do ig=1,ng + do jg=1,ng + if(jym_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,2*ng+jg)*savg(5*ll4f+ll4x+jym_m,jg),4) + if(jym /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,3*ng+jg)*savg(5*ll4f+ll4x+jym,jg),4) + if(jyp_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,5*ng+jg)*savg(5*ll4f+ll4x+jyp_m,jg),4) + if(jyp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,6*ng+jg)*savg(5*ll4f+ll4x+jyp,jg),4) + enddo + enddo + endif + call ALSB(ng,1,A,ier,ng) + if(ier /= 0) call XABORT('NSSANM2: singular matrix.(3)') + if(jxp /= 0) savg(5*ll4f+jxp,:ng)=A(:ng,ng+1) + enddo + !---- + ! one- and two-node relations along Y axis + !---- + do i=1,nx + nymin=1 + do j=1,ny + if(mat(i,j) > 0) exit + nymin=j+1 + enddo + if(nymin > ny) cycle + nymax=ny + do j=ny,1,-1 + if(mat(i,j) > 0) exit + nymax=j-1 + enddo + ! one-node relation at left + ind1=idl(i,nymin) + if(ind1 == 0) call XABORT('NSSANM2: invalid idl index.(5)') + iqf3=iqfr(3,i,nymin) + jxm=kn(1,i,nymin) ; jxp=kn(2,i,nymin) ; jym=kn(3,i,nymin) ; jyp=kn(4,i,nymin) + jxm_p=0 ; jxp_p=0 + if(nymin < ny) then + jxm_p=kn(1,i,nymin+1) ; jxp_p=kn(2,i,nymin+1) + endif + A(:ng,:ng+1)=0.0 + LLR(:ng,:8*ng)=0.0 + if((iqf3 > 0).or.(iqf3 == -1)) then + ! physical albedo + Lambda(:ng,:ng)=0.0 + do ig=1,ng + Lambda(ig,ig)=qfr(3,i,nymin,ig) + enddo + LLR(:ng,:8*ng)=matmul(Lambda(:ng,:ng),Ly(:ng,:8*ng,i,nymin)) + A(:ng,:ng)=-real(LLR(:ng,ng+1:2*ng),4) + do ig=1,ng + A(ig,ig)=-1.0+A(ig,ig) + enddo + else if(iqf3 == -2) then + ! zero net current + do ig=1,ng + A(ig,ig)=1.0 + enddo + else if(iqf3 == -3) then + ! zero flux + LLR(:ng,:8*ng)=Ly(:ng,:8*ng,i,nymin) + A(:ng,:ng)=real(-LLR(:ng,ng+1:2*ng),4) + else if(iqf3 == -4) then + call XABORT('NSSANM2: SYME boundary condition is not supported.(3)') + else + call XABORT('NSSANM2: illegal left Y-boundary condition.') + endif + if(iqf3 /= -2) then + A(:ng,ng+1)=real(matmul(LLR(:ng,:ng),savg(ind1,:ng)),4) + do ig=1,ng + do jg=1,ng + if(jxm /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,3*ng+jg)*savg(5*ll4f+jxm,jg),4) + if(jxm_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,4*ng+jg)*savg(5*ll4f+jxm_p,jg),4) + if(jxp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,6*ng+jg)*savg(5*ll4f+jxp,jg),4) + if(jxp_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,7*ng+jg)*savg(5*ll4f+jxp_p,jg),4) + enddo + enddo + endif + call ALSB(ng,1,A,ier,ng) + if(ier /= 0) call XABORT('NSSANM2: singular matrix.(4)') + if(jym /= 0) savg(5*ll4f+ll4x+jym,:ng)=A(:ng,ng+1) + ! + ! two-node relations + do j=nymin,nymax-1 + ind1=idl(i,j) + if(ind1 == 0) call XABORT('NSSANM2: invalid idl index.(6)') + ind2=idl(i,j+1) + if(ind2 == 0) call XABORT('NSSANM2: invalid idl index.(7)') + if(kn(3,i,j+1) /= kn(4,i,j)) call XABORT('NSSANM2: invalid kn index.(2)') + if(iqfr(4,i,j) /= 0) call XABORT('NSSANM2: invalid iqfr index.(3)') + if(iqfr(3,i,j+1) /= 0) call XABORT('NSSANM2: invalid iqfr index.(4)') + jxm=kn(1,i,j) ; jxp=kn(2,i,j) ; jym=kn(3,i,j) ; jyp=kn(4,i,j) + jxm_m=0 ; jxp_m=0 ; jxm_pp=0 ; jxp_pp=0 + if((j == 1).and.(iqfr(3,i,1) == -2)) then + jxm_m=kn(1,i,1) ; jxp_m=kn(2,i,1) + else if(j > 1) then + jxm_m=kn(1,i,j-1) ; jxp_m=kn(2,i,j-1) + endif + jxm_p=kn(1,i,j+1) ; jxp_p=kn(2,i,j+1) + if((j == ny-1).and.(iqfr(4,i,ny) == -2)) then + jxm_pp=kn(1,i,ny) ; jxp_pp=kn(2,i,ny) + else if(j < ny-1) then + jxm_pp=kn(1,i,j+2) ; jxp_pp=kn(2,i,j+2) + endif + ! + A(:ng,:ng+1)=0.0 + ! node j + LLR(:ng,:8*ng)=matmul(fd(mat(i,j),4,:ng,:ng),Ry(:ng,:8*ng,i,j)) + do ig=1,ng + A(:ng,ig)=A(:ng,ig)+real(LLR(:ng,ng+ig),4) + do jg=1,ng + A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,jg)*savg(ind1,jg),4) + if(jxm_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,2*ng+jg)*savg(5*ll4f+jxm_m,jg),4) + if(jxm /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,3*ng+jg)*savg(5*ll4f+jxm,jg),4) + if(jxm_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,4*ng+jg)*savg(5*ll4f+jxm_p,jg),4) + if(jxp_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,5*ng+jg)*savg(5*ll4f+jxp_m,jg),4) + if(jxp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,6*ng+jg)*savg(5*ll4f+jxp,jg),4) + if(jxp_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,7*ng+jg)*savg(5*ll4f+jxp_p,jg),4) + enddo + enddo + ! node j+1 + LLR(:ng,:8*ng)=matmul(fd(mat(i,j+1),3,:ng,:ng),Ly(:ng,:8*ng,i,j+1)) + do ig=1,ng + A(:ng,ig)=A(:ng,ig)+real(-LLR(:ng,ng+ig),4) + do jg=1,ng + A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,jg)*savg(ind2,jg),4) + if(jxm /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,2*ng+jg)*savg(5*ll4f+jxm,jg),4) + if(jxm_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,3*ng+jg)*savg(5*ll4f+jxm_p,jg),4) + if(jxm_pp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,4*ng+jg)*savg(5*ll4f+jxm_pp,jg),4) + if(jxp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,5*ng+jg)*savg(5*ll4f+jxp,jg),4) + if(jxp_p /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,6*ng+jg)*savg(5*ll4f+jxp_p,jg),4) + if(jxp_pp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(LLR(ig,7*ng+jg)*savg(5*ll4f+jxp_pp,jg),4) + enddo + enddo + call ALSB(ng,1,A,ier,ng) + if(ier /= 0) call XABORT('NSSANM2: singular matrix.(5)') + if(jyp /= 0) savg(5*ll4f+ll4x+jyp,:ng)=A(:ng,ng+1) + enddo + ! + ! one-node relation at right + ind1=idl(i,nymax) + if(ind1 == 0) call XABORT('NSSANM2: invalid idl index.(8)') + iqf4=iqfr(4,i,nymax) + jxm=kn(1,i,nymax) ; jxp=kn(2,i,nymax) ; jym=kn(3,i,nymax) ; jyp=kn(4,i,nymax) + jxm_m=0 ; jxp_m=0 + if(nymax > 1) then + jxm_m=kn(1,i,nymax-1) ; jxp_m=kn(2,i,nymax-1) + endif + A(:ng,:ng+1)=0.0 + LLR(:ng,:8*ng)=0.0 + if((iqf4 > 0).or.(iqf4 == -1)) then + ! physical albedo + Lambda(:ng,:ng)=0.0 + do ig=1,ng + Lambda(ig,ig)=qfr(4,i,nymax,ig) + enddo + LLR(:ng,:8*ng)=matmul(Lambda(:ng,:ng),Ry(:ng,:8*ng,i,nymax)) + A(:ng,:ng)=real(LLR(:ng,ng+1:2*ng),4) + do ig=1,ng + A(ig,ig)=-1.0+A(ig,ig) + enddo + else if(iqf4 == -2) then + ! zero net current + do ig=1,ng + A(ig,ig)=1.0 + enddo + else if(iqf4 == -3) then + ! zero flux + LLR(:ng,:8*ng)=Ry(:ng,:8*ng,i,nymax) + A(:ng,:ng)=real(LLR(:ng,ng+1:2*ng),4) + else if(iqf4 == -4) then + call XABORT('NSSANM2: SYME boundary condition is not supported.(4)') + else + call XABORT('NSSANM2: illegal right Y-boundary condition.') + endif + if(iqf4 /= -2) then + A(:ng,ng+1)=real(matmul(-LLR(:ng,:ng),savg(ind1,:ng)),4) + do ig=1,ng + do jg=1,ng + if(jxm_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,2*ng+jg)*savg(5*ll4f+jxm_m,jg),4) + if(jxm /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,3*ng+jg)*savg(5*ll4f+jxm,jg),4) + if(jxp_m /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,5*ng+jg)*savg(5*ll4f+jxp_m,jg),4) + if(jxp /= 0) A(ig,ng+1)=A(ig,ng+1)+real(-LLR(ig,6*ng+jg)*savg(5*ll4f+jxp,jg),4) + enddo + enddo + endif + call ALSB(ng,1,A,ier,ng) + if(ier /= 0) call XABORT('NSSANM2: singular matrix.(6)') + if(jyp /= 0) savg(5*ll4f+ll4x+jyp,:ng)=A(:ng,ng+1) + enddo + !---- + ! end of transverse current iterations + !---- + enddo + !---- + ! compute boundary fluxes + !---- + do j=1,ny + do i=1,nx + ind1=idl(i,j) + if(ind1 == 0) cycle + jxm=kn(1,i,j) ; jxp=kn(2,i,j) ; jym=kn(3,i,j) ; jyp=kn(4,i,j) + jym_m=0 ; jyp_m=0 ; jym_p=0 ; jyp_p=0 + if((i == 1).and.(iqfr(1,1,j) == -2)) then + jym_m=kn(3,1,j) ; jyp_m=kn(4,1,j) + else if(i > 1) then + jym_m=kn(3,i-1,j) ; jyp_m=kn(4,i-1,j) + endif + if((i == nx).and.(iqfr(2,nx,j) == -2)) then + jym_p=kn(3,nx,j) ; jyp_p=kn(4,nx,j) + else if(i < nx) then + jym_p=kn(3,i+1,j) ; jyp_p=kn(4,i+1,j) + endif + ! x- relations + savg(ll4f+ind1,:ng)=real(matmul(Lx(:ng,:ng,i,j),savg(ind1,:ng)),4) + if(jxm /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,ng+1:2*ng,i,j),savg(5*ll4f+jxm,:ng)),4) + if(jym_m /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,2*ng+1:3*ng,i,j),savg(5*ll4f+ll4x+jym_m,:ng)),4) + if(jym /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,3*ng+1:4*ng,i,j),savg(5*ll4f+ll4x+jym,:ng)),4) + if(jym_p /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,4*ng+1:5*ng,i,j),savg(5*ll4f+ll4x+jym_p,:ng)),4) + if(jyp_m /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,5*ng+1:6*ng,i,j),savg(5*ll4f+ll4x+jyp_m,:ng)),4) + if(jyp /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,6*ng+1:7*ng,i,j),savg(5*ll4f+ll4x+jyp,:ng)),4) + if(jyp_p /= 0) savg(ll4f+ind1,:ng)=savg(ll4f+ind1,:ng)+ & + & real(matmul(Lx(:ng,7*ng+1:8*ng,i,j),savg(5*ll4f+ll4x+jyp_p,:ng)),4) + ! + ! x+ relations + savg(2*ll4f+ind1,:ng)=real(matmul(Rx(:ng,:ng,i,j),savg(ind1,:ng)),4) + if(jxp /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,ng+1:2*ng,i,j),savg(5*ll4f+jxp,:ng)),4) + if(jym_m /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,2*ng+1:3*ng,i,j),savg(5*ll4f+ll4x+jym_m,:ng)),4) + if(jym /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,3*ng+1:4*ng,i,j),savg(5*ll4f+ll4x+jym,:ng)),4) + if(jym_p /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,4*ng+1:5*ng,i,j),savg(5*ll4f+ll4x+jym_p,:ng)),4) + if(jyp_m /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,5*ng+1:6*ng,i,j),savg(5*ll4f+ll4x+jyp_m,:ng)),4) + if(jyp /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,6*ng+1:7*ng,i,j),savg(5*ll4f+ll4x+jyp,:ng)),4) + if(jyp_p /= 0) savg(2*ll4f+ind1,:ng)=savg(2*ll4f+ind1,:ng)+ & + & real(matmul(Rx(:ng,7*ng+1:8*ng,i,j),savg(5*ll4f+ll4x+jyp_p,:ng)),4) + ! + jxm_m=0 ; jxp_m=0 ; jxm_p=0 ; jxp_p=0 + jxm=kn(1,i,j) ; jxp=kn(2,i,j) ; jym=kn(3,i,j) ; jyp=kn(4,i,j) + if((j == 1).and.(iqfr(3,i,1) == -2)) then + jxm_m=kn(1,i,1) ; jxp_m=kn(2,i,1) + else if(j > 1) then + jxm_m=kn(1,i,j-1) ; jxp_m=kn(2,i,j-1) + endif + if((j == ny).and.(iqfr(4,i,ny) == -2)) then + jxm_p=kn(1,i,ny) ; jxp_p=kn(2,i,ny) + else if(j < ny) then + jxm_p=kn(1,i,j+1) ; jxp_p=kn(2,i,j+1) + endif + ! y- relations + savg(3*ll4f+ind1,:ng)=real(matmul(Ly(:ng,:ng,i,j),savg(ind1,:ng)),4) + if(jym /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,ng+1:2*ng,i,j),savg(5*ll4f+ll4x+jym,:ng)),4) + if(jxm_m /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,2*ng+1:3*ng,i,j),savg(5*ll4f+jxm_m,:ng)),4) + if(jxm /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,3*ng+1:4*ng,i,j),savg(5*ll4f+jxm,:ng)),4) + if(jxm_p /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,4*ng+1:5*ng,i,j),savg(5*ll4f+jxm_p,:ng)),4) + if(jxp_m /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,5*ng+1:6*ng,i,j),savg(5*ll4f+jxp_m,:ng)),4) + if(jxp /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,6*ng+1:7*ng,i,j),savg(5*ll4f+jxp,:ng)),4) + if(jxp_p /= 0) savg(3*ll4f+ind1,:ng)=savg(3*ll4f+ind1,:ng)+ & + & real(matmul(Ly(:ng,7*ng+1:8*ng,i,j),savg(5*ll4f+jxp_p,:ng)),4) + ! + ! y+ relations + savg(4*ll4f+ind1,:ng)=real(matmul(Ry(:ng,:ng,i,j),savg(ind1,:ng)),4) + if(jyp /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,ng+1:2*ng,i,j),savg(5*ll4f+ll4x+jyp,:ng)),4) + if(jxm_m /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,2*ng+1:3*ng,i,j),savg(5*ll4f+jxm_m,:ng)),4) + if(jxm /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,3*ng+1:4*ng,i,j),savg(5*ll4f+jxm,:ng)),4) + if(jxm_p /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,4*ng+1:5*ng,i,j),savg(5*ll4f+jxm_p,:ng)),4) + if(jxp_m /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,5*ng+1:6*ng,i,j),savg(5*ll4f+jxp_m,:ng)),4) + if(jxp /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,6*ng+1:7*ng,i,j),savg(5*ll4f+jxp,:ng)),4) + if(jxp_p /= 0) savg(4*ll4f+ind1,:ng)=savg(4*ll4f+ind1,:ng)+ & + & real(matmul(Ry(:ng,7*ng+1:8*ng,i,j),savg(5*ll4f+jxp_p,:ng)),4) + enddo + enddo + !---- + ! scratch storage deallocation + !---- + deallocate(Ry,Ly,Rx,Lx) + deallocate(work5,work4,work3,work2,work1) + deallocate(LLR,Lambda,A) +end subroutine NSSANM2 -- cgit v1.2.3