summaryrefslogtreecommitdiff
path: root/Dragon/src/MUSACG.f90
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MUSACG.f90')
-rw-r--r--Dragon/src/MUSACG.f90723
1 files changed, 723 insertions, 0 deletions
diff --git a/Dragon/src/MUSACG.f90 b/Dragon/src/MUSACG.f90
new file mode 100644
index 0000000..4f68a52
--- /dev/null
+++ b/Dragon/src/MUSACG.f90
@@ -0,0 +1,723 @@
+!
+!---------------------------------------------------------------------
+!
+!Purpose:
+! To extract a macro geometry and construct its geometry basic
+! information.
+!
+!Copyright:
+! Copyright (C) 2025 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
+! ITRACK pointer to the tracking in IMACROth macro geometry.
+! IFTRK pointer to the TRACKING file in creation mode.
+! IPRINT print level.
+! IMACRO macro geometry index.
+! NBSLIN maximum number of segments in a single tracking line.
+! RCUTOF minimum distance between two surfacic elements.
+! GG geometry basic information.
+!
+!Parameters: output
+! LGINF leakage flag (=.true. if no leakage).
+! NBNODE_MACRO number of nodes in macro IMACRO.
+! NSURF_MACRO number of boundary surfaces in macro IMACRO.
+! GG geometry basic information.
+!
+!---------------------------------------------------------------------
+!
+SUBROUTINE MUSACG(ITRACK,IFTRK,IPRINT,IMACRO,NBSLIN,RCUTOF,GG,LGINF,NBNODE_MACRO,NSURF_MACRO)
+ USE GANLIB
+ USE PRECISION_AND_KINDS, ONLY : PDB, PI
+ USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC,NIPAR,NRPAR,G_BC_TYPE,ISPEC,TYPGEO,NBFOLD,ALLSUR
+ USE SAL_GEOMETRY_MOD, ONLY : SAL130_2,SAL130_4,SAL130_6,SAL130_8,SAL130_10,SAL131_2, &
+ & SAL140,SAL160_2,SAL170,SALFOLD_0
+ USE SAL_NUMERIC_MOD, ONLY : SAL141,FINDLC,DET_ROSETTA
+ USE SAL_TRACKING_TYPES, ONLY : PRTIND,NIPART,NRPART,EPS1
+ !----
+ ! Subroutine arguments
+ !----
+ TYPE(C_PTR) ITRACK
+ INTEGER IFTRK,IPRINT,IMACRO
+ REAL(PDB) RCUTOF
+ LOGICAL LGINF,LGBC
+ TYPE(T_G_BASIC) :: GG
+ !----
+ ! Local variables
+ !----
+ INTEGER, PARAMETER :: NSTATE=40
+ INTEGER, PARAMETER :: FOUT=6
+ INTEGER, PARAMETER :: NDIM=2 ! NUMBER OF DIMENSIONS
+ INTEGER ELEM, OK, TYPE
+ REAL(PDB) :: X1,X2,Y1,Y2,DET1,DET2
+ REAL(PDB) :: DGMESHX(2),DGMESHY(2)
+ INTEGER, DIMENSION(NSTATE) :: I_STATE,IEDIMG
+ CHARACTER(LEN=72) :: TEXT72
+ !----
+ ! Allocatable arrays
+ !----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: ELEM_LIST,NODE_LIST,NODE_MACRO, &
+ & ELEM_MACRO,PERIM_MACRO,AUX_ARR,SURF_MACRO,ICODE,IFOLD,IFOLD_GLOB, &
+ & PERIM_SURF
+ INTEGER, DIMENSION(:) , ALLOCATABLE :: ITAB ! LOCAL ARRAY
+ REAL(PDB), ALLOCATABLE, DIMENSION(:) :: ANGLE,ALBEDO
+ REAL, ALLOCATABLE, DIMENSION(:) :: VOLUME,GALBED
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: MATALB,KEYMRG,IBC
+ REAL(PDB), ALLOCATABLE, DIMENSION(:) :: VOLSUR
+ REAL(PDB), ALLOCATABLE, DIMENSION(:,:,:) :: ALIGN
+ TYPE(T_G_BASIC), ALLOCATABLE :: GG_MAC
+ LOGICAL, ALLOCATABLE, DIMENSION(:) :: LFOLD
+ !----
+ ! GG_MAC allocation
+ !----
+ ALLOCATE(GG_MAC)
+ !----
+ ! Set isotropic tracking
+ !----
+ ISPEC=0
+ !----
+ ! List of nodes and surface elements in IMACRO
+ !----
+ IF(IPRINT.GT.0) WRITE(FOUT,'(/A,I6,A)') ' ********** IMACRO=',IMACRO,' **********'
+ ALLOCATE(ELEM_LIST(GG%NB_ELEM),NODE_LIST(GG%NB_NODE))
+ ELEM_LIST(:GG%NB_ELEM) = 0
+ NODE_LIST(:GG%NB_NODE) = -9999
+ I=0
+ DO ELEM=1,GG%NB_ELEM
+ I2=GG%IPAR(2,ELEM)
+ I3=GG%IPAR(3,ELEM)
+ IF((I2.GT.0).AND.(ELEM_LIST(ELEM) == 0)) THEN
+ IF(GG%NUM_MACRO(GG%NUM_MERGE(I2)) == IMACRO) THEN
+ I=I+1
+ ELEM_LIST(ELEM)=I
+ NODE_LIST(I2)=1
+ ENDIF
+ ENDIF
+ IF((I3.GT.0).AND.(ELEM_LIST(ELEM) == 0)) THEN
+ IF(GG%NUM_MACRO(GG%NUM_MERGE(I3)) == IMACRO) THEN
+ I=I+1
+ ELEM_LIST(ELEM)=I
+ NODE_LIST(I3)=1
+ ENDIF
+ ENDIF
+ ENDDO ! ELEM
+ DO J=1,GG%NB_NODE
+ IF(NODE_LIST(J) /= -9999) THEN
+ NBNODE_MACRO = NBNODE_MACRO+1
+ NODE_LIST(J) = NBNODE_MACRO
+ ENDIF
+ ENDDO
+ NELEM_MACRO=I
+ ALLOCATE(NODE_MACRO(NBNODE_MACRO),ELEM_MACRO(NELEM_MACRO))
+ DO I=1,NELEM_MACRO
+ ELEM = FINDLC(ELEM_LIST,I)
+ ELEM_MACRO(I) = ELEM
+ ENDDO
+ DO I=1,NBNODE_MACRO
+ NODE_MACRO(I) = FINDLC(NODE_LIST,I)
+ ENDDO
+ DEALLOCATE(NODE_LIST,ELEM_LIST)
+ !----
+ ! Find perimeter. Three points are colinear if the determinant is zero.
+ !----
+ ALLOCATE(PERIM_MACRO(NELEM_MACRO),ANGLE(NELEM_MACRO),ALBEDO(NELEM_MACRO),LFOLD(NELEM_MACRO))
+ ALLOCATE(ALIGN(3,3,NELEM_MACRO))
+ ALIGN(:3,3,:NELEM_MACRO)=1.0_PDB
+ PERIM_MACRO(:NELEM_MACRO)=0
+ NPERIM=0
+ ITER0: DO I=1,NELEM_MACRO
+ ELEM = ELEM_MACRO(I)
+ DO J=1,NBNODE_MACRO
+ IF(GG%IPAR(2,ELEM).EQ.NODE_MACRO(J)) GO TO 10
+ ENDDO
+ GO TO 20
+ 10 DO J=1,NBNODE_MACRO
+ IF(GG%IPAR(3,ELEM).EQ.NODE_MACRO(J)) CYCLE ITER0
+ ENDDO
+ 20 IF(GG%IPAR(1,ELEM)==1) THEN
+ X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
+ X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
+ DO J=1,NPERIM
+ ALIGN(3,1,J)=X1; ALIGN(3,2,J)=Y1;
+ DET1 = DET_ROSETTA(ALIGN(1,1,J),3)
+ ALIGN(3,1,J)=X2; ALIGN(3,2,J)=Y2;
+ DET2 = DET_ROSETTA(ALIGN(1,1,J),3)
+ IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) THEN
+ PERIM_MACRO(I) = J
+ CYCLE ITER0
+ ENDIF
+ ENDDO
+ NPERIM=NPERIM+1
+ PERIM_MACRO(I) = NPERIM
+ ANGLE(NPERIM)=ATAN((Y2-Y1)/(X2-X1))
+ IF(ABS(ANGLE(NPERIM)).LE.1.0E-5) ANGLE(NPERIM)=0.0
+ ALIGN(1,1,NPERIM)=X1; ALIGN(1,2,NPERIM)=Y1
+ ALIGN(2,1,NPERIM)=X2; ALIGN(2,2,NPERIM)=Y2
+ ! Recover albedo from global geometry
+ ALBEDO(NPERIM)=1.0
+ DO IB=1,GG%NBBCDA
+ J = FINDLC(GG%BCDATAREAD(IB)%ELEMNB,ELEM)
+ IF(J.EQ.1) THEN
+ ALBEDO(NPERIM)=GG%BCDATAREAD(IB)%BCDATA(6)
+ EXIT
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO ITER0
+ !----
+ ! Printouts
+ !----
+ IF(IPRINT.GT.2) THEN
+ WRITE(FOUT,'(/39H MUSACG: IPAR VALUES IN GLOBAL GEOMETRY)')
+ WRITE(FOUT,'(5H ELEM,20I5/(5X,20I5))') ELEM_MACRO(:NELEM_MACRO)
+ WRITE(FOUT,'(5H TYPE,20I5/(5X,20I5))') GG%IPAR(1,ELEM_MACRO(:NELEM_MACRO))
+ WRITE(FOUT,'(5H -,20I5/(5X,20I5))') GG%IPAR(2,ELEM_MACRO(:NELEM_MACRO))
+ WRITE(FOUT,'(5H +,20I5/(5X,20I5))') GG%IPAR(3,ELEM_MACRO(:NELEM_MACRO))
+ ENDIF
+ !----
+ ! Create volume and merge information
+ !----
+ ALLOCATE(GG_MAC%VOL_NODE(NBNODE_MACRO),GG_MAC%MED(NBNODE_MACRO),GG_MAC%NAME_MACRO(1), &
+ & GG_MAC%NUM_MERGE(NBNODE_MACRO), STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory I,R')
+ GG_MAC%NAME_MACRO(1) = 'MACR000001'
+ DO I=1,NBNODE_MACRO
+ J = NODE_MACRO(I)
+ GG_MAC%VOL_NODE(I) = GG%VOL_NODE(J)
+ GG_MAC%MED(I) = GG%MED(J)
+ GG_MAC%NUM_MERGE(I) = GG%NUM_MERGE(J)
+ ENDDO
+ GG_MAC%NB_FLUX = GG%NB_FLUX
+ DO I=GG%NB_FLUX,1,-1
+ J = FINDLC(GG_MAC%NUM_MERGE,I)
+ IF(J.GT.0) CYCLE
+ DO J=1,NBNODE_MACRO
+ IF(GG_MAC%NUM_MERGE(J).GT.I) GG_MAC%NUM_MERGE(J)=GG_MAC%NUM_MERGE(J)-1
+ ENDDO
+ GG_MAC%NB_FLUX = GG_MAC%NB_FLUX-1
+ ENDDO
+ IF(IPRINT.GT.1) THEN
+ WRITE(FOUT,'(/32H MUSACG: number of flux in macro,I6.6,1H=,I5)') IMACRO,GG_MAC%NB_FLUX
+ WRITE(FOUT,'(5H NODE,20I5/(5X,20I5))') (I,I=1,NBNODE_MACRO)
+ WRITE(FOUT,'(5H MERG,20I5/(5X,20I5))') (GG_MAC%NUM_MERGE(I),I=1,NBNODE_MACRO)
+ ENDIF
+ ALLOCATE(GG_MAC%NUM_MACRO(GG_MAC%NB_FLUX), STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory NUM_MACRO')
+ GG_MAC%NUM_MACRO(:GG_MAC%NB_FLUX) = 1
+ !----
+ ! Fill the GG_MAC tracking structure
+ !----
+ GG_MAC%NBBCDA = NPERIM
+ GG_MAC%NB_ELEM = NELEM_MACRO
+ GG_MAC%NB_NODE = NBNODE_MACRO
+ GG_MAC%NB_MACRO = 1
+ GG_MAC%DEFAUL = GG%DEFAUL
+ ALLOCATE(GG_MAC%IPAR(NIPAR,NELEM_MACRO),GG_MAC%RPAR(NRPAR,NELEM_MACRO), STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory I,R')
+ GG_MAC%IPAR(:3,:GG_MAC%NB_ELEM)=0
+ GG_MAC%RPAR(:3,:GG_MAC%NB_ELEM)=0.0
+ DO ELEM=1, GG_MAC%NB_ELEM
+ GG_MAC%IPAR(1,ELEM) = GG%IPAR(1,ELEM_MACRO(ELEM)) ! elem type
+ I = GG%IPAR(2,ELEM_MACRO(ELEM)) ! node-
+ J = FINDLC(NODE_MACRO,I)
+ IF(J.GT.0) THEN
+ GG_MAC%IPAR(2,ELEM) = J
+ ELSE
+ GG_MAC%IPAR(3,ELEM) = MIN(I,0)
+ ENDIF
+ I = GG%IPAR(3,ELEM_MACRO(ELEM)) ! node+
+ J = FINDLC(NODE_MACRO,I)
+ IF(J.GT.0) THEN
+ GG_MAC%IPAR(3,ELEM) = J
+ ELSE
+ GG_MAC%IPAR(3,ELEM) = MIN(I,0)
+ ENDIF
+ GG_MAC%RPAR(:NRPAR,ELEM) = GG%RPAR(:NRPAR,ELEM_MACRO(ELEM))
+ ENDDO
+ IF(IPRINT.GT.1) THEN
+ WRITE(FOUT,'(/38H MUSACG: IPAR VALUES IN MACRO GEOMETRY)')
+ WRITE(FOUT,'(5H ELEM,20I5/(5X,20I5))') (I,I=1,GG_MAC%NB_ELEM)
+ WRITE(FOUT,'(5H TYPE,20I5/(5X,20I5))') GG_MAC%IPAR(1,:GG_MAC%NB_ELEM)
+ WRITE(FOUT,'(5H -,20I5/(5X,20I5))') GG_MAC%IPAR(2,:GG_MAC%NB_ELEM)
+ WRITE(FOUT,'(5H +,20I5/(5X,20I5))') GG_MAC%IPAR(3,:GG_MAC%NB_ELEM)
+ ENDIF
+ GG_MAC%ALBEDO=GG%ALBEDO
+ ALLOCATE(GG_MAC%BCDATAREAD(NPERIM))
+ DO IB=1, GG_MAC%NBBCDA
+ GG_MAC%BCDATAREAD(IB)%NBER = COUNT(PERIM_MACRO(:NELEM_MACRO) == IB)
+ ALLOCATE(GG_MAC%BCDATAREAD(IB)%ELEMNB(GG_MAC%BCDATAREAD(IB)%NBER))
+ J=0
+ DO I=1,NELEM_MACRO
+ IF(PERIM_MACRO(I) == IB) THEN
+ J=J+1
+ GG_MAC%BCDATAREAD(IB)%ELEMNB(J) = I
+ ENDIF
+ ENDDO
+ GG_MAC%BCDATAREAD(IB)%BCDATA(1) = ALIGN(1,1,IB)
+ GG_MAC%BCDATAREAD(IB)%BCDATA(2) = ALIGN(1,2,IB)
+ GG_MAC%BCDATAREAD(IB)%BCDATA(3) = COS(ANGLE(IB))
+ GG_MAC%BCDATAREAD(IB)%BCDATA(4) = SIN(ANGLE(IB))
+ GG_MAC%BCDATAREAD(IB)%BCDATA(5) = ANGLE(IB)
+ GG_MAC%BCDATAREAD(IB)%BCDATA(6) = ALBEDO(IB)
+ GG_MAC%BCDATAREAD(IB)%SALTYPE = 0
+ ENDDO
+ TYPGEO=0 ; NBFOLD=0
+ !----
+ ! Find if the perimeter contains an unfolding axis
+ !----
+ ALLOCATE(IFOLD_GLOB((2**NPERIM)*GG_MAC%NB_ELEM))
+ IFOLD_GLOB(:GG_MAC%NB_ELEM)=(/ (ELEM, ELEM=1,GG_MAC%NB_ELEM) /)
+ DO IPASS=1,2 ! two passes are required to get rid of unfolding axis
+ IF(IPRINT.GT.2) THEN
+ WRITE(FOUT,'(/45H MUSACG: PERIMETERS BEFORE UNFOLDING AT PASS=,I2)') IPASS
+ WRITE(FOUT,'(7H PERIM,10I12/(7X,10I12))') (IB,IB=1,GG_MAC%NBBCDA)
+ WRITE(FOUT,'(7H X,1P,10e12.4/(7X,10e12.4))') GG_MAC%BCDATAREAD(:GG_MAC%NBBCDA)%BCDATA(1)
+ WRITE(FOUT,'(7H Y,1P,10e12.4/(7X,10e12.4))') GG_MAC%BCDATAREAD(:GG_MAC%NBBCDA)%BCDATA(2)
+ WRITE(FOUT,'(7H ANGLE,1P,10e12.4/(7X,10e12.4))') GG_MAC%BCDATAREAD(:GG_MAC%NBBCDA)%BCDATA(5)*180._PDB/PI
+ WRITE(FOUT,'(7H ALBEDO,1P,10e12.4/(7X,10e12.4))') GG_MAC%BCDATAREAD(:GG_MAC%NBBCDA)%BCDATA(6)
+ ENDIF
+ OUT1: DO IB=1,GG_MAC%NBBCDA
+ ALIGN(:3,:3,IB)=1.0D0
+ DO I=1,GG_MAC%BCDATAREAD(IB)%NBER
+ INDBC=GG_MAC%BCDATAREAD(IB)%ELEMNB(I)
+ IF(INDBC==0) CYCLE
+ X1=GG_MAC%RPAR(1,INDBC); Y1=GG_MAC%RPAR(2,INDBC)
+ X2=X1+GG_MAC%RPAR(3,INDBC); Y2=Y1+GG_MAC%RPAR(4,INDBC)
+ ALIGN(1,1,IB)=X1; ALIGN(1,2,IB)=Y1
+ ALIGN(2,1,IB)=X2; ALIGN(2,2,IB)=Y2
+ EXIT
+ ENDDO
+ LFOLD(IB)=.FALSE.
+ DO ELEM=1,GG_MAC%NB_ELEM
+ IF(GG_MAC%IPAR(1,ELEM)==3) THEN
+ CALL SAL141(3,GG_MAC%RPAR(:,ELEM),X1,Y1,1)
+ CALL SAL141(3,GG_MAC%RPAR(:,ELEM),X2,Y2,2)
+ ALIGN(3,1,IB)=X1; ALIGN(3,2,IB)=Y1;
+ DET1 = DET_ROSETTA(ALIGN(1,1,IB),3)
+ ALIGN(3,1,IB)=X2; ALIGN(3,2,IB)=Y2;
+ DET2 = DET_ROSETTA(ALIGN(1,1,IB),3)
+ LFOLD(IB)=((ABS(DET1).LE.1.0E-4).OR.(ABS(DET2).LE.1.0E-4))
+ IF(LFOLD(IB)) CYCLE OUT1
+ ENDIF
+ ENDDO
+ ENDDO OUT1
+ IF(IPRINT.GT.2) THEN
+ WRITE(FOUT,'(7H UNFOLD,1P,10L12/(6X,10L12))') LFOLD(:GG_MAC%NBBCDA)
+ ENDIF
+ !----
+ ! Unfold macro geometry (many times, if required)
+ !----
+ DO IB=1,GG_MAC%NBBCDA
+ IF(LFOLD(IB)) THEN
+ ALLOCATE(IFOLD(2*GG_MAC%NB_ELEM))
+ CALL SALFOLD_0(GG_MAC,IPASS,IB,GG_MAC%NBBCDA,ALIGN,LFOLD,IFOLD)
+ DO ELEM=1,GG_MAC%NB_ELEM
+ IF(IFOLD(ELEM).GT.SIZE(IFOLD_GLOB,1)) CALL XABORT('MUSACG: IFOLD overflow')
+ IFOLD(ELEM)=IFOLD_GLOB(IFOLD(ELEM))
+ ENDDO
+ IFOLD_GLOB(:GG_MAC%NB_ELEM)=IFOLD(:GG_MAC%NB_ELEM)
+ DEALLOCATE(IFOLD)
+ ENDIF
+ ENDDO
+ ENDDO
+ IFOLD_GLOB(:GG_MAC%NB_ELEM)=ELEM_MACRO(IFOLD_GLOB(:GG_MAC%NB_ELEM))
+ IF(PRTIND>2) WRITE(6,'(/15H MUSACG: IFOLD=,20I5/(15X,20I5))') IFOLD_GLOB(:GG_MAC%NB_ELEM)
+ DEALLOCATE(ALIGN,LFOLD)
+ IF(IPRINT>0) WRITE(FOUT,*) 'MUSACG: after unfolding -- NB_ELEM=',GG_MAC%NB_ELEM,' NB_PERIM=',GG_MAC%NBBCDA
+ IF(PRTIND>5) THEN
+ !* print surfacic file
+ WRITE(FOUT,'(5H--cut,70(1H-),I5)') IMACRO
+ WRITE(FOUT,'(5HBEGIN)')
+ WRITE(FOUT,'(42H* typgeo nbfold nbnode nbelem nbmacr nbreg)')
+ WRITE(FOUT,'(6I7)') TYPGEO,NBFOLD,GG_MAC%NB_NODE,GG_MAC%NB_ELEM,GG_MAC%NB_MACRO,GG_MAC%NB_NODE
+ WRITE(FOUT,'(20H* index kndex prec)')
+ WRITE(FOUT,'(4I7)') 0,0,1
+ WRITE(FOUT,'(18H* eps eps0)')
+ WRITE(FOUT,'(1P,2E18.9)') 1.0E-03,1.0E-05
+ WRITE(FOUT,'(20H* num_of_region/mesh)')
+ WRITE(FOUT,'(10I7)') (GG_MAC%NUM_MERGE(I),I=1,GG_MAC%NB_NODE)
+ WRITE(FOUT,'(13H* macro names)')
+ WRITE(FOUT,'(4(3x,a10,2x))') (GG_MAC%NAME_MACRO(I),I=1,GG_MAC%NB_MACRO)
+ WRITE(FOUT,'(35H* macro_order_index_per_flux_region)')
+ WRITE(FOUT,'(10I7)') (GG_MAC%NUM_MACRO(I),I=1,GG_MAC%NB_FLUX)
+ DO ELEM=1,GG_MAC%NB_ELEM
+ TYPE=GG_MAC%IPAR(1,ELEM)
+ WRITE(FOUT,'(7h elem =,I6)') ELEM
+ WRITE(FOUT,'(22H*type node- node+)')
+ WRITE(FOUT,'(3I6)') (GG_MAC%IPAR(I,ELEM),I=1,3)
+ WRITE(FOUT,'(63H*cx cy ex_or_R ey_or_theta1 theta2)')
+ IF(TYPE<=2) THEN
+ WRITE(FOUT,'(1P,5E18.9)') (GG_MAC%RPAR(I,ELEM),I=1,5)
+ ELSE IF(TYPE==3) THEN
+ WRITE(FOUT,'(1P,5E18.9)') (GG_MAC%RPAR(I,ELEM),I=1,3),GG_MAC%RPAR(4,ELEM)*180._PDB/PI, &
+ (GG_MAC%RPAR(5,ELEM)-GG_MAC%RPAR(4,ELEM))*180._PDB/PI
+ ENDIF
+ ENDDO
+ WRITE(FOUT,'(40H*defaul nbbcda allsur divsur ndivsur)')
+ WRITE(FOUT,'(1P,5I8)') GG_MAC%DEFAUL,GG_MAC%NBBCDA,ALLSUR,0,0
+ WRITE(FOUT,'(17H*albedo deltasur)')
+ WRITE(FOUT,'(1P,2E18.9)') GG_MAC%ALBEDO,0.0
+ DO IB=1,GG_MAC%NBBCDA
+ WRITE(FOUT,'(37H particular boundary condition number,i12)') IB
+ WRITE(FOUT,'(13H*type nber)')
+ WRITE(FOUT,'(1P,2I8)') GG_MAC%BCDATAREAD(IB)%SALTYPE,GG_MAC%BCDATAREAD(IB)%NBER
+ WRITE(FOUT,'(14H*elems(1,nber))')
+ WRITE(FOUT,'(1P,10I8)') (GG_MAC%BCDATAREAD(IB)%ELEMNB(I),I=1,GG_MAC%BCDATAREAD(IB)%NBER)
+ IF(GG_MAC%BCDATAREAD(IB)%SALTYPE/=0) CALL XABORT('MUSACG: SALTYPE=0 expected')
+ WRITE(FOUT,'(7H*albedo)')
+ WRITE(FOUT,'(1P,E18.9)') GG_MAC%BCDATAREAD(IB)%BCDATA(6)
+ ENDDO
+ WRITE(FOUT,'(12H* mil(nbreg))')
+ WRITE(FOUT,'(10I7)') (GG_MAC%MED(I),I=1,GG_MAC%NB_NODE)
+ WRITE(FOUT,'(3HEND)')
+ WRITE(FOUT,'(5H--cut,70(1H-),I5)') IMACRO
+ ENDIF
+ IF(GG_MAC%NBBCDA.GT.6) CALL XABORT('MUSACG: The unfolded geometry has more than 6 perimeters')
+ !****
+ !* compute node perimeters for the macro
+ ALLOCATE (GG_MAC%PPERIM_NODE(GG_MAC%NB_NODE+1),STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory PPERIM_NODE')
+ ALLOCATE(AUX_ARR(2*GG_MAC%NB_ELEM))
+ CALL SAL130_2(GG_MAC%NB_ELEM,GG_MAC%NB_NODE,GG_MAC%IPAR,GG_MAC%PPERIM_NODE, &
+ GG_MAC%PERIM_NODE,AUX_ARR)
+ !
+ !* - compute number of bc's per 2D macro,
+ ! NB_BC2 counts total nber of 2D bc's
+ ! - compute IBC2_ELEM, keep relative 2D bc nber to elements
+ ! - allocation : 2D bc structures
+ ! 2D perimeter structure for a macro
+ ! - get list of elements in 2d macro perimeter
+ ALLOCATE(GG_MAC%IBC2_ELEM(GG_MAC%NB_ELEM),GG_MAC%ISURF2_ELEM(GG_MAC%NB_ELEM))
+ CALL SAL130_4(GG_MAC%NB_ELEM,NN,GG_MAC%IPAR,GG_MAC%IBC2_ELEM,AUX_ARR)
+ GG_MAC%NB_BC2 = NN
+ GG_MAC%NALBG = GG_MAC%NBBCDA
+ !* allocate bcdata
+ ALLOCATE (GG_MAC%BCDATA(6,GG_MAC%NALBG), STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory R')
+ DO IB=1,GG_MAC%NALBG
+ GG_MAC%BCDATA(:6,IB)=GG_MAC%BCDATAREAD(IB)%BCDATA(:6)
+ ENDDO
+ !
+ !* put default value in all bc elements:
+ ALLOCATE(GG_MAC%TYPE_BC2(NN),GG_MAC%IDATA_BC2(NN),GG_MAC%PERIM_MAC2(NN),STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory I,R')
+ GG_MAC%PERIM_MAC2(1:NN) = AUX_ARR(1:NN)
+ GG_MAC%NPERIM_MAC2 = NN
+ DEALLOCATE(AUX_ARR)
+ GG_MAC%TYPE_BC2(:NN) = 0
+ CALL SAL131_2(GG_MAC%NB_ELEM,GG_MAC%DEFAUL,GG_MAC%IPAR,GG_MAC%IBC2_ELEM,GG_MAC%TYPE_BC2, &
+ & GG_MAC%IDATA_BC2)
+ ITBC=0
+ IF(GG_MAC%NBBCDA>0)THEN
+ DO I=1,GG%NBBCDA
+ ITBC=ITBC+1
+ TYPE=GG_MAC%BCDATAREAD(I)%SALTYPE
+ IF(TYPE.NE.0) CALL XABORT('MUSACG: TYPE=0 EXPECTED.')
+ !
+ ! modify notation for boundary conditions
+ NBER=GG_MAC%BCDATAREAD(I)%NBER
+ DO J=1,NBER
+ ELEM=GG_MAC%BCDATAREAD(I)%ELEMNB(J)
+ IF(ELEM>GG_MAC%NB_ELEM.OR.ELEM<=0) CALL XABORT('MUSACG: unknown bc element')
+ ! get local surface nber
+ IB=GG_MAC%IBC2_ELEM(ELEM)
+ LGBC=GG_MAC%IPAR(2,ELEM)<=0
+ II=0
+ IF(LGBC)THEN
+ II=2
+ ELSE
+ LGBC=GG_MAC%IPAR(3,ELEM)<=0
+ IF(LGBC) II=3
+ ENDIF
+ IF(.NOT.LGBC) THEN
+ WRITE(*,*) 'elem :',ELEM
+ WRITE(*,*) 'GG_MAC%IPAR(:,ELEM) :',GG_MAC%IPAR(:,ELEM)
+ CALL XABORT('MUSACG: wrong bc element')
+ ENDIF
+ ! put bc type
+ GG_MAC%IPAR(II,ELEM)=G_BC_TYPE(TYPE)
+ GG_MAC%TYPE_BC2(IB)=G_BC_TYPE(TYPE)
+ ! put bc data position :
+ GG_MAC%IDATA_BC2(IB)=ITBC
+ ENDDO
+ ENDDO
+ ENDIF
+ !
+ !* - set BCDATA position for surfaces of type G_BC_TYPE(-1)
+ ! - compute the nber of surfaces (type -1,0,-12,-13,-14,-15) : nbsur2
+ ! - allocate structures for the surfaces
+ ! - compute surf_mac2
+ ALLOCATE(AUX_ARR(GG_MAC%NB_ELEM))
+ AUX_ARR(:GG_MAC%NB_ELEM)=0
+ GG_MAC%ISURF2_ELEM(:GG_MAC%NB_ELEM)=0
+ NELEM_MACRO=GG_MAC%NB_ELEM
+ NSURF_MACRO=0
+ DO IB=1,GG_MAC%NB_BC2
+ ! relative element nber
+ IELEM=GG_MAC%PERIM_MAC2(IB)
+ ! count 2D surfaces number
+ IF(GG_MAC%TYPE_BC2(IB)==G_BC_TYPE(-1) .OR. GG_MAC%TYPE_BC2(IB)==G_BC_TYPE(0) .OR. &
+ & GG_MAC%TYPE_BC2(IB)==G_BC_TYPE(1)) THEN
+ NSURF_MACRO=NSURF_MACRO+1
+ AUX_ARR(NSURF_MACRO)=IB
+ GG_MAC%ISURF2_ELEM(IELEM)=NSURF_MACRO
+ ELSE
+ GG_MAC%ISURF2_ELEM(IELEM)=0
+ ENDIF
+ ENDDO
+ GG_MAC%NB_SURF2 = NSURF_MACRO
+ ALLOCATE(SURF_MACRO(NSURF_MACRO))
+ IF(NSURF_MACRO>0) THEN
+ DO I=1,NSURF_MACRO
+ ELEM = FINDLC(GG_MAC%ISURF2_ELEM,I)
+ IF(ELEM.GT.NELEM_MACRO) CALL XABORT('MUSACG: ELEM_MACRO OVERFLOW.')
+ SURF_MACRO(I) = IFOLD_GLOB(ELEM)
+ ENDDO
+ ALLOCATE (GG_MAC%IBC2_SURF2(NSURF_MACRO),GG_MAC%IELEM_SURF2(NSURF_MACRO),STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: NOT ENOUGH MEMORY I,R')
+ GG_MAC%IBC2_SURF2(1:NSURF_MACRO)=AUX_ARR(1:NSURF_MACRO)
+ !
+ ! - define IELEM_SURF2 ???
+ ALLOCATE(GG_MAC%SURF2(GG_MAC%NB_SURF2),STAT = OK)
+ IF(OK /= 0) CALL XABORT('MUSACG: not enough memory I,R')
+ CALL SAL130_6(GG_MAC%NB_SURF2,GG_MAC%IBC2_SURF2,GG_MAC%PERIM_MAC2, &
+ & GG_MAC%IELEM_SURF2)
+ ELSE
+ NULLIFY(GG_MAC%IBC2_SURF2,GG_MAC%IELEM_SURF2)
+ ENDIF
+ DO I=1,NBNODE_MACRO
+ NODE_MACRO(I)=GG%NUM_MERGE(NODE_MACRO(I))
+ ENDDO
+ CALL LCMSIX(ITRACK,'SURFACIC_TMP',1)
+ IF(NSURF_MACRO>0) CALL LCMPUT(ITRACK,'SURF_MACRO',NSURF_MACRO,1,SURF_MACRO)
+ CALL LCMPUT(ITRACK,'MERGE_MACRO',NBNODE_MACRO,1,NODE_MACRO)
+ CALL LCMSIX(ITRACK,' ',2)
+ DEALLOCATE(IFOLD_GLOB,ELEM_MACRO,NODE_MACRO)
+ DEALLOCATE(AUX_ARR,ALBEDO,ANGLE,PERIM_MACRO,SURF_MACRO)
+ !
+ !* topological check
+ ALLOCATE (GG_MAC%VOL_NODE(GG_MAC%NB_NODE), STAT=OK)
+ IF(OK/=0) CALL XABORT('MUSACG: not enough memory VOL')
+ CALL SAL140(GG_MAC%NB_NODE,GG_MAC%RPAR,GG_MAC%IPAR,GG_MAC%PPERIM_NODE,GG_MAC%PERIM_NODE)
+ !
+ !* volumes, surfaces, put local nbers in node, and read media:
+ CALL SAL160_2(GG_MAC%NB_ELEM,GG_MAC%IPAR,GG_MAC%RPAR,GG_MAC%VOL_NODE,GG_MAC%ISURF2_ELEM, &
+ GG_MAC%NB_SURF2,GG_MAC%SURF2)
+ !
+ !* printout basic domain
+ CALL SAL170(GG_MAC)
+ !----
+ ! Save MUST tracking information on LCM
+ !----
+ CALL LCMSIX(ITRACK,'GEOMETRY ',1)
+ CALL LCMPUT(ITRACK,'NB_ELEM ',1,1,GG_MAC%NB_ELEM)
+ CALL LCMPUT(ITRACK,'NIPAR ',1,1,SIZE(GG_MAC%IPAR,1))
+ CALL LCMPUT(ITRACK,'IPAR ',SIZE(GG_MAC%IPAR),1,GG_MAC%IPAR)
+ CALL LCMPUT(ITRACK,'RPAR ',SIZE(GG_MAC%RPAR),4,GG_MAC%RPAR)
+ CALL LCMPUT(ITRACK,'ISURF2_ELEM ',SIZE(GG_MAC%ISURF2_ELEM),1,GG_MAC%ISURF2_ELEM)
+ CALL LCMPUT(ITRACK,'NB_NODE ',1,1,GG_MAC%NB_NODE)
+ CALL LCMPUT(ITRACK,'VOL_NOD ',GG_MAC%NB_NODE,4,GG_MAC%VOL_NODE)
+ CALL LCMPUT(ITRACK,'NB_SURF2 ',1,1,GG_MAC%NB_SURF2)
+ IF(GG_MAC%NBBCDA.GT.0) THEN
+ LGINF = .TRUE.
+ DO IB=1, GG_MAC%NBBCDA
+ LGINF = LGINF .AND. (GG_MAC%BCDATAREAD(IB)%BCDATA(6) == 1._PDB)
+ ENDDO
+ ELSE
+ LGINF = (GG_MAC%ALBEDO == 1._PDB)
+ ENDIF
+ IF(GG_MAC%NB_SURF2 > 0) THEN
+ CALL LCMPUT(ITRACK,'IBC2_SURF2 ',SIZE(GG_MAC%IBC2_SURF2),1,GG_MAC%IBC2_SURF2)
+ CALL LCMPUT(ITRACK,'IELEM_SURF2 ',SIZE(GG_MAC%IELEM_SURF2),1,GG_MAC%IELEM_SURF2)
+ CALL LCMPUT(ITRACK,'SURF2 ',SIZE(GG_MAC%SURF2),4,GG_MAC%SURF2)
+ ENDIF
+ CALL LCMPUT(ITRACK,'NPERIM_MAC2 ',1,1,GG_MAC%NPERIM_MAC2)
+ CALL LCMPUT(ITRACK,'PERIM_MAC2 ',SIZE(GG_MAC%PERIM_MAC2),1,GG_MAC%PERIM_MAC2)
+ CALL LCMPUT(ITRACK,'PERIM_NODE ',SIZE(GG_MAC%PERIM_NODE),1,GG_MAC%PERIM_NODE)
+ CALL LCMPUT(ITRACK,'PPERIM_NODE ',SIZE(GG_MAC%PPERIM_NODE),1,GG_MAC%PPERIM_NODE)
+ CALL LCMPUT(ITRACK,'BC_DATA_DIM2',1,1,SIZE(GG_MAC%BCDATA,2))
+ CALL LCMPUT(ITRACK,'NB_BC2 ',1,1,GG_MAC%NB_BC2)
+ CALL LCMSIX(ITRACK,' ',2) ! come back to father directory
+ !----
+ ! Print tracking object directory
+ !----
+ IF(IPRINT.GT.4) THEN
+ WRITE(FOUT,'(/14H MUSACG: MACRO,I6.6,20H GEOMETRY DIRECTORY:)') IMACRO
+ CALL LCMLIB(ITRACK)
+ CALL LCMSIX(ITRACK,'GEOMETRY',1)
+ CALL LCMLIB(ITRACK)
+ CALL LCMSIX(ITRACK,' ',2)
+ ENDIF
+ !----
+ ! store the STATE VECTOR
+ !----
+ NREG=MAXVAL(GG_MAC%NUM_MERGE)
+ LEAK=1
+ IF(.NOT.LGINF) LEAK=0 ! reset the leakage flag
+ CALL LCMGET(ITRACK,'STATE-VECTOR',I_STATE)
+ I_STATE(1) = NREG ! number of regions
+ I_STATE(2) = NREG ! number of unknowns in DRAGON
+ I_STATE(3) = LEAK ! 1 = absent leakage, 0 leakage
+ I_STATE(4) = MAXVAL(GG_MAC%MED(1:GG_MAC%NB_NODE)) ! maximum number of mixture
+ I_STATE(5) = GG_MAC%NB_SURF2 ! number of outer surface
+ I_STATE(9) = ISPEC
+ I_STATE(24)= 0
+ NSOUT=GG_MAC%NB_SURF2
+ CALL LCMPUT(ITRACK,'STATE-VECTOR',NSTATE,1,I_STATE)
+ !
+ ! fill-in medium number per region
+ ALLOCATE(ITAB(NREG),VOLUME(NREG), STAT =OK)
+ IF(OK /= 0) CALL XABORT('MUSACG: failure to allocate integer ITAB')
+ ! fill in MATCOD
+ DO J=1,GG_MAC%NB_NODE
+ ITAB(GG_MAC%NUM_MERGE(J)) = GG_MAC%MED(J)
+ ENDDO
+ CALL LCMPUT(ITRACK,'MATCOD',NREG,1,ITAB(1:NREG) )
+ ! fill-in KEYFLX per region
+ ITAB(1:NREG)=(/ (I,I=1,NREG) /)
+ CALL LCMPUT(ITRACK,'MERGE',NREG,1,ITAB)
+ CALL LCMPUT(ITRACK,'KEYFLX',NREG,1,ITAB)
+ ! fill-in volumes per region
+ VOLUME(:NREG) =0.
+ DO I=1,GG_MAC%NB_NODE
+ VOLUME(GG_MAC%NUM_MERGE(I)) = VOLUME(GG_MAC%NUM_MERGE(I)) + REAL(GG_MAC%VOL_NODE(I))
+ ENDDO
+ CALL LCMPUT(ITRACK,'VOLUME',NREG,2,VOLUME)
+ DEALLOCATE(VOLUME,ITAB)
+
+ ! useful values in SAL_TRACKING_TYPES module
+ NFREG=GG_MAC%NB_NODE
+ CALL LCMSIX(ITRACK,'NXTRecords',1)
+ DGMESHX=(/ 1.E10_PDB , -1.E10_PDB /)
+ DGMESHY=(/ 1.E10_PDB , -1.E10_PDB /)
+ DO ELEM=1,GG_MAC%NB_ELEM
+ DGMESHX(1)=MIN(DGMESHX(1),GG_MAC%RPAR(1,ELEM))
+ DGMESHX(2)=MAX(DGMESHX(2),GG_MAC%RPAR(1,ELEM))
+ DGMESHY(1)=MIN(DGMESHY(1),GG_MAC%RPAR(2,ELEM))
+ DGMESHY(2)=MAX(DGMESHY(2),GG_MAC%RPAR(2,ELEM))
+ ENDDO
+ CALL LCMPUT(ITRACK,'G00000001SMX',2,4,DGMESHX)
+ CALL LCMPUT(ITRACK,'G00000001SMY',2,4,DGMESHY)
+ IEDIMG(:NSTATE)=0
+ IEDIMG(1)=NDIM
+ IEDIMG(2)=0 ! Cartesian geometry
+ IF(TYPGEO.EQ.8) IEDIMG(2)=2 ! Isocel geometry with specular reflection
+ IF(TYPGEO.EQ.9) IEDIMG(2)=3 ! Hexagonal geometry with translation
+ IF(TYPGEO.EQ.10) IEDIMG(2)=4 ! Isocel geometry with RA60 symmetry
+ IF(TYPGEO.EQ.11) IEDIMG(2)=5 ! Lozenge geometry with R120 rotation
+ IF(TYPGEO.EQ.12) IEDIMG(2)=6 ! S30 geometry with specular reflection
+ IEDIMG(5)=1 ! 1 cellule
+ IEDIMG(13)=1 ! 1 cellule
+ IEDIMG(14)=1 ! 1 cellule
+ IEDIMG(22)=NSOUT ! number of external surfaces for this geometry
+ IEDIMG(23)=NFREG ! number of regions for this geometry
+ IEDIMG(25)=GG_MAC%NB_NODE
+ CALL LCMPUT(ITRACK,'G00000001DIM',NSTATE,1,IEDIMG)
+ CALL LCMSIX(ITRACK,' ',2) ! come back to father directory
+ !----
+ ! process boundary conditions
+ !----
+ IF(IPRINT.GT.0) WRITE(FOUT,*) 'number of merged regions,surfaces,nodes',NREG,NSOUT,NFREG
+ ALLOCATE(MATALB(-NSOUT:NFREG),VOLSUR(-NSOUT:NFREG),KEYMRG(-NSOUT:NFREG))
+ CALL LCMGET(ITRACK,'MATCOD',MATALB(1))
+ ALLOCATE(VOLUME(NREG))
+ CALL LCMGET(ITRACK,'VOLUME',VOLUME)
+ VOLSUR(1:NREG)=VOLUME(:NREG)
+ DEALLOCATE(VOLUME)
+ ! boundary conditions structures
+ ALLOCATE(ICODE(GG_MAC%NALBG),GALBED(GG_MAC%NALBG))
+ ICODE(1:GG_MAC%NALBG)=(/ (-I,I=1,GG_MAC%NALBG) /)
+ GALBED(:GG_MAC%NALBG)=REAL(GG_MAC%ALBEDO)
+ DO I=1,NSOUT
+ KEYMRG(-I)=-I
+ VOLSUR(-I)=GG_MAC%SURF2(I)
+ INDEX=GG_MAC%IDATA_BC2(GG_MAC%IBC2_SURF2(I))
+ IF(INDEX.EQ.0) THEN
+ ! Use the default albedo
+ MATALB(-I)=-1
+ GALBED(1)=REAL(GG_MAC%ALBEDO)
+ ELSE
+ IF(INDEX.GT.6) CALL XABORT('MUSACG: SDIRE overflow.')
+ IF(INDEX > GG_MAC%NALBG) THEN
+ CALL XABORT('MUSACG: Albedo array overflow(2).')
+ ENDIF
+ MATALB(-I)=-INDEX
+ IF(SIZE(GG_MAC%BCDATA) > 0) THEN
+ GALBED(INDEX)=REAL(GG_MAC%BCDATA(6,INDEX))
+ ELSE
+ GALBED(INDEX)=REAL(GG_MAC%ALBEDO)
+ ENDIF
+ ENDIF
+ ENDDO
+ MATALB(0)=0
+ KEYMRG(0)=0
+ VOLSUR(0)=0._PDB
+ KEYMRG(1:NREG)=(/ (I,I=1,NREG) /)
+ !
+ IF(IPRINT.GT.1) THEN
+ CALL PRINDM('VOLUME',VOLSUR(-NSOUT),NREG+NSOUT+1)
+ CALL PRINIM('MATALB',MATALB(-NSOUT),NREG+NSOUT+1)
+ CALL PRINIM('KEYMRG',KEYMRG(-NSOUT),NREG+NSOUT+1)
+ ENDIF
+ IF(IPRINT.GT.0) THEN
+ CALL PRINIM('ICODE ',ICODE(1),GG_MAC%NALBG)
+ CALL PRINAM('GALBED',GALBED(1),GG_MAC%NALBG)
+ ENDIF
+ !----
+ ! fill in tracking LCM object in excelt format
+ !----
+ TEXT72='SAL TRACKING'
+ CALL LCMPTC(ITRACK,'TITLE',72,TEXT72)
+ CALL LCMPUT(ITRACK,'ICODE',GG_MAC%NALBG,1,ICODE)
+ CALL LCMSIX(ITRACK,'NXTRecords',1)
+ CALL LCMPUT(ITRACK,'SAreaRvolume',NREG+NSOUT+1,4,VOLSUR(-NSOUT))
+ CALL LCMPUT(ITRACK,'MATALB',NREG+NSOUT+1,1,MATALB(-NSOUT))
+ CALL LCMPUT(ITRACK,'KEYMRG',NREG+NSOUT+1,1,KEYMRG(-NSOUT))
+ CALL LCMSIX(ITRACK,' ',2)
+ IF(NSOUT>0) THEN
+ ALLOCATE(IBC(NSOUT))
+ IBC(1:NSOUT)=(/ (I,I=1,NSOUT) /)
+ CALL LCMPUT(ITRACK,'BC-REFL+TRAN',NSOUT,1,IBC)
+ DEALLOCATE(IBC)
+ ENDIF
+ ALLOCATE(PERIM_SURF(GG_MAC%NB_ELEM))
+ PERIM_SURF(:GG_MAC%NB_ELEM)=0
+ DO IB=1,GG_MAC%NBBCDA
+ DO I=1,GG_MAC%BCDATAREAD(IB)%NBER
+ ELEM=GG_MAC%BCDATAREAD(IB)%ELEMNB(I)
+ IF(ELEM.GT.GG_MAC%NB_ELEM) CALL XABORT('MUSACG: inconsistent perimeter(1)')
+ ISURF=GG_MAC%ISURF2_ELEM(ELEM)
+ IF(ISURF.GT.NSURF_MACRO) CALL XABORT('MUSACG: inconsistent perimeter(2)')
+ PERIM_SURF(ISURF)=IB
+ ENDDO
+ ENDDO
+ CALL LCMPUT(ITRACK,'MATCOD',NREG,1,MATALB(1))
+ CALL LCMPUT(ITRACK,'ALBEDO',GG_MAC%NALBG,2,GALBED)
+ CALL LCMPUT(ITRACK,'PERIM_SURF',NSURF_MACRO,1,PERIM_SURF)
+ DEALLOCATE(PERIM_SURF,GALBED,ICODE)
+ DEALLOCATE(KEYMRG,VOLSUR,MATALB)
+ !----
+ ! Track macro geometry IMACR
+ !----
+ PRTIND=IPRINT
+ F_GEO=FGEO
+ EPS1=1.E-5_PDB
+ IF(RCUTOF>0._PDB) THEN
+ EPS1=RCUTOF
+ IF(PRTIND>0) WRITE(*,*) "MUSACG: set eps1 to ",EPS1
+ ENDIF
+ IGTRK=1
+ CALL SALTCG(ITRACK, IFTRK, IPRINT, IGTRK, NBSLIN, GG_MAC)
+ !----
+ ! Release allocated memory for macro IMACR
+ !----
+ CALL SALEND(GG_MAC)
+ DEALLOCATE(GG_MAC, STAT= OK)
+ IF(OK /= 0) CALL XABORT('MUSACG: failure to deallocate GG_MAC')
+END SUBROUTINE MUSACG
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:23:26 +0000