diff options
Diffstat (limited to 'Donjon/src/DREKOU.f')
| -rw-r--r-- | Donjon/src/DREKOU.f | 511 |
1 files changed, 511 insertions, 0 deletions
diff --git a/Donjon/src/DREKOU.f b/Donjon/src/DREKOU.f new file mode 100644 index 0000000..a9b1a97 --- /dev/null +++ b/Donjon/src/DREKOU.f @@ -0,0 +1,511 @@ +*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 |