! ! This file is copyrighted by Colorado ! State University. It may be distributed ! freely, as long as it is distributed ! completely, including this copyright ! message. ! MODULE kinds ! IMPLICIT NONE INTEGER, PARAMETER :: kindi = SELECTED_INT_KIND(9) INTEGER, PARAMETER :: kindr = SELECTED_REAL_KIND(p=11) END MODULE kinds ! ------------------------------------------------------------------- MODULE types USE kinds TYPE bound_plane INTEGER(kind=kindi) :: mat_num INTEGER(kind=kindi) :: int_level INTEGER(kind=kindi) :: iabs INTEGER(kind=kindi) :: np_abs INTEGER(kind=kindi) :: mod_num REAL(kind=kindr), DIMENSION(3) :: n REAL(kind=kindr) :: d REAL(kind=kindr), DIMENSION(4,3) :: V REAL(kind=kindr), DIMENSION(3) :: o REAL(kind=kindr), DIMENSION(2,3) :: M1,M2 REAL(kind=kindr), DIMENSION(3,2) :: W1,W2 REAL(kind=kindr) :: area REAL(kind=kindr), DIMENSION(3) :: max REAL(kind=kindr), DIMENSION(3) :: min END TYPE bound_plane TYPE circular_cyl INTEGER(kind=kindi) :: mat_num INTEGER(kind=kindi) :: int_level INTEGER(kind=kindi) :: iabs INTEGER(kind=kindi) :: np_abs INTEGER(kind=kindi) :: mod_num REAL(kind=kindr) :: r REAL(kind=kindr), DIMENSION(3) :: o REAL(kind=kindr) :: xmin REAL(kind=kindr) :: xmax REAL(kind=kindr) :: r_sqrd REAL(kind=kindr), DIMENSION(3) :: max REAL(kind=kindr), DIMENSION(3) :: min END TYPE circular_cyl TYPE opaque_material REAL(kind=kindr) :: rho_s REAL(kind=kindr) :: rho_ss REAL(kind=kindr) :: rho_d REAL(kind=kindr) :: r_ss REAL(kind=kindr) :: r_d REAL(kind=kindr), DIMENSION(3) :: pdf REAL(kind=kindr) :: rp1inv_ss REAL(kind=kindr) :: rp1inv_d END TYPE opaque_material TYPE transparent_material REAL(kind=kindr) :: n_glass REAL(kind=kindr) :: k_glass REAL(kind=kindr) :: t REAL(kind=kindr), DIMENSION(0:91) :: rho_glass REAL(kind=kindr), DIMENSION(0:91) :: tau_glass REAL(kind=kindr), DIMENSION(0:91) :: drho_glass REAL(kind=kindr), DIMENSION(0:91) :: dtau_glass END TYPE transparent_material END MODULE types ! ------------------------------------------------------------------- MODULE constants USE kinds ! REAL(kind=kindr) :: pi, p2i, pi2, pi4, pi2inv, pi4inv, dtor, rtod END MODULE constants ! ------------------------------------------------------------------- MODULE control USE kinds ! CHARACTER(len=1), DIMENSION(72) :: head INTEGER(kind=kindi) :: num_photons, num_ploops,max_int_level INTEGER(kind=kindi) :: seed0, seed_from_time, nploop INTEGER(kind=kindi) :: idata, iprint REAL(kind=kindr) :: tol ! END MODULE control ! ------------------------------------------------------------------- MODULE directions USE kinds USE types ! INTEGER(kind=kindi) :: n_long, n_tran, num_dirs, norm REAL(kind=kindr), ALLOCATABLE, DIMENSION(:) :: thetatd, thetald REAL(kind=kindr), ALLOCATABLE, DIMENSION(:) :: thetat, thetal REAL(kind=kindr), ALLOCATABLE, DIMENSION(:) :: altituded,phid TYPE(bound_plane), ALLOCATABLE, DIMENSION(:) :: emission END MODULE directions ! ------------------------------------------------------------------- MODULE materials USE kinds USE types ! INTEGER(kind=kindi) :: num_mats_opaq, num_mats_glass TYPE(opaque_material), ALLOCATABLE, DIMENSION(:) :: opaque TYPE(transparent_material), ALLOCATABLE, DIMENSION(:) :: transparent END MODULE materials ! ------------------------------------------------------------------- MODULE collector USE kinds USE types ! INTEGER(kind=kindi), PARAMETER :: max_of_type=2 INTEGER(kind=kindi), DIMENSION(-1:max_of_type) :: num_of_type INTEGER(kind=kindi) :: num_of_modules TYPE(bound_plane) :: aperture TYPE(bound_plane), ALLOCATABLE, DIMENSION(:) :: surf_type0 TYPE(bound_plane), ALLOCATABLE, DIMENSION(:) :: surf_type1 TYPE(circular_cyl), ALLOCATABLE, DIMENSION(:) :: surf_type2 REAL(kind=kindr), ALLOCATABLE, DIMENSION(:) :: y0_save END MODULE collector ! ------------------------------------------------------------------- MODULE photons USE kinds ! INTEGER(kind=kindi), ALLOCATABLE, DIMENSION(:) :: npabs_a, npabs_glass, & npabs_coll, npabs_surf, npabs_sum, npabs_tot, npabs_caps INTEGER(kind=kindi), ALLOCATABLE, DIMENSION(:) :: np_glass,np_coll, & np_surf,np_sum INTEGER(kind=kindi), ALLOCATABLE, DIMENSION(:) :: nplost,np_error REAL(kind=kindr), ALLOCATABLE, DIMENSION (:) :: iam, tau_alpha REAL(kind=kindr) :: tad,glass_frac END MODULE photons ! ------------------------------------------------------------------- MODULE dates USE kinds ! CHARACTER(len=8) :: adate, atime INTEGER :: valuesa(8,0:2) END MODULE dates ! ------------------------------------------------------------------- MODULE random_funcs USE kinds USE constants IMPLICIT NONE ! INTEGER(kind=kindi), PARAMETER :: modulus = 2147483647 INTEGER(kind=kindi) :: seed(17), it1 = 17, it2 = 12 ! REAL(kind=kindr), PARAMETER :: odenom = 1.0_kindr/(2.0_kindr**31) CONTAINS ! !*********************************************************************** SUBROUTINE raninit(iprnt,seed0) !*********************************************************************** !* subroutine to initialize pseudo-random number generator * !*********************************************************************** ! USE kinds ! IMPLICIT NONE ! INTEGER(kind=kindi) :: n, iprnt, seed0 INTEGER(kind=kindi), PARAMETER :: ia_mult = 16807 ! ! Generate initial seeds using multiplicative, congruential generator ! seed(1) = MOD(ia_mult*seed0,modulus) DO n = 2, 17 seed(n) = MOD(ia_mult*seed(n-1),modulus) IF (seed(n) < 0) THEN seed(n) = seed(n) + modulus END IF END DO ! ! Check for an odd value ! DO n=1, 17 IF (MOD(seed(n),2) == 1) THEN go to 250 END IF END DO seed(1) = IOR (seed(1), 1) ! ! Output of initial variables ! 250 CONTINUE IF (iprnt < 2 ) RETURN CALL header (9) WRITE (9,1000) n = 0 WRITE (9,1001) n,seed0,seed0 DO n=1, 17 WRITE (9,1001) n,seed(n),seed(n) END DO ! RETURN ! 1000 FORMAT (/10x,'I n i t i a l S e e d s'/ & / 7x,'n',12x,'seed',8x,'seed' & / 20x,'(dec)',8x,'(hex)') 1001 FORMAT (5x,i5,5x,i12,5x,z8) ! END SUBROUTINE raninit ! !*********************************************************************** FUNCTION ranf() !*********************************************************************** !* function to generate pseudo-random numbers * !*********************************************************************** ! IMPLICIT NONE ! REAL(kind=kindr) :: ranf ! ! Get new seed ! seed(it1) = seed(it1) + seed(it2) seed(it1) = IAND(seed(it1), modulus) ! ! Generate normalized random number ! ranf = seed(it1)*odenom ! ! Cyclic permutation of taps ! it1 = it1 - 1 it2 = it2 - 1 IF (it1 < 1) THEN it1 = 17 ELSEIF (it2 < 1) THEN it2 = 17 END IF ! RETURN ! END FUNCTION ranf END MODULE random_funcs ! ------------------------------------------------------------------- ! !********************************************************************* ! MODULE type1_funcs ! !********************************************************************* ! USE kinds USE constants IMPLICIT NONE CONTAINS ! !*********************************************************************** ! FUNCTION normal(V) ! !*********************************************************************** !* function to surface normal * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(3) :: normal REAL(kind=kindr), DIMENSION(4,3), INTENT(IN) :: V REAL(kind=kindr), DIMENSION(3) :: v1,v2 v1=V(3,:)-V(1,:) v2=V(2,:)-V(1,:) normal(1)=v1(2)*v2(3)-v2(2)*v1(3) normal(2)=v2(1)*v1(3)-v1(1)*v2(3) normal(3)=v1(1)*v2(2)-v2(1)*v1(2) IF (normal(3)<0) normal=-normal normal=normal/SQRT(DOT_PRODUCT(normal,normal)) END FUNCTION normal ! !*********************************************************************** ! FUNCTION origin(V) ! !*********************************************************************** !* function to set surface origin * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(3) :: origin REAL(kind=kindr), DIMENSION(4,3), INTENT(IN) :: V origin=V(1,:) END FUNCTION origin ! !*********************************************************************** ! FUNCTION matrix_M(V,i1,i2) ! !*********************************************************************** !* function to set matrix M * !*********************************************************************** ! IMPLICIT NONE INTEGER(kind=kindi), INTENT(IN) :: i1,i2 REAL(kind=kindr), DIMENSION(2,3) :: matrix_M REAL(kind=kindr), DIMENSION(4,3), INTENT(IN) :: V matrix_M(1,:)=V(i1,:)-V(1,:) matrix_M(2,:)=V(i2,:)-V(1,:) END FUNCTION matrix_M ! !*********************************************************************** ! FUNCTION matrix_W(M) ! !*********************************************************************** !* function to calculate conditional inverse of matrix M * !*********************************************************************** ! IMPLICIT NONE INTEGER(kind=kindi) :: i,j,ii,jj REAL(kind=kindr), DIMENSION(3,2) :: matrix_W REAL(kind=kindr), DIMENSION(2,3), INTENT(IN) :: M REAL(kind=kindr), DIMENSION(3) ::dets INTEGER(kind=kindi), DIMENSION(1) :: ip REAL(kind=kindr) :: det REAL(kind=kindr), DIMENSION(2) ::v REAL(kind=kindr), DIMENSION(2,2) ::Sub,Inv dets(1)=M(1,2)*M(2,3)-M(2,2)*M(1,3) dets(2)=M(1,1)*M(2,3)-M(2,1)*M(1,3) dets(3)=M(1,1)*M(2,2)-M(2,1)*M(1,2) ip=MAXLOC(ABS(dets)) det=dets(ip(1)) v=PACK((/1,2,3/),(/1,2,3/)/=ip(1)) DO j=1,2 jj=v(j) Sub(:,j)=M(:,jj) END DO Inv(1,1)=Sub(2,2) Inv(1,2)=-Sub(1,2) Inv(2,1)=-Sub(2,1) Inv(2,2)=Sub(1,1) matrix_W=0.0_kindr DO i=1,2 ii=v(i) matrix_W(ii,:)=Inv(i,:)/det END DO END FUNCTION matrix_W ! !*********************************************************************** ! FUNCTION b_coord(x,o,W) ! !*********************************************************************** !* function to calculate barycentric coordinates from cartesian * !* coordinates * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(2) :: b_coord REAL(kind=kindr), DIMENSION(3), INTENT(IN) :: x,o REAL(kind=kindr), DIMENSION(3,2), INTENT(IN) :: W b_coord=MATMUL((x-o),W) END FUNCTION b_coord ! !*********************************************************************** ! FUNCTION c_coord(p,o,M) ! !*********************************************************************** !* FUNCTION to calculate cartesian coordinates from barycentric * !* coordinates * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(3) :: c_coord REAL(kind=kindr), DIMENSION(2), INTENT(IN) :: p REAL(kind=kindr), DIMENSION(3), INTENT(IN) :: o REAL(kind=kindr), DIMENSION(2,3), INTENT(IN) :: M INTEGER(kind=kindi) i DO i=1,3 c_coord(i)=o(i)+p(1)*M(1,i)+p(2)*M(2,i) END DO END FUNCTION c_coord ! !*********************************************************************** ! ! FUNCTION area ! !*********************************************************************** !* FUNCTION to calculates the area of the parallelogram with sides * !* M(1,:) and M(2,:) * !*********************************************************************** ! FUNCTION area(M) IMPLICIT NONE REAL(kind=kindr) :: area REAL(kind=kindr), DIMENSION(2,3) , INTENT(IN) :: M area = (M(1,2)*M(2,3)-M(2,2)*M(1,3))**2 & +(M(1,1)*M(2,3)-M(2,1)*M(1,3))**2 & +(M(1,1)*M(2,2)-M(2,1)*M(1,2))**2 area=SQRT(area) END FUNCTION area END MODULE type1_funcs ! !********************************************************************* ! MODULE surface_interactions ! !********************************************************************* ! USE kinds USE constants USE random_funcs USE materials IMPLICIT NONE REAL(kind=kindr), DIMENSION(3) :: photon_orig,photon_dir CONTAINS ! !*********************************************************************** ! FUNCTION interaction(mat_num,i,n) ! !*********************************************************************** !* SUBROUTINE to determine photon-material interactions * !* output is interaction and new direction (not new X) * !* interaction = 2 for error * !* 1 for transmission * !* 0 for absorption * !* -1 for diffuse reflection * !* -2 for specular reflection * !* -3 for semi-specular reflection * !* Note: except for absorption, glass is treated as infinitely thin * !*********************************************************************** ! IMPLICIT NONE ! INTEGER (kind=kindi), INTENT(IN) :: mat_num ! material number REAL(kind=kindr), DIMENSION(:), INTENT(INOUT) :: i ! incident direction REAL(kind=kindr), DIMENSION(:), INTENT(IN) :: n ! surface normal direction ! INTEGER (kind=kindi) :: interaction ! returns type of interaction ! INTEGER(kind=kindi) :: mn INTEGER(kind=kindi) :: itheta_i REAL(kind=kindr) :: rn_in REAL(kind=kindr) :: theta_i, rhos, taus, frac ! mn=mat_num interaction = 2 rn_in = ranf() ! ! ***** Opaque materials ***** ! IF (mn>0) THEN IF (rn_in > opaque(mn)%pdf(3)) THEN ! Absorbed interaction = 0 ELSE IF (rn_in > opaque(mn)%pdf(2)) THEN ! Diffuse reflection interaction = -1 photon_dir=diffuse(n,opaque(mn)%rp1inv_d) ELSE IF (rn_in > opaque(mn)%pdf(1)) THEN ! Specular reflection interaction = -2 photon_dir=specular(i,n) ELSE ! Semi-specular reflection interaction = -3 photon_dir=semi_specular(i,n,opaque(mn)%rp1inv_ss) END IF ELSE mn=-mn ! ! ***** Tranparent materials ***** ! ! Interpolation on incident angle to determine material properties. ! theta_i = ACOS(-DOT_PRODUCT(i,n))*rtod itheta_i = theta_i frac = theta_i - itheta_i rhos = transparent(mn)%rho_glass(itheta_i) + frac*transparent(mn)%drho_glass(itheta_i) taus = transparent(mn)%tau_glass(itheta_i) + frac*transparent(mn)%dtau_glass(itheta_i) ! IF (rn_in < taus) THEN ! Transmitted (no change in direction) interaction = 1 ELSE IF (rn_in < taus + rhos) THEN ! Specular reflection interaction = -2 photon_dir=specular(i,n) ELSE ! Absorbed interaction = 0 END IF END IF ! END FUNCTION interaction ! !*********************************************************************** ! FUNCTION specular(i,n) ! !*********************************************************************** !* FUNCTION to calculate direction of specular reflection * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(:), INTENT(IN) :: i ! incident direction REAL(kind=kindr), DIMENSION(:), INTENT(IN) :: n ! surface normal direction REAL(kind=kindr), DIMENSION(3) :: specular ! returns reflected direction REAL(kind=kindr) :: q1 q1=-2.0_kindr*DOT_PRODUCT(i,n) specular=i+q1*n END FUNCTION specular ! !*********************************************************************** ! FUNCTION diffuse(n,rp1inv) ! !*********************************************************************** !* FUNCTION to calculate direction of diffuse reflection * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(:), INTENT(IN) :: n ! surface normal direction REAL(kind=kindr), INTENT(IN) :: rp1inv ! 1/(1+r) for the material REAL(kind=kindr), DIMENSION(3) :: diffuse ! returns reflected direction REAL(kind=kindr), DIMENSION(3) :: r REAL(kind=kindr) :: theta_o,phi_o,s_theta_o REAL(kind=kindr) :: q1,q2,q3,q4 ! ! Calculate the random angles about the (x',y',z') surface coordinate system. ! theta_o=ACOS(ranf()**rp1inv) phi_o=p2i*ranf() ! ! Calculate the direction of the reflected vector in (x',y',z') coordinate system. ! s_theta_o=SIN(theta_o) r(1)=COS(phi_o)*s_theta_o r(2)=SIN(phi_o)*s_theta_o r(3)=COS(theta_o) ! ! Transform the reflected vector into the world coordinate system. ! IF (n(1)==0) THEN diffuse(1)=r(1) diffuse(2)=n(3)*r(2)+n(2)*r(3) diffuse(3)=-n(2)*r(2)+n(3)*r(3) ELSE q1=SQRT(n(1)**2+n(2)**2) q2=n(1)/q1 q3=n(2)/q1 q4=n(3)*r(2) diffuse(1)=q3*r(1)+q2*q4+n(1)*r(3) diffuse(2)=-q2*r(1)+q3*q4+n(2)*r(3) diffuse(3)=-q1*r(2)+n(3)*r(3) END IF END FUNCTION diffuse ! !*********************************************************************** ! FUNCTION semi_specular(i,n,rp1inv) ! !*********************************************************************** !* FUNCTION to direction of semi specular reflection * !*********************************************************************** ! IMPLICIT NONE REAL(kind=kindr), DIMENSION(:), INTENT(IN) :: i ! incident direction REAL(kind=kindr), DIMENSION(:), INTENT(IN) :: n ! surface normal direction REAL(kind=kindr), INTENT(IN) :: rp1inv ! 1/(1+r) for the material REAL(kind=kindr), DIMENSION(3) :: semi_specular ! returns reflected direction REAL(kind=kindr), DIMENSION(3) :: s,r REAL(kind=kindr) :: theta_o,phi_o,s_theta_o REAL(kind=kindr) :: q1,q2,q3,q4 ! ! Calculate the direction for purely specular reflection. ! s=specular(i,n) ! ! Calculate the random angles about the (x',y',z') surface coordinate system. ! theta_o=(pi2-ACOS(-DOT_PRODUCT(i,n)))*ACOS(ranf()**rp1inv)/pi2 phi_o=p2i*ranf() ! ! Calculate the direction of the reflected vector in (x',y',z') coordinate system. ! s_theta_o=SIN(theta_o) r(1)=COS(phi_o)*s_theta_o r(2)=SIN(phi_o)*s_theta_o r(3)=COS(theta_o) ! ! Transform the reflected vector into the world coordinate system. ! IF (n(1)==0) THEN semi_specular(1)=r(1) semi_specular(2)=s(3)*r(2)+s(2)*r(3) semi_specular(3)=-s(2)*r(2)+s(3)*r(3) ELSE q1=SQRT(s(1)**2+s(2)**2) q2=s(1)/q1 q3=s(2)/q1 q4=s(3)*r(2) semi_specular(1)=q3*r(1)+q2*q4+s(1)*r(3) semi_specular(2)=-q2*r(1)+q3*q4+s(2)*r(3) semi_specular(3)=-q1*r(2)+s(3)*r(3) END IF END FUNCTION semi_specular END MODULE surface_interactions ! !********************************************************************* ! * PROGRAM main ! * !********************************************************************* !* program to compute IAM's for simple cylindrical geometries * !********************************************************************* ! ! - modules ! USE kinds USE random_funcs USE control USE constants USE directions USE materials USE collector USE photons USE dates USE type1_funcs ! IMPLICIT NONE ! ! Variable definitions ! CHARACTER(len=1), DIMENSION(40) :: headin CHARACTER(len=8) :: version CHARACTER(len=8) :: adate_compiled INTEGER(kind=kindi) :: nd INTEGER(kind=kindi) :: i, iu, nlines7, npdo, m, n INTEGER(kind=kindi) :: nd_converged, nd_not_converged REAL(kind=kindr) :: tol_check, zz = 1.96 REAL(kind=kindr) :: areaa ! !****************************************************************************** ! ! Initialization phase. ! !****************************************************************************** ! ! Get initial date and time. ! CALL getdt (0) CALL getdt (1) ! Just in case of a fatal error during input ! ! Version information. ! adate_compiled(1:8) = '06-01-98' version(1:8) = 'ver. 3.0' ! ! Initialize constants. ! pi = 4.0_kindr*ATAN(1.0_kindr) p2i = 2.0_kindr*pi pi2 = 0.5d0*pi pi4 = 0.25d0*pi pi4inv = 1.0_kindr/pi4 dtor = pi/180.0_kindr rtod = 1.0_kindr/dtor ! ! OPEN files ! ! - file 1: input file w/o comments (& in col. 1), containing: ! + control variables ! + directions ! + material properties ! + geometry ! - file 2: output file of iam's in tabular FORMAT ! - file 7: input file containing: ! + control variables ! + directions ! + material properties ! + geommetry ! - file 9: output file ! iu = 1 OPEN (iu, file = 'iam.tp1', status = 'unknown', err = 4000) iu = 2 OPEN (iu, file = 'iam.tab', status = 'unknown', err = 4000) iu = 7 OPEN (iu, file = 'iam.in', status = 'old', err = 4000) iu = 9 OPEN (iu, file = 'iam.out', status = 'unknown', err = 4000) ! WRITE (*,*) ' ' WRITE (*,*) ' Copyrighted by Colorado State University' WRITE (*,*) ' ' IF (iprint>0) THEN WRITE (9,*) ' Copyrighted by Colorado State University' WRITE (9,*) ' ' END IF ! ! Convert input file (strip comments w/ & in col. 1). ! CALL files (7, 1, nlines7) CLOSE (7) WRITE (*,*) ' ', nlines7,' - lines read from input file ' WRITE (*,2004) ! !****************************************************************************** ! ! Read data file. ! !****************************************************************************** ! ! Read file header. ! READ (1, 1000) (headin(i), i=1,40) ! ! Write title page. ! head = ' ' head(1:40) = headin(1:40) DO i = 1,8 head(54+i) = adate(i:i) head(64+i) = atime(i:i) END DO ! ! Read photon control information. ! READ (1, *) num_photons, num_ploops, seed0, tol ! ! Select initial seed, If 0, get from time. ! IF (seed0 == 0) THEN seed0 = seed_from_time ELSE IF (seed0 < 0) THEN WRITE (*,*) 'Error -- seed0 must be >= 0 -- stopping' WRITE (9,*) 'Error -- seed0 must be >= 0 -- stopping' CALL adios (2) END IF IF (MOD(seed0, 2) /= 1) THEN WRITE (*,*) 'Error -- seed0 must be odd - adding one and continuing' WRITE (9,*) 'Error -- seed0 must be odd - adding one and continuing' seed0 = seed0 + 1 END IF ! ! Read print control and debug information ! READ (1,*) iprint, idata IF (iprint>0) THEN CALL header (9) WRITE (9, 2010) head(1:40), adate, atime, version, adate_compiled WRITE (9, 2011) num_photons, num_ploops, tol, seed0 WRITE (9, 2014) iprint, idata END IF ! ! Initialize RNG ! CALL raninit (iprint, seed0) ! ! Read direction information and allocate arrays. ! READ (1, *) n_long, n_tran num_dirs = n_long*n_tran WRITE (*,*) 'Allocating direction arrays of length - ', num_dirs IF (num_dirs > 0) THEN ALLOCATE (thetald(n_long+1), thetatd(n_tran+1)) ALLOCATE (thetal(num_dirs), thetat(num_dirs)) ALLOCATE (altituded(num_dirs), phid(num_dirs)) ALLOCATE (emission(num_dirs)) END IF ! CALL read_directions ! ! Read material information and allocate arrays. ! READ (1, *) num_mats_opaq,num_mats_glass IF (num_mats_glass > 0) THEN WRITE (*,*) 'Allocating opaque material arrays of length - ', num_mats_opaq ALLOCATE (opaque(num_mats_opaq)) END IF IF (num_mats_glass > 0) THEN WRITE (*,*) 'Allocating glass material arrays of length - ', num_mats_glass ALLOCATE (transparent(num_mats_glass)) END IF CALL read_materials ! ! Read collector geometry information. ! CALL read_collector ! ! Absorption Arrays ! WRITE (*,*) 'Allocating absorption arrays' ! ! - 1-D arrays ! ALLOCATE (np_glass(0:num_of_modules)) ALLOCATE (np_coll(0:num_of_modules)) ALLOCATE (np_surf(0:num_of_modules)) ALLOCATE (np_sum(0:num_of_modules)) ALLOCATE (npabs_a(num_dirs)) ALLOCATE (npabs_glass(num_dirs), npabs_coll(num_dirs)) ALLOCATE (npabs_surf(num_dirs), npabs_sum(num_dirs)) ALLOCATE (npabs_tot(num_dirs), npabs_caps(num_dirs)) ALLOCATE (nplost(num_dirs), np_error(num_dirs)) ALLOCATE (iam(num_dirs), tau_alpha(num_dirs)) ! ! Finished input phase ! WRITE (*, *) ' Finished Data Check ' IF (iprint>1) THEN WRITE (9,*) WRITE (9, *) ' Finished Data Check ' END IF CALL getdt (1) ! ! data check termination ! IF (idata /= 0) THEN CALL adios(1) END IF ! ! Initializations for run. ! nd_converged = 0 nd_not_converged = 0 nplost = 0 np_error = 0 npabs_caps = 0 ! ! Write convergence header. ! WRITE (*,2003) ! !****************************************************************************** ! ! Start of incident directions loop. ! !****************************************************************************** ! DO nd=1,num_dirs ! ! Initialize arrays counting absorbed photons. ! aperture%np_abs=0 surf_type0%np_abs=0 surf_type1%np_abs=0 surf_type2%np_abs=0 ! ! Define the photon emission plane. ! CALL geteplane(nd) ! !****************************************************************************** ! ! Start of photon convergence loop - in blocks of num_photons. ! !****************************************************************************** ! npdo = 0 DO nploop=1,num_ploops npdo = npdo + 1 CALL photon_loop(nd) ! ! Photon total and IAM calculations. ! ! Total number of emitted photons which hit the aperture plane. ! npabs_a(nd)=aperture%np_abs ! ! Total number of photons collected on collector surfaces ! npabs_glass(nd)=0 npabs_coll(nd)=0 npabs_surf(nd)=0 IF (num_of_type(0) > 0) THEN npabs_glass(nd)=SUM(surf_type0(1:num_of_type(0))%np_abs, & 1,(surf_type0(1:num_of_type(0))%mat_num<0 .and. & surf_type0(1:num_of_type(0))%iabs/=1)) npabs_coll(nd)=SUM(surf_type0(1:num_of_type(0))%np_abs, & 1,surf_type0(1:num_of_type(0))%iabs==1) npabs_surf(nd)=SUM(surf_type0(1:num_of_type(0))%np_abs, & 1,(surf_type0(1:num_of_type(0))%mat_num>0 .and. & surf_type0(1:num_of_type(0))%iabs/=1)) END IF IF (num_of_type(1) > 0) THEN npabs_glass(nd)=npabs_glass(nd)+SUM(surf_type1(1:num_of_type(1))%np_abs, & 1,(surf_type1%mat_num<0 .and. surf_type1%iabs/=1)) npabs_coll(nd)=npabs_coll(nd)+SUM(surf_type1(1:num_of_type(1))%np_abs, & 1,surf_type1%iabs==1) npabs_surf(nd)=npabs_surf(nd)+SUM(surf_type1(1:num_of_type(1))%np_abs, & 1,(surf_type1%mat_num>0 .and. surf_type1%iabs/=1)) END IF IF (num_of_type(2) > 0) THEN npabs_glass(nd)=npabs_glass(nd)+SUM(surf_type2(1:num_of_type(2))%np_abs, & 1,(surf_type2%mat_num<0 .and. surf_type2%iabs/=1)) npabs_coll(nd)=npabs_coll(nd)+SUM(surf_type2(1:num_of_type(2))%np_abs, & 1,surf_type2%iabs==1) npabs_surf(nd)=npabs_surf(nd)+SUM(surf_type2(1:num_of_type(2))%np_abs, & 1,(surf_type2%mat_num>0 .and. surf_type2%iabs/=1)) END IF npabs_sum(nd)=npabs_glass(nd)+npabs_coll(nd)+npabs_surf(nd) ! ! Total number of photons collected by all surfaces. ! npabs_tot(nd)=npabs_sum(nd)+npabs_caps(nd)+nplost(nd) ! ! Print photon totals. ! IF (iprint>2) THEN CALL print_npabs (nd,9) END IF ! ! Calculate tau_alpha. ! IF (npabs_tot(nd) > 0) THEN tau_alpha(nd) = REAL(npabs_coll(nd))/REAL(npabs_a(nd)) ELSE WRITE (*,*) 'Error -- no absorbed photons: nd, nploop = ', nd, nploop WRITE (9,*) 'Error -- no absorbed photons: nd, nploop = ', nd, nploop CALL adios (2) END IF ! ! Check Convergence ! IF (tau_alpha(nd) > 0.) THEN tol_check = zz*SQRT(ABS((1.0_kindr-tau_alpha(nd))/(tau_alpha(nd)* & REAL(npabs_a(nd))))) ELSE IF (tau_alpha(nd) < 0.0_kindr) THEN WRITE (*,3070) nd, nploop, tau_alpha(nd) WRITE (9,3070) nd, nploop, tau_alpha(nd) CALL adios (2) ELSE tol_check = 1.0_kindr END IF IF (tol_check <= tol) THEN WRITE (*, 2000) nd, nploop, npabs_tot(nd), tol_check, tol nd_converged = nd_converged + 1 go to 600 ELSE WRITE (*, 2001) nd, nploop, npabs_tot(nd), tol_check, tol END IF END DO nd_not_converged = nd_not_converged + 1 600 WRITE (*,*) 'Finished Run for Direction No. - ', nd ! !****************************************************************************** ! ! End of photon convergence loop - Output for direction nd. ! !****************************************************************************** ! IF (iprint>0) CALL print_npabs(nd,9) ! ! Check for lost photons. ! IF (npabs_tot(nd) /= npdo*num_photons) THEN WRITE (*,5000) nd, npabs_tot(nd), np_error(nd), npdo*num_photons WRITE (9,5000) nd, npabs_tot(nd), np_error(nd), npdo*num_photons END IF END DO ! !****************************************************************************** ! ! End of incident directions loop - Finial output. ! !****************************************************************************** ! ! Calculate IAM's ! DO nd=1,num_dirs iam(nd) = tau_alpha(nd)/tau_alpha(norm) END DO ! ! Final photon tally ! CALL tau_alpha_d ! IF (iprint>0) THEN CALL header(9) WRITE (9,2020) DO nd=1,num_dirs WRITE (9,2021) nd,tau_alpha(nd),iam(nd),npabs_coll(nd),npabs_glass(nd) END DO END IF ! ! Write IAM table output ! DO n=0,n_long-1 WRITE (2,1001) (iam(m), m = n*n_tran + 1, (n + 1)*n_tran) END DO WRITE (2,*) ' ' WRITE (2,1002) 'tau_alpha normal = ', tau_alpha(norm) WRITE (2,1002) 'tau_alpha diffuse = ', tad areaa=aperture%area WRITE (2,1002) 'aperture area [m^2] = ', areaa/10000 WRITE (2,1002) 'tolerance = ', tol WRITE (2,1004) 'seed0 = ', seed0 WRITE (2,1004) '# of directions converged on = ', nd_converged ! ! Finished solution phase - end of run ! CALL adios (1) ! ! Trap errors opening files ! 4000 CONTINUE WRITE (*,4001) iu CALL adios (2) ! ! FORMAT statements ! 1000 FORMAT (40a1) 1001 FORMAT (18(f5.3,1x)) 1002 FORMAT (a38, f10.4) 1003 FORMAT (15(f6.2,1x)) 1004 FORMAT (a38, i10) 1005 FORMAT (a38, f7.2) 2000 FORMAT (1x,' Converged - ', 1x, 2(i8,1x), i9, 4x, f10.6, 2x, f9.3) 2001 FORMAT (1x,'Not Converged - ', 1x, 2(i8,1x), i9, 4x, f10.6, 2x, f9.3) 2003 FORMAT (/18x,'S o l u t i o n P h a s e'/ & / 20x,' direction photon tot. no. confidence tolerance' & / 20x,' no. loop of interval' & / 20x,' no. photons') 2004 FORMAT (/21x,'I n p u t P h a s e'/ ) 2010 FORMAT (/21x,'I A M P R O G R A M'/ & / 5x,'Problem heading = ',40a1 & / 5x,'Date of run = ',a8 & / 5x,'Time of Run = ',a8 & / 5x,'IAM program version = ',a8 & / 5x,'Compile date = ',a8 ) 2011 FORMAT (/12x,'P h o t o n C o n t r o l I n f o .'/ & / 5x,'Number of photons per photon loop = ',i9 & / 5x,'Maximum number of photon loops = ',i9 & / 5x,'Convergence tolerance for photons = ',f9.3 & / 5x,'Initial seed for RNG = ',i9 & / 8x,' If 0 entered, generated ''randomly'' from the time') 2014 FORMAT (/16x,'P r i n t i n g C o n t r o l' / & / 5x,'Print control (iprint) = ',i5, & / 8x,' = 0, Only errors are printed' & / 8x,' = 1, Modest printing of results in output file' & / 8x,' = 2, More detailed printing of results in output file' & / 8x,' = 3, Most detailed printing of results in output file' & / 5x,'Data check flag (idata) = ',i5 & / 8x,'/= 0, data check only' & / 8x,' = 0, normal execution') 2020 FORMAT(/'Direction No. tau-alpha __IAM__ npabs_coll npabs_glass') 2021 FORMAT(i8, 9x, f7.5, 7x, f7.5, i9, i9) ! 3070 FORMAT (/2x,'Error in tau_alpha for nd, nploop = ',2i5, & 2x,'tau_alpha = ',f10.5) 4001 FORMAT (1x,'Error opening file - ',i5,' - stopping') 5000 FORMAT (1x,'Warning -- truely lost photons for nd - ', i5 & / 5x,'Number absorbed = ',i10 & / 5x,'Number lost = ',i10 & / 5x,'Number emitted = ',i10 ) ! END program main ! !*********************************************************************** ! SUBROUTINE getdt (k) ! !*********************************************************************** !* Obtain ASCII date and time, and 'random' initial seed * !*********************************************************************** ! USE kinds USE random_funcs USE control USE dates ! IMPLICIT NONE ! CHARACTER(len=5) :: zone CHARACTER(len=10) :: time CHARACTER(len=8) :: date ! INTEGER(kind=kindi) :: iyr, imon, iday, ihr, imin, isec, imil INTEGER :: values(8), k ! CALL DATE_AND_TIME (date, time, zone, values) valuesa(1:8, k) = values IF (k == 0) THEN ! ! Get ASCII date and ASCII Time. ! adate(1:2) = date(5:6) adate(3:3) = '-' adate(4:5) = date(7:8) adate(6:6) = '-' adate(7:8) = date(3:4) ! atime(1:2) = time(1:2) atime(3:3) = ':' atime(4:5) = time(3:4) atime(6:6) = ':' atime(7:8) = time(5:6) ! ! Get integer values. ! iyr = values(1) IF (iyr > 99) THEN iyr = MOD(iyr, 100) END IF imon = values(2) iday = values(3) ihr = values(5) imin = values(6) isec = values(7) imil = values(8) ! ! get initial seed, seed_from_time ! seed_from_time=isec+60*(imin+60*(ihr+24*(iday+28*(imon+iyr)))) IF (MOD(seed_from_time,2) == 0) THEN seed_from_time = IOR(seed_from_time, 1) END IF ! IF (seed_from_time < 0) THEN WRITE (*,*) 'seed_from_time is negative = ',seed_from_time,' correcting' WRITE (9,*) 'seed_from_time is negative = ',seed_from_time,' correcting' seed_from_time = seed_from_time + modulus END IF END IF ! END SUBROUTINE getdt ! !********************************************************************* ! SUBROUTINE files (unit_in, unit_out, nlines) ! !********************************************************************* !* SUBROUTINE to strip comment (& in col. 1) lines from files * !********************************************************************* ! USE kinds ! IMPLICIT NONE ! INTEGER(kind=kindi) :: unit_out, unit_in, ierr, nlines, nchars, n, i INTEGER(kind=kindi), PARAMETER :: nlines_max = 10000 CHARACTER(len=1) :: aline(80) ! DO n = 1, nlines_max READ (unit_in, 1000, END = 100, err = 200) aline IF ( aline(1) /= '&') THEN DO i = 80, 1, -1 IF (aline(i) /= ' ') THEN go to 10 END IF END DO 10 nchars = i + 1 WRITE (unit_out, 1000) aline(1:nchars) END IF END DO ! WRITE (*,*) 'Error from files - too many lines in file no. - ', unit_in WRITE (*,*) ' --- Stopping ---' CALL adios (2) ! 100 nlines = n - 1 REWIND (unit_out) RETURN ! 200 WRITE (*,*) 'Error reading line ',n,' from file no. - ', unit_in WRITE (9,*) 'Error reading line ',n,' from file no. - ', unit_in ierr = 1 CALL adios (2) ! 1000 FORMAT (80a1) ! END SUBROUTINE files ! !********************************************************************* ! SUBROUTINE header (iu) ! !********************************************************************* !* SUBROUTINE to print header to file iu * !********************************************************************* ! USE kinds USE control ! IMPLICIT NONE ! INTEGER :: iu, n CHARACTER(len=8) :: stars8='********' CHARACTER(len=7) :: stars7='*******' ! WRITE (iu, 1000) (stars8, n=1,9), stars7 WRITE (iu, 1001) (head(n), n=1,72) WRITE (iu, 1002) (stars8, n=1,9),stars7 ! 1000 FORMAT (' ',/9a8,a7) 1001 FORMAT (4x, 72a1) 1002 FORMAT (' ',9a8,a7) ! END SUBROUTINE header ! !********************************************************************* ! SUBROUTINE adios (k) ! !********************************************************************* !* SUBROUTINE to terminate execution * !********************************************************************* ! USE kinds USE control USE dates USE photons USE directions ! IMPLICIT NONE ! INTEGER(kind=kindi) :: k, imon, idays, ihrs, imins, isecs INTEGER, DIMENSION(12) :: mondays = & (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/) REAL(kind=kindr) :: seconds, rate ! CALL getdt (2) ! CALL header (9) WRITE (9, 1000) ! IF (k == 1) THEN WRITE (*,*) ' Normal Termination' WRITE (9,*) ' ' WRITE (9,*) ' Normal Termination' ELSE WRITE (*,*) ' Error Termination - see output file/' WRITE (*,*) ' Hit Return to finish' READ (*,*) WRITE (9,*) ' ' WRITE (9,*) ' Error Termination' END IF ! ! Compute time and execution rate ! (for now, based on wall clock time - no seconds function yet) ! idays = valuesa(3,2) - valuesa(3,1) IF (idays < 0) THEN imon = valuesa(2,1) idays = mondays(imon) - valuesa(3,1) + valuesa(3,2) END IF ! IF (idays < 0) THEN WRITE (9, *) 'Error in days -- setting to zero' idays = 0 END IF ihrs = valuesa(5,2)-valuesa(5,1) + idays*24 imins = valuesa(6,2)-valuesa(6,1) isecs = valuesa(7,2)-valuesa(7,1) WRITE (9, 1001) ihrs, imins, isecs ! ! Compute rate, based on wall clock time ! seconds = ((idays*24. + ihrs)*60. + imins)*60. + isecs IF (seconds > 0.) THEN rate = REAL(SUM(npabs_tot))/seconds ELSE rate = 0. END IF WRITE (2, 1002) seconds, SUM(npabs_tot), SUM(np_error), rate WRITE (9, 1002) seconds, SUM(npabs_tot), SUM(np_error), rate ! STOP ! 1000 FORMAT (/5x,'R u n T e r m i n a t i o n I n f o r m a t i o n') 1001 FORMAT (/5x,'Elapsed Wall Clock Times' / & / 8x,'Hours = ', i5 & / 8x,'Minutes = ', i5 & / 8x,'Seconds = ', i5) 1002 FORMAT (/5x,'Execution Rate' / & / 8x,'Total seconds = ', f10.2, ' (seconds)' & / 8x,'Total No. of Photons Emitted = ', i10 & / 8x,'Total No. of Photons Lost = ', i10 & / 8x,'Execution Rate = ', f10.2, ' (photons/sec)') ! END SUBROUTINE adios ! !********************************************************************* ! SUBROUTINE read_directions ! !********************************************************************* !* SUBROUTINE to read and convert input directions * !********************************************************************* ! USE control USE directions USE constants ! IMPLICIT NONE ! INTEGER(kind=kindi) :: n, nl, nt, ierr, nn REAL(kind=kindr) :: long, tran REAL(kind=kindr) :: tan_thetal,tan_thetat INTEGER(kind=kindi), ALLOCATABLE, DIMENSION(:) :: ltest, ttest ! WRITE (*,*) ' - Reading directions' IF (iprint>1) THEN CALL header (9) WRITE (9,1000) END IF ! ierr = 0 ! ! Allocate test arrays ! ALLOCATE (ltest(n_long), ttest(n_tran)) ltest = 0 ttest = 0 ! ! Read Longitudinal Incident Angles ! DO n = 1, n_long READ (1, *, END = 3000, err = 4000) nn, long ! IF (nn < 1 .or. nn > n_long) THEN ierr = -1 WRITE (*,*) 'Error -- reading thetald on line - ',n WRITE (*,*) 'thetald no. -- ',nn,' out of range' WRITE (9,*) 'Error -- reading thetald on line - ',n WRITE (9,*) 'thetald no. -- ',nn,' out of range' ELSE thetald(nn) = long END IF ! IF (ltest(nn) /= 0) THEN ierr = -1 WRITE (*,*) 'Error -- reading thetald on direction - ',n WRITE (*,*) 'thetald no. -- ',nn,' already entered' WRITE (9,*) 'Error -- reading thetald on direction - ',n WRITE (9,*) 'thetald no. -- ',nn,' already entered' ELSE ltest(nn) = 1 END IF END DO ! DO n = 1, n_tran READ (1, *, END = 3000, err = 4000) nn, tran ! IF (nn < 1 .or. nn > n_tran) THEN ierr = -1 WRITE (*,*) 'Error -- reading thetatd on line - ',n WRITE (*,*) 'thetatd no. -- ',nn,' out of range' WRITE (9,*) 'Error -- reading thetatd on line - ',n WRITE (9,*) 'thetatd no. -- ',nn,' out of range' ELSE thetatd(nn) = tran END IF ! IF (ttest(nn) /= 0) THEN ierr = -1 WRITE (*,*) 'Error -- reading thetatd on direction - ',n WRITE (*,*) 'thetatd no. -- ',nn,' already entered' WRITE (9,*) 'Error -- reading thetatd on direction - ',n WRITE (9,*) 'thetatd no. -- ',nn,' already entered' ELSE ttest(nn) = 1 END IF END DO ! ! Test if all directions have been entered ! DO n = 1, n_long IF (ltest(n) == 0) THEN WRITE (*,*) 'Error -- in thetald entered' WRITE (*,*) 'thetald no. -- ',n,' is missing' WRITE (9,*) 'Error -- in thetald entered' WRITE (9,*) 'thetald no. -- ',n,' is missing' ierr = -1 END IF END DO ! DO n = 1, n_tran IF (ttest(n) == 0) THEN WRITE (*,*) 'Error -- in thetatd entered' WRITE (*,*) 'thetatd no. -- ',n,' is missing' WRITE (9,*) 'Error -- in thetatd entered' WRITE (9,*) 'thetatd no. -- ',n,' is missing' ierr = -1 END IF END DO ! IF (ierr == -1) THEN WRITE (*,*) 'Stopping in read_directions - errors detected' WRITE (9,*) 'Stopping in read_directions - errors detected' CALL adios (2) END IF ! ! Process and check input ! norm = 0 n = 0 DO nl = 1, n_long IF (thetald(nl) < 0. .or. thetald(nl) >= 90.) THEN WRITE (*,*) 'Error -- thetald(',nl,') out of range: [0, 90) deg.' WRITE (9,*) 'Error -- thetald(',nl,') out of range: [0, 90) deg.' ierr = -1 END IF ! tan_thetal = TAN(thetald(nl)*dtor) ! DO nt = 1, n_tran n = n + 1 IF (thetatd(nt) <= -90. .or. thetatd(nt) >= 90.) THEN WRITE (*,*) 'Error -- thetatd(',nt,') out of range: (-90, 90) deg.' WRITE (9,*) 'Error -- thetatd(',nt,') out of range: (-90, 90) deg.' ierr = -1 END IF ! IF (ierr == -1) THEN WRITE (*,*) 'Stopping in read_directions - errors detected' WRITE (9,*) 'Stopping in read_directions - errors detected' CALL adios (2) END IF ! ! Conversions ! IF (thetald(nl) == 0 .and. thetatd(nt) == 0) norm = n tan_thetat = TAN(thetatd(nt)*dtor) ! ! Calculate the direction cosines of the emission plane. ! emission(n)%n(3) = -1.0_kindr/SQRT(tan_thetal**2+tan_thetat**2+1.0_kindr) emission(n)%n(1) = emission(n)%n(3)*tan_thetal emission(n)%n(2) = -emission(n)%n(3)*tan_thetat ! ! Calculate the altitude angle and surface azimuth angle for information only. ! thetal(n) = thetald(nl) thetat(n) = thetatd(nt) altituded(n) = ASIN(-emission(n)%n(3))*rtod IF (altituded(n)<90.0_kindr) THEN phid(n)=-emission(n)%n(1)/COS(altituded(n)*dtor) IF (phid(n)>1.0_kindr) phid(n)=1.0_kindr phid(n) = SIGN(1.0_kindr,emission(n)%n(2))*ACOS(phid(n))*rtod ELSE phid(n) = 0.0_kindr END IF END DO END DO ! ! Make sure normal incidence was entered ! IF (norm == 0) THEN WRITE (*,*) 'Stopping in read_directions - normal incidence needed' WRITE (9,*) 'Stopping in read_directions - normal incidence needed' CALL adios (2) END IF ! ! Output ! IF (iprint>1) THEN DO n = 1, num_dirs WRITE (9,1001) n, altituded(n), phid(n), emission(n)%n(1), emission(n)%n(2), & emission(n)%n(3) END DO END IF WRITE (2, 1002) (thetald(n), n = 1, n_long) WRITE (2, 1002) (thetatd(n), n = 1, n_tran) ! ! Deallocate local arrays. ! DEALLOCATE (ltest, ttest) ! IF (ierr == -1) THEN WRITE (*,*) 'Errors detected in read_directions - stopping ' WRITE (9,*) 'Errors detected in read_directions - stopping ' CALL adios (2) END IF ! RETURN ! 3000 WRITE (*,*) 'Error in read_directions -- reached end of file -- Stopping' WRITE (9,*) 'Error in read_directions -- reached end of file -- Stopping' CALL adios (2) RETURN ! 4000 WRITE (*,*) 'Error from read_directions -- ', & ' error reading input on line - ',n,' -- stopping' WRITE (9,*) 'Error from read_directions -- ', & ' error reading input on line - ',n,' -- CALL stopping' CALL adios (2) ! 1000 FORMAT (/22x,'D i r e c t i o n s'/ & / 5x,'No.',9x,'Altitude',8x,'Azimuth',6x,' EX',10x,'EY',10x,'EZ' & / 18x,'(deg.)',10x,'(deg.)'/) 1001 FORMAT (2x,i5,'.',6x,f10.3,5x,f10.3,1x,3(1x,f11.4)) 1002 FORMAT (15(f5.0, 1x)) ! END SUBROUTINE read_directions ! !********************************************************************* ! SUBROUTINE read_materials ! !********************************************************************* !* SUBROUTINE to read materials from input file * !********************************************************************* ! USE kinds USE constants USE control USE materials ! IMPLICIT NONE ! INTEGER(kind=kindi) :: n, ierr, nn, it INTEGER(kind=kindi), ALLOCATABLE, DIMENSION (:) :: itest REAL(kind=kindr) :: ng, kg, rhos, rhoss, rhod, rss, rd, tn, l REAL(kind=kindr) :: th1, th2, ct1, ct2, st1, st2, tau REAL(kind=kindr) :: rho_par, rho_perp, tau_par, tau_perp, alpha_glass ! WRITE (*,*) ' - Reading opaque materials' IF (iprint>1) THEN CALL header (9) WRITE (9,1000) END IF ! ierr = 0 ! ! Initialization of PDF array ! opaque%pdf(1) = 0.0_kindr opaque%pdf(2) = 0.0_kindr opaque%pdf(3) = 0.0_kindr ! ! ALLOCATE itest array ! ALLOCATE (itest(num_mats_opaq)) itest = 0 ! DO n = 1, num_mats_opaq READ (1, *, END = 3000, err = 4000) nn, rhos, rhoss, rhod, rss, rd IF (nn < 1 .or. nn > num_mats_opaq) THEN ierr = -1 WRITE (*,*) 'Error -- reading opaque materials on line - ',n WRITE (*,*) 'Opaque material no. -- ',nn,' out of range' WRITE (9,*) 'Error -- reading opaque materials on line - ',n WRITE (9,*) 'Opaque material no. -- ',nn,' out of range' END IF opaque(nn)%rho_s = rhos opaque(nn)%rho_ss = rhoss opaque(nn)%rho_d = rhod opaque(nn)%r_ss = rss opaque(nn)%r_d = rd ! ! Test if already entered ! IF (itest(nn) /= 0) THEN ierr = 1 WRITE (*,*) 'Error -- reading opaque materials on material - ',n WRITE (*,*) 'Opaque material no. -- ',nn,' already entered' WRITE (9,*) 'Error -- reading opaque materials on material - ',n WRITE (9,*) 'Opaque material no. -- ',nn,' already entered' ELSE itest(nn) = 1 END IF END DO ! ! Test if all surfaces have been entered ! DO n = 1, num_mats_opaq IF (itest(n) == 0) THEN WRITE (*,*) 'Error -- in opaque materials entered' WRITE (*,*) 'Opaque material no. -- ',n,' is missing' WRITE (9,*) 'Error -- in opaque materials entered' WRITE (9,*) 'Opaque material no. -- ',n,' is missing' ierr = -1 END IF END DO ! ! Check Values ! DO n = 1, num_mats_opaq IF (opaque(n)%rho_s < 0. .or. opaque(n)%rho_s > 1.) THEN WRITE (9, 3001) n WRITE (9, 3003) opaque(n)%rho_s ierr = -1 END IF IF (opaque(n)%rho_ss < 0. .or. opaque(n)%rho_ss > 1.) THEN WRITE (9, 3001) n WRITE (9, 3004) opaque(n)%rho_ss ierr = -1 END IF IF (opaque(n)%rho_d < 0. .or. opaque(n)%rho_d > 1.) THEN WRITE (9, 3001) n WRITE (9, 3005) opaque(n)%rho_d ierr = -1 END IF IF (opaque(n)%r_ss < 0. .or. opaque(n)%r_ss > 10.) THEN WRITE (9, 3001) n WRITE (9, 3006) opaque(n)%r_ss ierr = -1 END IF IF (opaque(n)%r_d < 0. .or. opaque(n)%r_d > 10.) THEN WRITE (9, 3001) n WRITE (9, 3007) opaque(n)%r_d ierr = -1 END IF END DO ! ! Output of opaque materials ! DO n = 1, num_mats_opaq opaque(n)%pdf(1) = opaque(n)%rho_ss opaque(n)%pdf(2) = opaque(n)%pdf(1) + opaque(n)%rho_s opaque(n)%pdf(3) = opaque(n)%pdf(2) + opaque(n)%rho_d ! opaque(n)%rp1inv_ss = 1.0_kindr/(1.0_kindr+opaque(n)%r_ss) opaque(n)%rp1inv_d = 1.0_kindr/(1.0_kindr+opaque(n)%r_d) ! IF (iprint>1) THEN WRITE (9, 1001) n, opaque(n)%rho_s, opaque(n)%rho_ss, opaque(n)%rho_d, & opaque(n)%r_ss, opaque(n)%r_d END IF ! ! Test range ! IF (opaque(n)%pdf(3) > 1.) THEN WRITE (*,3001) n WRITE (*,3010) opaque(n)%pdf(3) ierr = -1 END IF END DO ! DEALLOCATE (itest) ! ! Transparent Materials ! IF (num_mats_glass > 0) THEN ALLOCATE (itest(num_mats_glass)) itest = 0 DO n = 1, num_mats_glass READ (1, *, END = 3000, err = 4000) nn, ng, kg, tn IF (nn < 1 .or. nn > num_mats_glass) THEN ierr = -1 WRITE (*,*) 'Error -- reading transparent materials on line - ',n WRITE (*,*) 'Transparent material no. -- ',nn,' out of range' WRITE (9,*) 'Error -- reading transparent materials on line - ',n WRITE (9,*) 'Transparent material no. -- ',nn,' out of range' END IF transparent(nn)%n_glass = ng transparent(nn)%k_glass = kg transparent(nn)%t = tn ! ! Test if already entered ! IF (itest(nn) /= 0) THEN ierr = 1 WRITE (*,*) 'Error -- reading transparent materials on material - ',n WRITE (*,*) 'Transparent material no. -- ',nn,' already entered' WRITE (9,*) 'Error -- reading transparent materials on material - ',n WRITE (9,*) 'Transparent material no. -- ',nn,' already entered' ELSE itest(nn) = 1 END IF END DO ! ! Check Values ! DO n = 1, num_mats_glass IF (transparent(n)%n_glass < 1. .or. transparent(n)%n_glass > 10.) THEN WRITE (9, 3011) n WRITE (9, 3012) transparent(n)%n_glass ierr = -1 END IF IF (transparent(n)%k_glass <= 0. .or. transparent(n)%k_glass > 1.) THEN WRITE (9, 3011) n WRITE (9, 3013) transparent(n)%k_glass ierr = -1 END IF IF (transparent(n)%t < 1.e-6) THEN WRITE (*,*) 'Error -- t(',n,') < 1.e-6' WRITE (9,*) 'Error -- t(',n,') < 1.e-6' ierr = -1 END IF END DO ! ! Output of transparent materials ! IF (iprint>1) WRITE (9, 1010) DO n = 1, num_mats_glass DO it = 0, 89 th1 =it*dtor ct1 = COS(th1) st1 = SIN(th1) st2 = st1/transparent(n)%n_glass th2 = ASIN(st2) ct2 = COS(th2) l = transparent(n)%t/ct2 tau = EXP(-transparent(n)%k_glass*l) rho_par = ((ct2 - transparent(n)%n_glass*ct1)**2)/((ct2 + transparent(n)%n_glass*ct1)**2) rho_perp = ((ct1 - transparent(n)%n_glass*ct2)**2)/((ct1 + transparent(n)%n_glass*ct2)**2) rho_par = MAX(rho_par, 0.0_kindr) rho_perp = MAX(rho_perp, 0.0_kindr) tau_par = ((1.-rho_par)**2)*tau/(1.-tau**2*rho_par**2) tau_perp = ((1.-rho_perp)**2)*tau/(1.-tau**2*rho_perp**2) transparent(n)%tau_glass(it) = 0.5*(tau_par + tau_perp) rho_par = rho_par*(1. + (1.-rho_par)**2*tau**2/(1.-rho_par**2*tau**2)) rho_perp =rho_perp*(1.+(1.-rho_perp)**2*tau**2/(1.-rho_perp**2*tau**2)) transparent(n)%rho_glass(it) = 0.5*(rho_par + rho_perp) END DO transparent(n)%tau_glass(90) = 0.0_kindr transparent(n)%rho_glass(90) = 1.0_kindr IF (iprint>1) WRITE (9, 1011) n, transparent(n)%n_glass, transparent(n)%k_glass ! ! compute increment in transparent properties ! DO it = 0, 89 transparent(n)%drho_glass(it) = transparent(n)%rho_glass(it+1) - transparent(n)%rho_glass(it) transparent(n)%dtau_glass(it) = transparent(n)%tau_glass(it+1) - transparent(n)%tau_glass(it) END DO transparent(n)%drho_glass(90) = 0.0_kindr transparent(n)%dtau_glass(90) = 0.0_kindr ! IF (iprint>1) THEN WRITE (9,1012) DO it = 0, 90 alpha_glass = 1.0_kindr - transparent(n)%tau_glass(it) - transparent(n)%rho_glass(it) WRITE (9,1013) it, transparent(n)%tau_glass(it), transparent(n)%rho_glass(it), alpha_glass END DO END IF END DO ! DEALLOCATE (itest) END IF ! ! Error termination ! IF (ierr == -1) THEN WRITE (*,*) 'Fatal errors detected in read_materials -- stopping' WRITE (9,*) 'Fatal errors detected in read_materials -- stopping' CALL adios (2) END IF ! RETURN ! 3000 WRITE (*,*) 'Error from read_materials --', & ' reached end of file at material sequence - ',n WRITE (9,*) 'Error from read_materials --', & ' reached end of file at material sequence - ',n CALL adios (2) ! 4000 WRITE (*,*) 'Error from read_materials -- ', & ' error reading input on line - ',n WRITE (9,*) 'Error from read_materials-- ', & ' error reading input on line - ',n CALL adios (2) ! 1000 FORMAT (/13x,' O u t p u t o f M a t e r i a l s'/ & / 10x,'O p a q u e M a t e r i a l s :' & / 5x,'No. rho_s rho_ss rho_d r_ss r_d') 1001 FORMAT (2x,i5,'.',3(2x,f8.3),5x,2(2x,f6.3)) ! 1010 FORMAT (/ 10x,'T r a n s p a r e n t M a t e r i a l s :' & / 5x,'No. n_glass k_glass ') 1011 FORMAT (2x,i5,'.',2x,g12.5,2x,g14.6) 1012 FORMAT (/5x,'Transparent Properties' & / 7x,'theta',9x,'tau',12x,'rho',11x,'alpha') 1013 FORMAT (5x,i5,'.',3(3x, f12.5)) ! 3001 FORMAT (1x,'Error from read_materials, surf. mat. no. - ',i5) 3003 FORMAT (1x,' -- specular reflectance out of range; [0,1], rho_s = ',g10.5) 3004 FORMAT (1x,' -- semi-specular reflectance out of range; [0,1], rho_ss = ', & g10.5) 3005 FORMAT (1x,' -- diffuse reflectance out of range; [0,1], rho_d = ',g10.5) 3006 FORMAT (1x,' -- semi-specular power out of range; [0,10], r_ss = ',g10.5) 3007 FORMAT (1x,' -- diffuse power out of range; [0,10], r_d = ',g10.5) 3010 FORMAT (1x,' -- material property sum out of range; [0,1], sum = ',g10.5) 3011 FORMAT (1x,'Error from read_materials, transparent . mat. no. - ',i5) 3012 FORMAT (1x,' -- index of refraction out of range [0,5], n_g = ',g10.5) 3013 FORMAT (1x,' -- extinction coefficient out of range [0,1], k_g = ',g10.5) ! END SUBROUTINE read_materials ! !********************************************************************* ! SUBROUTINE read_collector ! !********************************************************************* !* SUBROUTINE to read module definitions from input file * !********************************************************************* ! USE kinds USE control USE directions USE collector USE photons USE type1_funcs ! IMPLICIT NONE ! CHARACTER(LEN=80) :: line CHARACTER(LEN=80), ALLOCATABLE, DIMENSION(:) :: lines INTEGER(kind=kindi) :: i,j,n,mn,ln ! loop counters INTEGER(kind=kindi) :: ierr INTEGER(kind=kindi) :: nn INTEGER(kind=kindi) :: num_mod_defs,num_elements,surf_type INTEGER(kind=kindi) :: num_modules INTEGER(kind=kindi), ALLOCATABLE, DIMENSION(:) :: num_lines INTEGER(kind=kindi), DIMENSION(0:max_of_type) :: type_count INTEGER(kind=kindi) :: mat_num,int_level,iabs REAL(kind=kindr) :: x0,y0,z0 REAL(kind=kindr) :: xmin,xmax REAL(kind=kindr) :: x1,y1,z1,r REAL(kind=kindr), DIMENSION(3) :: max_bound,min_bound ! ! Open temporary file to store collector definition. ! OPEN (20, file = 'iam.tp2',status = 'unknown',err = 4000) ! ierr = 0 ! ! Read and process aperture plane definition. ! ALLOCATE(lines(4)) READ (1,4001) (lines(j),j=1,4) CALL test_quads(lines,-1,1) DO j=1,4 READ (lines(j),*) nn,x1,y1,z1 aperture%V(j,:)=(/x1,y1,z1/) END DO aperture%n=normal(aperture%V) aperture%o=origin(aperture%V) aperture%d=DOT_PRODUCT(aperture%n,aperture%o) aperture%M1=matrix_M(aperture%V,3,2) aperture%M2=matrix_M(aperture%V,3,4) aperture%W1=matrix_W(aperture%M1) aperture%W2=matrix_W(aperture%M2) aperture%area=0.5_kindr*(area(aperture%M1)+area(aperture%M2)) aperture%np_abs=0 ! ! Read collector definitions to determine array sizes. *********************** ! ! Read module element definitions. ! READ (1,4001) line WRITE (20,4002) line READ (line,*) num_mod_defs ALLOCATE(num_lines(num_mod_defs)) num_lines=0 num_of_type=0 DO mn=1,num_mod_defs type_count=0 READ (1,4001) line WRITE (20,4002) line READ (line,*) num_elements DO n=1,num_elements READ (1,4001) line WRITE (20,4002) line READ (line,*) nn,surf_type SELECT CASE(surf_type) CASE(1) type_count(1)=type_count(1)+1 num_lines(mn)=num_lines(mn)+5 READ (1,4001) (lines(j),j=1,4) CALL test_quads(lines,mn,n) WRITE (20,4002) (lines(j),j=1,4) CASE(2) type_count(2)=type_count(2)+1 num_lines(mn)=num_lines(mn)+2 READ (1,4001) line WRITE (20,4002) line CASE DEFAULT ierr=-1 END SELECT END DO ! ! Read module placement information. ! READ (1,4001) line WRITE (20,4002) line READ (line,*) num_modules IF (mn==1) ALLOCATE(y0_save(num_modules)) DO n=1,num_modules READ (1,4001) line WRITE (20,4002) line END DO num_of_type(1:max_of_type)=num_of_type(1:max_of_type)+type_count(1:max_of_type)*num_modules END DO ! ! Read additional plane surface information. ! READ (1,4001) line WRITE (20,4002) line READ (line,*) num_of_type(0) DO n=1,num_of_type(0) READ (1,4001) line WRITE (20,4002) line READ (1,4001) (lines(j),j=1,4) CALL test_quads(lines,0,n) WRITE (20,4002) (lines(j),j=1,4) END DO DEALLOCATE(lines) ! ! Allocate and initialize arrays for each of the surface types. ! WRITE (*,*) 'Allocating external surfaces array of length - ', num_of_type(0) ALLOCATE(surf_type0(-1:num_of_type(0))) IF (num_of_type(1)>0) THEN WRITE (*,*) 'Allocating bounded planes array of length - ', num_of_type(1) ALLOCATE(surf_type1(num_of_type(1))) END IF IF (num_of_type(2)>0) THEN WRITE (*,*) 'Allocating circular cylinder array of length - ', num_of_type(2) ALLOCATE(surf_type2(num_of_type(2))) END IF ! ! Rewind temporary file. ! REWIND(20) ! ! Read collector definitions and process ************************************** ! WRITE (*,*) ' - Reading collector definitions' IF (iprint>1) THEN CALL header (9) WRITE (9,1000) WRITE (9,1001) num_of_type(0:max_of_type) WRITE (9,1005) WRITE(9,1022) DO j=1,4 WRITE(9,1023) aperture%V(j,:) END DO WRITE (9,1010) END IF ! ! Read and save (locally) the module element definitions. ! num_of_modules=0 READ (20,*) num_mod_defs DO mn=1,num_mod_defs READ (20,*) num_elements IF (num_lines(mn)>0) THEN ALLOCATE(lines(num_lines(mn))) END IF DO n=1,num_lines(mn) READ (20,4001) lines(n) END DO ! ! Read module positions, make copies, and perform coordinate translation. ! type_count=0 READ (20,*) num_modules DO n=1,num_modules num_of_modules=num_of_modules+1 IF (iprint>1) WRITE (9,1011) mn,n READ (20,*) nn,x0,y0,z0 IF (mn==1) y0_save(n)=y0 ln=1 DO i=1,num_elements READ (lines(ln),*) nn,surf_type,mat_num,int_level,iabs ln=ln+1 SELECT CASE(surf_type) ! ! Type 1 surfaces - bounded planes. ! CASE(1) type_count(1)=type_count(1)+1 DO j=1,4 READ (lines(ln),*) nn,x1,y1,z1 ln=ln+1 x1=x1+x0 y1=y1+y0 z1=z1+z0 surf_type1(type_count(1))%V(j,:)=(/x1,y1,z1/) END DO surf_type1(type_count(1))%mat_num=mat_num surf_type1(type_count(1))%int_level=int_level surf_type1(type_count(1))%iabs=iabs surf_type1(type_count(1))%mod_num=num_of_modules IF (iprint>1) THEN WRITE(9,1020) type_count(1) WRITE(9,1021) mat_num,int_level,iabs WRITE(9,1022) DO j=1,4 WRITE(9,1023) surf_type1(type_count(1))%V(j,:) END DO END IF ! ! Type 2 surfaces - circular cylinders. ! CASE(2) type_count(2)=type_count(2)+1 READ (lines(ln),*) x1,y1,z1,r,xmin,xmax ln=ln+1 x1=x1+x0 y1=y1+y0 z1=z1+z0 xmin=xmin+x0 xmax=xmax+x0 surf_type2(type_count(2))%o=(/x1,y1,z1/) surf_type2(type_count(2))%r=r surf_type2(type_count(2))%xmin=xmin surf_type2(type_count(2))%xmax=xmax surf_type2(type_count(2))%mat_num=mat_num surf_type2(type_count(2))%int_level=int_level surf_type2(type_count(2))%iabs=iabs surf_type2(type_count(2))%mod_num=num_of_modules IF (iprint>1) THEN WRITE(9,1030) type_count(2) WRITE(9,1021) mat_num,int_level,iabs WRITE(9,1031) surf_type2(type_count(2))%r WRITE(9,1032) surf_type2(type_count(2))%o END IF CASE DEFAULT ierr=-1 END SELECT END DO IF (ln /= num_lines(mn)+1) ierr=-1 END DO IF (ALLOCATED(lines)) DEALLOCATE(lines) END DO IF (ALLOCATED(num_lines)) DEALLOCATE(num_lines) ! ! Find maximum intersection level. ! max_int_level=MAXVAL(surf_type0%int_level) max_int_level=MAX(max_int_level,MAXVAL(surf_type1%int_level)) max_int_level=MAX(max_int_level,MAXVAL(surf_type2%int_level)) ! ! Read additional plane surface information. ! IF (iprint>1) WRITE (9,1090) READ (20,*) type_count(0) DO n=1,type_count(0) READ (20,*) nn,mat_num,int_level,iabs surf_type0(n)%mat_num=mat_num surf_type0(n)%int_level=int_level surf_type0(n)%iabs=iabs surf_type0(n)%mod_num=0 DO j=1,4 READ (20,*) nn,x1,y1,z1 surf_type0(n)%V(j,:)=(/x1,y1,z1/) END DO IF (iprint>1) THEN WRITE(9,1091) n WRITE(9,1021) mat_num,int_level,iabs WRITE(9,1022) DO j=1,4 WRITE(9,1023) surf_type0(n)%V(j,:) END DO END IF END DO ! ! Close temporary file. ! CLOSE(20) ! ! Preprocess surface information. ! surf_type0%np_abs=0 DO n=1,num_of_type(0) surf_type0(n)%n=normal(surf_type0(n)%V) surf_type0(n)%o=origin(surf_type0(n)%V) surf_type0(n)%d=DOT_PRODUCT(surf_type0(n)%n,surf_type0(n)%o) surf_type0(n)%M1=matrix_M(surf_type0(n)%V,3,2) surf_type0(n)%M2=matrix_M(surf_type0(n)%V,3,4) surf_type0(n)%W1=matrix_W(surf_type0(n)%M1) surf_type0(n)%W2=matrix_W(surf_type0(n)%M2) surf_type0(n)%area=0.5_kindr*(area(surf_type0(n)%M1)+area(surf_type0(n)%M2)) surf_type0(n)%max=MAXVAL(surf_type0(n)%V,1) surf_type0(n)%min=MINVAL(surf_type0(n)%V,1) END DO ! surf_type1%np_abs=0 DO n=1,num_of_type(1) surf_type1(n)%n=normal(surf_type1(n)%V) surf_type1(n)%o=origin(surf_type1(n)%V) surf_type1(n)%d=DOT_PRODUCT(surf_type1(n)%n,surf_type1(n)%o) surf_type1(n)%M1=matrix_M(surf_type1(n)%V,3,2) surf_type1(n)%M2=matrix_M(surf_type1(n)%V,3,4) surf_type1(n)%W1=matrix_W(surf_type1(n)%M1) surf_type1(n)%W2=matrix_W(surf_type1(n)%M2) surf_type1(n)%area=0.5_kindr*(area(surf_type1(n)%M1)+area(surf_type1(n)%M2)) surf_type1(n)%max=MAXVAL(surf_type1(n)%V,1) surf_type1(n)%min=MINVAL(surf_type1(n)%V,1) END DO ! surf_type2%np_abs=0 surf_type2%r_sqrd=surf_type2%r*surf_type2%r DO n=1,num_of_type(2) surf_type2(n)%max(1)=surf_type2(n)%xmax surf_type2(n)%min(1)=surf_type2(n)%xmin surf_type2(n)%max(2)=surf_type2(n)%o(2)+surf_type2(n)%r surf_type2(n)%min(2)=surf_type2(n)%o(2)-surf_type2(n)%r surf_type2(n)%max(3)=surf_type2(n)%o(3)+surf_type2(n)%r surf_type2(n)%min(3)=surf_type2(n)%o(3)-surf_type2(n)%r END DO ! ! Create a bounding box around all collector surfaces (including aperture plane). ! max_bound=MAXVAL(aperture%V,1) min_bound=MINVAL(aperture%V,1) DO n=1,num_of_type(1) max_bound=MAX(max_bound,surf_type1(n)%max) min_bound=MIN(min_bound,surf_type1(n)%min) END DO DO n=1,num_of_type(2) max_bound=MAX(max_bound,surf_type2(n)%max) min_bound=MIN(min_bound,surf_type2(n)%min) END DO surf_type0(-1)%V(1,:)=(/max_bound(1),min_bound(2),min_bound(3)/) surf_type0(-1)%V(2,:)=(/min_bound(1),min_bound(2),min_bound(3)/) surf_type0(-1)%V(3,:)=(/min_bound(1),max_bound(2),min_bound(3)/) surf_type0(-1)%V(4,:)=(/max_bound(1),max_bound(2),min_bound(3)/) surf_type0(0)%V(1,:)=(/max_bound(1),min_bound(2),max_bound(3)/) surf_type0(0)%V(2,:)=(/min_bound(1),min_bound(2),max_bound(3)/) surf_type0(0)%V(3,:)=(/min_bound(1),max_bound(2),max_bound(3)/) surf_type0(0)%V(4,:)=(/max_bound(1),max_bound(2),max_bound(3)/) IF (iprint>1) THEN WRITE (9,2000) WRITE (9,2001) DO j=1,4 WRITE (9,2003) surf_type0(-1)%V(j,:) END DO WRITE (9,2002) DO j=1,4 WRITE (9,2003) surf_type0(0)%V(j,:) END DO END IF ! ! Error termination ! IF (ierr == -1) THEN WRITE (*,*) 'Fatal errors detected in read_collector -- stopping' WRITE (9,*) 'Fatal errors detected in read_collector -- stopping' CALL adios (2) END IF ! RETURN ! 3000 WRITE (*,*) 'Error from read_collector --', & ' reached end of file at material sequence - ',n WRITE (9,*) 'Error from read_collector --', & ' reached end of file at material sequence - ',n CALL adios (2) ! 4000 WRITE (*,*) 'Error from read_collector -- ', & ' error reading input on line - ',n WRITE (9,*) 'Error from read_collector-- ', & ' error reading input on line - ',n CALL adios (2) ! 1000 FORMAT (/13x,' O u t p u t o f C o l l e c t o r'/ & /5x,'Collector Surfaces') 1001 FORMAT (10x,'Number of External Surfaces = ',i3, & /10x,'Number of Bound Planes = ',i3, & /10x,'Number of Circular Cylinders = ',i3) ! 1005 FORMAT (/5x,'Aperture Plane') ! 1010 FORMAT (/5x,'Module Details') 1011 FORMAT (/7x,'Module Type ',i3,' Number ',i3) ! 1020 FORMAT (/10x,'Bound Plane ',i3) 1021 FORMAT (/15x,'Material Number = ',i3, & /15x,'Intersection Level = ',i3, & /15x,'Absorber = ',i3) 1022 FORMAT (15x,'Vertice',9x,'x',9x,'y',9x,'z') 1023 FORMAT (25x,3f10.3) ! 1030 FORMAT (/10x,'Circular Cylinder ',i3) 1031 FORMAT (15x,'Radius = ',f7.3) 1032 FORMAT (15x,'Center',10x,'x',9x,'y',9x,'z',/,25x,3f10.3) ! 1090 FORMAT (/5x,'External Surface Details') 1091 FORMAT (/7x,'External Surface ',i3) ! 2000 FORMAT (/5x,'Bounding Box for the Collector Surfaces') 2001 FORMAT (/10x,'The Lower Bounding Plane is:'/) 2002 FORMAT (/10x,'The Upper Bounding Plane is:'/) 2003 FORMAT (15x,'x = ',f10.3,' y = ',f10.3,' z = ',f10.3) ! 4001 FORMAT (A80) 4002 FORMAT (1x,A80) ! END SUBROUTINE read_collector ! !********************************************************************* ! SUBROUTINE test_quads(lines,mod_num,elem_num) ! !********************************************************************* !* SUBROUTINE to test plane surfaces * !********************************************************************* ! USE kinds USE type1_funcs ! IMPLICIT NONE ! CHARACTER(LEN=80), DIMENSION(4), INTENT(INOUT) :: lines CHARACTER(LEN=80) :: line INTEGER(kind=kindi), INTENT(IN) :: mod_num,elem_num INTEGER(kind=kindi) :: i,nn,n_try ! loop counters REAL(kind=kindr) :: p2 REAL(kind=kindr), DIMENSION(3) :: n REAL(kind=kindr), DIMENSION(4,3) :: V REAL(kind=kindr), DIMENSION(2,3) :: M REAL(kind=kindr), DIMENSION(3,2) :: W REAL(kind=kindr) :: tolerance=1.0e-06 ! DO n_try=1,2 DO i=1,4 READ(lines(i),*) nn,V(i,1),V(i,2),V(i,3) END DO n=normal(V) M(1,:)=V(3,:)-V(1,:) M(2,:)=V(2,:)-V(1,:) ! ! Vertices must be coplainer. ! IF (DOT_PRODUCT((V(4,:)-V(1,:)),n)>tolerance) THEN WRITE(*,1000) mod_num,elem_num WRITE(9,1000) mod_num,elem_num CALL adios (2) END IF ! ! Special case. If V(4,:) is a repeat of any other point ==> triangle. ! IF (V(4,1)==V(1,1) .AND. V(4,2)==V(1,2) .AND. V(4,3)==V(1,3)) THEN lines(4)=lines(2) RETURN END IF IF (V(4,1)==V(2,1) .AND. V(4,2)==V(2,2) .AND. V(4,3)==V(2,3)) THEN lines(4)=lines(2) RETURN END IF IF (V(4,1)==V(3,1) .AND. V(4,2)==V(3,2) .AND. V(4,3)==V(3,3)) THEN lines(4)=lines(2) RETURN END IF ! ! Vertices must be entered in order. If not, fix. ! W=matrix_W(M) p2=DOT_PRODUCT((V(4,:)-V(1,:)),W(:,2)) IF (p2>0.0_kindr) THEN ! vertices are out of order WRITE(*,1001) mod_num,elem_num WRITE(9,1001) mod_num,elem_num line=lines(3) IF (p2<1.0_kindr) THEN ! switch lines 2 & 3 lines(3)=lines(2) lines(2)=line ELSE IF (p2>1.0_kindr) THEN ! switch lines 3 & 4 lines(3)=lines(4) lines(4)=line END IF ELSE RETURN END IF END DO ! ! Atempt to fix order of points failed. ! WRITE(*,1002) mod_num,elem_num WRITE(9,1002) mod_num,elem_num CALL adios (2) ! 1000 FORMAT (1x,'ERROR from read_collector, module definition - ',i5, & ' surface - ',i5, & /1x,' -- coordinates given for vertices do not lie in a plane') 1001 FORMAT (1x,'WARNING from read_collector, module definition - ',i5, & ' surface - ',i5, & /1x,' -- vertices entered in wrong of order, will correct') 1002 FORMAT (1x,'ERROR from read_collector, module definition - ',i5, & ' surface - ',i5, & /1x,' -- vertices entered in wrong of order, could not fix') ! END SUBROUTINE test_quads ! !*********************************************************************** ! SUBROUTINE photon_loop(nd) ! !*********************************************************************** !* Inner photon loop. Generates photon, searches for intersection * !* and determines material interaction. * !*********************************************************************** ! USE kinds USE control USE directions USE collector USE random_funcs USE type1_funcs USE surface_interactions USE photons ! IMPLICIT NONE ! INTEGER(kind=kindi), INTENT(IN) :: nd ! INTEGER(kind=kindi) :: np,ns INTEGER(kind=kindi) :: intersect_level INTEGER(kind=kindi) :: surf_type,surf_num,material,int_level INTEGER(kind=kindi), DIMENSION(3) :: intersect_array INTEGER(kind=kindi) :: interaction_case INTEGER(kind=kindi) :: endcap REAL(kind=kindr), DIMENSION(2) :: random_numbers REAL(kind=kindr) :: min_dist,dist,denominator REAL(kind=kindr), DIMENSION(3) :: point,npoint,surf_normal REAL(kind=kindr), DIMENSION(2) :: p REAL(kind=kindr), DIMENSION(2) :: oc REAL(kind=kindr) :: d,factor,oc_sqrd,d_sqrd,ov REAL(kind=kindr) :: d1,d2,dc1,dc2 REAL(kind=kindr), PARAMETER :: tolerance=1.0e-6 ! ! Photon loop. ! DO np=1,num_photons ! ! Generate a new photon release from the emission plane. ! intersect_level=0 random_numbers=(/ranf(),ranf()/) photon_orig=c_coord(random_numbers,emission(nd)%o,emission(nd)%M1) photon_dir=emission(nd)%n ! ! Check (virtual) intersection with aperture plane. ! denominator=-DOT_PRODUCT(aperture%n,photon_dir) IF (denominator/=0.0_kindr) THEN dist=(DOT_PRODUCT(aperture%n,photon_orig)-aperture%d)/denominator npoint=(photon_orig+photon_dir*dist)-aperture%o p(2)=DOT_PRODUCT(npoint,aperture%W1(:,2)) IF (p(2)>0.0_kindr) THEN IF (p(2)>1.0_kindr) GOTO 50 ! outside of bounds p(1)=DOT_PRODUCT(npoint,aperture%W1(:,1)) IF (p(1)<0.0_kindr) GOTO 50 ! outside of bounds IF (p(1)+p(2)>1.0_kindr) GOTO 50 ! outside of bounds ELSE p(1)=DOT_PRODUCT(npoint,aperture%W2(:,1)) IF (p(1)<0.0_kindr) GOTO 50 ! outside of bounds p(2)=DOT_PRODUCT(npoint,aperture%W2(:,2)) IF (p(1)+p(2)>1.0_kindr) GOTO 50 ! outside of bounds IF (p(2)<0.0_kindr) CYCLE END IF aperture%np_abs=aperture%np_abs+1 END IF 50 CONTINUE ! ! Search for intersection. ! 100 CONTINUE ! upon reflection/transmission, continue searching ! ! Initialize local intersection variables. ! min_dist=HUGE(dist) endcap=0 surf_type=-1 surf_num=-1 intersect_array=(/-1,intersect_level,intersect_level+1/) ! ! Type 0 surfaces ! DO ns=1,num_of_type(0) IF (ANY(intersect_array==surf_type0(ns)%int_level)) THEN denominator=-DOT_PRODUCT(surf_type0(ns)%n,photon_dir) IF (denominator==0.0_kindr) THEN CYCLE ! parallel to surface ELSE dist=(DOT_PRODUCT(surf_type0(ns)%n,photon_orig)-surf_type0(ns)%d)/denominator END IF IF (dist>min_dist) CYCLE ! not closest intersection IF (dist0.0_kindr) THEN p(1)=DOT_PRODUCT(npoint,surf_type0(ns)%W1(:,1)) IF (p(1)<0.0_kindr) CYCLE ! outside of bounds IF (p(1)+p(2)>1.0_kindr) CYCLE ! outside of bounds ELSE p(1)=DOT_PRODUCT(npoint,surf_type0(ns)%W2(:,1)) IF (p(1)<0.0_kindr) CYCLE ! outside of bounds p(2)=DOT_PRODUCT(npoint,surf_type0(ns)%W2(:,2)) IF (p(1)+p(2)>1.0_kindr) CYCLE ! outside of bounds IF (p(2)<0.0_kindr) CYCLE END IF min_dist=dist surf_type=0 surf_num=ns END IF END DO ! ! Type 1 surfaces ! DO ns=1,num_of_type(1) ! ! Fast rejection of surface. ! IF (photon_dir(2) >= 0.0_kindr) THEN IF (photon_orig(2) > surf_type1(ns)%max(2)) CYCLE ELSE IF (photon_orig(2) < surf_type1(ns)%min(2)) CYCLE END IF ! IF (ANY(intersect_array==surf_type1(ns)%int_level)) THEN denominator=-DOT_PRODUCT(surf_type1(ns)%n,photon_dir) IF (denominator==0.0_kindr) THEN CYCLE ! parallel to surface ELSE dist=(DOT_PRODUCT(surf_type1(ns)%n,photon_orig)-surf_type1(ns)%d)/denominator END IF IF (dist>min_dist) CYCLE ! not closest intersection IF (dist0.0_kindr) THEN p(1)=DOT_PRODUCT(npoint,surf_type1(ns)%W1(:,1)) IF (p(1)<0.0_kindr) CYCLE ! outside of bounds IF (p(1)+p(2)>1.0_kindr) CYCLE ! outside of bounds ELSE ! outside of bounds p(1)=DOT_PRODUCT(npoint,surf_type1(ns)%W2(:,1)) IF (p(1)<0.0_kindr) CYCLE ! outside of bounds p(2)=DOT_PRODUCT(npoint,surf_type1(ns)%W2(:,2)) IF (p(1)+p(2)>1.0_kindr) CYCLE ! outside of bounds IF (p(2)<0.0_kindr) CYCLE END IF min_dist=dist surf_type=1 surf_num=ns END IF END DO ! ! Type 2 surfaces ! denominator=SQRT(1.0_kindr-photon_dir(1)**2) DO ns=1,num_of_type(2) ! ! Fast rejection of surface. ! IF (photon_dir(2) >= 0.0_kindr) THEN IF (photon_orig(2) > surf_type2(ns)%max(2)) CYCLE ELSE IF (photon_orig(2) < surf_type2(ns)%min(2)) CYCLE END IF ! IF (ANY(intersect_array==surf_type2(ns)%int_level)) THEN ! ! Special case: Photon direction is parallel to y-z plane. ! IF (denominator==1.0_kindr) THEN IF (photon_orig(1) > surf_type2(ns)%max(1)) CYCLE ! beyond end IF (photon_orig(1) < surf_type2(ns)%min(1)) CYCLE ! beyond end oc=surf_type2(ns)%o(2:3)-photon_orig(2:3) oc_sqrd=DOT_PRODUCT(oc,oc) ov=DOT_PRODUCT(oc,photon_dir(2:3)) d_sqrd=surf_type2(ns)%r_sqrd-(oc_sqrd-ov**2) IF (d_sqrd<0.0_kindr) CYCLE ! no intersection d=SQRT(d_sqrd) IF (oc_sqrd>surf_type2(ns)%r_sqrd+tolerance) THEN ! outside dist=ov-d IF (distmin_dist) CYCLE ! not closest intersection ELSE ! inside dist=ov+d IF (distmin_dist) CYCLE ! not closest intersection END IF endcap=0 ! ! Special case: Photon direction is parallel to x-axis. ! ELSE IF (denominator==0.0_kindr) THEN oc_sqrd=(photon_orig(2)-surf_type2(ns)%o(2))**2 & +(photon_orig(3)-surf_type2(ns)%o(3))**2 IF (oc_sqrd>=surf_type2(ns)%r_sqrd) CYCLE ! no intersection IF (photon_dir(1)<0.0_kindr) THEN dc2=surf_type2(ns)%xmin-photon_orig(1) IF (dc2min_dist) CYCLE ! not closest intersection dist=dc2 endcap=2 ! intersect far end cap ELSE IF (dc1>min_dist) CYCLE ! not closest intersection dist=dc1 endcap=1 ! intersect near end cap END IF ! ! General case: ! ELSE IF (photon_dir(1)<0.0_kindr) THEN dc2=(surf_type2(ns)%xmin-photon_orig(1))/photon_dir(1) IF (dc2min_dist) CYCLE ! not closest intersection dist=d2 ! intersect inside wall endcap=0 ELSE IF (dc2>min_dist) CYCLE ! not closest intersection dist=dc2 ! intersect far end cap endcap=2 END IF ELSE ! before cylinder near end IF (dc1>d2) CYCLE IF (dc1>min_dist) CYCLE ! not closest intersection dist=dc1 ! intersect near end cap endcap=1 END IF ELSE ! outside (infinite) cylinder IF (d_sqrd<0.0_kindr) CYCLE ! no intersection d=SQRT(d_sqrd) d1=factor*(ov-d) IF (d1dc2) CYCLE ! no intersection IF (dc1min_dist) CYCLE ! not closest intersection dist=d1 endcap=0 ELSE ! before cylinder near end d2=factor*(ov+d) IF (d2dc1) THEN IF (d1>min_dist) CYCLE ! not closest intersection dist=d1 ! intersect outside wall endcap=0 ELSE IF (dc1>min_dist) CYCLE ! not closest intersection dist=dc1 ! intersect near end cap endcap=1 END IF END IF END IF END IF min_dist=dist surf_type=2 surf_num=ns END IF END DO ! ! No intersection if front of photon. ! IF (surf_type==-1) THEN nplost(nd)=nplost(nd)+1 IF (intersect_level > 0) THEN np_error(nd)=np_error(nd)+1 WRITE(*,1000) nd,np,intersect_level END IF CYCLE END IF ! ! Absorption on end caps. ! IF (endcap>0) THEN npabs_caps(nd)=npabs_caps(nd)+1 CYCLE END IF ! ! Intersection. Get surface information. ! point=photon_orig+photon_dir*min_dist SELECT CASE(surf_type) CASE (0) surf_normal=surf_type0(surf_num)%n material=surf_type0(surf_num)%mat_num int_level=surf_type0(surf_num)%int_level CASE (1) surf_normal=surf_type1(surf_num)%n material=surf_type1(surf_num)%mat_num int_level=surf_type1(surf_num)%int_level CASE (2) IF (endcap>0) THEN surf_normal(1)=1.0_kindr surf_normal(2)=0.0_kindr surf_normal(3)=0.0_kindr ELSE surf_normal(1)=0.0_kindr surf_normal(2)=(point(2)-surf_type2(surf_num)%o(2))/surf_type2(surf_num)%r surf_normal(3)=(point(3)-surf_type2(surf_num)%o(3))/surf_type2(surf_num)%r END IF material=surf_type2(surf_num)%mat_num int_level=surf_type2(surf_num)%int_level CASE DEFAULT WRITE(9,*) 'Intersection routine did not work' CALL adios (2) END SELECT ! ! Check to see if intersection is on the oposite side. ! IF (DOT_PRODUCT(photon_dir,surf_normal)>0) surf_normal=-surf_normal ! ! Determine material interaction. ! interaction_case=interaction(material,photon_dir,surf_normal) SELECT CASE (interaction_case) CASE (-3:-1) ! reflected photon_orig=point ! new photon origin GO TO 100 ! search for next intersection CASE (0) ! absorbed SELECT CASE (surf_type) CASE (0) surf_type0(surf_num)%np_abs=surf_type0(surf_num)%np_abs+1 CASE (1) surf_type1(surf_num)%np_abs=surf_type1(surf_num)%np_abs+1 CASE (2) surf_type2(surf_num)%np_abs=surf_type2(surf_num)%np_abs+1 CASE DEFAULT CALL adios (2) END SELECT CASE (1) ! transmitted IF(int_level==intersect_level) THEN intersect_level=intersect_level-1 ! transmitted thru outer cover ELSE intersect_level=intersect_level+1 ! transmitted thru inner cover END IF photon_orig=point ! new photon origin GO TO 100 ! search for next intersection CASE (2) ! ERROR CALL adios (2) CASE DEFAULT ! ERROR WRITE(9,*) 'ERROR - interaction is wrong - ',interaction_case CALL adios (2) END SELECT ! ! End of inner photon loop. ! END DO ! 999 FORMAT(' ',3i10,3f10.3) 1000 FORMAT(/'WARNING: Lost photon: direction, ',i2,' photon, ',i5,' level',i2/) END SUBROUTINE photon_loop ! !********************************************************************* ! SUBROUTINE geteplane(nd) ! !********************************************************************* !* SUBROUTINE to define emission plane * !********************************************************************* ! USE kinds USE control USE directions USE collector USE type1_funcs ! IMPLICIT NONE ! INTEGER(kind=kindi), INTENT(IN) :: nd INTEGER(kind=kindi) :: i,j,n_surf,n_surfs REAL(kind=kindr), ALLOCATABLE, DIMENSION(:,:) :: displ REAL(kind=kindr), ALLOCATABLE, DIMENSION(:,:,:) :: proj REAL(kind=kindr), ALLOCATABLE, DIMENSION(:,:,:) :: p REAL(kind=kindr), DIMENSION(2) :: pmin,pmax ! IF (iprint>0) THEN CALL header (9) WRITE (9,1000) nd, phid(nd), altituded(nd),thetal(nd),thetat(nd) END IF ! ! Allocate arrays. ! n_surfs=num_of_type(0) ALLOCATE(displ(4,-1:n_surfs)) ALLOCATE(proj(4,3,-1:n_surfs)) ALLOCATE(p(4,2,-1:n_surfs)) ! ! Calculate d for the emission plane. ! DO n_surf=-1,n_surfs displ(:,n_surf) = MATMUL(surf_type0(n_surf)%V,emission(nd)%n) END DO emission(nd)%d = MINVAL(displ)-0.5_kindr ! ! Project the surfaces into the emission plane. ! DO n_surf=-1,n_surfs displ(:,n_surf) = displ(:,n_surf)-emission(nd)%d END DO DO n_surf=-1,n_surfs proj(:,:,n_surf) = surf_type0(n_surf)%V & -MATMUL(RESHAPE(displ(:,n_surf),(/4,1/)),RESHAPE(emission(nd)%n,(/1,3/))) END DO ! ! Define the barycentric coordinate system (use projection of top of bounding box). ! emission(nd)%o=origin(proj(:,:,0)) emission(nd)%M1=matrix_M(proj(:,:,0),2,4) emission(nd)%W1=matrix_W(emission(nd)%M1) ! ! Now calculate the barycentric coordinates of the vertices of all the surfaces. ! DO n_surf=-1,n_surfs DO i=1,4 p(i,:,n_surf)=b_coord(proj(i,:,n_surf),emission(nd)%o,emission(nd)%W1) END DO END DO ! ! Find the minimum and maximum parameters. ! pmin(1)=MINVAL(p(:,1,:)) pmin(2)=MINVAL(p(:,2,:)) pmax(1)=MAXVAL(p(:,1,:)) pmax(2)=MAXVAL(p(:,2,:)) ! ! Define the vertices of the emission plane. ! emission(nd)%V(1,:)=c_coord(pmin,emission(nd)%o,emission(nd)%M1) emission(nd)%V(2,:)=c_coord((/pmax(1),pmin(2)/),emission(nd)%o,emission(nd)%M1) emission(nd)%V(3,:)=c_coord(pmax,emission(nd)%o,emission(nd)%M1) emission(nd)%V(4,:)=c_coord((/pmin(1),pmax(2)/),emission(nd)%o,emission(nd)%M1) ! ! Redefine the barycentric coordinate system for the emission plane. ! emission(nd)%o=origin(emission(nd)%V) emission(nd)%M1=matrix_M(emission(nd)%V,2,4) emission(nd)%W1=matrix_W(emission(nd)%M1) ! ! Deallocate arrays. ! DEALLOCATE(displ) DEALLOCATE(proj) DEALLOCATE(p) ! ! Print Out Points ! IF (iprint>1) THEN WRITE (9,2000) nd DO i = 1, 4 WRITE (9,2001) i, (emission(nd)%V(i,j), j = 1,3) END DO END IF ! 1000 FORMAT (/5x,'O u t p u t F o r D i r e c t i o n' ,' -', i5,/ & / 10x, 'azimuthal angle, phi (degrees) = ',f10.2 & / 10x, 'altitude angle, altitude (degrees) = ',f10.2 & / 10x, 'longitudinal incident angle, (degrees) = ',f10.2 & / 10x, 'lateral incident angle, altitude (degrees) = ',f10.2) 2000 FORMAT (/5x,'E m i s s i o n P l a n e F o r D i r e c t i o n' ,' -', i5,/ & / 27x,'x',12x,'y',12x,'z') 2001 FORMAT (10x,'Point -', i5, '.', 4(1x,f12.5)) ! ! That's all. ! 100 FORMAT(' ',2f10.5) 110 FORMAT(' ',3f10.5) END SUBROUTINE geteplane ! !*********************************************************************** ! SUBROUTINE print_npabs (nd,iu) ! !*********************************************************************** !* SUBROUTINE to print photon totals * !*********************************************************************** ! ! - modules ! USE kinds USE control USE collector USE photons ! IMPLICIT NONE ! INTEGER(kind=kindi), INTENT(IN) :: nd,iu INTEGER(kind=kindi) :: i,n,mn ! ! ! Print npabs output. ! np_glass=0 np_coll=0 np_surf=0 IF (num_of_type(0) > 0) THEN np_glass(0)=SUM(surf_type0(1:num_of_type(0))%np_abs, & 1,(surf_type0(1:num_of_type(0))%mat_num<0 .and. & surf_type0(1:num_of_type(0))%iabs/=1)) np_coll(0)=SUM(surf_type0(1:num_of_type(0))%np_abs, & 1,surf_type0(1:num_of_type(0))%iabs==1) np_surf(0)=SUM(surf_type0(1:num_of_type(0))%np_abs, & 1,(surf_type0(1:num_of_type(0))%mat_num>0 .and. & surf_type0(1:num_of_type(0))%iabs/=1)) END IF DO mn=1,num_of_modules IF (num_of_type(1) > 0) THEN np_glass(mn)=SUM(surf_type1(1:num_of_type(1))%np_abs, & 1,(surf_type1%mat_num<0 .and. surf_type1%iabs/=1 .and. & surf_type1%mod_num==mn)) np_coll(mn)=SUM(surf_type1(1:num_of_type(1))%np_abs, & 1,(surf_type1%iabs==1 .and. surf_type1%mod_num==mn)) np_surf(mn)=SUM(surf_type1(1:num_of_type(1))%np_abs, & 1,(surf_type1%mat_num>0 .and. surf_type1%iabs/=1 .and. & surf_type1%mod_num==mn)) END IF IF (num_of_type(2) > 0) THEN np_glass(mn)=np_glass(mn)+SUM(surf_type2(1:num_of_type(2))%np_abs, & 1,(surf_type2%mat_num<0 .and. surf_type2%iabs/=1 .and. & surf_type2%mod_num==mn)) np_coll(mn)=np_coll(mn)+SUM(surf_type2(1:num_of_type(2))%np_abs, & 1,(surf_type2%iabs==1 .and. surf_type2%mod_num==mn)) np_surf(mn)=np_surf(mn)+SUM(surf_type2(1:num_of_type(2))%np_abs, & 1,(surf_type2%mat_num>0 .and. surf_type2%iabs/=1 .and. & surf_type2%mod_num==mn)) END IF END DO np_sum=np_glass+np_coll+np_surf WRITE (iu,1000) WRITE (iu,1001) npabs_a(nd) WRITE (iu,1010) WRITE (iu,1011) np_glass(0),np_coll(0),np_surf(0),np_sum(0) DO mn=1,num_of_modules WRITE (iu,1012) mn,np_glass(mn),np_coll(mn),np_surf(mn),np_sum(mn) END DO WRITE (iu,1013) npabs_glass(nd),npabs_coll(nd),npabs_surf(nd),npabs_sum(nd) WRITE (iu,1014) npabs_sum(nd),npabs_caps(nd),nplost(nd),npabs_tot(nd) glass_frac = REAL(npabs_glass(nd))/REAL(npabs_tot(nd)) WRITE (iu,1015) glass_frac,tau_alpha(nd) ! ! Print surface photon totals by intersection level. ! IF (iprint<2) RETURN WRITE(iu,2001) DO i=-1,max_int_level WRITE (iu,2002) i WRITE (iu,2003) DO n=1,num_of_type(0) IF (surf_type0(n)%int_level==i) WRITE (iu, 2004) 0,n,surf_type0(n)%np_abs END DO DO n=1,num_of_type(1) IF (surf_type1(n)%int_level==i) WRITE (iu, 2004) 1,n,surf_type1(n)%np_abs END DO DO n=1,num_of_type(2) IF (surf_type2(n)%int_level==i) WRITE (iu, 2004) 2,n,surf_type2(n)%np_abs END DO END DO ! 1000 FORMAT(/5x,'S u r f a c e A b s o r p t i o n T o t a l s') 1001 FORMAT(/10x,'Aperture Plane -',i10) 1010 FORMAT(/10x,'Collector Surfaces' & /30x,'Glass',2x,'Absorber',5x,'Other',5x'Total') 1011 FORMAT(12x,'External ',4i10) 1012 FORMAT(12x,'Module',i3,1x,4i10) 1013 FORMAT(12x,'Total ',4i10) 1014 FORMAT(/10x,'Collector Surfaces - ',i10 & /10x,'Cylinder End Caps - ',i10 & /10x,'Lost from Volume - ',i10 & //10x,'Total for direction - ',i10) 1015 FORMAT(/10x,'Glass Fraction - ',f10.5 & /10x,'Tau-Alpha - ',f10.5) ! 2001 FORMAT (/5x,'P h o t o n T o t a l s b y S u r f a c e') 2002 FORMAT (/10x,'Surfaces with Intersection Level = ',i5) 2003 FORMAT (12x, 'Surface Type', 3x, ' Number', 3x, & & 'Photons Absorbed') 2004 FORMAT (15x, i5, 8x, i5, 5x, i10) ! END SUBROUTINE print_npabs ! !********************************************************************* ! SUBROUTINE tau_alpha_d ! !********************************************************************* !* SUBROUTINE to calculate tau-alpha diffuse * !********************************************************************* ! USE control USE photons USE directions USE constants ! IMPLICIT NONE ! INTEGER(kind=kindi) :: i, j, k, l, m, q, t, ndelt REAL(kind=kindr) :: tad1, tad2, talp, ddelt, dthet, dthett, dthetl, sintt, & costt, tantt, azm REAL(kind=kindr), ALLOCATABLE, DIMENSION(:,:) :: taul REAL(kind=kindr) :: a,b ! ALLOCATE (taul(n_long+1,n_tran+1)) ! ! ! Make a 2-D array of tau-alpha, with a 90 degree buffer value: ! q = 0 DO i = 1, n_long DO j = 1, n_tran q = q + 1 taul(i,j) = tau_alpha(q) END DO taul(i, n_tran+1) = taul(i, n_tran) END DO DO j = 1, n_tran taul(n_long+1,j) = taul(n_long,j) END DO taul(n_long+1,n_tran+1) = taul(n_long,n_tran) thetatd(n_tran+1) = 90 thetald(n_long+1) = 90 ! tad2 = 0 ndelt = 30 ! ! Integrate tau_alpha over the azimuth (azm) ! DO i = 1, ndelt tad1 = 0 ddelt = pi2/ndelt azm = (i - 1)*ddelt + ddelt/2 ! ! Integrate tau_alpha over the zenith (dthet) ! DO j = 1, ndelt dthet = (j-1)*ddelt + ddelt/2 sintt = SIN(dthet) costt = COS(dthet) tantt = TAN(dthet) ! ! Convert zenith and azimuth to transverse and longitudinal angles ! dthett = ATAN(tantt*SIN(azm))*rtod dthetl = ATAN(tantt*COS(azm))*rtod ! ! Determine the interval dthett & dtehtl fall into ! t = 0 l = 0 DO k = 1, n_tran IF (dthett >= thetatd(k) .and. dthett < thetatd(k+1)) THEN t = k END IF END DO IF (t == 0) THEN t = n_tran END IF DO m = 1, n_long IF (dthetl >= thetald(m) .and. dthetl < thetald(m+1)) THEN l = m END IF END DO IF (l == 0) THEN l = n_long END IF ! ! Use bilinear interpolation to determine tau-alpha(dthetl,dthett) ! a = (dthetl - thetald(l))/(thetald(l+1) - thetald(l)) b = (dthett - thetatd(t))/(thetatd(t+1) - thetatd(t)) talp = (1-a)*(1-b)*taul(l,t) + a*(1-b)*taul(l+1,t) + a*b*taul(l+1,t+1) & + (1-a)*b*taul(l,t+1) tad1 = tad1 + talp*costt*sintt*ddelt END DO tad2 = tad2 + tad1*ddelt END DO tad = 4*tad2/pi ! RETURN ! END SUBROUTINE tau_alpha_d