*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