*DECK DREKOU SUBROUTINE DREKOU(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,RMSD) * *----------------------------------------------------------------------- * *Purpose: * Compute the GPT sources corresponding to the gradient of the RMS * absorption distribution. Case with direct effect. * *Copyright: * Copyright (C) 2017 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_SOURCE 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 * RMSD RMS error on rate distribution. * *Parameters: * ITYPE * IELEM * *----------------------------------------------------------------------- * 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) 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,GRATOT,SOUTOT,AB1TOT,AB2TOT,FI1TOT, > FI2TOT,SUM1,DSUM,DELTA,OUT,SA,SF,SUNGAR,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,GRAD,RHS,CONST DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: SIGA,SIGF *---- * 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('DREKOU: 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 DREKOU: MACRO B2=,1P,E12.4)') B21 ENDIF IF((ILEAK2.EQ.1).AND.(IPRINT.GT.0)) THEN WRITE(6,'(/22H DREKOU: 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('DREKOU: MISSING ACTUAL FLUX.') CALL LCMLEN(KPMAC2,'FLUX-INTG',ILG,ITYLCM) IF(ILG.NE.NMIL) CALL XABORT('DREKOU: 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.EQ.1) OUTG2(IG)=OUTG2R(IG,1)-OUTG2R(IG,2) IF((NALBP.GT.0).AND.(OUTG2(IG).LT.-1.0E-6)) THEN WRITE(HSMG,'(44HDREKOU: 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) WRITE(6,130) IG,OUTG1(IG)/REAL(AB1TOT), 1 OUTG2(IG)/REAL(AB2TOT) 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 DREKOU: MACRO KEFF=,1P,E12.5)') FI1TOT/DEN1 WRITE(6,'(/24H DREKOU: 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('DREKOU: NO DIRECT EFFECT WITH ' > //'THIS TYPE OF PERTURBATION.') IF(ISTATE(7).NE.1) CALL XABORT('DREKOU: IBM1=1 EXPECTED.') IF(ISTATE(8).NE.NMIL) CALL XABORT('DREKOU: IBM2=NMIL EXPECTED.') IMC=ISTATE(4)-2 NGR1=ISTATE(5) NGR2=ISTATE(6) IF(IMC.LE.2) THEN NPERT=(NMIL+NALBP)*(NGR2-NGR1+1) ELSE NPERT=NALBP*(NGR2-NGR1+1) ENDIF ALLOCATE(VARV(NPERT)) CALL LCMGET(IPGRAD,'VAR-VALUE',VARV) *---- * COMPUTE THE RMS FUNCTIONAL AND CONSTRAINTS *---- IREL(:NGR2-NGR1+1)=0 RHS(:NGR2-NGR1+1)=0.0D0 WEI=REAL(NMIL) RMSD=0.0D0 IF(IMC.LE.2) THEN IPERT=0 DO IG=NGR1,NGR2 SUM1=0.0D0 DSUM=0.0D0 DO IBM=1,NMIL IPERT=IPERT+1 ABS2M=MAX(EPS*AB2TOT,DBLE(ABS2(IBM,IG))) DELTA=ABS1(IBM,IG)*AB2TOT/(ABS2M*AB1TOT)-ABS2(IBM,IG)/ABS2M RMSD=RMSD+DELTA**2 SUM1=SUM1+PHI2(IBM,IG)/VARV(IPERT) DSUM=DSUM+PHI2(IBM,IG) ENDDO DELTA=SUM1/DSUM-1.0D0 RMSD=RMSD+DELTA**2 CONST(IG-NGR1+1)=DELTA IPERT=IPERT+NALBP ENDDO ENDIF IF(NALBP.GT.0) THEN DO IG=1,NG OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(IG))) DELTA=OUTG1(IG)*FI2TOT/(OUT2M*FI1TOT)-OUTG2(IG)/OUT2M RMSD=RMSD+WEI*DELTA**2 ENDDO ENDIF IF(IPRINT.GT.0) THEN WRITE(6,100) RMSD IF(IMC.LE.2) THEN DO IG=NGR1,NGR2 WRITE(6,110) IG,CONST(IG-NGR1+1) ENDDO ENDIF ENDIF IF((IPRINT.GT.2).AND.(IMC.LE.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',1,4,RMSD) IF(LNO) GO TO 20 *---- * COMPUTE THE GRADIENT OF THE RMS FUNCTIONAL *---- ALLOCATE(SUNK(NUN)) JPGPT=LCMLID(IPGPT,'ASOUR',1) KPGPT=LCMLIL(JPGPT,1,NG) DO IG=1,NG SUNK(:NUN)=0.0 DO IR=1,NREG IUNK=KEYFLX(IR) IF(IUNK.EQ.0) CYCLE IBM=MATCOD(IR) IF(IBM.EQ.0) CYCLE SA=SIGA(IBM,IG) SF=SIGF(IBM,IG) SOUT1=0.0D0 SOUT2=0.0D0 SUNGAR=0.0D0 IF(IMC.LE.2) THEN DO JG=1,NG DO JBM=1,NMIL ABS2M=MAX(EPS*AB2TOT,DBLE(ABS2(JBM,JG))) DELTA=ABS1(JBM,JG)*AB2TOT/(ABS2M*AB1TOT)-ABS2(JBM,JG)/ 1 ABS2M IF((IG.EQ.JG).AND.(IBM.EQ.JBM)) THEN SOUT1=SOUT1+DELTA/ABS2M ENDIF SOUT2=SOUT2+(ABS1(JBM,JG)/AB1TOT)*DELTA/ABS2M ENDDO ENDDO SUNGAR=2.0D0*VOL(IR)*SA*AB2TOT*(SOUT1-SOUT2)/AB1TOT ENDIF IF(NALBP.GT.0) THEN SOUT1=0.0D0 SOUT2=0.0D0 DO JG=1,NG OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(JG))) DELTA=OUTG1(JG)*FI2TOT/(OUT2M*FI1TOT)-OUTG2(JG)/OUT2M IF(IG.EQ.JG) SOUT1=SOUT1+DELTA*SA/OUT2M SOUT2=SOUT2+(OUTG1(JG)/FI1TOT)*DELTA*SF/OUT2M ENDDO SUNGAR=SUNGAR+2.0D0*VOL(IR)*WEI*FI2TOT*(SOUT1-SOUT2)/FI1TOT ENDIF SUNK(IUNK)=REAL(SUNGAR) ENDDO IF(NALBP.GT.0) THEN OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(IG))) DELTA=OUTG1(IG)*FI2TOT/(OUT2M*FI1TOT)-OUTG2(IG)/OUT2M DO IUNK=1,NUN SOUT1=DELTA*GAMMA(IUNK,IG)/OUT2M SUNK(IUNK)=SUNK(IUNK)+2.0*WEI*REAL(FI2TOT*SOUT1/FI1TOT) ENDDO ENDIF CALL LCMPDL(KPGPT,IG,NUN,2,SUNK) ENDDO *---- * CHECK SOURCE ORTHOGONALITY *---- ALLOCATE(FLUX(NUN)) JPFLX=LCMGID(IPFLX,'FLUX') 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.0) THEN WRITE(6,'(/21H DREKOU: DOT PRODUCT=,1P,E11.4)') DSUM ENDIF IF(ABS(DSUM).GT.1.0E-4) THEN WRITE(HSMG,'(36HDREKOU: NON ORTHOGONAL SOURCE (DSUM=,1P,E11.3, 1 2H).)') DSUM CALL XABORT(HSMG) ENDIF DEALLOCATE(FLUX,SUNK) *---- * COMPUTE THE DIRECT GRADIENTS *---- ALLOCATE(GRAD(NPERT)) GRAD(:NPERT)=0.0D0 IF(IMC.GT.2) GO TO 10 IPERT=0 DO IG=NGR1,NGR2 DSUM=0.0D0 DO IBM=1,NMIL DSUM=DSUM+PHI2(IBM,IG) ENDDO DO IBM=1,NMIL IPERT=IPERT+1 GRATOT=0.0D0 DO JG=1,NG DO JBM=1,NMIL SOUTOT=0.0D0 IF((IG.EQ.JG).AND.(IBM.EQ.JBM)) SOUTOT=1.0 SOUTOT=SOUTOT-ABS1(IBM,IG)/AB1TOT ABS2M=MAX(EPS*AB2TOT,DBLE(ABS2(JBM,JG))) DELTA=ABS1(JBM,JG)*AB2TOT/(ABS2M*AB1TOT)-ABS2(JBM,JG)/ABS2M GRATOT=GRATOT+SOUTOT*ABS1(JBM,JG)*DELTA*AB2TOT/ABS2M ENDDO ENDDO GRAD(IPERT)=2.0D0*GRATOT/AB1TOT/VARV(IPERT) IF(NALBP.GT.0) THEN SOUT1=0.0D0 SOUT2=0.0D0 DO JG=1,NG OUT2M=MAX(EPSL*FI2TOT,DBLE(OUTG2(JG))) DELTA=OUTG1(JG)*FI2TOT/(OUT2M*FI1TOT)-OUTG2(JG)/OUT2M IF(IG.EQ.JG) SOUT1=SOUT1+ABS1(IBM,IG)*DELTA/OUT2M SOUT2=SOUT2+(OUTG1(JG)/FI1TOT)*SUM(LHS1(IBM,:NG,IG))* 1 DELTA/OUT2M ENDDO GRAD(IPERT)=GRAD(IPERT)+2.0D0*WEI*FI2TOT*(SOUT1-SOUT2)/ 1 FI1TOT/VARV(IPERT) ENDIF * equality constraints GRAD(IPERT)=GRAD(IPERT)-2.0D0*CONST(IG-NGR1+1)*PHI2(IBM,IG)/ 1 (DSUM*VARV(IPERT)**2) ENDDO IPERT=IPERT+NALBP ENDDO 10 CALL LCMPUT(IPGRAD,'GRADIENT-DIR',NPERT,4,GRAD) DEALLOCATE(GRAD) *---- * SCRATCH STORAGE DEALLOCATION *---- 20 DEALLOCATE(VARV,SIGF,SIGA,OUTG2R,OUTG2,OUTG1,GAMMA,RHS,IREL,CONST, 1 LHS2,RHS2,LHS1,RHS1,ABS2,ABS1,PHI2,PHI1) RETURN * 100 FORMAT(/40H DREKOU: RMS ERROR ON RATE DISTRIBUTION=,1P,E11.4) 110 FORMAT(23H DREKOU: 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