Commit 2d4f2cf0 authored by Femi Kadri's avatar Femi Kadri

original files

parents
SUBROUTINE INIT_GEOMETRY_backwards_step(&
& global_x_min_real,&
& global_x_max_real,&
& global_x_min_ghost,&
& global_x_max_ghost,&
& global_y_min_real,&
& global_y_max_real,&
& global_y_min_ghost,&
& global_y_max_ghost,&
& global_z_min_real,&
& global_z_max_real,&
& My_MPI_Process_ID,&
& local_x_min_ghost,&
& local_x_max_ghost,&
& local_y_min_ghost,&
& local_y_max_ghost,&
& local_z_min_ghost,&
& local_z_max_ghost,&
& local_walls)
!THIS SUBROUTINE CREATES A BACKWARDS STEP (SINCE THE CODE ASSUMES PERIODIC BC, THE LENGTH OF THE STEP IS CONTROLLED)
USE debug_module
IMPLICIT NONE
!-----------INTENT(IN)-----------------------------------------------------------------------------------------------
INTEGER,INTENT(IN) :: global_x_min_real,global_x_max_real,global_x_min_ghost,global_x_max_ghost
INTEGER,INTENT(IN) :: global_y_min_real,global_y_max_real,global_y_min_ghost,global_y_max_ghost
INTEGER,INTENT(IN) :: global_z_min_real,global_z_max_real,My_MPI_Process_ID
INTEGER,INTENT(IN) :: local_x_min_ghost,local_x_max_ghost
INTEGER,INTENT(IN) :: local_y_min_ghost,local_y_max_ghost
INTEGER,INTENT(IN) :: local_z_min_ghost,local_z_max_ghost
LOGICAL*4,PARAMETER :: debug = .TRUE. !GLOBAL_DEBUG
!-----------INTENT(INOUT)-----------------------------------------------------------------------------------------------
LOGICAL*4,INTENT(INOUT),DIMENSION(local_x_min_ghost:local_x_max_ghost, &
& local_y_min_ghost:local_y_max_ghost,local_z_min_ghost:local_z_max_ghost) &
& :: local_walls
!-----------NO INTENT-----------------------------------------------------------------------------------------------
INTEGER :: global_STEP_x_max_real,global_STEP_y_max_real
INTEGER :: I,J,K
IF(debug) WRITE (9,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: entered subroutine'
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: entered subroutine'
global_STEP_x_max_real = INT( 0.1 * REAL(global_x_max_real) ) !USER INPUT (CONTROLS STEP WIDTH RELATIVE TO DOMAIN WIDTH)
global_STEP_y_max_real = INT( 0.5 * REAL(global_y_max_real) ) !USER INPUT (CONTROLS STEP HEIGHT RELATIVE TO DOMAIN HEIGHT)
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_STEP_x_max_real',global_STEP_x_max_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_STEP_y_max_real',global_STEP_y_max_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: local_z_min_ghost',local_z_min_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: local_z_max_ghost',local_z_max_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: local_y_min_ghost',local_y_min_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: local_y_max_ghost',local_y_max_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: local_x_min_ghost',local_x_min_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: local_x_max_ghost',local_x_max_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_z_min_real',global_z_min_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_z_max_real',global_z_max_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_y_min_real',global_y_min_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_y_max_real',global_y_max_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_x_min_real',global_x_min_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: global_x_max_real',global_x_max_real
DO K=local_z_min_ghost,local_z_max_ghost
DO J=local_y_min_ghost,local_y_max_ghost
DO I=local_x_min_ghost,local_x_max_ghost
IF( (J.LE.global_STEP_y_max_real).AND.(I.LE.global_STEP_x_max_real) ) THEN
local_walls(I,J,K) = .TRUE.
ENDIF !!(J.LE.global_STEP_y_max_real).AND.(I.LE.global_STEP_x_max_real)
IF( (J.EQ.global_y_min_real).OR.(J.EQ.global_y_max_real) ) THEN
local_walls(I,J,K) = .TRUE.
ENDIF !!J.EQ.global_y_min_real).AND.(J.EQ.global_y_max_real)
END DO !k = 1, nz
END DO !j = 1, ny
END DO !i = 1, nx
IF(debug) WRITE (9,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: about to return from subroutine'
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_backwards_step: about to return from subroutine'
RETURN
END SUBROUTINE INIT_GEOMETRY_backwards_step
\ No newline at end of file
SUBROUTINE INIT_GEOMETRY_comb(&
& global_x_min_real,&
& global_x_max_real,&
& global_x_min_ghost,&
& global_x_max_ghost,&
& global_y_min_real,&
& global_y_max_real,&
& global_y_min_ghost,&
& global_y_max_ghost,&
& global_z_min_real,&
& global_z_max_real,&
& My_MPI_Process_ID,&
& local_x_min_ghost,&
& local_x_max_ghost,&
& local_y_min_ghost,&
& local_y_max_ghost,&
& local_z_min_ghost,&
& local_z_max_ghost,&
& local_walls)
!THIS SUBROUTINE CREATES A SLIT GEOMETRY (TOP WALL AND BOTTOM WALL)
USE debug_module
IMPLICIT NONE
!-----------INTENT(IN)-----------------------------------------------------------------------------------------------
INTEGER,INTENT(IN) :: global_x_min_real,global_x_max_real,global_x_min_ghost,global_x_max_ghost
INTEGER,INTENT(IN) :: global_y_min_real,global_y_max_real,global_y_min_ghost,global_y_max_ghost
INTEGER,INTENT(IN) :: global_z_min_real,global_z_max_real,My_MPI_Process_ID
INTEGER,INTENT(IN) :: local_x_min_ghost,local_x_max_ghost
INTEGER,INTENT(IN) :: local_y_min_ghost,local_y_max_ghost
INTEGER,INTENT(IN) :: local_z_min_ghost,local_z_max_ghost
LOGICAL*4,PARAMETER :: debug = GLOBAL_DEBUG
!-----------INTENT(INOUT)-----------------------------------------------------------------------------------------------
LOGICAL*4,INTENT(INOUT),DIMENSION(local_x_min_ghost:local_x_max_ghost, &
& local_y_min_ghost:local_y_max_ghost,local_z_min_ghost:local_z_max_ghost) &
& :: local_walls
!-----------NO INTENT-----------------------------------------------------------------------------------------------
INTEGER :: I,J,K
!USER DEFINED SETUP PARAMETERS ARE ENTERED HERE(IN LATTICE UNITS):
LOGICAL,PARAMETER :: INIT_GEOMETRY_comb_exists_in_XY = .TRUE.
LOGICAL,PARAMETER :: INIT_GEOMETRY_comb_exists_in_YZ = .TRUE.
LOGICAL,PARAMETER :: INIT_GEOMETRY_comb_exists_in_ZX = .TRUE.
INTEGER,PARAMETER :: INIT_GEOMETRY_XY_comb_start_X = 1
INTEGER,PARAMETER :: INIT_GEOMETRY_YZ_comb_start_Y = 1
INTEGER,PARAMETER :: INIT_GEOMETRY_ZX_comb_start_Z = 1
INTEGER,PARAMETER :: INIT_GEOMETRY_XY_comb_end_X = 10
INTEGER,PARAMETER :: INIT_GEOMETRY_YZ_comb_end_Y = 10
INTEGER,PARAMETER :: INIT_GEOMETRY_ZX_comb_end_Z = 10
INTEGER,PARAMETER :: INIT_GEOMETRY_XY_comb_Y_spacing = 3
INTEGER,PARAMETER :: INIT_GEOMETRY_YZ_comb_Z_spacing = 3
INTEGER,PARAMETER :: INIT_GEOMETRY_ZX_comb_X_spacing = 3
INTEGER,PARAMETER :: INIT_GEOMETRY_XY_comb_Z_position = 10
INTEGER,PARAMETER :: INIT_GEOMETRY_YZ_comb_X_position = 10
INTEGER,PARAMETER :: INIT_GEOMETRY_ZX_comb_Y_position = 10
!------------------------------------------------------------------
IF(debug) WRITE (9,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: entered subroutine'
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: entered subroutine'
IF(debug) WRITE (0,*) My_MPI_Process_ID,'entered SUBROUTINE INIT_GEOMETRY_comb'
! IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_z_min_ghost',global_z_min_ghost
! IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_z_max_ghost',global_z_max_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_y_min_ghost',global_y_min_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_y_max_ghost',global_y_max_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_x_min_ghost',global_x_min_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_x_max_ghost',global_x_max_ghost
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_z_min_real',global_z_min_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_z_max_real',global_z_max_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_y_min_real',global_y_min_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_y_max_real',global_y_max_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_x_min_real',global_x_min_real
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: global_x_max_real',global_x_max_real
DO K=local_z_min_ghost,local_z_max_ghost
DO J=local_y_min_ghost,local_y_max_ghost
DO I=local_x_min_ghost,local_x_max_ghost
!!INITIALIZE THE COMB GEOMETRY IN THE XY PLANE
IF(INIT_GEOMETRY_comb_exists_in_XY) THEN
IF (K.EQ.INIT_GEOMETRY_XY_comb_Z_position) THEN
IF( (MOD(J,INIT_GEOMETRY_XY_comb_Y_spacing).EQ.1).AND.&
&( I.GE.INIT_GEOMETRY_XY_comb_start_X ).AND.&
&( I.LE.INIT_GEOMETRY_XY_comb_end_X ) ) THEN
local_walls(I,J,K) = .TRUE.
ENDIF !!(MOD(J,INIT_GEOMETRY_XY_comb_Y_spacing).EQ.1).AND. ...
ENDIF !! K.EQ.INIT_GEOMETRY_XY_comb_Z_position
ENDIF !!INIT_GEOMETRY_comb_exists_in_XY
!!INITIALIZE THE COMB GEOMETRY IN THE YZ PLANE
IF(INIT_GEOMETRY_comb_exists_in_YZ) THEN
IF (I.EQ.INIT_GEOMETRY_YZ_comb_X_position) THEN
IF( (MOD(K,INIT_GEOMETRY_YZ_comb_Z_spacing).EQ.1).AND.&
&( J.GE.INIT_GEOMETRY_YZ_comb_start_Y ).AND.&
&( J.LE.INIT_GEOMETRY_YZ_comb_end_Y ) ) THEN
local_walls(I,J,K) = .TRUE.
ENDIF !!(MOD(K,INIT_GEOMETRY_YZ_comb_Z_spacing).EQ.1).AND. ...
ENDIF !! I.EQ.INIT_GEOMETRY_YZ_comb_X_position
ENDIF !!INIT_GEOMETRY_comb_exists_in_YZ
!!INITIALIZE THE COMB GEOMETRY IN THE ZX PLANE
IF(INIT_GEOMETRY_comb_exists_in_ZX) THEN
IF (J.EQ.INIT_GEOMETRY_ZX_comb_Y_position) THEN
IF( (MOD(I,INIT_GEOMETRY_ZX_comb_X_spacing).EQ.1).AND.&
&( K.GE.INIT_GEOMETRY_ZX_comb_start_Z ).AND.&
&( K.LE.INIT_GEOMETRY_ZX_comb_end_Z ) ) THEN
local_walls(I,J,K) = .TRUE.
ENDIF !!(MOD(I,INIT_GEOMETRY_ZX_comb_X_spacing).EQ.1).AND. ...
ENDIF !! J.EQ.INIT_GEOMETRY_ZX_comb_Y_position
ENDIF !!INIT_GEOMETRY_comb_exists_in_ZX
END DO !K=local_z_min_ghost,local_z_max_ghost
END DO !J=local_y_min_ghost,local_y_max_ghost
END DO !I=local_x_min_ghost,local_x_max_ghost
IF(debug) WRITE (9,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: about to return from subroutine'
IF(debug) WRITE (0,*) My_MPI_Process_ID,'INIT_GEOMETRY_comb: about to return from subroutine'
RETURN
END SUBROUTINE INIT_GEOMETRY_comb
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
SUBROUTINE INIT_GEOMETRY_fiber_layers_check(&
& My_MPI_Process_ID,&
& xx,yy,zz,bundle,xc_bundle,yc_bundle,zc_bundle,temp,&
& flag_bundle,bundle_radius_lu,layer_length_in_lattice_units,&
& layer_width_in_lattice_units,layer_thickness_in_lattice_units,&
& number_of_deltal_bundle,number_of_bundles_per_layer)
!
! ---------------------------------------------------------------------
! Create function ggnqt on the convex with normal distribution
! ---------------------------------------------------------------------
!
! REAL FUNCTION randomnormal(random_seed)
! IMPLICIT NONE
! REAL, EXTERNAL :: ranf
! REAL :: ssumm, zzz, random_seed
! INTEGER :: i
! ssumm=0.d0
! DO i=1,48
! zzz=ranf(random_seed)
! ssumm=ssumm+zzz
! END DO
! randomnormal=(ssumm-24.0d0)*0.5d0
! RETURN
! END FUNCTION randomnormal
!
!--------------------------------------------------------------------------------
! Create random numbers from 0 to 1 on the convex with uniform distribution
!--------------------------------------------------------------------------------
! REAL FUNCTION ranf(random_seed)
! IMPLICIT NONE
! REAL :: random_seed, oldseed, d2p31m, d2p31
! d2p31m=2147483647.0d0
! d2p31=2147483648.0d0
! oldseed=random_seed
! random_seed=mod(16807.d0*random_seed,d2p31m)
! ranf=random_seed/d2p31
! END FUNCTION ranf
!
!----------------------------------------------------------------------------------------------------------------------------
! Calculation of whether a computational cell contains a whether one of its corners
! is inside a cylinder. First, calculate the angle of this point of interrest and the 2 ends of the cylinders
! then, calculate the distance between this point and the cylinder axis to see if it is smaller that the radius or not.
!----------------------------------------------------------------------------------------------------------------------------
!
USE debug_module
IMPLICIT NONE
INCLUDE "mpif.h"
!-----------INTENT(IN)------------------------------------------
INTEGER,INTENT(IN) :: My_MPI_Process_ID
INTEGER,INTENT(IN) :: xx,yy,zz,bundle
INTEGER,INTENT(IN) :: layer_length_in_lattice_units,layer_width_in_lattice_units,layer_thickness_in_lattice_units
INTEGER,INTENT(IN) :: number_of_deltal_bundle,number_of_bundles_per_layer
REAL,INTENT(IN) :: bundle_radius_lu
REAL,DIMENSION(number_of_deltal_bundle,number_of_bundles_per_layer),&
& INTENT(IN) :: xc_bundle,yc_bundle,zc_bundle
LOGICAL*4,PARAMETER :: debug = GLOBAL_DEBUG
!-----------INTENT(OUT)------------------------------------------
!-----------INTENT(INOUT)------------------------------------------
INTEGER, INTENT(INOUT) :: temp
INTEGER,DIMENSION(0:layer_length_in_lattice_units,0:layer_width_in_lattice_units,0:layer_thickness_in_lattice_units),&
& INTENT(INOUT)::flag_bundle
!-----------NO INTENT------------------------------------------
INTEGER :: x6,y6,z6,count
REAL :: distance_xc1_xc2,distance_walker_xc1
REAL :: distance_walker_xc2,distance_walker_bundle,angel1,angel2
IF(debug) WRITE (0,*) My_MPI_Process_ID,'entered SUBROUTINE INIT_GEOMETRY_fiber_layers_check'
IF(debug) WRITE (9,*) My_MPI_Process_ID,'entered SUBROUTINE INIT_GEOMETRY_fiber_layers_check'
count=0
temp = 0
!
DO x6 = xx,xx+1
DO y6 = yy,yy+1
DO z6 = zz,zz+1
distance_xc1_xc2 = ABS(SQRT((xc_bundle(number_of_deltal_bundle,bundle) &
& -xc_bundle(1,bundle))**2 +(yc_bundle(number_of_deltal_bundle,bundle) &
& -yc_bundle(1,bundle))**2 +(zc_bundle(number_of_deltal_bundle,bundle) &
& -zc_bundle(1,bundle))**2))
!
distance_walker_xc1 = ABS(SQRT((REAL(x6) - xc_bundle(1,bundle))**2 &
& +(REAL(y6)-yc_bundle(1,bundle))**2+(REAL(z6)-zc_bundle(1,bundle))**2))
!
distance_walker_xc2 = ABS(SQRT((REAL(x6) - xc_bundle(number_of_deltal_bundle,bundle)) &
& **2+(REAL(y6) - yc_bundle(number_of_deltal_bundle,bundle))**2+(REAL(z6) &
& - zc_bundle(number_of_deltal_bundle,bundle))**2))
!
angel1= ((distance_walker_xc1)**2+(distance_xc1_xc2)**2 &
& -(distance_walker_xc2)**2) &
& /(2.0d0*SQRT(distance_walker_xc1)*SQRT(distance_xc1_xc2))
angel2= ((distance_walker_xc2)**2+(distance_xc1_xc2)**2 &
& -(distance_walker_xc1)**2) &
& /(2.0d0*SQRT(distance_walker_xc2)*SQRT(distance_xc1_xc2))
!
IF ((angel1.GE.0.0d0).AND.(angel2.GE.0.0d0)) THEN
distance_walker_bundle = distance_walker_xc1 &
& * SQRT(1.0d0- (((distance_xc1_xc2)**2+ &
& (distance_walker_xc1)**2-(distance_walker_xc2) &
& **2)/(2.0d0*distance_xc1_xc2 &
& *distance_walker_xc1))**2)
IF(distance_walker_bundle.LE.bundle_radius_lu) THEN
count=count+1
END IF !!distance_walker_bundle.LE.bundle_radius_lu
END IF !!((angel1.GE.0.0d0).AND.(angel2.GE.0.0d0))
END DO !z6
END DO !y6
END DO !x6
IF (count.EQ.0) THEN
temp = 0
ELSE !count>0)
IF ((flag_bundle(xx,yy,zz).EQ.0).OR.(flag_bundle(xx,yy,zz).EQ.bundle)) THEN
temp = 1
END IF !!((flag_bundle(xx,yy,zz).EQ.0).OR.(flag_bundle(xx,yy,zz).EQ.bundle))
IF ((flag_bundle(xx,yy,zz).NE.bundle).AND.(flag_bundle(xx,yy,zz).NE.0)) THEN
temp = 2
END IF !((flag_bundle(xx,yy,zz).NE.bundle).AND.(flag_bundle(xx,yy,zz).NE.0))
END IF !count.EQ. !
IF(debug) WRITE (0,*) My_MPI_Process_ID,'about to END SUBROUTINE INIT_GEOMETRY_fiber_layers_check'
IF(debug) WRITE (9,*) My_MPI_Process_ID,'about to END SUBROUTINE INIT_GEOMETRY_fiber_layers_check'
RETURN
END SUBROUTINE INIT_GEOMETRY_fiber_layers_check
SUBROUTINE INIT_GEOMETRY_fiber_layers_check_tube(&
& My_MPI_Process_ID,&
& xx,yy,zz,tube,xc_tube,yc_tube,zc_tube,temp,&
& flag_tube,tube_radius_lu,layer_length_in_lattice_units,&
& layer_width_in_lattice_units,layer_thickness_in_lattice_units,&
& number_of_deltal_tube,number_of_tubes_per_layer)
!
! ---------------------------------------------------------------------
! Create function ggnqt on the convex with normal distribution
! ---------------------------------------------------------------------
!
! REAL FUNCTION randomnormal(random_seed)
! IMPLICIT NONE
! REAL, EXTERNAL :: ranf
! REAL :: ssumm, zzz, random_seed
! INTEGER :: i
! ssumm=0.d0
! DO i=1,48
! zzz=ranf(random_seed)
! ssumm=ssumm+zzz
! END DO
! randomnormal=(ssumm-24.0d0)*0.5d0
! RETURN
! END FUNCTION randomnormal
!
!--------------------------------------------------------------------------------
! Create random numbers from 0 to 1 on the convex with uniform distribution
!--------------------------------------------------------------------------------
! REAL FUNCTION ranf(random_seed)
! IMPLICIT NONE
! REAL :: random_seed, oldseed, d2p31m, d2p31
! d2p31m=2147483647.0d0
! d2p31=2147483648.0d0
! oldseed=random_seed
! random_seed=mod(16807.d0*random_seed,d2p31m)
! ranf=random_seed/d2p31
! END FUNCTION ranf
!
!----------------------------------------------------------------------------------------------------------------------------
! Calculation of whether a computational cell contains a whether one of its corners
! is inside a cylinder. First, calculate the angle of this point of interrest and the 2 ends of the cylinders
! then, calculate the distance between this point and the cylinder axis to see if it is smaller that the radius or not.
!----------------------------------------------------------------------------------------------------------------------------
!
USE debug_module
IMPLICIT NONE
INCLUDE "mpif.h"
!-----------INTENT(IN)------------------------------------------
INTEGER,INTENT(IN) :: My_MPI_Process_ID
INTEGER,INTENT(IN) :: xx,yy,zz,tube
INTEGER,INTENT(IN) :: layer_length_in_lattice_units,layer_width_in_lattice_units,layer_thickness_in_lattice_units
INTEGER,INTENT(IN) :: number_of_deltal_tube,number_of_tubes_per_layer
REAL,INTENT(IN) :: tube_radius_lu
REAL,DIMENSION(number_of_deltal_tube,number_of_tubes_per_layer),&
& INTENT(IN) :: xc_tube,yc_tube,zc_tube
LOGICAL*4,PARAMETER :: debug = GLOBAL_DEBUG
!-----------INTENT(OUT)------------------------------------------
!-----------INTENT(INOUT)------------------------------------------
INTEGER, INTENT(INOUT) :: temp
INTEGER,DIMENSION(0:layer_length_in_lattice_units,0:layer_width_in_lattice_units,0:layer_thickness_in_lattice_units),&
& INTENT(INOUT)::flag_tube
!-----------NO INTENT------------------------------------------
INTEGER :: x6,y6,z6,count
REAL :: distance_xc1_xc2,distance_walker_xc1
REAL :: distance_walker_xc2,distance_walker_tube,angel1,angel2
IF(debug) WRITE (0,*) My_MPI_Process_ID,'entered SUBROUTINE INIT_GEOMETRY_fiber_layers_check_tube'
IF(debug) WRITE (9,*) My_MPI_Process_ID,'entered SUBROUTINE INIT_GEOMETRY_fiber_layers_check_tube'
count=0
temp = 0
!
DO x6 = xx,xx+1
DO y6 = yy,yy+1
DO z6 = zz,zz+1
! write(9999,*), 'tube=', tube
distance_xc1_xc2 = ABS(SQRT((xc_tube(number_of_deltal_tube,tube) &
& -xc_tube(1,tube))**2 +(yc_tube(number_of_deltal_tube,tube) &
& -yc_tube(1,tube))**2 +(zc_tube(number_of_deltal_tube,tube) &
& -zc_tube(1,tube))**2))
!
distance_walker_xc1 = ABS(SQRT((REAL(x6) - xc_tube(1,tube))**2 &
& +(REAL(y6)-yc_tube(1,tube))**2+(REAL(z6)-zc_tube(1,tube))**2))
!
distance_walker_xc2 = ABS(SQRT((REAL(x6) - xc_tube(number_of_deltal_tube,tube)) &
& **2+(REAL(y6) - yc_tube(number_of_deltal_tube,tube))**2+(REAL(z6) &
& - zc_tube(number_of_deltal_tube,tube))**2))
!
angel1= ((distance_walker_xc1)**2+(distance_xc1_xc2)**2 &
& -(distance_walker_xc2)**2) &
& /(2.0d0*SQRT(distance_walker_xc1)*SQRT(distance_xc1_xc2))
angel2= ((distance_walker_xc2)**2+(distance_xc1_xc2)**2 &
& -(distance_walker_xc1)**2) &
& /(2.0d0*SQRT(distance_walker_xc2)*SQRT(distance_xc1_xc2))
!
IF ((angel1.GE.0.0d0).AND.(angel2.GE.0.0d0)) THEN
distance_walker_tube = distance_walker_xc1 &
& * SQRT(1.0d0- (((distance_xc1_xc2)**2+ &
& (distance_walker_xc1)**2-(distance_walker_xc2) &
& **2)/(2.0d0*distance_xc1_xc2 &
& *distance_walker_xc1))**2)
IF(distance_walker_tube.LE.tube_radius_lu) THEN
count=count+1
END IF !!distance_walker_tube.LE.tube_radius_lu
END IF !!((angel1.GE.0.0d0).AND.(angel2.GE.0.0d0))
END DO !z6
END DO !y6
END DO !x6