summaryrefslogtreecommitdiff
path: root/Dragon/src/DMASOU.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/DMASOU.f')
-rw-r--r--Dragon/src/DMASOU.f555
1 files changed, 555 insertions, 0 deletions
diff --git a/Dragon/src/DMASOU.f b/Dragon/src/DMASOU.f
new file mode 100644
index 0000000..495d32d
--- /dev/null
+++ b/Dragon/src/DMASOU.f
@@ -0,0 +1,555 @@
+*DECK DMASOU
+ SUBROUTINE DMASOU(IPRINT,IPDMA,IPMAC,IPFLX,NG,NREG,NMIL,NL,
+ 1 NDEL,NED,NAMEAD,NUN,NMERGE,NGCOND,NCST,IMERGE,INDGRP,MATCOD,
+ 2 KEYFLX,VOL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the GPT sources corresponding to the gradient of a macrolib.
+*
+*Copyright:
+* Copyright (C) 2008 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
+* IPRINT print parameter
+* IPDMA pointer to the DMA data structure.
+* IPMAC pointer to the macrolib structure.
+* IPFLX pointer to the multigroup flux.
+* NG number of energy groups.
+* NREG number of regions.
+* NMIL number of material mixtures.
+* NL number of Legendre orders.
+* NDEL number of delayed precursors.
+* NED number of extra edit vectors.
+* NAMEAD names of these extra edits.
+* NUN number of unknowns per energy group.
+* NMERGE number of merged regions.
+* NGCOND number of condensed energy groups.
+* NCST number of DMA fixed sources.
+* IMERGE merging indices.
+* INDGRP condensation indices.
+* MATCOD material mixture indices per region.
+* KEYFLX position of averaged fluxes in unknown vector.
+* VOL volumes.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPDMA,IPMAC,IPFLX
+ INTEGER IPRINT,NG,NREG,NMIL,NL,NDEL,NED,NAMEAD(2,NED),NMERGE,
+ 1 NGCOND,NCST,IMERGE(NREG),INDGRP(NG),MATCOD(NREG),KEYFLX(NREG)
+ REAL VOL(NREG)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (NSECT=4,EPSMAX=1.0E-7)
+ TYPE(C_PTR) JPFLX,JPDMA,KPDMA,JPMAC,KPMAC
+ CHARACTER TEXT12*12,TEXB12*12,HSECT(NSECT)*12,CM*2
+ DOUBLE PRECISION WW,SUM,FUNC,ZN
+*----
+* DATA STATEMENTS
+*----
+ DATA HSECT/'NTOT0','SIGS00','N2N','N3N'/
+*----
+* ALLOCATABLE ARRAYS
+*----
+ TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) ::IKEP
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJS00,NJJS00,IPOS00
+ REAL, ALLOCATABLE, DIMENSION(:) :: FLUX,SIGT,CHI,EPS,SCAT
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XSSNN
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: GAR1,GAR2
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) ::GAR3
+ REAL, POINTER, DIMENSION(:) :: SUNK
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(IKEP(NG))
+ ALLOCATE(FLUX(NUN),SIGT(NMIL),CHI(NMIL),XSSNN(NMIL,NG,NG),
+ 1 EPS(NCST))
+ ALLOCATE(GAR1(NMERGE,NGCOND),GAR2(NMERGE,NGCOND),
+ 1 GAR3(NMERGE,NGCOND,NGCOND))
+*
+ IOF=0
+ EPS(:NCST)=0.0
+ JPFLX=LCMGID(IPFLX,'FLUX')
+ JPDMA=LCMLID(IPDMA,'ASOUR',NCST)
+*----
+* NWT0 INFORMATION
+*----
+ SUM=0.0D0
+ GAR2(:NMERGE,:NGCOND)=0.0D0
+ DO IG=1,NG
+ IGCND=INDGRP(IG)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IMERG=IMERGE(IR)
+ WW=FLUX(KEYFLX(IR))*VOL(IR)
+ SUM=SUM+WW
+ IF((IGCND.NE.0).AND.(IMERG.NE.0)) THEN
+ GAR2(IMERG,IGCND)=GAR2(IMERG,IGCND)+WW
+ ENDIF
+ ENDDO
+ ENDDO
+ DO IGCND=1,NGCOND
+ DO IMERG=1,NMERGE
+ IOF=IOF+1
+ IF(IOF.GT.NCST) CALL XABORT('DMASOU: NCST OVERFLOW(1).')
+ DO IG=1,NG
+ IKEP(IG)=LCMARA(NUN)
+ CALL C_F_POINTER(IKEP(IG),SUNK,(/ NUN /))
+ SUNK(:NUN)=0.0
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IUNK=KEYFLX(IR)
+ FUNC=VOL(IR)*GAR2(IMERG,IGCND)/SUM
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND))
+ 1 THEN
+ SUNK(IUNK)=REAL(FUNC/GAR2(IMERG,IGCND))
+ ENDIF
+ SUNK(IUNK)=REAL(SUNK(IUNK)-FUNC/SUM)
+ ZN=SUM/FUNC
+ EPS(IOF)=MAX(EPS(IOF),ABS(SUNK(IUNK)*REAL(ZN)))
+ ENDDO
+ ENDDO
+ IF(EPS(IOF).GT.EPSMAX) THEN
+ KPDMA=LCMLIL(JPDMA,IOF,NG)
+ DO IG=1,NG
+ CALL LCMPPL(KPDMA,IG,NUN,2,IKEP(IG))
+ ENDDO
+ ELSE
+ DO IG=1,NG
+ CALL LCMDRD(IKEP(IG))
+ ENDDO
+ ENDIF
+ ENDDO
+ ENDDO
+*----
+* SET OF NSECT BASIC CROSS SECTIONS
+*----
+ JPMAC=LCMGID(IPMAC,'GROUP')
+ DO ISECT=1,NSECT
+ TEXT12=HSECT(ISECT)
+ GAR1(:NMERGE,:NGCOND)=0.0D0
+ GAR2(:NMERGE,:NGCOND)=0.0D0
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ IGCND=INDGRP(IG)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IMERG=IMERGE(IR)
+ WW=FLUX(KEYFLX(IR))*VOL(IR)
+ IF((IGCND.NE.0).AND.(IMERG.NE.0)) THEN
+ IBM=MATCOD(IR)
+ GAR1(IMERG,IGCND)=GAR1(IMERG,IGCND)+WW
+ IF(IBM.GT.0) THEN
+ GAR2(IMERG,IGCND)=GAR2(IMERG,IGCND)+SIGT(IBM)*WW
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ DO IGCND=1,NGCOND
+ DO IMERG=1,NMERGE
+ IOF=IOF+1
+ IF(IOF.GT.NCST) CALL XABORT('DMASOU: NCST OVERFLOW(2).')
+ IF(GAR2(IMERG,IGCND).NE.0.0) THEN
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IKEP(IG)=LCMARA(NUN)
+ CALL C_F_POINTER(IKEP(IG),SUNK,(/ NUN /))
+ SUNK(:NUN)=0.0
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND))
+ 1 THEN
+ IUNK=KEYFLX(IR)
+ IBM=MATCOD(IR)
+ FUNC=VOL(IR)*GAR2(IMERG,IGCND)/GAR1(IMERG,IGCND)
+ IF(IBM.EQ.0) THEN
+ SUNK(IUNK)=REAL(-FUNC/GAR1(IMERG,IGCND))
+ ELSE
+ SUNK(IUNK)=REAL(FUNC*(SIGT(IBM)/
+ 1 GAR2(IMERG,IGCND)-1.0D0/GAR1(IMERG,IGCND)))
+ ENDIF
+ ZN=SUM/FUNC
+ EPS(IOF)=MAX(EPS(IOF),ABS(SUNK(IUNK)*REAL(ZN)))
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ IF(EPS(IOF).GT.EPSMAX) THEN
+ KPDMA=LCMLIL(JPDMA,IOF,NG)
+ DO IG=1,NG
+ CALL LCMPPL(KPDMA,IG,NUN,2,IKEP(IG))
+ ENDDO
+ ELSE
+ DO IG=1,NG
+ CALL LCMDRD(IKEP(IG))
+ ENDDO
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+*----
+* SCATTERING CROSS SECTION INFORMATION
+*----
+ ALLOCATE(IJJS00(NMIL),NJJS00(NMIL),IPOS00(NMIL))
+ DO IL=1,NL
+ WRITE(CM,'(I2.2)') IL-1
+ XSSNN(:NMIL,:NG,:NG)=0.0
+ DO JG=1,NG
+ KPMAC=LCMGIL(JPMAC,JG)
+ CALL LCMGET(KPMAC,'IJJS'//CM,IJJS00)
+ CALL LCMGET(KPMAC,'NJJS'//CM,NJJS00)
+ CALL LCMGET(KPMAC,'IPOS'//CM,IPOS00)
+ IMAX=0
+ DO IBM=1,NMIL
+ IMAX=IMAX+NJJS00(IBM)
+ ENDDO
+ ALLOCATE(SCAT(IMAX))
+ CALL LCMGET(KPMAC,'SCAT'//CM,SCAT)
+ DO IBM=1,NMIL
+ IPOS=IPOS00(IBM)
+ IG=IJJS00(IBM)
+ IENBR=NJJS00(IBM)
+ DO WHILE (IENBR.GE.1)
+ XSSNN(IBM,JG,IG)=SCAT(IPOS) ! JG <-- IG
+ IPOS=IPOS+1
+ IENBR=IENBR-1
+ IG=IG-1
+ ENDDO
+ ENDDO
+ DEALLOCATE(SCAT)
+ ENDDO
+ GAR1(:NMERGE,:NGCOND)=0.0D0
+ GAR3(:NMERGE,:NGCOND,:NGCOND)=0.0D0
+ DO IG=1,NG
+ IGCND=INDGRP(IG)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IMERG=IMERGE(IR)
+ WW=FLUX(KEYFLX(IR))*VOL(IR)
+ IF((IGCND.NE.0).AND.(IMERG.NE.0)) THEN
+ IBM=MATCOD(IR)
+ GAR1(IMERG,IGCND)=GAR1(IMERG,IGCND)+WW
+ IF(IBM.GT.0) THEN
+ DO JG=1,NG
+ JGCND=INDGRP(JG)
+ GAR3(IMERG,JGCND,IGCND)=GAR3(IMERG,JGCND,IGCND)
+ 1 +XSSNN(IBM,JG,IG)*WW
+ ENDDO
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ DO JGCND=1,NGCOND
+ DO IGCND=1,NGCOND
+ DO IMERG=1,NMERGE
+ IOF=IOF+1
+ IF(IOF.GT.NCST) CALL XABORT('DMASOU: NCST OVERFLOW(3).')
+ IF(GAR3(IMERG,JGCND,IGCND).NE.0.0) THEN
+ DO IG=1,NG
+ IKEP(IG)=LCMARA(NUN)
+ CALL C_F_POINTER(IKEP(IG),SUNK,(/ NUN /))
+ SUNK(:NUN)=0.0
+ DO JG=1,NG
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND)
+ 1 .AND.(INDGRP(JG).EQ.JGCND)) THEN
+ IBM=MATCOD(IR)
+ IF(IBM.NE.0) THEN
+ IUNK=KEYFLX(IR)
+ FUNC=VOL(IR)*GAR3(IMERG,JGCND,IGCND)/
+ 1 GAR1(IMERG,IGCND)
+ SUNK(IUNK)=REAL(SUNK(IUNK)+FUNC*
+ 1 XSSNN(IBM,JG,IG)/GAR3(IMERG,JGCND,IGCND))
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND))
+ 1 THEN
+ IUNK=KEYFLX(IR)
+ FUNC=VOL(IR)*GAR3(IMERG,JGCND,IGCND)/
+ 1 GAR1(IMERG,IGCND)
+ SUNK(IUNK)=SUNK(IUNK)-REAL(FUNC/GAR1(IMERG,IGCND))
+ ZN=SUM/FUNC
+ EPS(IOF)=MAX(EPS(IOF),ABS(SUNK(IUNK)*REAL(ZN)))
+ ENDIF
+ ENDDO
+ ENDDO
+ IF(EPS(IOF).GT.EPSMAX) THEN
+ KPDMA=LCMLIL(JPDMA,IOF,NG)
+ DO IG=1,NG
+ CALL LCMPPL(KPDMA,IG,NUN,2,IKEP(IG))
+ ENDDO
+ ELSE
+ DO IG=1,NG
+ CALL LCMDRD(IKEP(IG))
+ ENDDO
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ DEALLOCATE(IPOS00,NJJS00,IJJS00)
+*----
+* FISSION INFORMATION
+*----
+ DO IDEL=1,1+NDEL
+ IF(IDEL.EQ.1) THEN
+ TEXT12='NUSIGF'
+ TEXB12='CHI'
+ ELSE
+ WRITE(TEXT12,'(6HNUSIGF,I2.2)') IDEL-1
+ WRITE(TEXB12,'(3HCHI,I2.2)') IDEL-1
+ ENDIF
+ GAR1(:NMERGE,:NGCOND)=0.0D0
+ GAR2(:NMERGE,:NGCOND)=0.0D0
+ GAR3(:NMERGE,:NGCOND,1)=0.0D0
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ CALL LCMGET(KPMAC,TEXB12,CHI)
+ IGCND=INDGRP(IG)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IMERG=IMERGE(IR)
+ WW=FLUX(KEYFLX(IR))*VOL(IR)
+ IF((IGCND.NE.0).AND.(IMERG.NE.0)) THEN
+ IBM=MATCOD(IR)
+ GAR1(IMERG,IGCND)=GAR1(IMERG,IGCND)+WW
+ IF(IBM.GT.0) THEN
+ GAR2(IMERG,IGCND)=GAR2(IMERG,IGCND)+SIGT(IBM)*WW
+ GAR3(IMERG,IGCND,1)=GAR3(IMERG,IGCND,1)+CHI(IBM)*
+ 1 SIGT(IBM)*WW
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ DO IGCND=1,NGCOND
+ DO IMERG=1,NMERGE
+ IOF=IOF+1
+ IF(IOF.GT.NCST) CALL XABORT('DMASOU: NCST OVERFLOW(4).')
+ IF(GAR2(IMERG,IGCND).NE.0.0) THEN
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IKEP(IG)=LCMARA(NUN)
+ CALL C_F_POINTER(IKEP(IG),SUNK,(/ NUN /))
+ SUNK(:NUN)=0.0
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND))
+ 1 THEN
+ IUNK=KEYFLX(IR)
+ IBM=MATCOD(IR)
+ FUNC=VOL(IR)*GAR2(IMERG,IGCND)/GAR1(IMERG,IGCND)
+ IF(IBM.EQ.0) THEN
+ SUNK(IUNK)=REAL(-FUNC/GAR1(IMERG,IGCND))
+ ELSE
+ SUNK(IUNK)=REAL(FUNC*(SIGT(IBM)/
+ 1 GAR2(IMERG,IGCND)-1.0D0/GAR1(IMERG,IGCND)))
+ ENDIF
+ ZN=SUM/FUNC
+ EPS(IOF)=MAX(EPS(IOF),ABS(SUNK(IUNK)*REAL(ZN)))
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ IF(EPS(IOF).GT.EPSMAX) THEN
+ KPDMA=LCMLIL(JPDMA,IOF,NG)
+ DO IG=1,NG
+ CALL LCMPPL(KPDMA,IG,NUN,2,IKEP(IG))
+ ENDDO
+ ELSE
+ DO IG=1,NG
+ CALL LCMDRD(IKEP(IG))
+ ENDDO
+ ENDIF
+ ENDIF
+ IOF=IOF+1
+ IF(IOF.GT.NCST) CALL XABORT('DMASOU: NCST OVERFLOW(5).')
+ IF(GAR3(IMERG,IGCND,1).NE.0.0) THEN
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IKEP(IG)=LCMARA(NUN)
+ CALL C_F_POINTER(IKEP(IG),SUNK,(/ NUN /))
+ SUNK(:NUN)=0.0
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ CALL LCMGET(KPMAC,TEXB12,CHI)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND))
+ 1 THEN
+ IBM=MATCOD(IR)
+ IF(IBM.NE.0) THEN
+ IUNK=KEYFLX(IR)
+ FUNC=VOL(IR)*SIGT(IBM)*GAR3(IMERG,IGCND,1)/
+ 1 GAR2(IMERG,IGCND)
+ SUNK(IUNK)=REAL(FUNC*(CHI(IBM)/
+ 1 GAR3(IMERG,IGCND,1)-1.0D0/GAR2(IMERG,IGCND)))
+ ZN=SUM/FUNC
+ EPS(IOF)=MAX(EPS(IOF),ABS(SUNK(IUNK)*REAL(ZN)))
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ IF(EPS(IOF).GT.EPSMAX) THEN
+ KPDMA=LCMLIL(JPDMA,IOF,NG)
+ DO IG=1,NG
+ CALL LCMPPL(KPDMA,IG,NUN,2,IKEP(IG))
+ ENDDO
+ ELSE
+ DO IG=1,NG
+ CALL LCMDRD(IKEP(IG))
+ ENDDO
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+*----
+* ADDITIONAL CROSS SECTION INFORMATION
+*----
+ DO IED=1,NED
+ WRITE(TEXT12,'(2A4)') NAMEAD(1,IED),NAMEAD(2,IED)
+ GAR1(:NMERGE,:NGCOND)=0.0D0
+ GAR2(:NMERGE,:NGCOND)=0.0D0
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ IGCND=INDGRP(IG)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IMERG=IMERGE(IR)
+ WW=FLUX(KEYFLX(IR))*VOL(IR)
+ IF((IGCND.NE.0).AND.(IMERG.NE.0)) THEN
+ IBM=MATCOD(IR)
+ GAR1(IMERG,IGCND)=GAR1(IMERG,IGCND)+WW
+ IF(IBM.GT.0) THEN
+ GAR2(IMERG,IGCND)=GAR2(IMERG,IGCND)+SIGT(IBM)*WW
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ DO IGCND=1,NGCOND
+ DO IMERG=1,NMERGE
+ IOF=IOF+1
+ IF(IOF.GT.NCST) CALL XABORT('DMASOU: NCST OVERFLOW(6).')
+ IF(GAR2(IMERG,IGCND).NE.0.0) THEN
+ DO IG=1,NG
+ KPMAC=LCMGIL(JPMAC,IG)
+ CALL LCMLEN(KPMAC,TEXT12,ILONG,ITYLCM)
+ IKEP(IG)=LCMARA(NUN)
+ CALL C_F_POINTER(IKEP(IG),SUNK,(/ NUN /))
+ SUNK(:NUN)=0.0
+ IF(ILONG.GT.0) THEN
+ CALL LCMGET(KPMAC,TEXT12,SIGT)
+ DO IR=1,NREG
+ IF(KEYFLX(IR).EQ.0) CYCLE
+ IF((IMERGE(IR).EQ.IMERG).AND.(INDGRP(IG).EQ.IGCND))
+ 1 THEN
+ IUNK=KEYFLX(IR)
+ IBM=MATCOD(IR)
+ FUNC=VOL(IR)*GAR2(IMERG,IGCND)/GAR1(IMERG,IGCND)
+ IF(IBM.EQ.0) THEN
+ SUNK(IUNK)=REAL(-FUNC/GAR1(IMERG,IGCND))
+ ELSE
+ SUNK(IUNK)=REAL(FUNC*(SIGT(IBM)/
+ 1 GAR2(IMERG,IGCND)-1.0D0/GAR1(IMERG,IGCND)))
+ ENDIF
+ ZN=SUM/FUNC
+ EPS(IOF)=MAX(EPS(IOF),ABS(SUNK(IUNK)*REAL(ZN)))
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ IF(EPS(IOF).GT.EPSMAX) THEN
+ KPDMA=LCMLIL(JPDMA,IOF,NG)
+ DO IG=1,NG
+ CALL LCMPPL(KPDMA,IG,NUN,2,IKEP(IG))
+ ENDDO
+ ELSE
+ DO IG=1,NG
+ CALL LCMDRD(IKEP(IG))
+ ENDDO
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+*----
+* CHECK SOURCE ORTHOGONALITY
+*----
+ ALLOCATE(SUNK(NUN))
+ DO IOF=1,NCST
+ CALL LCMLEL(JPDMA,IOF,ILONG,ITYLCM)
+ IF(ILONG.NE.0) THEN
+ KPDMA=LCMGIL(JPDMA,IOF)
+ SUM=0.0D0
+ DO IG=1,NG
+ CALL LCMGDL(KPDMA,IG,SUNK)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IR=1,NREG
+ IUNK=KEYFLX(IR)
+ IF(IUNK.GT.0) SUM=SUM+SUNK(IUNK)*FLUX(IUNK)
+ ENDDO
+ ENDDO
+ IF(IPRINT.GT.0) THEN
+ WRITE(6,'(14H SOURCE INDEX=,I10,14H DOT PRODUCT=,1P,E11.4,
+ 1 19H SOURCE INTENSITY=,E11.4)') IOF,ABS(SUM),EPS(IOF)
+ ENDIF
+ IF(ABS(SUM).GT.1.0E-5) THEN
+ WRITE(TEXT12,'(I10,2X)') IOF
+ CALL XABORT('DMASOU: NON ORTHOGONAL SOURCE (IOF='//
+ 1 TEXT12(:10)//').')
+ ENDIF
+ ENDIF
+ ENDDO
+ DEALLOCATE(SUNK)
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(GAR3,GAR2,GAR1)
+ DEALLOCATE(EPS,XSSNN,CHI,SIGT,FLUX)
+ DEALLOCATE(IKEP)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:15:36 +0000