From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/SALTCG.f | 606 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 606 insertions(+) create mode 100644 Dragon/src/SALTCG.f (limited to 'Dragon/src/SALTCG.f') diff --git a/Dragon/src/SALTCG.f b/Dragon/src/SALTCG.f new file mode 100644 index 0000000..123a288 --- /dev/null +++ b/Dragon/src/SALTCG.f @@ -0,0 +1,606 @@ +*DECK SALTCG + SUBROUTINE SALTCG(IPTRK ,IFTRK ,IPRINT,IGTRK ,NBSLIN, GG) +* +*----------------------------------------------------------------------- +* +*Purpose: +* To track an assembly of cells containing clusters using the new SALT +* tracking procedure (based on NXTTCG.f). +* +*Copyright: +* Copyright (C) 2014 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 and G. Marleau +* +*Parameters: input +* IPTRK pointer to the TRACKING data structure in +* update or creation mode. +* IFTRK pointer to the TRACKING file in creation mode. +* IPRINT print level. +* IGTRK flag to generate the tracking file. In the case where +* IGTRK=1, the tracking is performed and +* used to evaluate the track normalisation factor and the +* tracking file is generated. When IGTRK=0, the tracking is +* still performed and used to evaluate the track normalisation +* factor but the tracking file is not generated. +* NBSLIN maximum number of segments in a single tracking line. +* computed by default in SALTCG but limited to 100000 +* elements. This default value can be bypassed using +* keyword NBSLIN. +* GG geometry basic information. +* +*---------- +* + USE GANLIB + USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC,ISPEC + USE SAL_TRACKING_TYPES, ONLY : NMAX2,ITRAC2,RTRAC2,IPART,RPART, + > NIPART,NRPART + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + TYPE(C_PTR) IPTRK + INTEGER IFTRK + INTEGER IPRINT,IGTRK + INTEGER NBSLIN + TYPE(T_G_BASIC) GG +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='SALTCG') + INTEGER NSTATE + PARAMETER (NSTATE=40) + INTEGER NMAX0 + DOUBLE PRECISION PI,DZERO,DONE,DTWO,DSUM + PARAMETER (PI=3.14159265358979, DZERO=0.0D0,DONE=1.0D0, + > DTWO=2.0D0,NMAX0=100000) +*---- +* Functions +*---- + INTEGER KDROPN,IFTEMP,KDRCLS,ICLS +*---- +* Local variables +*---- + INTEGER ISTATE(NSTATE),IEDIMG(NSTATE),ICODE(6) + REAL RSTATT(NSTATE),ALBEDO(6) + INTEGER RENO,LTRK,AZMOAQ,ISYMM,POLQUA,POLOAQ,AZMQUA, + > AZMNBA,OK + DOUBLE PRECISION DENUSR,RCUTOF,DENLIN,SPACLN,WEIGHT + DOUBLE PRECISION RADIUS,CENTER(3) + INTEGER NDIM,ITYPBC,IDIRG,NBOCEL,NBUCEL,IDIAG, + > ISAXIS(3),NOCELL(3),NUCELL(3),MXMSH,MAXMSH, + > MAXREG,NBTCLS,MAXMSP,MAXRSP,NFSUR,NFREG, + > MXGSUR,NUNK,NPLANE,NPOINT,NTLINE,NBTDIR, + > MAXSUB,MAXSGL,NBDR,ILONG,ITYLCM,IPER(3) + INTEGER JJ,KK,NCOR,NQUAD,NANGL,NBANGL,LINMAX + DOUBLE PRECISION DQUAD(4),ABSC(3,2),RCIRC,SIDEH,ANGLE + CHARACTER CTRK*4,COMENT*80 + INTEGER IFMT,NEREG,NESUR +*---- +* Allocatable arrays +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYMRG,MATALB + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NBSANG + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: SURVOL,DGMESH, + > DNSANG,DDANG,DVNOR,DSNOR + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DDENWT + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: DANGLT +*---- +* Processing starts: +* print routine opening header if required +* and initialize various parameters. +*---- + IF(IPRINT .GE. 1) WRITE(IOUT,6000) NAMSBR +*---- +* Open temporary tracking file if required +*---- + IF(IGTRK .EQ. 1) THEN + IFTEMP= KDROPN('DUMMYSQ',0,2,0) + IF(IFTEMP .LE. 0) WRITE(IOUT,9010) NAMSBR + ENDIF +*---- +* Get state vectors +*---- + ISTATE(:NSTATE)=0 + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + CALL LCMGET(IPTRK,'EXCELTRACKOP',RSTATT) + NEREG=ISTATE(1) + NESUR=ISTATE(5) + NUNK=NEREG+NESUR+1 + RENO =ISTATE(8) + LTRK =ISTATE(9) + AZMOAQ=ISTATE(11) + ISYMM =ISTATE(12) + POLQUA=ISTATE(13) + POLOAQ=ISTATE(14) + AZMQUA=ISTATE(15) + AZMNBA=ISTATE(16) + IFMT=ISTATE(21) + DENUSR=DBLE(RSTATT(2)) + RCUTOF=DBLE(RSTATT(3)) + DENLIN=DBLE(RSTATT(4)) + SPACLN=DBLE(RSTATT(5)) + WEIGHT=RSTATT(6) +*---- +* Get main tracking records +*---- + CALL LCMLEN(IPTRK,'ICODE ',ILONG,ITYLCM) + IF(ILONG.GT.6) CALL XABORT('SALTCG: ALBEDO OVERFLOW.') + CALL LCMGET(IPTRK,'ICODE ',ICODE ) + CALL LCMGET(IPTRK,'ALBEDO ',ALBEDO) + CALL LCMSIX(IPTRK,'NXTRecords ',1) +*---- +* Get general dimensioning vector for geometry tracking +*---- + IEDIMG(:NSTATE)=0 + CALL LCMGET(IPTRK,'G00000001DIM',IEDIMG) + NDIM =IEDIMG( 1) + ITYPBC =IEDIMG( 2) + IDIRG =IEDIMG( 3) + NBOCEL =IEDIMG( 4) + NBUCEL =IEDIMG( 5) + IDIAG =IEDIMG( 6) + ISAXIS(1)=IEDIMG( 7) + ISAXIS(2)=IEDIMG( 8) + ISAXIS(3)=IEDIMG( 9) + NOCELL(1)=IEDIMG(10) + NOCELL(2)=IEDIMG(11) + NOCELL(3)=IEDIMG(12) + NUCELL(1)=IEDIMG(13) + NUCELL(2)=IEDIMG(14) + NUCELL(3)=IEDIMG(15) + MXMSH =IEDIMG(16) + MAXREG =IEDIMG(17) + NBTCLS =IEDIMG(18) + MAXMSP =IEDIMG(20) + MAXRSP =IEDIMG(21) + NFSUR =IEDIMG(22) + NFREG =IEDIMG(23) + MXGSUR =IEDIMG(24) + MAXMSH=MAX(1,MXMSH,MAXMSP,MAXREG) + IF(IPRINT .GE. 1) THEN + WRITE(IOUT,6011) NFREG,NEREG,NFSUR,NESUR + ENDIF + IF((ITYPBC .EQ. 0).OR.(ITYPBC .EQ. 2)) THEN +*---- +* Define Cell for periodicity +* Cartesian Boundary +*---- + IPER(1)=2 + IF(ABS(ISAXIS(1)) .EQ. 3) IPER(1)=1 + IPER(2)=2 + IF(ABS(ISAXIS(2)) .EQ. 3) IPER(2)=1 + IPER(3)=2 + IF(ABS(ISAXIS(3)) .EQ. 3) IPER(3)=1 +*---- +* Use intrinsic geometry symmetries +* to simplify tracking unless +* NOSY tracking option activated +*---- + IF(ISYMM .NE. 0) THEN + ISYMM=0 + IF(ABS(ISAXIS(1)) .EQ. 1 .OR. ABS(ISAXIS(1)) .EQ. 2) THEN +*---- +* X SYMMETRY +*---- + ISYMM=2 + ENDIF + IF(ABS(ISAXIS(2)) .EQ. 1 .OR. ABS(ISAXIS(2)) .EQ. 2) THEN +*---- +* Y SYMMETRY +*---- + ISYMM=4+2*ISYMM + ENDIF + IF(NDIM .EQ. 3) THEN + IF(ABS(ISAXIS(3)) .EQ. 1 .OR. ABS(ISAXIS(3)) .EQ. 2) THEN +*---- +* Z SYMMETRY +*---- + ISYMM=16+ISYMM + ENDIF + ENDIF + IF(ISYMM .EQ. 0) ISYMM=1 + ENDIF + ENDIF +*---- +* Read global mesh for geometry +*---- + ALLOCATE(DGMESH((MAXMSH+2)*4)) + IF(ITYPBC.EQ.0) THEN + CALL NXTXYZ(IPTRK ,IPRINT,NDIM ,ITYPBC,MAXMSH,NUCELL, + > ABSC,DGMESH) + ELSE IF(ITYPBC.GE.2) THEN + ! hexagonal geometry + CALL LCMGET(IPTRK,'G00000001SMX',DGMESH) + SIDEH=DGMESH(2)-DGMESH(1) + ABSC(:2,1)=SIDEH*SQRT(3.0) + ABSC(:2,2)=ABSC(:2,1) + ELSE + CALL XABORT(NAMSBR//': geometry not implemented') + ENDIF +*---- +* Verify tracking parameters and compute number of angles +* associated with angular order and spatial quadrature parameters +* 1. Isotropic tracking +*---- + NCOR= 1 + NPLANE=1 + IF(LTRK .EQ. 0) THEN + IF(NDIM .EQ. 3) THEN + CALL XABORT(NAMSBR//': 3-D geometry is not allowed') + ENDIF + NQUAD=2 + DQUAD(1)=DONE + DQUAD(2)=DONE + NANGL=AZMOAQ + NBANGL=NANGL + IF(ISYMM .EQ. 2 .OR. ISYMM .EQ. 8) THEN + DQUAD(1)=DONE/DTWO + DQUAD(2)=DZERO + ENDIF + ALLOCATE(DANGLT(NDIM,NBANGL,NQUAD),DDENWT(NBANGL,NQUAD)) + CALL NXTQAS(IPRINT,NDIM ,AZMQUA,NANGL ,NQUAD ,NBANGL, + > DQUAD ,DANGLT,DDENWT) + LINMAX=NBUCEL*(4+MXGSUR+16) +*---- +* Select standard spatial tracking parameters +*---- + CALL NXTQSS(IPRINT,NDIM ,ITYPBC,MAXMSH,NUCELL,DENUSR, + > DGMESH,NPLANE,NPOINT,DENLIN,SPACLN,WEIGHT, + > RADIUS,CENTER) + CALL LCMPUT(IPTRK,'TrackingDirc',NDIM*NBANGL*NQUAD,4,DANGLT) + ELSE +*---- +* 2. Specular tracking +*---- + NPOINT=0 + NQUAD=2 + IF(NDIM .EQ. 3) CALL XABORT(NAMSBR// + > ': TSPC option not valid for 3-D geometries') + IF(ITYPBC .EQ. 0) THEN + ! Cartesian geometries + IF(AZMOAQ .GT. 24) THEN + WRITE(IOUT,9002) NAMSBR,AZMOAQ,24,30 + AZMOAQ=30 + ELSE IF(AZMOAQ .GT. 20) THEN + IF(AZMOAQ .NE. 24) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,20,24,24 + AZMOAQ=24 + ENDIF + ELSE IF(AZMOAQ .GT. 18) THEN + IF(AZMOAQ .NE. 20) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,18,20,20 + AZMOAQ=20 + ENDIF + ELSE IF(AZMOAQ .GT. 14) THEN + IF(AZMOAQ .NE. 18) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,14,18,18 + AZMOAQ=18 + ENDIF + ELSE IF(AZMOAQ .GT. 12) THEN + IF(AZMOAQ .NE. 14) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,12,14,14 + AZMOAQ=14 + ENDIF + ELSE IF(AZMOAQ .GT. 8) THEN + IF(AZMOAQ .NE. 12) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,8,12,12 + AZMOAQ=12 + ENDIF + ELSE IF(AZMOAQ .GT. 6) THEN + IF(AZMOAQ .NE. 8) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,6,8,8 + AZMOAQ=8 + ENDIF + ELSE IF(AZMOAQ .GT. 2) THEN + IF(AZMOAQ .NE. 6) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,2,6,6 + AZMOAQ=6 + ENDIF + ELSE IF(AZMOAQ .GE. 0) THEN + IF(AZMOAQ .NE. 2) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,-1,2,2 + AZMOAQ=2 + ENDIF + ENDIF + ELSE IF(ITYPBC .GE.2) THEN + ! hexagonal geometries + IF(AZMOAQ .GT. 12) THEN + WRITE(IOUT,9002) NAMSBR,AZMOAQ,12,18 + AZMOAQ=18 + ELSE IF(AZMOAQ .GT. 6) THEN + IF(AZMOAQ .NE. 12) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,6,12,12 + AZMOAQ=12 + ENDIF + ELSE IF(AZMOAQ .GT. 3) THEN + IF(AZMOAQ .NE. 6) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,3,6,6 + AZMOAQ=6 + ENDIF + ELSE IF(AZMOAQ .GT. 1) THEN + IF(AZMOAQ .NE. 3) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,1,3,3 + AZMOAQ=3 + ENDIF + ELSE IF(AZMOAQ .GE. 0) THEN + IF(AZMOAQ .NE. 1) THEN + WRITE(IOUT,9003) NAMSBR,AZMOAQ,-1,1,1 + AZMOAQ=1 + ENDIF + ENDIF + ENDIF + NBANGL=AZMOAQ + NANGL =AZMOAQ + ALLOCATE(NBSANG(5,NBANGL)) + ALLOCATE(DANGLT(NDIM,NBANGL,4),DDENWT(NBANGL,4),DNSANG(NBANGL), + > DDANG(NBANGL)) + LINMAX=8*NANGL*NBUCEL*(4+MXGSUR+16) + RCIRC=SQRT(ABSC(1,1)**2+ABSC(2,1)**2) + ABSC(1,1)= ABSC(1,1)/RCIRC + ABSC(2,1)= ABSC(2,1)/RCIRC + CALL NXTQAC(IPRINT,NDIM ,NANGL ,NBANGL,ITYPBC,DENUSR, + > ABSC ,RCIRC ,AZMQUA,IPER ,DANGLT,DDENWT, + > DNSANG,NBSANG,DDANG) + DEALLOCATE(DDANG) + DO JJ=1,NBANGL + DANGLT(1,NBANGL-JJ+1,2)=-DANGLT(1,JJ,1) + DANGLT(2,NBANGL-JJ+1,2)=DANGLT(2,JJ,1) + DDENWT(NBANGL-JJ+1,2)=DDENWT(JJ,1) + ENDDO + DO JJ=1,NBANGL + DANGLT(1,NBANGL-JJ+1,4)=DANGLT(1,JJ,1) + DANGLT(2,NBANGL-JJ+1,4)=-DANGLT(2,JJ,1) + DDENWT(NBANGL-JJ+1,4)=DDENWT(JJ,1) + DANGLT(1,NBANGL-JJ+1,3)=DANGLT(1,JJ,2) + DANGLT(2,NBANGL-JJ+1,3)=-DANGLT(2,JJ,2) + DDENWT(NBANGL-JJ+1,3)=DDENWT(JJ,2) + ENDDO + IF(IPRINT.GT.1) THEN + WRITE(IOUT,'(/34H SALTCG: CYCLIC ANGULAR QUADRATURE/8X, + > 5HANGLE,8X,7HCOSINES,16(1h-),3X,6HWEIGHT,11X, + > 17HCYCLIC PARAMETERS)') + DSUM=0.D0 + DO KK=1,4 + DO JJ=1,NANGL + ANGLE=SIGN(ACOS(DANGLT(1,JJ,KK))/PI*180.0,DANGLT(2,JJ,KK)) + WRITE(IOUT,'(1X,I4,1P,4E13.4,5X,2I4)') (KK-1)*NANGL+JJ, + > ANGLE,DANGLT(:2,JJ,KK),0.5D0/DDENWT(JJ,KK),NBSANG(:2,JJ) + DSUM=DSUM+0.5D0/DDENWT(JJ,KK) + ENDDO + ENDDO + WRITE(IOUT,'(39X,5HDSUM=,1P,E13.4)') DSUM + ENDIF + CALL LCMPUT(IPTRK,'TrackingDirc',NDIM*NBANGL*4,4,DANGLT) + CALL LCMPUT(IPTRK,'TrackingTrkW',NBANGL*4,4,DDENWT) + CALL LCMPUT(IPTRK,'TrackingSpaD',NBANGL,4,DNSANG) + CALL LCMPUT(IPTRK,'TrackingNbST',5*NBANGL,1,NBSANG) + ENDIF + RSTATT(4)=REAL(DENLIN) + RSTATT(5)=REAL(SPACLN) + RSTATT(6)=REAL(WEIGHT) + RSTATT(7)=REAL(RADIUS) + RSTATT(8)=REAL(CENTER(1)) + RSTATT(9)=REAL(CENTER(2)) +*---- +* Allocate memory to hold tracking data +*---- +** tracking data buffer: +* integers +* ITRAC2(NMAX OR 2*NMAX) = integer tracking array +* +* *integer descriptors in itrac2: +* 1 = address of last data +* 2 = total number of sub-trajectories +* 3 = +* 4 = phi for trajectory (2D) +* +* reals +* RTRAC2(nmax or 2*nmax) = real tracking array +* +* *real descriptors: +* 1 = +* 2 = cos phi entering basic +* 3 = sin phi left surface +* 4 = sin phi right surface +* 5 = cos phi left surface +* 6 = cos phi right surface +* 7 = total weight (DELR*WPHI) +* 8 = radial weight (DELR) +* +* IPART(NIPART,MXELEM) = to store integer intersection data +* RPART(NRPART,MXELEM) = to store real intersection data +* + IF(NBSLIN <= 0)THEN + NMAX2=NMAX0 + ELSE + NMAX2=NBSLIN + ENDIF + IF(ISPEC == 1) NMAX2=NMAX2*100 + ALLOCATE(ITRAC2(2*NMAX2),IPART(NIPART,GG%NB_ELEM),RTRAC2(NMAX2), + 1 RPART(NRPART,GG%NB_ELEM),STAT=OK) + IF(OK/=0) CALL XABORT('SALTCG: not enough memory IRD') +*---- +* Track +*---- + LINMAX=MAX(LINMAX,NBSLIN) + IF(IPRINT .GE. 10) WRITE(IOUT,6010) LINMAX + NBDR=1 + IF(RENO .EQ. -1) THEN + IF(LTRK .EQ. 0) THEN + NBDR=NQUAD*NBANGL+1 + ELSE IF(LTRK .EQ. 1) THEN + NBDR=4*NBANGL+1 + ENDIF + ENDIF + ALLOCATE(DVNOR(NFREG*NBDR),DSNOR(NFSUR*NQUAD*NBANGL)) + IF(LTRK .EQ. 0) THEN +*---- +* Standard (isotropic) tracking (white boundary conditions) +*---- + IF(IPRINT .GE. 1) WRITE(IOUT,6030) NBANGL*NQUAD,NPOINT + MAXSUB=1 + CALL SALTLS(IFTEMP,IPRINT,IGTRK ,NFREG ,NBANGL, NQUAD , + > RENO ,NBDR ,IFMT ,DENUSR,DANGLT, DDENWT, + > GG ,NBTDIR,MAXSGL,NTLINE,DVNOR ) + ELSE +*---- +* Cyclic (specular) tracking (mirror like boundary conditions) +*---- + NBTDIR=0 + CALL SALTLC(IFTEMP,IPRINT,IGTRK,NDIM,NFREG,NBANGL,RENO,NBDR, + > IFMT,DENUSR,DANGLT,DDENWT,NBSANG,GG,MAXSUB,MAXSGL, + > NTLINE,DVNOR ) + ENDIF +*---- +* Release allocated memory for SALT tracking +*---- + DEALLOCATE(RPART,RTRAC2,IPART,ITRAC2,STAT =OK) + IF(OK /= 0) CALL XABORT('SALTCG: failure to deallocate storage') +*---- +* Save track normalisation vector +*---- + CALL LCMPUT(IPTRK,'NumMerge ',NFREG,1,GG%NUM_MERGE) + CALL LCMPUT(IPTRK,'VolMerge ',NFREG,4,GG%VOL_NODE) + CALL LCMPUT(IPTRK,'VTNormalize ',NFREG,4,DVNOR) + IF(NBDR .GT. 1) THEN + CALL LCMPUT(IPTRK,'VTNormalizeD',NFREG*(NBDR-1),4, + > DVNOR(NFREG+1)) + ENDIF +*---- +* Get cell description of geometry +*---- + ALLOCATE(KEYMRG(-NFSUR:NFREG),MATALB(-NFSUR:NFREG), + > SURVOL(-NFSUR:NFREG)) + CALL LCMLEN(IPTRK,'KEYMRG ',ILONG,ITYLCM) + IF(ILONG>NFREG+NFSUR+1) CALL XABORT('SALTCG: KEYMRG OVERLOW.') + CALL LCMGET(IPTRK,'KEYMRG ',KEYMRG) + CALL LCMGET(IPTRK,'MATALB ',MATALB) + CALL LCMGET(IPTRK,'SAreaRvolume',SURVOL) +*---- +* Build NXTRecords directory +*---- + CALL LCMSIX(IPTRK,'NXTRecords ',2) + ISTATE(12)=ISYMM + ISTATE(14)=POLOAQ + ISTATE(17)=NPOINT + ISTATE(18)=LINMAX + ISTATE(19)=NTLINE + ISTATE(20)=NBTDIR + IF(LTRK .EQ. 0) THEN + ISTATE(21)=NQUAD*NBANGL + ELSE IF(LTRK .EQ. 1) THEN + ISTATE(21)=4*NBANGL + ENDIF + ISTATE(22)=NPLANE + CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,ISTATE) + CALL LCMPUT(IPTRK,'EXCELTRACKOP',NSTATE,2,RSTATT) +*---- +* Renormalize tracks if required and transfer to final tracking file +*---- + IF(IGTRK .EQ. 1) THEN + CTRK = '$TRK' + WRITE(IFTRK) CTRK,5,NTLINE,IFMT + COMENT='CREATOR : DRAGON' + WRITE(IFTRK) COMENT + COMENT='MODULE : SALTCG' + WRITE(IFTRK) COMENT + COMENT='TYPE : CARTESIAN' + WRITE(IFTRK) COMENT + IF(RENO .EQ. -1) THEN + COMENT='TRKNOR : Directional ' + ELSE IF(RENO .EQ. 0) THEN + COMENT='TRKNOR : Global ' + ELSE + COMENT='TRKNOR : Off ' + ENDIF + WRITE(IFTRK) COMENT + IF(IFMT .EQ. 1) THEN + COMENT='OPTION : Extended ' + WRITE(IFTRK) COMENT + ELSE + COMENT='OPTION : Short ' + WRITE(IFTRK) COMENT + ENDIF +*---- +* Compress VOLSUR and MATALB according to KEYMRG and save on IFTRK +*---- + IF(LTRK .EQ. 0) THEN + WRITE(IFTRK) NDIM,LTRK,NEREG,NESUR,6,NCOR,NQUAD*NBANGL,MAXSUB, + > MAXSGL + ELSE IF(LTRK .EQ. 1) THEN + WRITE(IFTRK) NDIM,LTRK,NEREG,NESUR,6,NCOR,4*NBANGL,MAXSUB, + > MAXSGL + ENDIF + KEYMRG(1:NFREG)=GG%NUM_MERGE(:NFREG) + SURVOL(1:NFREG)=GG%VOL_NODE(:NFREG) + MATALB(1:NFREG)=GG%MED(:NFREG) + CALL NXTCVM(IFTRK,IPRINT,NFREG,NFSUR,NEREG,NESUR,MATALB,SURVOL, + > KEYMRG) + WRITE(IFTRK) ( ICODE(JJ),JJ=1,6) + WRITE(IFTRK) (ALBEDO(JJ),JJ=1,6) + IF(LTRK .EQ. 0) THEN + CALL NXTSQD(IFTRK,IPRINT,NDIM,NQUAD,NBANGL,DANGLT,DDENWT) + ELSE IF(LTRK .EQ. 1) THEN + WRITE(IFTRK) ((DANGLT(1,JJ,KK),DANGLT(2,JJ,KK),JJ=1,NBANGL), + > KK=1,4) + WRITE(IFTRK) ((DDENWT(JJ,KK),JJ=1,NBANGL),KK=1,4) + ENDIF + REWIND IFTEMP + CALL NXTTNS(IFTRK ,IFTEMP,IPRINT,RENO ,NFSUR ,NFREG ,NDIM , + > MAXSUB,MAXSGL,NTLINE,NBDR ,IFMT ,KEYMRG,DVNOR) +*---- +* Close temporary tracking file if required +*---- + ICLS=KDRCLS(IFTEMP,2) + IF(ICLS .NE. 0) WRITE(IOUT,9011) NAMSBR + ENDIF +*---- +* Deallocate memory +*---- + DEALLOCATE(DSNOR,DVNOR) + IF(LTRK .EQ. 1) DEALLOCATE(NBSANG,DNSANG) + DEALLOCATE(DDENWT,DANGLT) + DEALLOCATE(DGMESH,SURVOL,MATALB,KEYMRG) +*---- +* Processing finished: +* print routine closing output header if required +* and return +*---- + IF(IPRINT .GE. 1) THEN + WRITE(IOUT,6012) + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6010 FORMAT(' Maximum length of a line =',I10) + 6011 FORMAT(' Tracking of geometry begins:'/ + > ' Number of regions before merge =',I10/ + > ' Number of regions after merge =',I10/ + > ' Number of surfaces before merge=',I10/ + > ' Number of surfaces after merge =',I10) + 6012 FORMAT(' Tracking of geometry completed') + 6030 FORMAT(' Number of directions for tracking = ',I10/ + > ' Number of lines per direction = ',I10) + 9002 FORMAT(' ***** Warning in ',A6,' *****'/ + > ' Number of specular angles requested :',I10/ + > ' For values > ',I10,' use ',I10) + 9003 FORMAT(' ***** Warning in ',A6,' *****'/ + > ' Number of specular angles requested :',I10/ + > ' For values > ',I10,' and < ',I10,' use ',I10) + 9010 FORMAT(' ***** Warning in ',A6,' *****'/ + > ' Impossible to open temporary tracking file ') + 9011 FORMAT(' ***** Warning in ',A6,' *****'/ + > ' Impossible to close temporary tracking file ') + END -- cgit v1.2.3