Initial commit from Polytechnique Montreal

1 files changed, 106 insertions, 0 deletions

diff --git a/Trivac/src/TRIALB.f b/Trivac/src/TRIALB.f
new file mode 100755
index 0000000..9b26ca2
--- /dev/null
+++ b/Trivac/src/TRIALB.f
@@ -0,0 +1,106 @@
+*DECK TRIALB
+      SUBROUTINE TRIALB(IPTRK,IPMACR,IPMACP,IPSYS,NGRP,NALBP,IPR,GAMMA)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Process physical albedo information and calculation of multigroup
+* albedo functions.
+*
+*Copyright:
+* Copyright (C) 2018 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
+* IPTRK   L_TRACK pointer to the TRIVAC tracking information.
+* IPMACR  L_MACROLIB pointer to the unperturbed cross sections.
+* IPMACP  L_MACROLIB pointer to the perturbed cross sections if
+*         IPR.gt.0. Equal to IPMACR if IPR=0.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* NGRP    number of energy groups.
+* NALBP   number of physical albedos per energy group.
+* IPR     type of assembly:
+*         =0: calculation of the system matrices;
+*         =1: calculation of the derivative of these matrices;
+*         =2: calculation of the first variation of these matrices;
+*         =3: identical to IPR=2, but these variation are added to
+*         unperturbed system matrices.
+*
+*Parameters: output
+* GAMMA   albedo functions
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPMACR,IPMACP,IPSYS
+      INTEGER NGRP,NALBP,IPR
+      REAL GAMMA(NALBP,NGRP)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (NSTATE=40)
+      INTEGER ISTATE(NSTATE)
+      CHARACTER TEXT12*12
+      REAL, DIMENSION(:,:), ALLOCATABLE :: ALBP,DALBP
+*
+      ALB(X)=0.5*(1.0-X)/(1.0+X)
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(ALBP(NALBP,NGRP),DALBP(NALBP,NGRP))
+*----
+*  RECOVER PHYSICAL ALBEDOS
+*----
+      IF(NALBP.EQ.0) CALL XABORT('TRIALB: NO PHYSICAL ALBEDOS.')
+      CALL LCMGET(IPMACR,'ALBEDO',ALBP)
+      IF(IPR.GT.0) CALL LCMGET(IPMACP,'ALBEDO',DALBP)
+*----
+*  COMPUTE ALBEDO FUNCTIONS
+*----
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      ICHX=ISTATE(12)
+      DO IGR=1,NGRP
+        GAMMA(:NALBP,IGR)=0.0
+        DO IALB=1,NALBP
+          IF(ICHX.NE.2) THEN
+            IF(IPR.EQ.0) THEN
+              GAMMA(IALB,IGR)=ALB(ALBP(IALB,IGR))
+            ELSE
+              GAMMA(IALB,IGR)=ALB(DALBP(IALB,IGR))
+            ENDIF
+          ELSE IF((ICHX.EQ.2).AND.(ALBP(IALB,IGR).NE.1.0)) THEN
+            IF(IPR.EQ.0) THEN
+              GAMMA(IALB,IGR)=1.0/ALB(ALBP(IALB,IGR))
+            ELSE IF(IPR.EQ.1) THEN
+              GG=ALB(ALBP(IALB,IGR))
+              DGG=ALB(DALBP(IALB,IGR))
+              GAMMA(IALB,IGR)=-DGG/(GG**2)
+            ELSE
+              GG=ALB(ALBP(IALB,IGR))
+              DGG=ALB(ALBP(IALB,IGR))+ALB(DALBP(IALB,IGR))
+              GAMMA(IALB,IGR)=1.0/DGG-1.0/GG
+            ENDIF
+          ELSE IF((ICHX.EQ.2).AND.(ALBP(IALB,IGR).EQ.1.0)) THEN
+            GAMMA(IALB,IGR)=1.0E20
+          ENDIF
+        ENDDO
+*----
+*  SAVE ALBEDO FUNCTIONS ON IPSYS
+*----
+         WRITE(TEXT12,'(9HALBEDO-FU,I3.3)') IGR
+         CALL LCMPUT(IPSYS,TEXT12,NALBP,2,GAMMA(1,IGR))
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(DALBP,ALBP)
+      RETURN
+      END