summaryrefslogtreecommitdiff
path: root/Donjon/src/DRENOU.f
diff options
context:
space:
mode:
Diffstat (limited to 'Donjon/src/DRENOU.f')
-rw-r--r--Donjon/src/DRENOU.f549
1 files changed, 549 insertions, 0 deletions
diff --git a/Donjon/src/DRENOU.f b/Donjon/src/DRENOU.f
new file mode 100644
index 0000000..ea5941f
--- /dev/null
+++ b/Donjon/src/DRENOU.f
@@ -0,0 +1,549 @@
+*DECK DRENOU
+ SUBROUTINE DRENOU(IPRINT,IPGPT,IPMAC1,IPMAC2,IPFLX,IPTRK,IPGRAD,
+ 1 NG,NREG,ITYPE,IELEM,NMIL,NALBP,NUN,NFIS1,NFIS2,ILEAK1,ILEAK2,
+ 2 IDF2,MATCOD,KEYFLX,VOL,LNO,NFUNC,RMSD)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the GPT sources corresponding to the gradient of the RMS
+* absorption distribution. Case with NFUNC individual components to be
+* used with a Newtonian method.
+*
+*Copyright:
+* Copyright (C) 2019 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
+* IPGPT pointer to the L_GPT data structure.
+* IPMAC1 pointer to the actual macrolib structure.
+* IPMAC2 pointer to the reference macrolib structure.
+* IPFLX pointer to the multigroup flux.
+* IPTRK pointer to the tracking object.
+* IPGRAD pointer to the L_OPTIMIZE object.
+* NG number of energy groups.
+* NREG number of regions.
+* NMIL number of material mixtures.
+* NALBP number of physical albedos.
+* NUN number of unknowns per energy group.
+* NFIS1 number of fissile isotopes in actual macrolib.
+* NFIS2 number of fissile isotopes in reference macrolib.
+* ILEAK1 type of leakage calculation in actual macrolib
+* =0: no leakage; =1: homogeneous leakage (Diffon).
+* ILEAK2 type of leakage calculation in reference macrolib.
+* IDF2 ADF type, 0 = none, 1 = Albedo, 2 = FD_B/FD_C/..., 3 = ADF.
+* MATCOD material mixture indices per region.
+* KEYFLX position of averaged fluxes in unknown vector.
+* VOL volumes.
+* LNO flag set to .true. to exit after calculation of RMS.
+*
+*Parameters: output
+* NFUNC number of individual components in the gradient terms.
+* RMSD RMS error on rate distribution.
+*
+*Parameters:
+* ITYPE
+* IELEM
+*
+* Reference:
+* A. Hebert,"Developpement de la methode SPH: Homogeneisation de
+* cellules dans un reseau non uniforme et calcul des parametres de
+* reflecteur," Note CEA-N-2209, Sect. 3.5.1, 1981.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPGPT,IPMAC1,IPMAC2,IPFLX,IPTRK,IPGRAD
+ INTEGER IPRINT,NG,NREG,ITYPE,IELEM,NMIL,NALBP,NUN,NFIS1,NFIS2,
+ > ILEAK1,ILEAK2,IDF2,MATCOD(NREG),KEYFLX(NREG),NFUNC
+ REAL VOL(NREG)
+ DOUBLE PRECISION RMSD
+ LOGICAL LNO
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (NSTATE=40)
+ TYPE(C_PTR) JPMAC1,JPMAC2,KPMAC1,KPMAC2,JPFLX,JPGPT,KPGPT
+ INTEGER ISTATE(NSTATE)
+ DOUBLE PRECISION SOUT1,SOUT2,SOUTOT,AB1TOT,AB2TOT,FI1TOT,FI2TOT,
+ > SUM1,DSUM,DELTA,OUT,SA,SF,ABS2M,OUT2M,AIL,BIL,DEN1,DEN2
+ CHARACTER HSMG*131
+ DOUBLE PRECISION, PARAMETER :: EPS=1.0E-4,EPSL=1.0E-4
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS,IREL,KN,IQFR
+ REAL, ALLOCATABLE, DIMENSION(:) :: GAR,WORK,SUNK,FLUX,QFR,OUTG1,
+ > OUTG2,DIFHOM,DIFF
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: PHI1,PHI2,ABS1,ABS2,NUF1,
+ > NUF2,GAMMA,OUTG2R
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: CHI1,CHI2,RHS1,LHS1,RHS2,
+ > LHS2
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VARV,RHS,CONST,FF
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: SIGA,SIGF,DFF
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(PHI1(NMIL,NG),PHI2(NMIL,NG),ABS1(NMIL,NG),ABS2(NMIL,NG),
+ 1 RHS1(NMIL,NG,NG),LHS1(NMIL,NG,NG),RHS2(NMIL,NG,NG),
+ 2 LHS2(NMIL,NG,NG),CONST(NG),IREL(NG),RHS(NG),GAMMA(NUN,NG),
+ 3 OUTG1(NG),OUTG2(NG),OUTG2R(NG,2),SIGA(NMIL,NG),SIGF(NMIL,NG))
+*----
+* COMPUTE THE ACTUAL AND REFERENCE REACTION RATE MATRICES
+*----
+ CALL LCMGET(IPMAC1,'K-EFFECTIVE',ZKEFF1)
+ CALL LCMGET(IPMAC2,'K-EFFECTIVE',ZKEFF2)
+ IF(IDF2.EQ.1) THEN
+ CALL LCMSIX(IPMAC2,'ADF',1)
+ CALL LCMLEN(IPMAC2,'ALBS00',ILCMLN,ITYLCM)
+ IF(ILCMLN.NE.2*NG) CALL XABORT('DRENOU: WRONG ALBS00 LENGTH.')
+ CALL LCMGET(IPMAC2,'ALBS00',OUTG2R)
+ CALL LCMSIX(IPMAC2,' ',2)
+ ENDIF
+ CALL LCMLEN(IPMAC1,'B2 B1HOM',ILCMLN,ITYLCM)
+ IF(ILCMLN.EQ.1) THEN
+ CALL LCMGET(IPMAC1,'B2 B1HOM',B21)
+ ELSE
+ B21=0.0
+ ENDIF
+ CALL LCMLEN(IPMAC2,'B2 B1HOM',ILCMLN,ITYLCM)
+ IF(ILCMLN.EQ.1) THEN
+ CALL LCMGET(IPMAC2,'B2 B1HOM',B22)
+ ELSE
+ B22=0.0
+ ENDIF
+ IF((ILEAK1.EQ.1).AND.(IPRINT.GT.0)) THEN
+ WRITE(6,'(/22H DRENOU: MACRO B2=,1P,E12.4)') B21
+ ENDIF
+ IF((ILEAK2.EQ.1).AND.(IPRINT.GT.0)) THEN
+ WRITE(6,'(/22H DRENOU: REFERENCE B2=,1P,E12.4)') B22
+ ENDIF
+ RHS1(:NMIL,:NG,:NG)=0.0
+ LHS1(:NMIL,:NG,:NG)=0.0
+ RHS2(:NMIL,:NG,:NG)=0.0
+ LHS2(:NMIL,:NG,:NG)=0.0
+ SIGA(:NMIL,:NG)=0.0D0
+ SIGF(:NMIL,:NG)=0.0D0
+ JPMAC1=LCMGID(IPMAC1,'GROUP')
+ JPMAC2=LCMGID(IPMAC2,'GROUP')
+ ALLOCATE(IJJ(NMIL),NJJ(NMIL),IPOS(NMIL),GAR(NMIL),WORK(NMIL*NG),
+ 1 CHI1(NMIL,NFIS1,NG),NUF1(NMIL,NFIS1),CHI2(NMIL,NFIS2,NG),
+ 2 NUF2(NMIL,NFIS2),DIFHOM(NG),DIFF(NMIL))
+ DO IG=1,NG
+ KPMAC1=LCMGIL(JPMAC1,IG)
+ CALL LCMGET(KPMAC1,'CHI',CHI1(1,1,IG))
+ KPMAC2=LCMGIL(JPMAC2,IG)
+ CALL LCMGET(KPMAC2,'CHI',CHI2(1,1,IG))
+ CALL LCMLEN(KPMAC1,'FLUX-INTG',ILG,ITYLCM)
+ IF(ILG.NE.NMIL) CALL XABORT('DRENOU: MISSING ACTUAL FLUX.')
+ CALL LCMLEN(KPMAC2,'FLUX-INTG',ILG,ITYLCM)
+ IF(ILG.NE.NMIL) CALL XABORT('DRENOU: MISSING REFERENCE FLUX.')
+ CALL LCMGET(KPMAC1,'FLUX-INTG',PHI1(1,IG))
+ CALL LCMGET(KPMAC2,'FLUX-INTG',PHI2(1,IG))
+ ENDDO
+ DO IG=1,NG
+ KPMAC1=LCMGIL(JPMAC1,IG)
+ KPMAC2=LCMGIL(JPMAC2,IG)
+ IF(ILEAK1.EQ.1) THEN
+ CALL LCMLEN(KPMAC1,'DIFF',ILCMLN,ITYLCM)
+ IF(ILCMLN.GT.0) THEN
+ CALL LCMGET(KPMAC1,'DIFF',DIFF)
+ ELSE
+ CALL LCMGET(IPMAC1,'DIFHOMB1HOM',DIFHOM)
+ DO IBM=1,NMIL
+ DIFF(IBM)=DIFHOM(IG)
+ ENDDO
+ ENDIF
+ ELSE
+ DIFF(:NMIL)=0.0
+ ENDIF
+ CALL LCMGET(KPMAC1,'NTOT0',GAR)
+ CALL LCMGET(KPMAC1,'SCAT00',WORK)
+ CALL LCMGET(KPMAC1,'NJJS00',NJJ)
+ CALL LCMGET(KPMAC1,'IJJS00',IJJ)
+ CALL LCMGET(KPMAC1,'IPOS00',IPOS)
+ DO IBM=1,NMIL
+ SIGA(IBM,IG)=SIGA(IBM,IG)+GAR(IBM)
+ IPOSDE=IPOS(IBM)
+ DO JG=IJJ(IBM),IJJ(IBM)-NJJ(IBM)+1,-1
+* IG <-- JG
+ RHS1(IBM,IG,JG)=RHS1(IBM,IG,JG)-WORK(IPOSDE)*PHI1(IBM,JG)
+ SIGA(IBM,JG)=SIGA(IBM,JG)-WORK(IPOSDE)
+ IPOSDE=IPOSDE+1
+ ENDDO
+ RHS1(IBM,IG,IG)=RHS1(IBM,IG,IG)+(GAR(IBM)+B21*DIFF(IBM))*
+ > PHI1(IBM,IG)
+ ENDDO
+ CALL LCMGET(KPMAC1,'NUSIGF',NUF1)
+ DO IBM=1,NMIL
+ DO IFIS=1,NFIS1
+ DO JG=1,NG
+ LHS1(IBM,JG,IG)=LHS1(IBM,JG,IG)+CHI1(IBM,IFIS,JG)*
+ > NUF1(IBM,IFIS)*PHI1(IBM,IG)
+ SIGF(IBM,IG)=SIGF(IBM,IG)+CHI1(IBM,IFIS,JG)*NUF1(IBM,IFIS)
+ ENDDO
+ ENDDO
+ ENDDO
+*
+ IF(ILEAK2.EQ.1) THEN
+ CALL LCMLEN(KPMAC2,'DIFF',ILCMLN,ITYLCM)
+ IF(ILCMLN.GT.0) THEN
+ CALL LCMGET(KPMAC2,'DIFF',DIFF)
+ ELSE
+ CALL LCMGET(IPMAC2,'DIFHOMB1HOM',DIFHOM)
+ DO IBM=1,NMIL
+ DIFF(IBM)=DIFHOM(IG)
+ ENDDO
+ ENDIF
+ ELSE
+ DIFF(:NMIL)=0.0
+ ENDIF
+ CALL LCMGET(KPMAC2,'NTOT0',GAR)
+ CALL LCMGET(KPMAC2,'SCAT00',WORK)
+ CALL LCMGET(KPMAC2,'NJJS00',NJJ)
+ CALL LCMGET(KPMAC2,'IJJS00',IJJ)
+ CALL LCMGET(KPMAC2,'IPOS00',IPOS)
+ DO IBM=1,NMIL
+ IPOSDE=IPOS(IBM)
+ DO JG=IJJ(IBM),IJJ(IBM)-NJJ(IBM)+1,-1
+* IG <-- JG
+ RHS2(IBM,IG,JG)=RHS2(IBM,IG,JG)-WORK(IPOSDE)*PHI2(IBM,JG)
+ IPOSDE=IPOSDE+1
+ ENDDO
+ RHS2(IBM,IG,IG)=RHS2(IBM,IG,IG)+(GAR(IBM)+B22*DIFF(IBM))*
+ > PHI2(IBM,IG)
+ ENDDO
+ CALL LCMGET(KPMAC2,'NUSIGF',NUF2)
+ DO IBM=1,NMIL
+ DO IFIS=1,NFIS2
+ DO JG=1,NG
+ LHS2(IBM,JG,IG)=LHS2(IBM,JG,IG)+CHI2(IBM,IFIS,JG)*
+ > NUF2(IBM,IFIS)*PHI2(IBM,IG)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ DEALLOCATE(DIFF,DIFHOM,NUF2,CHI2,NUF1,CHI1,WORK,GAR,IPOS,NJJ,IJJ)
+*----
+* COMPUTE THE ACTUAL AND REFERENCE ABSORPTION AND FISSION RATES
+*----
+ AB1TOT=0.0D0
+ AB2TOT=0.0D0
+ FI1TOT=0.0D0
+ FI2TOT=0.0D0
+ DO IG=1,NG
+ OUTG1(IG)=0.0
+ OUTG2(IG)=0.0
+ DO IBM=1,NMIL
+ OUTG1(IG)=OUTG1(IG)+SUM(LHS1(IBM,IG,:NG))/ZKEFF1-
+ 1 SUM(RHS1(IBM,IG,:NG))
+ OUTG2(IG)=OUTG2(IG)+SUM(LHS2(IBM,IG,:NG))/ZKEFF2-
+ 1 SUM(RHS2(IBM,IG,:NG))
+ ABS1(IBM,IG)=SUM(RHS1(IBM,:NG,IG))
+ ABS2(IBM,IG)=SUM(RHS2(IBM,:NG,IG))
+ AB1TOT=AB1TOT+ABS1(IBM,IG)
+ AB2TOT=AB2TOT+ABS2(IBM,IG)
+ FI1TOT=FI1TOT+SUM(LHS1(IBM,:NG,IG))
+ FI2TOT=FI2TOT+SUM(LHS2(IBM,:NG,IG))
+ ENDDO
+ IF(IDF2.GT.0) OUTG2(IG)=OUTG2R(IG,1)-OUTG2R(IG,2)
+ IF((NALBP.GT.0).AND.(OUTG2(IG).LT.-1.0E-6)) THEN
+ WRITE(HSMG,'(44HDRENOU: INCONSISTENT REFERENCE LEAKAGE IN GR,
+ 1 3HOUP,I4,7H. LEAK=,1P,E13.4)') IG,OUTG2(IG)
+ CALL XABORT(HSMG)
+ ENDIF
+ ENDDO
+*----
+* COMPUTE THE ACTUAL LEAKAGE FROM OUT-CURRENTS
+*----
+ OUT=0.0D0
+ GAMMA(:NUN,:NG)=0.0
+ IF(NALBP.GT.0) THEN
+ CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
+ CALL LCMLEN(IPTRK,'QFR',MAXQF,ITYLCM)
+ ALLOCATE(KN(MAXKN),QFR(MAXQF),IQFR(MAXQF),FLUX(NUN))
+ CALL LCMGET(IPTRK,'KN',KN)
+ CALL LCMGET(IPTRK,'QFR',QFR)
+ CALL LCMGET(IPTRK,'IQFR',IQFR)
+ JPFLX=LCMGID(IPFLX,'FLUX')
+ DO IG=1,NG
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ CALL DREJ02(ITYPE,IELEM,NREG,NUN,MAXKN,MAXQF,MATCOD,KN,QFR,
+ 1 IQFR,VOL,FLUX,OUTG1(IG),GAMMA(1,IG))
+ OUT=OUT+OUTG1(IG)
+ IF(IPRINT.GT.0) THEN
+ WRITE(6,130) IG,OUTG1(IG)/REAL(AB1TOT),OUTG2(IG)/REAL(AB2TOT)
+ ENDIF
+ ENDDO
+ DEALLOCATE(FLUX,IQFR,QFR,KN)
+ ENDIF
+*----
+* COMPUTE MACRO AND REFERENCE K-EFFECTIVE
+*----
+ DEN1=0.0D0
+ DEN2=0.0D0
+ DO IG=1,NG
+ OUTG1(IG)=OUTG1(IG)+SUM(ABS1(:NMIL,IG))
+ OUTG2(IG)=OUTG2(IG)+SUM(ABS2(:NMIL,IG))
+ DEN1=DEN1+OUTG1(IG)
+ DEN2=DEN2+OUTG2(IG)
+ ENDDO
+ IF(IPRINT.GT.0) THEN
+ WRITE(6,'(/24H DRENOU: MACRO KEFF=,1P,E12.5)') FI1TOT/DEN1
+ WRITE(6,'(/24H DRENOU: REFERENCE KEFF=,1P,E12.5)') FI2TOT/DEN2
+ ENDIF
+*----
+* GET INFORMATION FROM L_OPTIMIZE OBJECT
+*----
+ CALL LCMGET(IPGRAD,'DEL-STATE',ISTATE)
+ IF(ISTATE(4).LE.2) CALL XABORT('DRENOU: NO DIRECT EFFECT WITH '
+ > //'THIS TYPE OF PERTURBATION.')
+ IF(ISTATE(7).NE.1) CALL XABORT('DRENOU: IBM1=1 EXPECTED.')
+ IF(ISTATE(8).NE.NMIL) CALL XABORT('DRENOU: IBM2=NMIL EXPECTED.')
+ NGR1=ISTATE(5)
+ NGR2=ISTATE(6)
+ NPERT=(NMIL+NALBP)*(NGR2-NGR1+1)
+ NFUNC=(NMIL+NALBP+1)*(NGR2-NGR1+1)
+ ALLOCATE(VARV(NPERT))
+ ALLOCATE(FF(NFUNC),DFF(NPERT,NFUNC))
+ CALL LCMGET(IPGRAD,'VAR-VALUE',VARV)
+*----
+* COMPUTE THE RMS FUNCTIONAL AND CONSTRAINTS
+*----
+ IREL(:NGR2-NGR1+1)=0
+ RHS(:NGR2-NGR1+1)=0.0D0
+ FF(:NFUNC)=0.0D0
+ WEI=REAL(NMIL)
+ RMSD=0.0D0
+ IPERT=0
+ IFUNC=0
+ DO IG=NGR1,NGR2
+ SUM1=0.0D0
+ DSUM=0.0D0
+ DO IBM=1,NMIL
+ IPERT=IPERT+1
+ IFUNC=IFUNC+1
+ ABS2M=MAX(EPS*AB2TOT,DBLE(ABS2(IBM,IG)))
+ DELTA=ABS1(IBM,IG)*AB2TOT/(ABS2M*AB1TOT)-ABS2(IBM,IG)/ABS2M
+ FF(IFUNC)=DELTA
+ RMSD=RMSD+DELTA**2
+ SUM1=SUM1+PHI2(IBM,IG)/VARV(IPERT)
+ DSUM=DSUM+PHI2(IBM,IG)
+ ENDDO
+ IF(NALBP.GT.0) THEN
+ OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(IG)))
+ DELTA=OUTG1(IG)*FI2TOT/(OUT2M*FI1TOT)-OUTG2(IG)/OUT2M
+ IFUNC=IFUNC+1
+ FF(IFUNC)=SQRT(WEI)*DELTA
+ RMSD=RMSD+WEI*DELTA**2
+ ENDIF
+ DELTA=SUM1/DSUM-1.0D0
+ IFUNC=IFUNC+1
+ FF(IFUNC)=DELTA
+ RMSD=RMSD+DELTA**2
+ CONST(IG-NGR1+1)=DELTA
+ IPERT=IPERT+NALBP
+ ENDDO
+ IF(IPRINT.GT.0) THEN
+ WRITE(6,100) RMSD,DOT_PRODUCT(FF,FF)
+ DO IG=NGR1,NGR2
+ WRITE(6,110) IG,CONST(IG-NGR1+1)
+ ENDDO
+ ENDIF
+ IF(IPRINT.GT.2) THEN
+ DO IG=1,NG
+ WRITE(6,'(7H GROUP=,I4)') IG
+ DO IBM=1,NMIL
+ WRITE(6,120) IBM,ABS1(IBM,IG)/REAL(AB1TOT),
+ 1 ABS2(IBM,IG)/REAL(AB2TOT)
+ ENDDO
+ ENDDO
+ ENDIF
+*----
+* STORE INFORMATION ON L_OPTIMIZE OBJECT
+*----
+ CALL LCMPUT(IPGRAD,'FOBJ-CST-VAL',NFUNC,4,FF)
+ IF(LNO) GO TO 10
+*----
+* COMPUTE THE GRADIENT MATRIX OF THE RMS FUNCTIONAL
+*----
+ ALLOCATE(SUNK(NUN))
+ JPGPT=LCMLID(IPGPT,'ASOUR',NFUNC)
+ IFUNC=0
+ DO IG=NGR1,NGR2
+ DO IBM=1,NMIL
+ ABS2M=MAX(EPS*AB2TOT,DBLE(ABS2(IBM,IG)))
+ SOUTOT=AB2TOT/AB1TOT/ABS2M
+ SOUT2=ABS1(IBM,IG)/AB1TOT
+ IFUNC=IFUNC+1
+ KPGPT=LCMLIL(JPGPT,IFUNC,NG)
+ DO JG=1,NG
+ SUNK(:NUN)=0.0
+ DO IR=1,NREG
+ IUNK=KEYFLX(IR)
+ IF(IUNK.EQ.0) CYCLE
+ JBM=MATCOD(IR)
+ IF(JBM.EQ.0) CYCLE
+ SA=SIGA(JBM,JG)
+ SOUT1=0.0D0
+ IF((IG.EQ.JG).AND.(IBM.EQ.JBM)) SOUT1=1.0D0
+ SUNK(IUNK)=REAL(SOUTOT*VOL(IR)*SA*(SOUT1-SOUT2))
+ ENDDO
+ CALL LCMPDL(KPGPT,JG,NUN,2,SUNK)
+ ENDDO
+ ENDDO
+ IF(NALBP.GT.0) THEN
+ IFUNC=IFUNC+1
+ OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(IG)))
+ SOUTOT=SQRT(WEI)*FI2TOT/FI1TOT/OUT2M
+ SOUT2=OUTG1(IG)/FI1TOT
+ KPGPT=LCMLIL(JPGPT,IFUNC,NG)
+ DO JG=1,NG
+ SUNK(:NUN)=0.0
+ DO IR=1,NREG
+ IUNK=KEYFLX(IR)
+ IF(IUNK.EQ.0) CYCLE
+ JBM=MATCOD(IR)
+ IF(JBM.EQ.0) CYCLE
+ SA=SIGA(JBM,JG)
+ SF=SIGF(JBM,JG)
+ SOUT1=0.0D0
+ IF(IG.EQ.JG) SOUT1=1.0D0
+ SUNK(IUNK)=REAL(SOUTOT*VOL(IR)*(SA*SOUT1-SF*SOUT2))
+ ENDDO
+ IF(IG.EQ.JG) THEN
+ DO IUNK=1,NUN
+ SOUT1=GAMMA(IUNK,IG)
+ SUNK(IUNK)=SUNK(IUNK)+REAL(SOUTOT*SOUT1)
+ ENDDO
+ ENDIF
+ CALL LCMPDL(KPGPT,JG,NUN,2,SUNK)
+ ENDDO
+ ENDIF
+ IFUNC=IFUNC+1
+ KPGPT=LCMLIL(JPGPT,IFUNC,NG)
+ SUNK(:NUN)=0.0
+ DO JG=1,NG
+ CALL LCMPDL(KPGPT,JG,NUN,2,SUNK)
+ ENDDO
+ ENDDO
+*----
+* CHECK SOURCE ORTHOGONALITY
+*----
+ ALLOCATE(FLUX(NUN))
+ JPFLX=LCMGID(IPFLX,'FLUX')
+ DO IFUNC=1,NFUNC
+ KPGPT=LCMGIL(JPGPT,IFUNC)
+ AIL=0.0D0
+ BIL=0.0D0
+ DO IG=1,NG
+ CALL LCMGDL(KPGPT,IG,SUNK)
+ CALL LCMGDL(JPFLX,IG,FLUX)
+ DO IUNK=1,NUN
+ GAZ=FLUX(IUNK)*SUNK(IUNK)
+ DAZ=FLUX(IUNK)**2
+ AIL=AIL+GAZ
+ BIL=BIL+DAZ
+ ENDDO
+ ENDDO
+ DSUM=ABS(AIL)/ABS(BIL)/REAL(NUN)
+ IF(IPRINT.GT.3) THEN
+ WRITE(6,'(/21H DRENOU: DOT PRODUCT=,1P,E11.4,11H COMPONENT=,
+ 1 I5)') DSUM,IFUNC
+ ENDIF
+ IF(ABS(DSUM).GT.1.0E-4) THEN
+ WRITE(HSMG,'(36HDRENOU: NON ORTHOGONAL SOURCE (DSUM=,1P,E11.3,
+ 1 26H) FOR INDIVIDUAL COMPONENT,I5,1H.)') DSUM,IFUNC
+ CALL XABORT(HSMG)
+ ENDIF
+ ENDDO
+ DEALLOCATE(FLUX,SUNK)
+*----
+* COMPUTE THE DIRECT GRADIENT MATRIX
+*----
+ IFUNC=0
+ DFF(:NPERT,:NFUNC)=0.0D0
+ DO IG=NGR1,NGR2
+ DSUM=0.0D0
+ DO IBM=1,NMIL
+ DSUM=DSUM+PHI2(IBM,IG)
+ ENDDO
+ DO IBM=1,NMIL
+ ABS2M=MAX(EPS*AB2TOT,DBLE(ABS2(IBM,IG)))
+ SOUTOT=ABS1(IBM,IG)*AB2TOT/AB1TOT/ABS2M
+ IFUNC=IFUNC+1
+ IPERT=0
+ DO JG=NGR1,NGR2
+ DO JBM=1,NMIL
+ IPERT=IPERT+1
+ IF((IG.EQ.JG).AND.(IBM.EQ.JBM)) THEN
+ DFF(IPERT,IFUNC)=DFF(IPERT,IFUNC)+SOUTOT/VARV(IPERT)
+ ENDIF
+ DFF(IPERT,IFUNC)=DFF(IPERT,IFUNC)-SOUTOT*ABS1(JBM,JG)/
+ > AB1TOT/VARV(IPERT)
+ ENDDO
+ IF(NALBP.GT.0) IPERT=IPERT+1
+ ENDDO
+ ENDDO
+ IF(NALBP.GT.0) THEN
+ IFUNC=IFUNC+1
+ IPERT=0
+ OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(IG)))
+ SOUTOT=SQRT(WEI)*FI2TOT/FI1TOT/OUT2M
+ DO JG=NGR1,NGR2
+ DO JBM=1,NMIL
+ IPERT=IPERT+1
+ IF(IG.EQ.JG) THEN
+ DFF(IPERT,IFUNC)=DFF(IPERT,IFUNC)+SOUTOT*ABS1(JBM,JG)/
+ > VARV(IPERT)
+ ENDIF
+ DFF(IPERT,IFUNC)=DFF(IPERT,IFUNC)-SOUTOT*OUTG1(IG)*
+ > SUM(LHS1(JBM,:NG,JG))/FI1TOT/VARV(IPERT)
+ ENDDO
+ IPERT=IPERT+1
+ ENDDO
+ ENDIF
+ IFUNC=IFUNC+1
+ IPERT=0
+ DO JG=NGR1,NGR2
+ DO JBM=1,NMIL
+ IPERT=IPERT+1
+ IF(IG.EQ.JG) THEN
+ DFF(IPERT,IFUNC)=DFF(IPERT,IFUNC)-PHI2(JBM,IG)/(DSUM*
+ > VARV(IPERT)**2)
+ ENDIF
+ ENDDO
+ IF(NALBP.GT.0) IPERT=IPERT+1
+ ENDDO
+ ENDDO
+ CALL LCMPUT(IPGRAD,'GRADIENT-DIR',NPERT*NFUNC,4,DFF)
+*----
+* MODIFY STATE VECTOR OF OPTIMIZE OBJECT
+*----
+ 10 CALL LCMGET(IPGRAD,'STATE-VECTOR',ISTATE)
+ ISTATE(2)=NFUNC-1
+ ISTATE(8)=4
+ CALL LCMPUT(IPGRAD,'STATE-VECTOR',NSTATE,1,ISTATE)
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(DFF,FF)
+ DEALLOCATE(VARV,SIGF,SIGA,OUTG2R,OUTG2,OUTG1,GAMMA,RHS,IREL,CONST,
+ 1 LHS2,RHS2,LHS1,RHS1,ABS2,ABS1,PHI2,PHI1)
+ RETURN
+*
+ 100 FORMAT(/40H DRENOU: RMS ERROR ON RATE DISTRIBUTION=,1P,2E11.4)
+ 110 FORMAT(23H DRENOU: CONSTRAINT(,I4,2H)=,1P,E11.4)
+ 120 FORMAT(5X,16HABSORPTION RATE(,I4,2H)=,1P,2E12.4)
+ 130 FORMAT(5X,6HGROUP=,I4,9H LEAKAGE=,1P,2E12.4)
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 14:55:01 +0000