Initial commit from Polytechnique Montreal

1 files changed, 383 insertions, 0 deletions

diff --git a/Dragon/src/COMCAT.f b/Dragon/src/COMCAT.f
new file mode 100644
index 0000000..e885642
--- /dev/null
+++ b/Dragon/src/COMCAT.f
@@ -0,0 +1,383 @@
+*DECK COMCAT
+      SUBROUTINE COMCAT(IPCPO,IPRHS,NORIG,NPARN,MUPCPO,LGNCPO,LWARN)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Catenate a RHS compo into the output multicompo.
+*
+*Copyright:
+* Copyright (C) 2002 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
+* IPCPO   pointer to the output multicompo.
+* IPRHS   pointer to the rhs multicompo (contains the new calculations).
+* NORIG   index of the elementary calculation associated to the
+*         father node in the parameter tree.
+* NPARN   number of global parameters in the output multicompo.
+* MUPCPO  tuple of the new global parameters in the output multicompo.
+* LGNCPO  LGNEW value of the new global parameters in the output
+*         multicompo.
+* LWARN   logical used in case if an elementary calculation in the RHS 
+*         is already present in CPO. If LWARN=.true. a warning is send 
+*         and the CPO values are kept, otherwise XABORT is called
+*         (default).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPCPO,IPRHS
+      INTEGER NORIG,NPARN,MUPCPO(NPARN)
+      LOGICAL LGNCPO(NPARN),LWARN
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPCPO,KPCPO,JPRHS,KPRHS,LPCPO,LPRHS,MPCPO,MPRHS
+      PARAMETER (NSTATE=40,MAXPAR=50,MAXVAL=1000)
+      INTEGER ISTATE(NSTATE),NVPO(2),NVALUE(2*MAXPAR),MUPLET(2*MAXPAR),
+     1 MUPRHS(2*MAXPAR)
+      CHARACTER HSMG*131,RECNAM*12,TEXT4*4,TEXT12*12,PARFMT(MAXPAR)*8,
+     1 VCHAR(MAXVAL)*12,PARKEY(MAXPAR)*12,PARCPO(MAXPAR)*12
+      LOGICAL COMTRE,LGERR,LGNEW(MAXPAR)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: ILCALR,MUOLD,IORRHS,IDEBAR,
+     1 IARBVA,JDEBAR,JARBVA,IORIGI,VINTE
+      REAL, ALLOCATABLE, DIMENSION(:) :: VREAL
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: LGOLD
+*
+      NIDEM=0
+      CALL LCMGET(IPRHS,'STATE-VECTOR',ISTATE)
+      NMIL=ISTATE(1)
+      NG=ISTATE(2)
+      NCALR=ISTATE(3)
+      NPAR=ISTATE(5)
+      NLOC=ISTATE(6)
+      IF(NCALR.EQ.0) CALL XABORT('COMCAT: NO CALCULATION IN RHS COMPO.')
+      ALLOCATE(ILCALR(NCALR))
+*
+      CALL LCMGET(IPCPO,'STATE-VECTOR',ISTATE)
+      NCAL=ISTATE(3)
+      IF(NPARN+NLOC.GT.2*MAXPAR) CALL XABORT('COMCAT: MAXPAR OVERFLOW.')
+      IF(NCAL.EQ.0) THEN
+*        COMPLETE STATE-VECTOR.
+         IF(ISTATE(1).EQ.0) THEN
+            ISTATE(1)=NMIL
+         ELSE IF(NMIL.NE.ISTATE(1)) THEN
+            WRITE(HSMG,'(42HCOMCAT: ELEMENTARY CALCULATION WITH AN INV,
+     1      22HALIB NB. OF MIXTURES =,I7,3H NE,I7,1H.)') NMIL,ISTATE(1)
+            CALL XABORT(HSMG)
+         ENDIF
+         ISTATE(2)=NG
+      ELSE
+         IF(NMIL.NE.ISTATE(1)) THEN
+            WRITE(HSMG,'(42HCOMCAT: ELEMENTARY CALCULATION WITH AN INV,
+     1      22HALIB NB. OF MIXTURES =,I7,3H NE,I7,1H.)') NMIL,ISTATE(1)
+            CALL XABORT(HSMG)
+         ELSE IF(NG.NE.ISTATE(2)) THEN
+            WRITE(HSMG,'(42HCOMCAT: ELEMENTARY CALCULATION WITH AN INV,
+     1      20HALIB NB. OF GROUPS =,I7,3H NE,I7,1H.)') NG,ISTATE(2)
+            CALL XABORT(HSMG)
+         ENDIF
+      ENDIF
+      IF(NPAR.GT.NPARN) THEN
+         WRITE(HSMG,'(42HCOMCAT: ELEMENTARY CALCULATION WITH AN INV,
+     1   31HALIB NB. OF GLOBAL PARAMETERS =,I7,3H NE,I7,1H.)') NPAR,
+     2   NPARN
+         CALL XABORT(HSMG)
+      ELSE IF(NLOC.NE.ISTATE(6)) THEN
+         WRITE(HSMG,'(42HCOMCAT: ELEMENTARY CALCULATION WITH AN INV,
+     1   30HALIB NB. OF LOCAL PARAMETERS =,I7,3H NE,I7,1H.)') NLOC,
+     2   ISTATE(6)
+         CALL XABORT(HSMG)
+      ENDIF
+*----
+*  ADJUST THE SIZE OF THE OUTPUT COMPO
+*----
+      ISTATE(3)=ISTATE(3)+NCALR
+      IF(ISTATE(3).GT.ISTATE(4)) THEN
+         ISTATE(4)=ISTATE(4)+NCALR+9
+         JPCPO=LCMLID(IPCPO,'MIXTURES',NMIL)
+         DO 10 IBM=1,NMIL
+         KPCPO=LCMDIL(JPCPO,IBM)
+         LPCPO=LCMLID(KPCPO,'CALCULATIONS',ISTATE(4))
+   10    CONTINUE
+      ENDIF
+      MAXCAL=ISTATE(4)
+*----
+*  UPDATE THE STATE-VECTOR
+*----
+      CALL LCMPUT(IPCPO,'STATE-VECTOR',NSTATE,1,ISTATE)
+*----
+*  MAIN LOOP OVER THE HOMOGENEOUS MIXTURES *********************
+*----
+      JPCPO=LCMLID(IPCPO,'MIXTURES',NMIL)
+      JPRHS=LCMGID(IPRHS,'MIXTURES')
+      ALLOCATE(MUOLD(NCALR*NPARN),LGOLD(NCALR*NPARN))
+      DO 190 IBM=1,NMIL
+      KPCPO=LCMDIL(JPCPO,IBM)
+      KPRHS=LCMGIL(JPRHS,IBM)
+*----
+*  MAIN LOOP OVER THE NCALR ELEMENTARY CALCULATIONS OF THE RHS COMPO
+*----
+      ILCALR(:NCALR)=1
+      NIDEM=0
+      DO 170 ICAL=1,NCALR
+*----
+*  COMPUTE THE MUPLET VECTOR FROM THE RHS COMPO
+*----
+      CALL LCMSIX(KPRHS,'TREE',1)
+      CALL LCMLEN(KPRHS,'ARBVAL',MAXNVP,ITYLCM)
+      CALL LCMLEN(KPRHS,'ORIGIN',MAXNCA,ITYLCM)
+      ALLOCATE(IORRHS(MAXNCA))
+      CALL LCMGET(KPRHS,'ORIGIN',IORRHS)
+      ALLOCATE(JDEBAR(MAXNVP+1),JARBVA(MAXNVP))
+      CALL LCMGET(KPRHS,'NCALS',NCALS)
+      CALL LCMGET(KPRHS,'DEBARB',JDEBAR)
+      CALL LCMGET(KPRHS,'ARBVAL',JARBVA)
+      CALL LCMGET(KPRHS,'NVP',NVPO)
+      DO 30 I=NVPO(1)-NCALS+1,NVPO(1)
+      IF(JDEBAR(I+1).EQ.ICAL) THEN
+         I0=I
+         GO TO 40
+      ENDIF
+   30 CONTINUE
+      CALL XABORT('COMCAT: MUPLET ALGORITHM FAILURE 1.')
+   40 MUPRHS(NPAR+NLOC)=JARBVA(I0)
+      DO 65 IPAR=NPAR+NLOC-1,1,-1
+      DO 50 I=1,NVPO(1)-NCALS
+      IF(JDEBAR(I+1).GT.I0) THEN
+         I0=I
+         GO TO 60
+      ENDIF
+   50 CONTINUE
+      CALL XABORT('COMCAT: MUPLET ALGORITHM FAILURE 2.')
+   60 MUPRHS(IPAR)=JARBVA(I0)
+   65 CONTINUE
+      DEALLOCATE(JARBVA,JDEBAR)
+*----
+*  RECOVER THE GLOBAL PARAMETERS
+*----
+      IF(IBM.EQ.1) THEN
+         DO 70 I=1,NPARN
+         MUPLET(I)=MUPCPO(I)
+         LGNEW(I)=LGNCPO(I)
+   70    CONTINUE
+         CALL LCMSIX(IPCPO,'GLOBAL',1)
+         CALL LCMSIX(IPRHS,'GLOBAL',1)
+         CALL LCMGTC(IPCPO,'PARKEY',12,NPARN,PARCPO)
+         CALL LCMGTC(IPRHS,'PARKEY',12,NPAR,PARKEY)
+         CALL LCMGTC(IPRHS,'PARFMT',8,NPAR,PARFMT)
+         CALL LCMGET(IPRHS,'NVALUE',NVALUE)
+         DO 100 IPAR=1,NPAR
+            DO 80 I0=1,NPARN
+            IF(PARKEY(IPAR).EQ.PARCPO(I0)) THEN
+               IPARN=I0
+               GO TO 90
+            ENDIF
+   80       CONTINUE
+            CALL XABORT('COMCAT: UNABLE TO FIND '//PARKEY(IPAR)//'.')
+   90       WRITE(RECNAM,'(''pval'',I8.8)') IPAR
+            IVAL=MUPRHS(IPAR)
+            IF(PARFMT(IPAR).EQ.'REAL') THEN
+               ALLOCATE(VREAL(NVALUE(IPAR)))
+               CALL LCMGET(IPRHS,RECNAM,VREAL)
+               FLOTT=VREAL(IVAL)
+               DEALLOCATE(VREAL)
+            ELSE IF(PARFMT(IPAR).EQ.'INTEGER') THEN
+               ALLOCATE(VINTE(NVALUE(IPAR)))
+               CALL LCMGET(IPRHS,RECNAM,VINTE)
+               NITMA=VINTE(IVAL)
+               DEALLOCATE(VINTE)
+            ELSE IF(PARFMT(IPAR).EQ.'STRING') THEN
+               IF(NVALUE(IPAR).GT.MAXVAL) CALL XABORT('COMCAT: MAXVAL '
+     1         //'OVERFLOW.')
+               CALL LCMGTC(IPRHS,RECNAM,12,NVALUE(IPAR),VCHAR)
+               TEXT12=VCHAR(IVAL)
+            ENDIF
+            CALL COMPAV(IPCPO,IPARN,NPARN,PARFMT(IPAR),FLOTT,NITMA,
+     1      TEXT12,MUPLET(IPARN),LGNEW(IPARN))
+  100    CONTINUE
+         DO 110 IPARN=1,NPARN
+         MUOLD((ICAL-1)*NPARN+IPARN)=MUPLET(IPARN)
+         LGOLD((ICAL-1)*NPARN+IPARN)=LGNEW(IPARN)
+  110    CONTINUE
+         CALL LCMSIX(IPRHS,' ',2)
+         CALL LCMSIX(IPCPO,' ',2)
+      ELSE
+         DO 120 IPARN=1,NPARN
+         MUPLET(IPARN)=MUOLD((ICAL-1)*NPARN+IPARN)
+         LGNEW(IPARN)=LGOLD((ICAL-1)*NPARN+IPARN)
+  120    CONTINUE
+      ENDIF
+*----
+*  RECOVER THE LOCAL PARAMETERS
+*----
+      CALL LCMSIX(KPCPO,'TREE',1)
+      DO 130 ILOC=1,NLOC
+         WRITE(RECNAM,'(''pval'',I8.8)') ILOC
+         IVAL=MUPRHS(NPAR+ILOC)
+         CALL LCMLEN(KPRHS,RECNAM,ILONG,ITYLCM)
+         ALLOCATE(VREAL(ILONG))
+         CALL LCMGET(KPRHS,RECNAM,VREAL)
+         FLOTT=VREAL(IVAL)
+         DEALLOCATE(VREAL)
+         CALL COMPAV(KPCPO,ILOC,NLOC,PARFMT(ILOC),FLOTT,NITMA,TEXT12,
+     1   MUPLET(NPARN+ILOC),LGNEW(NPARN+ILOC))
+  130 CONTINUE
+      CALL LCMSIX(KPRHS,' ',2)
+*----
+*  UPDATE THE PARAMETER TREE IN THE OUTPUT COMPO
+*----
+      CALL LCMLEN(KPCPO,'NVP',ILONG,ITYLCM)
+      IF(ILONG.EQ.0) THEN
+         MAXNVP=20*(NPARN+NLOC+1)
+         ALLOCATE(IDEBAR(MAXNVP+1),IARBVA(MAXNVP))
+         IDEBAR(:MAXNVP+1)=0
+         IARBVA(:MAXNVP)=0
+         IARBVA=0
+         DO 140 I=1,NPARN+NLOC
+         IDEBAR(I)=I+1
+         IARBVA(I+1)=1
+  140    CONTINUE
+         IDEBAR(NPARN+NLOC+1)=NPARN+NLOC+2
+         IDEBAR(NPARN+NLOC+2)=1
+         NCALS=1
+         NVPNEW=NPARN+NLOC+1
+      ELSE
+         CALL LCMLEN(KPCPO,'ARBVAL',JLONG,ITYLCM)
+         ALLOCATE(JDEBAR(JLONG+1),JARBVA(JLONG))
+         CALL LCMGET(KPCPO,'NCALS',NCALS)
+         CALL LCMGET(KPCPO,'DEBARB',JDEBAR)
+         CALL LCMGET(KPCPO,'ARBVAL',JARBVA)
+         CALL LCMGET(KPCPO,'NVP',NVPO)
+         DO 150 IPAR=1,NPARN+NLOC
+         IF(LGNEW(IPAR)) THEN
+            II=IPAR
+            GO TO 160
+         ENDIF
+  150    CONTINUE
+         II=NPARN+NLOC+1
+  160    LGERR=COMTRE(NPARN+NLOC,NVPO(1),JARBVA,JDEBAR,MUPLET,KK,I0,
+     1   IORD,JJ,LAST)
+         IF((II.GT.NPARN+NLOC).AND.LGERR) THEN
+            WRITE(TEXT4,'(I4)') IORD
+            IF(LWARN) THEN
+               WRITE(6,*)'COMCAT: ELEMENTARY CALCULATION HAS THE ',
+     1         'SAME PARAMETERS AS ELEMENTARY CALCULATION NB ',TEXT4
+               DEALLOCATE(JARBVA,JDEBAR,IORRHS)
+               CALL LCMSIX(KPCPO,' ',2)
+               ILCALR(ICAL)=0
+               NIDEM=NIDEM+1
+               GOTO 170
+            ELSE
+               CALL XABORT('COMCAT: ELEMENTARY CALCULATION HAS THE '//
+     1         'SAME PARAMETERS AS ELEMENTARY CALCULATION NB '//TEXT4)
+            ENDIF
+         ENDIF
+*
+*        Size of the new tree.
+*
+         NVPNEW=NVPO(1)+NPARN+NLOC+1-MIN(II,KK)
+         IF(NVPNEW.GT.MAXNVP) MAXNVP=NVPNEW+MAXNVP
+         ALLOCATE(IDEBAR(MAXNVP+1),IARBVA(MAXNVP))
+         IDEBAR(NVPNEW+2:MAXNVP+1)=0
+         IARBVA(NVPNEW+1:MAXNVP)=0
+*
+*        Update values and suppress old PARBRE.
+*
+         CALL COMARB(NPARN+NLOC,NVPO(1),NVPNEW,JDEBAR,JARBVA,LGNEW,
+     1   MUPLET,NCALS,IDEBAR,IARBVA)
+         DEALLOCATE(JARBVA,JDEBAR)
+      ENDIF
+      IF(NCALS.NE.NCAL+ICAL-NIDEM) CALL XABORT('COMCAT: INVALID NCALS.')
+      NVPO(1)=NVPNEW
+      NVPO(2)=MAXNVP
+      CALL LCMPUT(KPCPO,'NVP',2,1,NVPO)
+      CALL LCMPUT(KPCPO,'NCALS',1,1,NCALS)
+      CALL LCMPUT(KPCPO,'DEBARB',NVPO(1)+1,1,IDEBAR)
+      CALL LCMPUT(KPCPO,'ARBVAL',NVPO(1),1,IARBVA)
+      DEALLOCATE(IARBVA,IDEBAR)
+      IF(NCALS.EQ.1) THEN
+         ALLOCATE(IORIGI(MAXCAL))
+         IORIGI(:MAXCAL)=0
+      ELSE
+         CALL LCMLEN(KPCPO,'ORIGIN',MAXOLD,ITYLCM)
+         IF(MAXOLD.GT.MAXCAL) CALL XABORT('COMCAT: ORIGIN OVERFLOW(1).')
+         ALLOCATE(IORIGI(MAXCAL))
+         IORIGI(:MAXCAL)=0
+         CALL LCMGET(KPCPO,'ORIGIN',IORIGI)
+      ENDIF
+      IF(NCALS.GT.MAXCAL) CALL XABORT('COMCAT: ORIGIN OVERFLOW(2).')
+      IF(IORRHS(ICAL).EQ.0) THEN
+         IORIGI(NCALS)=NORIG
+      ELSE
+         IORIGI(NCALS)=NCAL+IORRHS(ICAL)
+      ENDIF
+      CALL LCMPUT(KPCPO,'ORIGIN',NCALS,1,IORIGI)
+      DEALLOCATE(IORIGI)
+      CALL LCMSIX(KPCPO,' ',2)
+      DEALLOCATE(IORRHS)
+  170 CONTINUE
+* END OF LOOP ON MIXTURES. ***********************************
+*----
+*  RECOVER THE MICROLIBS
+*----
+      LPCPO=LCMLID(KPCPO,'CALCULATIONS',ISTATE(4))
+      LPRHS=LCMGID(KPRHS,'CALCULATIONS')
+      NIDEM=0
+      DO 180 ICAL=1,NCALR
+         IF(ILCALR(ICAL).EQ.1)THEN
+            MPCPO=LCMDIL(LPCPO,NCAL+ICAL-NIDEM)
+            MPRHS=LCMGIL(LPRHS,ICAL)
+            CALL LCMEQU(MPRHS,MPCPO)
+         ELSE
+            NIDEM=NIDEM+1
+         ENDIF
+  180 CONTINUE
+  190 CONTINUE
+      DEALLOCATE(LGOLD,MUOLD)
+*----
+*  RECOVER THE DEPLETION CHAIN
+*----
+      CALL LCMLEN(IPRHS,'DEPL-CHAIN',ILONG,ITYLCM)
+      IF(ILONG.EQ.-1) THEN
+         CALL LCMSIX(IPCPO,'DEPL-CHAIN',1)
+         CALL LCMSIX(IPRHS,'DEPL-CHAIN',1)
+         CALL LCMEQU(IPRHS,IPCPO)
+         CALL LCMSIX(IPRHS,' ',2)
+         CALL LCMSIX(IPCPO,' ',2)
+      ENDIF
+*----
+*  RECOVER THE MACRO-GEOMETRIES
+*----
+      CALL LCMLEN(IPRHS,'GEOMETRIES',ILONG,ITYLCM)
+      IF(ILONG.GT.0) THEN
+         JPCPO=LCMLID(IPCPO,'GEOMETRIES',MAXCAL)
+         JPRHS=LCMGID(IPRHS,'GEOMETRIES')
+         DO 200 ICAL=1,NCALR
+           IF(ILCALR(ICAL).EQ.1)THEN
+             KPCPO=LCMDIL(JPCPO,NCAL+ICAL)
+             KPRHS=LCMGIL(JPRHS,ICAL)
+             CALL LCMEQU(KPRHS,KPCPO)
+           ENDIF
+  200    CONTINUE
+         ISTATE(11)=1
+      ENDIF
+      DEALLOCATE(ILCALR)
+      ISTATE(3)=ISTATE(3)-NIDEM
+      IF(ISTATE(15).EQ.-1) ISTATE(15)=0
+      CALL LCMPUT(IPCPO,'STATE-VECTOR',NSTATE,1,ISTATE)
+* END OF LOOP ON ELEMENTARY CALCULATIONS. ********************
+      RETURN
+      END