summaryrefslogtreecommitdiff
path: root/Dragon/src/PSOMON.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/PSOMON.f')
-rw-r--r--Dragon/src/PSOMON.f523
1 files changed, 523 insertions, 0 deletions
diff --git a/Dragon/src/PSOMON.f b/Dragon/src/PSOMON.f
new file mode 100644
index 0000000..8f4ea9e
--- /dev/null
+++ b/Dragon/src/PSOMON.f
@@ -0,0 +1,523 @@
+*DECK PSOMON
+ SUBROUTINE PSOMON(IPTRK,IPGEOM,NREG,LX,LY,LZ,NG,NDIM,
+ 1 NSOUR,ISOUR,ISISOUR,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XXX,YYY,ZZZ,
+ 2 MESHL,BSINFO,BS,MAXL,NORM,DIR,MONOP,NBST,NBS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute moments of fixed boundary monodirectionnal sources.
+*
+*Copyright:
+* Copyright (C) 2022 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): C. Bienvenue
+*
+*Parameters: input
+* IPTRK pointer to the tracking LCM object.
+* IPGEOM porinter to the geometry LCM object.
+* NREG number of regions.
+* LX number of meshes along X axis.
+* LY number of meshes along Y axis.
+* LZ number of meshes along Z axis.
+* NG number of energy groups.
+* NDIM geometry dimension.
+* NSOUR number of sources defined.
+* ISOUR intensity of the sources.
+* ISISOUR array with 0/1 values to indicate if energy group contain
+* sources
+* XMIN lower boundaries of the source along X axis.
+* XMAX upper boundaries of the source along X axis.
+* YMIN lower boundaries of the source along Y axis.
+* YMAX upper boundaries of the source along Y axis.
+* ZMIN lower boundaries of the source along Z axis.
+* ZMAX upper boundaries of the source along Z axis.
+* XXX regions boundaries along X axis.
+* YYY regions boundaries along Y axis.
+* ZZZ regions boundaries along Z axis.
+* MESHL number of regions along X-, Y- and Z-axis
+* MAXL number of intensity values needed to fully described each
+* boundary source.
+* DIR direction of the source particles.
+* MONOP value describing the boundary source location.
+* NBST total number of sources.
+*
+*Parameters: output
+* NORM normalization factor.
+* BS intensity values for each sources.
+* BSINFO energy group, discrete ordinate and location corresponding to
+* each sources.
+* NBS number of sources in each energy group.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK,IPGEOM
+ INTEGER NREG,NBST,LX,LY,LZ,NG,NDIM,NSOUR,MESHL(3),MONOP,
+ 1 ISISOUR(NG),BSINFO(2,NG,NBST),NBS(NG)
+ REAL XMIN(NSOUR),XMAX(NSOUR),YMIN(NSOUR),YMAX(NSOUR),ZMIN(NSOUR),
+ 1 ZMAX(NSOUR),ISOUR(NG),XXX(MESHL(1)),
+ 2 YYY(MESHL(2)),ZZZ(MESHL(3)),NORM,DIR(3),BS(MAXL,NG,NBST)
+*----
+* LOCAL VARIABLES
+*----
+ INTEGER SPLIT_LEN(3),XP(NREG),YP(NREG),ZP(NREG),
+ 1 IPN,MONOP2
+ REAL X(LX),Y(LY),Z(LZ),IPVAL,S
+*----
+* ALLOCATABLE ARRAYS
+*----
+ TYPE(C_PTR) U_PTR,DU_PTR,DE_PTR,DZ_PTR,W_PTR
+ INTEGER,ALLOCATABLE,DIMENSION(:) :: SPLITX,SPLITY,SPLITZ,MN
+ REAL, POINTER, DIMENSION(:) :: U,DU,DE,DZ,W
+ REAL,ALLOCATABLE,DIMENSION(:) :: P
+*----
+* RECOVER GEOMETRY INFORMATION
+*----
+ IF(ABS(MONOP).NE.1) THEN
+ CALL LCMLEN(IPGEOM,'SPLITX',SPLIT_LEN(1),ITYLCM)
+ ALLOCATE(SPLITX(SPLIT_LEN(1)))
+ CALL LCMGET(IPGEOM,'SPLITX',SPLITX)
+ ENDIF
+ IF(ABS(MONOP).NE.2.AND.NDIM.GE.2) THEN
+ CALL LCMLEN(IPGEOM,'SPLITY',SPLIT_LEN(2),ITYLCM)
+ ALLOCATE(SPLITY(SPLIT_LEN(2)))
+ CALL LCMGET(IPGEOM,'SPLITY',SPLITY)
+ ENDIF
+ IF(ABS(MONOP).NE.3.AND.NDIM.GE.3) THEN
+ CALL LCMLEN(IPGEOM,'SPLITZ',SPLIT_LEN(3),ITYLCM)
+ ALLOCATE(SPLITZ(SPLIT_LEN(3)))
+ CALL LCMGET(IPGEOM,'SPLITZ',SPLITZ)
+ ENDIF
+
+ ! CALCULATE X- OR Y- COORDINATES OF EACH VOXELS
+ IF(ABS(MONOP).NE.1) THEN
+ K=1
+ DO I=1,SPLIT_LEN(1)
+ DO J=1,SPLITX(I)
+ XP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+ DO 10 IX=1,LX
+ STEPX=(XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))
+ IF(XP(IX).EQ.1) THEN
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-1)
+ ELSE
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-SUM(SPLITX(1:XP(IX)-1))-1)
+ ENDIF
+ 10 CONTINUE
+ ENDIF
+ IF(ABS(MONOP).NE.2.AND.NDIM.GE.2) THEN
+ K=1
+ DO I=1,SPLIT_LEN(2)
+ DO J=1,SPLITY(I)
+ YP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+ DO 110 IY=1,LY
+ STEPY=(YYY(YP(IY)+1)-YYY(YP(IY)))/SPLITY(YP(IY))
+ IF(YP(IY).EQ.1) THEN
+ Y(IY)=YYY(YP(IY))+0.5*STEPY+STEPY*(IY-1)
+ ELSE
+ Y(IY)=YYY(YP(IY))+0.5*STEPY+STEPY*(IY-SUM(SPLITY(1:YP(IY)-1))-1)
+ ENDIF
+ 110 CONTINUE
+ ENDIF
+ IF(ABS(MONOP).NE.3.AND.NDIM.GE.3) THEN
+ K=1
+ DO I=1,SPLIT_LEN(3)
+ DO J=1,SPLITZ(I)
+ ZP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+ DO 220 IZ=1,LZ
+ STEPZ=(ZZZ(ZP(IZ)+1)-ZZZ(ZP(IZ)))/SPLITZ(ZP(IZ))
+ IF(ZP(IZ).EQ.1) THEN
+ Z(IZ)=ZZZ(ZP(IZ))+0.5*STEPZ+STEPZ*(IZ-1)
+ ELSE
+ Z(IZ)=ZZZ(ZP(IZ))+0.5*STEPZ+STEPZ*(IZ-SUM(SPLITZ(1:ZP(IZ)-1))-1)
+ ENDIF
+ 220 CONTINUE
+ ENDIF
+
+ IF(MONOP.EQ.-1) THEN
+ MONOP2=1
+ ELSE IF(MONOP.EQ.1) THEN
+ MONOP2=2
+ ELSE IF(MONOP.EQ.-2) THEN
+ MONOP2=3
+ ELSE IF(MONOP.EQ.2) THEN
+ MONOP2=4
+ ELSE IF(MONOP.EQ.-3) THEN
+ MONOP2=5
+ ELSE IF(MONOP.EQ.3) THEN
+ MONOP2=6
+ ELSE
+ MONOP2=0
+ ENDIF
+
+*----
+* 1D CARTESIAN CASE
+*----
+
+ IF(NDIM.EQ.1) THEN
+
+ ! Define the corresponding discrete ordinate
+
+ CALL LCMLEN(IPTRK,'U',NLF,ITYLCM)
+ CALL LCMGPD(IPTRK,'U',U_PTR)
+ CALL LCMGPD(IPTRK,'W',W_PTR)
+ CALL C_F_POINTER(U_PTR,U,(/ NLF /))
+ CALL C_F_POINTER(W_PTR,W,(/ NLF /))
+
+ IPN=1
+ IPVAL=(U(1)-DIR(1))**2
+ DO IP=2,NLF
+ IF((U(IP)-DIR(1))**2.LT.IPVAL) THEN
+ IPN=IP
+ IPVAL=(U(IP)-DIR(1))**2
+ ENDIF
+ ENDDO
+
+ IF(U(IPN).LT.0.0.AND.MONOP.EQ.-1) THEN
+ CALL XABORT('PSOMON: X- AND BACKWARD ORIENTED SOURCE.')
+ ELSEIF(U(IPN).GT.0.0.AND.MONOP.EQ.1) THEN
+ CALL XABORT('PSOMON: X+ AND FOWARD ORIENTED SOURCE.')
+ ENDIF
+
+ ! Normalization
+ NORM=SUM(ISOUR)
+
+ ! Save source information
+ DO IG=1,NG
+ IND=1
+ IF(ISISOUR(IG).EQ.1) THEN
+ BSINFO(1,IG,IND)=MONOP2
+ BSINFO(2,IG,IND)=IPN
+ BS(1,IG,IND)=ISOUR(IG)/W(IPN)
+ IND=IND+1
+ ENDIF
+ NBS(IG)=IND-1
+ ENDDO
+
+*----
+* 2D CARTESIAN CASE
+*----
+
+ ELSE IF(NDIM.EQ.2) THEN
+
+ ! Define the corresponding discrete ordinate
+
+ CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM)
+ CALL LCMGPD(IPTRK,'DU',DU_PTR)
+ CALL LCMGPD(IPTRK,'DE',DE_PTR)
+ CALL LCMGPD(IPTRK,'W',W_PTR)
+ CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /))
+ CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /))
+ CALL C_F_POINTER(W_PTR,W,(/ NPQ /))
+
+ ALLOCATE(MN(2),P(2))
+ INIT1=0
+ VAL1=8.0
+ VAL2=8.0
+ DO 300 M=1,NPQ
+ IF(W(M).EQ.0.0) GO TO 300
+ IF(INIT1.EQ.0) THEN
+ MN(1)=M
+ VAL1=(DU(M)-DIR(1))**2+(DE(M)-DIR(2))**2
+ INIT1=1
+ ELSE
+ VAL=(DU(M)-DIR(1))**2+(DE(M)-DIR(2))**2
+ IF(VAL1.GT.VAL) THEN
+ MN(2)=MN(1)
+ MN(1)=M
+ VAL2=VAL1
+ VAL1=VAL
+ ELSEIF(VAL2.GT.VAL) THEN
+ MN(2)=M
+ VAL2=VAL
+ ENDIF
+ ENDIF
+ 300 CONTINUE
+
+ VALTOT=VAL1+VAL2
+ P(1)=VAL1/VALTOT
+ P(2)=VAL2/VALTOT
+
+ IF(DU(MN(1)).LT.0.0.AND.DU(MN(2)).LT.0.0.AND.MONOP.EQ.-1) THEN
+ CALL XABORT('PSOMON: X- AND BACKWARD ORIENTED SOURCE.')
+ ELSEIF(DU(MN(1)).GT.0.0.AND.DU(MN(2)).GT.0.0.AND.MONOP.EQ.1) THEN
+ CALL XABORT('PSOMON: X+ AND FOWARD ORIENTED SOURCE.')
+ ELSEIF(DE(MN(1)).LT.0.0.AND.DE(MN(2)).LT.0.0.AND.MONOP.EQ.-2) THEN
+ CALL XABORT('PSOMON: Y- AND BACKWARD ORIENTED SOURCE.')
+ ELSEIF(DE(MN(1)).GT.0.0.AND.DE(MN(2)).GT.0.0.AND.MONOP.EQ.2) THEN
+ CALL XABORT('PSOMON: Y+ AND FOWARD ORIENTED SOURCE.')
+ ENDIF
+
+ ! Normalization
+ S=0.0
+ IF(ABS(MONOP).EQ.1) THEN
+ DO NS=1,NSOUR
+ S=S+(YMAX(NS)-YMIN(NS))
+ ENDDO
+ ENDIF
+ IF(ABS(MONOP).EQ.2) THEN
+ DO NS=1,NSOUR
+ S=S+(XMAX(NS)-XMIN(NS))
+ ENDDO
+ ENDIF
+ NORM=SUM(ISOUR)*S
+
+ ! Save source information
+ DO IG=1,NG
+ IND=1
+ DO IDIR=1,2
+ IF(ISISOUR(IG).EQ.1) THEN
+ BSINFO(1,IG,IND)=MONOP2
+ BSINFO(2,IG,IND)=MN(IDIR)
+ IF(ABS(MONOP).EQ.1) THEN
+ DO IY=1,LY
+ DO NS=1,NSOUR
+ IF(YMIN(NS).LE.Y(IY).AND.YMAX(NS).GE.Y(IY)) THEN
+ BS(IY,IG,IND)=ISOUR(IG)/W(MN(IDIR))*P(IDIR)
+ ELSE
+ BS(IY,IG,IND)=0.0
+ ENDIF
+ ENDDO
+ ENDDO
+ ELSEIF(ABS(MONOP).EQ.2) THEN
+ DO IX=1,LX
+ DO NS=1,NSOUR
+ IF(XMIN(NS).LE.X(IX).AND.XMAX(NS).GE.X(IX)) THEN
+ BS(IX,IG,IND)=ISOUR(IG)/W(MN(IDIR))*P(IDIR)
+ ELSE
+ BS(IX,IG,IND)=0.0
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDIF
+ IND=IND+1
+ ENDIF
+ ENDDO
+ NBS(IG)=IND-1
+ ENDDO
+
+ DEALLOCATE(MN,P)
+ IF(ALLOCATED(SPLITX)) DEALLOCATE(SPLITX)
+ IF(ALLOCATED(SPLITY)) DEALLOCATE(SPLITY)
+
+*----
+* 3D CARTESIAN CASE
+*----
+
+ ELSE IF(NDIM.EQ.3) THEN
+
+ ! Define the corresponding discrete ordinate
+
+ CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM)
+ CALL LCMGPD(IPTRK,'DU',DU_PTR)
+ CALL LCMGPD(IPTRK,'DE',DE_PTR)
+ CALL LCMGPD(IPTRK,'DZ',DZ_PTR)
+ CALL LCMGPD(IPTRK,'W',W_PTR)
+ CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /))
+ CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /))
+ CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /))
+ CALL C_F_POINTER(W_PTR,W,(/ NPQ /))
+
+ ALLOCATE(MN(4),P(4))
+ INIT1=0
+ INIT2=0
+ INIT3=0
+ VAL1=12.0
+ VAL2=12.0
+ VAL3=12.0
+ VAL4=12.0
+ DO 400 M=1,NPQ
+ IF(W(M).EQ.0.0) GO TO 400
+ IF(INIT1.EQ.0) THEN
+ MN(1)=M
+ VAL1=(DU(M)-DIR(1))**2+(DE(M)-DIR(2))**2+(DZ(M)-DIR(3))**2
+ INIT1=1
+ ELSEIF(INIT2.EQ.0) THEN
+ VAL=(DU(M)-DIR(1))**2+(DE(M)-DIR(2))**2+(DZ(M)-DIR(3))**2
+ IF(VAL1.GT.VAL) THEN
+ MN(2)=MN(1)
+ MN(1)=M
+ VAL2=VAL1
+ VAL1=VAL
+ ELSE
+ MN(2)=M
+ VAL2=VAL
+ ENDIF
+ INIT2=1
+ ELSEIF(INIT3.EQ.0) THEN
+ VAL=(DU(M)-DIR(1))**2+(DE(M)-DIR(2))**2+(DZ(M)-DIR(3))**2
+ IF(VAL1.GT.VAL) THEN
+ MN(3)=MN(2)
+ MN(2)=MN(1)
+ MN(1)=M
+ VAL3=VAL2
+ VAL2=VAL1
+ VAL1=VAL
+ ELSEIF(VAL2.GT.VAL) THEN
+ MN(3)=MN(2)
+ MN(2)=M
+ VAL3=VAL2
+ VAL2=VAL
+ ELSE
+ MN(3)=M
+ VAL3=VAL
+ ENDIF
+ INIT3=1
+ ELSE
+ VAL=(DU(M)-DIR(1))**2+(DE(M)-DIR(2))**2+(DZ(M)-DIR(3))**2
+ IF(VAL1.GT.VAL) THEN
+ MN(4)=MN(3)
+ MN(3)=MN(2)
+ MN(2)=MN(1)
+ MN(1)=M
+ VAL4=VAL3
+ VAL3=VAL2
+ VAL2=VAL1
+ VAL1=VAL
+ ELSEIF(VAL2.GT.VAL) THEN
+ MN(4)=MN(3)
+ MN(3)=MN(2)
+ MN(2)=M
+ VAL4=VAL3
+ VAL3=VAL2
+ VAL2=VAL
+ ELSEIF(VAL3.GT.VAL) THEN
+ MN(4)=MN(3)
+ MN(3)=M
+ VAL4=VAL3
+ VAL3=VAL
+ ELSEIF(VAL4.GT.VAL) THEN
+ MN(4)=M
+ VAL4=VAL
+ ENDIF
+ ENDIF
+ 400 CONTINUE
+
+ VALTOT=VAL1+VAL2+VAL3+VAL4
+ P(1)=VAL1/VALTOT
+ P(2)=VAL2/VALTOT
+ P(3)=VAL3/VALTOT
+ P(4)=VAL4/VALTOT
+
+ IF(DU(MN(1)).LT.0.0.AND.DU(MN(2)).LT.0.0.AND.
+ 1 DU(MN(3)).LT.0.0.AND.DU(MN(4)).LT.0.0.AND.MONOP.EQ.-1) THEN
+ CALL XABORT('PSOMON: X- AND BACKWARD X-ORIENTED SOURCE.')
+ ELSEIF(DU(MN(1)).GT.0.0.AND.DU(MN(2)).GT.0.0.AND.
+ 1 DU(MN(3)).GT.0.0.AND.DU(MN(4)).GT.0.0.AND.MONOP.EQ.1) THEN
+ CALL XABORT('PSOMON: X+ AND FOWARD X-ORIENTED SOURCE.')
+ ELSEIF(DE(MN(1)).LT.0.0.AND.DE(MN(2)).LT.0.0.AND.
+ 1 DE(MN(3)).LT.0.0.AND.DE(MN(4)).LT.0.0.AND.MONOP.EQ.-2) THEN
+ CALL XABORT('PSOMON: Y- AND BACKWARD Y-ORIENTED SOURCE.')
+ ELSEIF(DE(MN(1)).GT.0.0.AND.DE(MN(2)).GT.0.0.AND.
+ 1 DE(MN(3)).GT.0.0.AND.DE(MN(4)).GT.0.0.AND.MONOP.EQ.2) THEN
+ CALL XABORT('PSOMON: Y+ AND FOWARD Y-ORIENTED SOURCE.')
+ ELSEIF(DZ(MN(1)).LT.0.0.AND.DZ(MN(2)).LT.0.0.AND.
+ 1 DZ(MN(3)).LT.0.0.AND.DZ(MN(4)).LT.0.0.AND.MONOP.EQ.-3) THEN
+ CALL XABORT('PSOMON: Z- AND BACKWARD Z-ORIENTED SOURCE.')
+ ELSEIF(DZ(MN(1)).GT.0.0.AND.DZ(MN(2)).GT.0.0.AND.
+ 1 DZ(MN(3)).GT.0.0.AND.DZ(MN(4)).GT.0.0.AND.MONOP.EQ.3) THEN
+ CALL XABORT('PSOMON: Z+ AND FOWARD Z-ORIENTED SOURCE.')
+ ENDIF
+
+ ! Normalization
+ S=0.0
+ IF(ABS(MONOP).EQ.1) THEN
+ DO NS=1,NSOUR
+ S=S+(YMAX(NS)-YMIN(NS))*(ZMAX(NS)-ZMIN(NS))
+ ENDDO
+ ENDIF
+ IF(ABS(MONOP).EQ.2) THEN
+ DO NS=1,NSOUR
+ S=S+(XMAX(NS)-XMIN(NS))*(ZMAX(NS)-ZMIN(NS))
+ ENDDO
+ ENDIF
+ IF(ABS(MONOP).EQ.3) THEN
+ DO NS=1,NSOUR
+ S=S+(XMAX(NS)-XMIN(NS))*(YMAX(NS)-YMIN(NS))
+ ENDDO
+ ENDIF
+ NORM=SUM(ISOUR)*S
+
+ ! Save source information
+ DO IG=1,NG
+ IND=1
+ DO IDIR=1,4
+ IF(ISISOUR(IG).EQ.1) THEN
+ BSINFO(1,IG,IND)=MONOP2
+ BSINFO(2,IG,IND)=MN(IDIR)
+ IF(ABS(MONOP).EQ.1) THEN
+ DO IY=1,LY
+ DO IZ=1,LZ
+ IPOS=(IY-1)*LZ+IZ
+ DO NS=1,NSOUR
+ IF(YMIN(NS).LE.Y(IY).AND.YMAX(NS).GE.Y(IY).AND.
+ 1 ZMIN(NS).LE.Z(IZ).AND.ZMAX(NS).GE.Z(IZ)) THEN
+ BS(IPOS,IG,IND)=ISOUR(IG)/W(MN(IDIR))*P(IDIR)
+ ELSE
+ BS(IPOS,IG,IND)=0.0
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ELSEIF(ABS(MONOP).EQ.2) THEN
+ DO IX=1,LX
+ DO IZ=1,LZ
+ IPOS=(IX-1)*LZ+IZ
+ DO NS=1,NSOUR
+ IF(XMIN(NS).LE.X(IX).AND.XMAX(NS).GE.X(IX).AND.
+ 1 ZMIN(NS).LE.Z(IZ).AND.ZMAX(NS).GE.Z(IZ)) THEN
+ BS(IPOS,IG,IND)=ISOUR(IG)/W(MN(IDIR))*P(IDIR)
+ ELSE
+ BS(IPOS,IG,IND)=0.0
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ELSEIF(ABS(MONOP).EQ.3) THEN
+ DO IX=1,LX
+ DO IY=1,LY
+ IPOS=(IX-1)*LY+IY
+ DO NS=1,NSOUR
+ IF(XMIN(NS).LE.X(IX).AND.XMAX(NS).GE.X(IX).AND.
+ 1 YMIN(NS).LE.Y(IY).AND.YMAX(NS).GE.Y(IY)) THEN
+ BS(IPOS,IG,IND)=ISOUR(IG)/W(MN(IDIR))*P(IDIR)
+ ELSE
+ BS(IPOS,IG,IND)=0.0
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ IND=IND+1
+ ENDIF
+ ENDDO
+ NBS(IG)=IND-1
+ ENDDO
+
+ DEALLOCATE(MN,P)
+ IF(ALLOCATED(SPLITX)) DEALLOCATE(SPLITX)
+ IF(ALLOCATED(SPLITY)) DEALLOCATE(SPLITY)
+ IF(ALLOCATED(SPLITZ)) DEALLOCATE(SPLITZ)
+
+ ELSE
+ CALL XABORT('SOUR: INVALID GEOMETRY, ONLY 1D, 2D AND 3D CARTESIAN'
+ 1 //' GEOMETRY ARE ACTUALLY IMPLEMENTED.')
+ ENDIF
+ RETURN
+ END