1 files changed, 107 insertions, 0 deletions
diff --git a/Donjon/src/GRA001.f b/Donjon/src/GRA001.f
new file mode 100644
index 0000000..90569f5
--- /dev/null
+++ b/Donjon/src/GRA001.f
@@ -0,0 +1,107 @@
+*DECK GRA001
+      SUBROUTINE GRA001(IPFLX,IPGPT,NVAR,NCST,DERIV)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute new gradients of system characteristics (part 2).
+*
+*Copyright:
+* Copyright (C) 2012 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
+* IPFLX   pointer of a LCM object containing a set solutions of
+*         fixed-source eigenvalue problems.
+* IPGPT   pointer of a LCM object containing a set of fixed sources.
+* NVAR    number of control variables.
+* NCST    number of constraints with indirect effects (can be zero).
+*
+*Parameters: output
+* DERIV   gradient matrix.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPFLX,IPGPT
+      INTEGER NVAR,NCST
+      DOUBLE PRECISION DERIV(NVAR,NCST+1)
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPFLX,KPFLX,JPGPT,KPGPT
+      PARAMETER (NSTATE=40)
+      INTEGER ISTATE(NSTATE)
+      DOUBLE PRECISION SUM
+      REAL, ALLOCATABLE, DIMENSION(:) :: DFLUX,SOUR
+*
+      CALL LCMGET(IPFLX,'STATE-VECTOR',ISTATE)
+      NGRP=ISTATE(1)
+      NUN=ISTATE(2)
+      ITYPE=ISTATE(3)
+      NGPT=ISTATE(5)
+      CALL LCMGET(IPGPT,'STATE-VECTOR',ISTATE)
+      IF(ISTATE(1).NE.NGRP) CALL XABORT('GRA001: INVALID NGRP')
+      IF(ISTATE(2).NE.NUN) CALL XABORT('GRA001: INVALID NUN')
+      ND=ISTATE(3)
+      NA=ISTATE(4)
+*
+      ALLOCATE(SOUR(NUN),DFLUX(NUN))
+      DERIV(:NVAR,:NCST+1)=0.0D0
+      IF(ITYPE.EQ.100) THEN
+*       EXPLICIT APPROACH
+        IF(NVAR.NE.NGPT) CALL XABORT('GRA001: INVALID NGPT(1)')
+        IF(NCST+1.NE.NA) CALL XABORT('GRA001: INVALID NA(1)')
+        JPFLX=LCMGID(IPFLX,'DFLUX')
+        JPGPT=LCMGID(IPGPT,'ASOUR')
+        DO 25 IVAR=1,NVAR
+        KPFLX=LCMGIL(JPFLX,IVAR)
+        DO 20 ICST=1,NCST+1
+        KPGPT=LCMGIL(JPGPT,ICST)
+        SUM=0.0D0
+        DO 15 IGR=1,NGRP
+        CALL LCMGDL(KPGPT,IGR,SOUR)
+        CALL LCMGDL(KPFLX,IGR,DFLUX)
+        DO 10 IUN=1,NUN
+        SUM=SUM+SOUR(IUN)*DFLUX(IUN)
+   10   CONTINUE
+   15   CONTINUE
+        DERIV(IVAR,ICST)=SUM
+   20   CONTINUE
+   25   CONTINUE
+      ELSE IF(ITYPE.EQ.1000) THEN
+*       IMPLICIT APPROACH
+        IF(NVAR.NE.ND) CALL XABORT('GRA001: INVALID ND(2)')
+        IF(NCST+1.NE.NGPT) CALL XABORT('GRA001: INVALID NGPT(2)')
+        JPFLX=LCMGID(IPFLX,'ADFLUX')
+        JPGPT=LCMGID(IPGPT,'DSOUR')
+        DO 45 ICST=1,NCST+1
+        KPFLX=LCMGIL(JPFLX,ICST)
+        DO 40 IVAR=1,NVAR
+        KPGPT=LCMGIL(JPGPT,IVAR)
+        SUM=0.0D0
+        DO 35 IGR=1,NGRP
+        CALL LCMGDL(KPGPT,IGR,SOUR)
+        CALL LCMGDL(KPFLX,IGR,DFLUX)
+        DO 30 IUN=1,NUN
+        SUM=SUM+SOUR(IUN)*DFLUX(IUN)
+   30   CONTINUE
+   35   CONTINUE
+        DERIV(IVAR,ICST)=SUM
+   40   CONTINUE
+   45   CONTINUE
+      ELSE
+        CALL XABORT('GRA001: INVALID FLUX OBJECT')
+      ENDIF
+      DEALLOCATE(DFLUX,SOUR)
+      RETURN
+      END