PROGRAM DDSCAT IMPLICIT NONE c----------------------------------------------------------------------- c To enable MPI use, uncomment INCLUDE statement c comment out INTEGER MPI_COMM_WORLD statement c For non-MPI version, comment out INCLUDE statement c uncomment INTEGER MPI_COMM_WORLD statement c INCLUDE 'mpif.h' INTEGER MPI_COMM_WORLD c----------------------------------------------------------------------- C********************************************************************** C C DDSCAT is a program to use the Discrete Dipole Approximation (DDA) C to calculate the absorption and scattering properties of targets C of arbitrary geometry and dielectric properties. The only C approximation is that the target may be approximated by an array C of point dipoles on a cubic lattice. C DDA theory and validity criteria are presented by B. T. Draine C (1988; Astrophys. J. 333, pp. 848-872). C The original version of DDSCAT was developed by B. T. Draine, C Princeton University Observatory. C C Subsequent versions of DDSCAT were jointly developed by C B. T. Draine (Princeton University Observatory) and C P. J. Flatau (University of California, San Diego, Scripps Inst. C of Oceanography, La Jolla) C C DDSCAT.6.0 includes the following features: C 1. Optional generation of target arrays for to approximate several C standard target shapes, including C a. ellipsoid (sphere is special case) ("ELLIPS") C b. cylindrical prism ("CYLNDR") C c. rectangular prism ("RCTNGL") C d. hexagonal prism ("HEXGON") C e. regular tetrahedron ("TETRAH") C f. cylindrical prism with anisotropic diel. tensor ("UNICYL") C (uniaxial material with crystal axis = cylinder axis) C g. ellipsoid with anisotropic diel. tensor ("ANIELL") C h. two touching ellipsoids with either isotropic or aniostropic C dielectric tensors ("TWOELL" or "TWOAEL") C i. three thouching ellipsoids with either isotropic or aniostropic C dielectric constants ("THRELL" or "THRAEL") C j. two concentric ellipsoids ("CONELL") C k. target consisting of cubic blocks ("BLOCKS") C l. target consisting of union of N spheres ("NSPHER") C m. triangular prism ("PRISM3") C n. layered slab ("LYRSLB") C 2. Alternatively, program may read in target array properties from C ascii file ("FRMFIL") C 3. User may easily specify target orientation with respect to C incident plane wave. C 4. Solution is found iteratively using Complex Conjugate Gradient C methods. C 5. User may select FFT implementation: C a. GPFA code of Temperton (recommended for general use). C b. FFTW code of Frigo and Johnson C c. Convex vector library implementation (if on Convex); C 6. User may select prescription for determining dipole C polarizabilities: C a. Lattice Dispersion Relation ("LATTDR") is only option supported C in DDSCAT.6.0 C 7. Automatic computation of cross sections for: C a. absorption; C b. extinction; C c. scattering; C d. vector radiation pressure; C e. vector torque. C 8. User may easily specify scattering directions for which C scattering matrixs is to be computed. C 9. Automatic computation of orientational averages, with C orientation weights provided by simple subroutine. C 10. Multicomponent targets are allowed. C 11. Anisotropic dielectric tensors are allowed, although with C the restriction that the principal axes must be aligned with C the xyz directions in the "target frame". C 12. Wavelength-dependent dielectric functions are read from user- C provided files, with automatic interpolation if necessary. C 13. Automatic looping over radii and wavelengths if desired. C 14. User may select up to nine Mueller matrix elements for C computation. C 15. DDSCAT.6.0 now supports MPI (Message Passing Interface) for C parallel computations of different scattering orientations. C History: C Note: Overall history of significant changes to the DDSCAT package may C be found in the comments in subroutine version.f . C Here we limit comments to changes to this program module, DDSCAT.f: C 90.11.30 (BTD): Removed LSC3 portion of CXSC complex scratch space. C 90.12.03 (BTD): Reordered vectors. C Removed LSC2 and LSC1 portions of CXSC scratch space. C 90.12.21 (BTD): Added code to create output files C 'qtable' (summary of Q values) C 'mtable' (summary of diel. func. for material 1) C 91.01.03 (BTD): Added MXBETA,MXPHI,MXTHET to argument list for REAPAR C 91.01.05 (BTD): Changed I4 -> I7 in format statements 9010,9011. C Set IOSHP=-1 to suppress printing of target.out file. C 91.05.08 (BTD): Added CALPHA to argument lists for REAPAR and ALPHA C and added CALPHA to various WRITE statements. C 91.05.23 (BTD): Added IWRKSC to arg. list for REAPAR and use it to C control creation of wAArBBkCCC.sca files C Move lines calculating G. C Move code writing to qtable and mtable. C Added IDVOUT to argument list for TARGET C 91.08.14 (BTD): Add QBKSCA to arg. list for GETFML C Modify code to print out Q_bk C 91.08.15 (BTD): Provide different headings for qtable depending on C whether IORTH=1 or 2 (add statement 9043). C 91.09.17 (BTD): Remove calls to ALPHA and EVALA (needed to move these C to subroutine GETFML since alpha from Lattice C Dispersion Relation depends on propagation direction C and polarization state) C Add CALPHA,CXEPS,ICOMP,MXCOMP to argument list for C GETFML (information required by ALPHA) C 91.11.12 (BTD): Added variable IDVSHP to argument list for TARGET C (device number for REASHP to use in reading shape file) C 93.01.15 (BTD): Divided output file qtable into 2 files: C qtable (containing Q_ext,Q_abs,Q_sca,g,Q_bk) and C qtable2 (containing Q_pha, Q_pol, Q_cpol) C 93.01.16 (BTD): Modify so that qtable, qtable2, and mtable are closed C after writing, and reopened for each new write C 93.01.20 (BTD): Add MXNX, MXNY, MXNZ to argument list for subroutine C EXTEND to permit checking of target size against C maximum allowed dimensions C 93.03.11 (BTD): Deleted all code associated with unused variable C CMACHN (originally included to identify machine/OS) C 93.03.12 (BTD): Changed CDESCR*60 -> CDESCR*67 C Moved WRITE(IDVOUT,9010) to subr. TARGET C 93.06.02 (BTD): Corrected error in FORMAT statement 9045 C 93.06.03 (BTD): Corrected error in computation of angle-averaged C backscattering cross section QBKSUM(JO) (had neglected C to reset sum to zero for each new target). C 93.09.28 (BTD): Add "fix" to work around Sun compiler/OS bug (see C description below. C 93.12.15 (BTD): Corrected calculation of PPOL. Added comments on C correspondence between our notation and elements C of 2x2 amplitude scattering matrix and 4x4 Mueller C scattering matrix. C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) to C allow up to 6 shape parameters C 94.06.20 (PJF): Comment about NEWTMP FFT and corresponding changes C in reapar, eself, extend. Add FFT timing code C TSTFFT to distribution and changes to cxfft3 and C timeit routines. C 94.12.20 (BTD): Introduce subroutine VERSION to set value of C variable CSTAMP. version.f will now serve as C location to maintain log of significant code C changes. C 95.06.15 (BTD): Changed CSTAMP from character*11 to character*22 C to include a date. C 95.06.14 (BTD): Create version 5x.1 for internal use. C Version 5x.1 computes vector force and torque C on grain due to incident radiation. C Changed QSCATG(2) to QSCAG(3,2) to allow for C asymmetric scattering by target. C Added variable QTRQSC(1-3,1-2) and added to C argument list of GETFML C Replaced variable GSUM(1-2) to QSCGSUM(1-3,1-2) C Added variables QTRQAB(1-3,1-2) and C QTRQABSUM(1-3,1-2) C 95.06.20 (PJF+: Introduce changes into DDSCAT in order to support C BTD): use of modular iterative solvers (e.g., CCG). C : add CMDSOL to argument list to specify method C of solution C : add CMDTRQ to argument list of REAPAR and GETFML C 95.07.10 (BTD): deleted redundant declarations of NAT,NAT0,NAT3, C NX,NY,NZ (already declared in COMMON/M6/ C 95.07.20 (BTD): added new scratch array SCRRS2(MXNAT) for use in C torque calculations in SCAT C 95.07.27 (BTD): Added new scratch variable CXSCR1 for GETFML to C use for storing first polarization solution while C solving for second. C 96.01.26 (BTD): Replaced IX0,IY0,IZ0 by IXYZ0 C (with simultaneous change in TARGET) C 96.02.23 (BTD): Modified format statements to get additional digits C for a_1, a_2, k vector, and pol. vectors. C 96.11.05 (BTD): Modified to remove various computations to C subroutines GETMUELLER and WRITESCA C 96.11.14 (PJF): Add TIMERS, MXTIMERS, NTIMERS C Add cbinflag, cbinfile, IOBIN (binary file) C Add cbinflag to formal parameters of REAPAR routine C Remove GETSET routine calls C 96.11.15 (PJF): Add cnetflag, cnetfile to reapar call. Add netCDF close C 96.11.21 (BTD): Declared IDNC,NTIMERS,etc. C Removed call to NCCLOS C Added call to WRITESCA in order to call WRITENET in C order to call NCCLOS in order to close netCDF file C 96.12.02 (BTD): initialized smori to zero before looping over C orientations C eliminated variables S1212,S222,CX1112,CX1122, C CX1221,CX1222,CX2122 C 97.07.24 (BTD): Changed argument list for GETMUELLER, replacing C PHI by PHIN, in order to correct error in computation C of Mueller matrix elements and amplitude scattering C matrix elements (see GETMUELLER). C 97.11.01 (BTD): Modified to allow noncubic rectangular lattice. C Information on lattice anisotropy is contained C in variable DX(1-3)=(dx/d,dy/d,dz/d) where C dx,dy,dz=lattice spacings in x,y,z directions C d=(dx*dy*dz)**(1/3) C Note that by definition, DX(1)*DX(2)*DX(3)=1 C Pass DX through COMMON/M1/ for use by subroutines C MATVEC and CMATVEC C 97.12.25 (BTD): Added CXAOFF to COMMON/M7/ C 98.01.13 (BTD): Added DX to argument list in last 2 calls to WRITESCA C 98.05.01 (BTD): Removed DATA/IDVOUT/0 since IDVOUT is in COMMON C Reduced number of continuation lines in CALL WRITESCA C statements C 98.08.09 (BTD): Slight change to FORMAT statement 9020 C 98.12.12 (BTD): Add NSMELTS,SMIND1,SMIND2 to argument lists of REAPAR C and WRITESCA, to allow selection of Muller matrix C elements C 98.12.21 (BTD): changed dimension of CFLPAR from CHARACTER*40 to C CHARACTER*60 to allow longer file names C (also changed in reapar.f and dielec.f) C 99.01.26 (BTD): modify so we can restart calculation for specified C IWAV0,IRAD0,IORI0 C now need to input IWAV0,IRAD0,IORI0 through ddscat.par C 00.06.12 (BTD): added CFLSHP to argument list of REAPAR to support C option NSPHER (need to pass CFLSPH from REAPAR to C TARGET) C 03.02.13 (BTD): added 1 more digit of precision to qtable and qtable2 C outputs; modified header line in stmnts 9044,9046 C 03.04.13 (BTD): added character variable CFLLOG to contain name of C running output file C added call to NAMID to generate a unique name for output C file, to support MPI use. C Unfortunately, when this is done, output to device 0 C is no longer unbuffered. C If it is necessary to reenable unbuffered output for C debugging purposes, C OPEN(UNIT=IDVOUT,FILE=CFLLOG) C should be commented out C Added variables ITNUM(2) and MXITER to argument lists C for GETFML and WRITESCA C 03.07.13 (BTD): Added new variable ITNUMMX(2) to store maximum number C of iterations taken for any of the orientations for C each size and wavelength, for output in orientational C average files waarbbkccc.sca C end history C Copyright (C) 1993,1994,1995,1996,1997,1998,1999,2000,2003 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C Adjustable Parameters: C C MXNX = max. extent of target in x direction C MXNY = y C MXNZ = z C MXCOMP = max. number of different dielectric functions C MXTHET = max. number of target rotation angles THETA C MXBETA = max. number of target rotation angles BETA C MXPHI = max. number of target rotation angles PHI C MXSCA = max. number of scattered directions C MXRAD = max. number of radii C MXWAV = max. number of wavelengths C MXWAVT = max. number of wavelengths in dielectric tables C C*** Set parameters MXNX,MXNY,MXNZ ************************************* INTEGER MXNX,MXNY,MXNZ C PARAMETER(MXNX=8,MXNY=8,MXNZ=8) C PARAMETER(MXNX=12,MXNY=12,MXNZ=12) C PARAMETER(MXNX=16,MXNY=16,MXNZ=16) C PARAMETER(MXNX=12,MXNY=18,MXNZ=18) C PARAMETER(MXNX=16,MXNY=24,MXNZ=24) C PARAMETER(MXNX=20,MXNY=20,MXNZ=20) C PARAMETER(MXNX=24,MXNY=24,MXNZ=24) C PARAMETER(MXNX=25,MXNY=25,MXNZ=25) PARAMETER(MXNX=32,MXNY=32,MXNZ=32) C PARAMETER(MXNX=36,MXNY=36,MXNZ=36) C PARAMETER(MXNX=24,MXNY=48,MXNZ=42) C PARAMETER(MXNX=40,MXNY=40,MXNZ=40) C PARAMETER(MXNX=48,MXNY=16,MXNZ=16) C PARAMETER(MXNX=48,MXNY=48,MXNZ=48) C PARAMETER(MXNX=64,MXNY=64,MXNZ=64) C PARAMETER(MXNX=72,MXNY=48,MXNZ=32) C PARAMETER(MXNX=72,MXNY=72,MXNZ=72) C PARAMETER(MXNX=81,MXNY=54,MXNZ=36) C PARAMETER(MXNX=90,MXNY=60,MXNZ=40) C PARAMETER(MXNX=96,MXNY=48,MXNZ=24) C PARAMETER(MXNX=72,MXNY=72,MXNZ=72) C PARAMETER(MXNX=96,MXNY=32,MXNZ=32) C Experiment for maximum possible size on Intel-32 Linux C PARAMETER(MXNX=135,MXNY=135,MXNZ=135) C PARAMETER(MXNX=144,MXNY=144,MXNZ=144) C 150x150x150 is too large for 32-bit Redhat Linux and pgf77 C PARAMETER(MXNX=150,MXNY=150,MXNZ=150) C C C*** Set parameters MXCOMP,MXTHET,MXBETA,MXPHI,MXRAD,MXWAV,MXWAVT,MXSCA INTEGER MXCOMP PARAMETER(MXCOMP=9) INTEGER MXTHET,MXBETA,MXPHI PARAMETER(MXTHET=100,MXBETA=100,MXPHI=100) INTEGER MXRAD,MXWAV,MXWAVT PARAMETER(MXRAD=100,MXWAV=100,MXWAVT=1500) INTEGER MXSCA PARAMETER(MXSCA=1000) INTEGER MXTIMERS,NTIMERS PARAMETER(MXTIMERS=20,NTIMERS=12) C C*** Derived Parameters: ********************************************** INTEGER MXNAT,MXN3 PARAMETER(MXNAT=MXNX*MXNY*MXNZ,MXN3=3*MXNAT) INTEGER MXCXSC C MXCXSC specifies required complex scratch array CXSC C increase from 7*MXN3 required by PETR C to 10*MXN3 required by PIMCBICGSTAB method C (PIM Complex BiConjugate Gradient with Stabilization) PARAMETER(MXCXSC=10*MXN3) C C Common variables (required for communication with the C conjugate gradient package): C CHARACTER CMDFFT*6 INTEGER IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT3,NAT0,NX,NY,NZ INTEGER*2 IOCC(MXNAT) REAL AKR(3),DX(3) COMPLEX & CXADIA(MXNAT,3,3), & CXAOFF(MXNAT,3), & CXZC(MXNX+1,MXNY+1,MXNZ+1,6), & CXZW(2*MXNX,2*MXNY,2*MXNZ,3) C----------------------------------------------------------------------- COMMON/M1/AKR,DX COMMON/M2/CXADIA COMMON/M3/CXZC COMMON/M4/CXZW COMMON/M5/IOCC COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3, & NAT0,NX,NY,NZ,MXN3F,IDVOUT,CMDFFT COMMON/M7/CXAOFF C----------------------------------------------------------------------- C C Local variables: C CHARACTER CMDSOL*6,CMDTRQ*6,CMSGNM*70 C COMPLEX CXEN C REAL & AEFF,AK1,AK3,BETAD,BETAMI,BETAMX,BETMID,BETMXD, & DEGRAD,FREQ, & PHID,PHIMAX,PHIMID,PHIMIN,PHIMXD,PI,PIA2, & THETAD,THETMI,THETMX,THTMID,THTMXD,TOL, & WAVE,XX C C Target Properties: C AEFF=aeff=effective radius of target (physical units) C PIA2=pi*(aeff/d)**2=pi*(3*NAT/4*pi)**(2/3) C C Incident Wave Properties: C AK1=(k*d), where k=2*pi/wave=propagation vector and d=dipole spacing C AK3=(k*d)**3 C AKE2=(k*d \dot E0), where E0 is incident polarization vector C WAVE=wavelength (in vacuo) in physical units C XX=k*aeff=2*pi*aeff/wave C C Scattering Properties: C G= for scattered radiation C QSCGSUM(1-3,1-2)=running (weighted) sum of g(1-3)*Q_sca C over target orientations, for incident polarizations 1-2, C where C g(1)= for scattered radiation C g(2)= for scattered radiation C g(3)= for scattered radiation C theta is measured relative to incident direction C phi is measured relative to x,y plane (Lab Frame) C QPHA(1-2)=(phaseshift cross section)/(pi*aeff**2) C QEXSUM(1-2)=running sum of QEXT over target orientations C QABSUM(1-2)=running sum of QABS over target orientations C QBKSUM(1-2)=running sum of QBKSCA over target orientations C QSCSUM(1-2)=running sum of QSCA over target orientations C QTRQSCSUM(1-3,1-2)=running sum of QTRQSC(1-3,1-2) over target orientations C THETND=scattering angle theta in degrees C PHIND=scattering angle phi in degrees C SINPHI=sin(phi) C DEPOLR=depolarization ratio for scattered of initially polarized C light into direction theta,phi C PPOL=percentage polarization for scattering of initially unpolarized C light into direction theta,phi C INTEGER & IBETA,IBETA1,IBETH,IBETH1,ICTHM,IDIR,IDNC,IDVERR,IDVSHP,IERR, & INIT,ILIN10,ILIN12, & IOBIN,IOPAR,IORI,IORI0,IORI1,IORTH,IOSHP,IPHI,IPHI1,IPHIM,IRAD, & IRAD0,IRAD1,ITASK,ITHETA,ITHETA1,IWAV,IWAV0,IWAV1, & IWRKSC, & J,JO,LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0,MYID, & NBETA,NBETH,NCOMP,NORI,NPHI,NRAD,NSCAT,NSMELTS,NTHETA, & NUMPROCS,NWAV INTEGER MXITER INTEGER ITNUM(2),ITNUMMX(2) INTEGER SMIND1(9),SMIND2(9) C C MYID = MPI process identifier that runs from 0-(NUMPROCS-1). 0 is the C master process. C NUMPROCS = number of MPI parallel processes (including the master) C IERR = MPI error code, which we ignore C NAT=number of dipoles in target C NAT3=3*NAT C NBETA=number of different beta values for target orientation C NBETH=NBETA*NTHETA=number of target orientations for outer orientation loop C NTHETA=number of different theta values for target orientation C NPHI=number of different phi values for target orientation C NORI=NBETA*NTHETA*NPHI=total number of target orientations C NSCAT=number of different scattering directions for each orientation C NWAV=number of different wavelengths C NRAD=number of different target radii C NSMELTS=number of scattering matrix elements to print out C INIT=0,1,2 for choice of initial vector |x0> in complex conjugate gradient C method C IORTH=1 or 2 to do 1 or 2 incident polarizations C IWRKSC=0 or 1 to suppress or generate "wAArBBkCCC.sca" files for each C target orientation C NX=x-length of extended target C NY=y-length of extended target C NZ=z-length of extended target C C SMIND1(1-NSMELTS) = index 1 of scattering matrix elements to C be printed out (e.g., 1 for element S_{13}) C SMIND2(1-NSMELTS) = index 2 of scattering matrix elements to C be printed out (e.g., 3 for element S_{13}) CHARACTER & CALPHA*6,CDIEL*6,CSHAPE*6,CSTAMP*22, & CFLAVG*13,CFLPAR*13,CFLSCA*14,CFLSHP*13,CDESCR*67, & CBINFLAG*6, CBINFILE*14, CNETFLAG*6, CNETFILE*14 CHARACTER CFLLOG*14 CHARACTER*60 CFLEPS(MXCOMP) C .. C .. Local Arrays .. COMPLEX & CX1121(MXSCA), & CXALOF(MXN3),CXALOS(MXN3), & CXALPH(MXN3),CXALDS(MXN3), & CXE(MXN3),CXE01(3),CXE02(3), & CXE01R(3),CXE02R(3),CXEPS(MXCOMP), & CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),CXF22(MXSCA), & CXRLOC(MXCOMP+1,3,3), & CXRFR(MXCOMP),CXS1(MXSCA),CXS2(MXSCA),CXS3(MXSCA),CXS4(MXSCA), & CXSC(MXCXSC),CXSCR1(MXN3),CXXI(MXN3) C C Arrays describing target properties: C CXADIA(1-3*NAT)=diagonal elements of A matrix C CXEPS(1-MXCOMP)=dielectric const. for each of MXCOMP materials C CXREFR(1-3)=complex refractive index at current wavelength C CXALPH(1-3*NAT)=diagonal elements of dipole polarizability tensor alpha C CXALOF(1-3*NAT)=off-diagonal elements of dipole polarizability tensor alpha C CXALDS(1-3*NAT)=diagonal elements of dipole polarizability tensor alpha C calculated (for CALPHA.EQ.'SCLDR') by Rahmani et al. (2002) C method, and stored to save computation time and memory C CXALOS(1-3*NAT)=off-diagonal elements of dipole polarizability tensor C corresponding to CXALDS C C Arrays describing incident radiation: C CXE01(1-3)=incident polarization vector e01 prior to target rotation C CXE02(1-3)=incident polarization vector e02 prior to target rotation C CXE01R(1-3)=rotated incident polarization vector e01 in Target Frame C CXE02R(1-3)=rotated incident polarization vector e02 in Target Frame C CXE0R(1-3)=incident polarization vector (=CXE01R or CXE02R) being used C CXE(1-3*NAT)=incident electric field vector in Target Frame C C Array describing polarization of target for one incident wave: C CXXI(1-3*NAT)=polarization vector at each dipole C C Arrays describing scattering properties of target: C CXFF11(1-NSCAT)=scattering matrix element f11 for each of NSCAT directions C CXFF12(1-NSCAT)= f12 C CXFF21(1-NSCAT)= f21 C CXFF22(1-NSCAT)= f22 C CX1121(1-NSCAT)=sum of f11* \times f21 over orientations C C Scratch arrays: C CXSC(1-MXCXSC)=scratch array used by DDACCG C INTEGER*2 ICOMP(MXN3),ISCR1(MXN3),IXYZ0(MXNAT,3) C C IXYZ0(1-NAT0,1-3)=x/d,y/d,z/d for dipoles 1-NAT0 in real target C C Scratch arrays: C ISCR1,ISCR2,ISCR3,ISCR4 (used by EXTEND) C C REAL & A1(3),A2(3),AEFFA(MXRAD),AKSR(3,MXSCA), & BETA(MXBETA), & E1A(MXWAVT),E2A(MXWAVT),EM1(3,MXSCA),EM1R(3,MXSCA), & EM2(3,MXSCA),EM2R(3,MXSCA),EN0R(3),ENSC(3,MXSCA), & PHI(MXPHI),PHIN(MXSCA), & QABS(2),QABSUM(2),QBKSCA(2),QBKSUM(2),QEXT(2), & QEXSUM(2),QPHA(2),QPHSUM(2),QSCAT(2),QSCAG(3,2),QSCGSUM(3,2), & QSCSUM(2),QTRQAB(3,2),QTRQSC(3,2), & QTRQABSUM(3,2),QTRQSCSUM(3,2), & S1111(MXSCA),S2121(MXSCA), & SCRRS1(MXN3),SCRRS2(MXNAT),SHPAR(6),SM(MXSCA,4,4), & SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA), & WAVEA(MXWAV),WGTA(MXTHET,MXPHI),WGTB(MXBETA),WVA(MXWAVT), & TIMERS(MXTIMERS) C C Arrays describing target properties: C A1(1-3)=initial direction of target axis a1 in Lab Frame C A2(1-3)=initial direction of target axis a2 in Lab Frame C AEFFA(1-NRAD)=values of effective radius of target (physical units) C SHPAR(1-6)=up to 6 parameters for defining target geometry C C Arrays describing target rotation: C THETA(1-NTHETA)=desired values of target rotation angle theta C BETA(1-NBETA)=desired values of target rotation angle beta C PHI(1-NPHI)=desired values of target rotation angle phi C WGTA(1-NTHETA,1-NPHI)=weight factor for each orientation of target C axis a1 in Lab Frame C WGTB(1-NBETA)=weight factor for rotation of target around a1 C Arrays describing incident wave: C WAVEA(1-NWAV)=values of wavelength (physical units) C EN0R(3)=unit vector in direction of incidence in (rotated) Target C Frame C AKR(3)=k*d vector in (rotated) Target Frame C C Arrays describing scattering properties: C ENSC(1-3,1-NSCAT)=scattering vectors in Lab Frame C AKSR(1-3,1-NSCAT)=scattering k-vectors in Target Frame C EM1(1-3,1-NSCAT)=scattering polarization state 1 in Lab Frame C EM2(1-3,1-NSCAT)=scattering polarization state 2 in Lab Frame C EM1R(1-3,1-NSCAT)=scattering polarization state 1 in Target Frame C EM1R(1-3,1-NSCAT)=scattering polarization state 2 in Target Frame C THETAN(1-NSCAT)=values of theta for NSCAT scattering directions nhat C at which scattering matrix fml is to be calculated C PHIN(1-NSCAT)=values of phi for scattering NSCAT directions nhat at C which scattering matrix fml is to be calculated C S1111(1-NSCAT)=sum of |f11|**2 over target orientations C S2121(1-NSCAT)= |f21|**2 C QABS(1-2)=(absorption cross section)/(pi*aeff**2) for 2 polarizations C QBKSCA(1-2)=4*pi*(diff.scatt.cross section for theta=pi)/(pi*aeff**2) C for 2 polarizations C QEXT(1-2)=(extinction cross section)/(pi*aeff**2) C QSCA(1-2)=(scattering cross section)/(pi*aeff**2)=QEXT-QABS C QSCAT(1-2)=(scattering cross section)/(pi*aeff**2) estimated by SCAT C C Scratch array: C SCRRS1(1-3*NAT)=scratch array used by DDACCG and SCAT C SCRRS2(1-NAT)=scratch array used by SCAT C C External subroutines: EXTERNAL DIELEC,EVALA,NAMER,REAPAR C C Intrinsic functions: INTRINSIC ATAN,SQRT,REAL C C VMS version: C C DATA CFLPAR,CFLSHP/'DDSCAT.PAR','SHAPE.DAT'/ C C Unix version: C DATA CFLPAR,CFLSHP/'ddscat.par','shape.dat'/ C C Set device numbers: C IDVERR = device number for error messages C IOPAR = device number for reading ddscat.par file C IOSHP = device number for writing target.out (IOSHP<0 to suppress) C IDVSHP = device number for reading shape file C DATA IDVERR/0/ DATA IDVSHP/10/ DATA IOPAR/1/ DATA IOSHP/-1/ C C IDVOUT = device number for running I/O C on Solaris IDVOUT=0 produces unbuffered output to standard C output C IDVOUT=6 produces buffered output to standard C output C Initialize IDVOUT (cannot be done in DATA statement because IDVOUT C is in COMMON/M6/ C IDVOUT=0 C C*********************************************************************** C MPI environment setup: CALL MPI_INIT(IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NUMPROCS, IERR) C C MPI_INIT is called by each process to initialize the MPI environment C MPI_COMM_RANK is called by each process to determine its identifier C which is stored in MYID C MPI_COMM_SIZE is called by each process to determine the number of C parallel processes, stored in NUMPROCS C If not using MPI, these calls can be turned into dummies by including C subroutines in mpi_fake.f (also need to comment out INCLUDE 'mpif.h' C and uncomment declaration of integer variable MPI_COMM_WORLD as noted C above) C Call subroutine NAMID to obtain unique filename for output file C written by this process. CALL NAMID(MYID,CFLLOG) C 2003.04.12 (BTD) When following statement is enabled, output is C buffered. If it is necessary to have unbuffered output for debugging C code, comment out the following OPEN statement: OPEN(UNIT=IDVOUT,FILE=CFLLOG) CALL VERSION(CSTAMP) IF(MYID.EQ.0)THEN WRITE(IDVOUT,FMT=9000)CSTAMP,MXNX,MXNY,MXNZ ENDIF C C Initialize various variables in common blocks which are derived C from parameters: C MXNATF=MXNAT MXNXF=MXNX MXNYF=MXNY MXNZF=MXNZ MXN3F=MXN3 C .. C Pseudo names of scratch array positions: LACE=1 LAXI=1*MXN3+1 LCLM=2*MXN3+1 LGI=3*MXN3+1 LPI=4*MXN3+1 LQI=5*MXN3+1 LSC0=6*MXN3+1 C PI=4.*ATAN(1.) DEGRAD=180./PI C C C**** Read program control parameters C CMDFFT=choice of solution methods C CSHAPE=description of shape C SHPAR(1-6)= up to 6 parameters for defining target geometry C INIT=selection of |x0> in CCG method C TOL=error tolerance in CCG method C ICTHM=number of theta values for computation of CSCA and G by SCAT C IPHIM=number of phi values for computation of CSCA and G by SCAT C MXWAV=dimensioning information for array WAVEA C NWAV=number of wavelengths to be used C WAVEA(1-NWAV)=wavelengths (physical units) C MXRAD=dimensioning information for array AEFFA C NRAD=number of radii to be used C AEFFA(1-NRAD)=effective radii (physical units) C IORTH=1 or 2 for 1 or 2 incident polarization states to be used C NBETA=number of beta values to be used for target orientation C BETAMI=minimum value of beta (radians) C BETAMX=maximum value of beta (radians) C NTHETA=number of theta values to be used for target orientation C THETMI=minimum value of theta (radians) C THETMX=maximum value of theta (radians) C NPHI=number of phi values to be used for target orientation C PHIMIN=minimum value of phi (radians) C PHIMAX=maximum value of phi (radians) C DEGRAD=degrees/radian (conversion factor) C MXSCA=dimensioning information for arrays PHIN, THETAN C NSCAT=number of scattered directions for computation of fml C THETAN(1-NSCAT)=theta values (radians) for scattered directions C PHIN(1-NSCAT)=phi values (radians) for scattered directions C IF(MYID.EQ.0)THEN CALL REAPAR(CBINFLAG,CNETFLAG,AEFFA, & BETAMI,BETAMX, & CALPHA,CDIEL,CFLEPS,CFLPAR,CFLSHP,CMDSOL,CMDFFT, & CSHAPE,CMDTRQ,CXE01,CXE02, & DEGRAD, & ICTHM,IDVERR,IDVOUT,INIT,IOPAR,IORTH,IPHIM,IWRKSC, & MXBETA,MXCOMP,MXPHI,MXRAD,MXSCA,MXTHET,MXWAV, & NBETA,NCOMP,NPHI,NRAD,NSCAT,NTHETA,NWAV, & IWAV0,IRAD0,IORI0,NSMELTS, & PHIMIN,PHIMAX,PHIN, & THETAN,THETMI,THETMX,TOL, & WAVEA,SHPAR,SMIND1,SMIND2,DX) C NORI=NTHETA*NBETA*NPHI NBETH=NTHETA*NBETA C Compute angles in degrees for printout BETMID=DEGRAD*BETAMI BETMXD=DEGRAD*BETAMX THTMID=DEGRAD*THETMI THTMXD=DEGRAD*THETMX PHIMID=DEGRAD*PHIMIN PHIMXD=DEGRAD*PHIMAX C C*** Obtain scattering vectors and scattering pol. vectors in Lab Frame C CALL SCAVEC(MXSCA,NSCAT,THETAN,PHIN,CXE01,ENSC,EM1,EM2) C C*** Specify target properties C C A1(1-3)=target axis 1 in Target Frame C A2(1-3)=target axis 2 in Target Frame C CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT0,IXYZ0,ICOMP,IDVOUT,MXN3,NAT3) C C NAT0=number of dipoles in "real" target C*** Extend target to rectangular volume suitable for FFT C CALL EXTEND(CMDFFT,ICOMP,ISCR1,IDVERR,IDVOUT,IOCC, & IXYZ0(1,1),IXYZ0(1,2),IXYZ0(1,3), & MXNX,MXNY,MXNZ,MXNAT,MXN3,NAT,NAT0,NAT3,NX,NY,NZ) C C Note: IXYZ0(1-NAT0,1-3) contain coordinates of occupied lattice sites. C WRITE(IDVOUT,FMT=9011)NAT,NX,NY,NZ C C**** Specify target orientations C C MXBETA=dimensioning information for array BETA of beta values C MXTHET=dimensioning information for array THETA of theta values C MXPHI=dimensioning information for array PHI of phi values C BETAMI,BETAMX=minimum,maximum beta values (radians) C THETMI,THETMX=minimum,maximum theta values (radians) C PHIMIN,PHIMAX=minimum,maximum phi values (radians) C BETA(1-NBETA)=beta values for target orientations (radians) C THETA(1-NTHETA)=theta values for target orientations (radians) C PHI(1-NPHI)=phi values for target orientations (radians) C WGTA(1-NTHETA,1-NPHI)=weight for each orientation of axis a1 in LF C WGTB(1-NBETA)=weight for each rotation of target around a1 C CALL ORIENT(BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX, & MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI, & BETA,THETA,PHI,WGTA,WGTB) C C*** Open file 'qtable' for summary of Q values and C 'mtable' for summary of refractive indices C OPEN(UNIT=10,FILE='qtable',STATUS='UNKNOWN') OPEN(UNIT=11,FILE='mtable',STATUS='UNKNOWN') OPEN(UNIT=12,FILE='qtable2',STATUS='UNKNOWN') C why not use cbinfile ? and open it in writesca ? IF(CBINFLAG.EQ.'ALLBIN'.OR.CBINFLAG.EQ.'ORIBIN')THEN IOBIN=51 OPEN(UNIT=IOBIN,file='dd.bin',form='unformatted') ENDIF C C*** Write out header for 'qtable', 'qtable2', and 'mtable' files WRITE(10,FMT=9020)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0,CDIEL ILIN10=7 WRITE(10,FMT=9021)(CFLEPS(J),J=1,NCOMP) ILIN10=ILIN10+NCOMP WRITE(12,FMT=9020)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0,CDIEL ILIN12=7 WRITE(12,FMT=9021)(CFLEPS(J),J=1,NCOMP) ILIN12=ILIN12+NCOMP IF(IORTH.EQ.1)THEN WRITE(10,FMT=9043)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN10=ILIN10+6 WRITE(12,FMT=9045)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN12=ILIN12+6 ELSE WRITE(10,FMT=9044)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN10=ILIN10+6 WRITE(12,FMT=9046)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN12=ILIN12+6 ENDIF CLOSE(10) CLOSE(12) WRITE(11,FMT=9020)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0,CDIEL WRITE(11,FMT=9021)(CFLEPS(J),J=1,NCOMP) WRITE(11,FMT=9048) ENDIF C C*** The master process needs to share information from the above section C with the slave processes. All processes call this subroutine, which C does nothing if MPI is not being used. C CALL SHARE1(A1,A2,AEFFA,BETA, & BETMID,BETMXD, & CALPHA,CBINFLAG,CDESCR,CDIEL, & CFLEPS,CMDFFT,CMDSOL,CMDTRQ, & CNETFLAG,CSHAPE,CXE01,CXE02, & DX,ENSC,EM1,EM2, & ICOMP,ICTHM,IDVERR,IDVOUT,ILIN10, & ILIN12,INIT,IOBIN,IOCC,IORI0,IORTH, & IPHIM,IRAD0,IWAV0,IWRKSC,IXYZ0, & MXBETA,MXCOMP,MXN3,MXNAT, & MXPHI,MXRAD,MXSCA,MXTHET,MXWAV, & NAT,NAT0,NAT3,NBETA,NBETH,NCOMP,NORI,NPHI, & NRAD,NSCAT,NSMELTS,NTHETA,NWAV,NX,NY,NZ, & PHI,PHIMID,PHIMXD,PHIN, & SHPAR,SMIND1,SMIND2,THETA,THETAN, & THTMID,THTMXD,TOL, & WAVEA,WGTA,WGTB) C C*** Loop over wavelengths: C IWAV1=IWAV0+1 IRAD1=IRAD0+1 IORI1=IORI0+1 ITHETA1=1+IORI0/(NBETA*NPHI) IBETA1=1+IORI0/NPHI-(ITHETA1-1)*NBETA IBETH1=1+IORI0/NPHI IPHI1=IORI1-(ITHETA1-1)*NBETA*NPHI-(IBETA1-1)*NPHI DO IWAV=IWAV1,NWAV WAVE=WAVEA(IWAV) C C*** Determine dielectric properties of target C IF(MYID.EQ.0)THEN CALL DIELEC(CDIEL,WAVE,IDVOUT,CFLEPS,CXEPS, & MXCOMP,MXWAVT,NCOMP,E1A,E2A,WVA) ENDIF C C*** Again the master process needs to share information with the slaves. C CALL SHARE2(MXCOMP,MXWAVT, & CXEPS,E1A,E2A,WVA) C C*** Obtain complex polarizabilities C DO J=1,NCOMP CXRFR(J)=SQRT(CXEPS(J)) IF(MYID.EQ.0)THEN WRITE(IDVOUT,FMT=9031)CXRFR(J),CXEPS(J),J ENDIF ENDDO C C**** Loop over grain sizes: C DO IRAD=IRAD1,NRAD AEFF=AEFFA(IRAD) XX=2.*PI*AEFF/WAVE IF(MYID.EQ.0)THEN WRITE(IDVOUT,FMT=9032)AEFF,WAVE,XX ENDIF C C Compute AK1=length of k vector in "natural" units=k*d C (natural unit of length = d = lattice spacing) C Remember that NAT = number of dipoles including "vacuum" sites in C extended target C NAT0= number of dipoles in "real" target AK1=XX*(4.*PI/(3.*REAL(NAT0)))**(1./3.) AK3=AK1**3 C PIA2=pi*(aeff/d)**2 PIA2=PI*(.75*NAT0/PI)**(2./3.) C C*** Initialize various sums over target orientation. DO JO=1,IORTH QSCSUM(JO)=0. QABSUM(JO)=0. QEXSUM(JO)=0. QBKSUM(JO)=0. QPHSUM(JO)=0. DO J=1,3 QSCGSUM(J,JO)=0. QTRQABSUM(J,JO)=0 QTRQSCSUM(J,JO)=0. ENDDO ENDDO DO IDIR=1,NSCAT S1111(IDIR)=0. S2121(IDIR)=0. CX1121(IDIR)=(0.,0.) ENDDO DO J=1,4 DO JO=1,4 DO IDIR=1,NSCAT SMORI(IDIR,JO,J)=0. ENDDO ENDDO ENDDO C C Loop over target orientations (target rotations in Lab Frame) C (actually, loop over lab rotations in Target Frame) C C**** The outer orientation loop over IBETH=1,NBETH is divided up among C the parallel processes. C IORI=IORI0 c ITNUMMX will store the maximum number of iterations used for any of c the orientations for current size and wavelength, so that this can c be written to the waarbbkccc.sca output file DO J=1,2 ITNUMMX(J)=0 ENDDO DO IBETH=MYID+IBETH1,NBETH,NUMPROCS ITHETA=INT((IBETH-1)/NBETA)+1 THETAD=THETA(ITHETA)*DEGRAD IBETA=MOD(IBETH-1,NBETA)+1 BETAD=BETA(IBETA)*DEGRAD C C*** collinge: for testing C IF(MYID.EQ.0.AND.IORI.EQ.IORI0)THEN IF(ITHETA.NE.ITHETA1)THEN WRITE(CMSGNM,FMT='(A,A)')'Got wrong ', & 'value for ITHETA1' CALL ERRMSG('FATAL','DDSCAT',CMSGNM) ENDIF IF(IBETA.NE.IBETA1)THEN WRITE(CMSGNM,FMT='(A,A)')'Got wrong ', & 'value for IBETA1' CALL ERRMSG('FATAL','DDSCAT',CMSGNM) ENDIF ENDIF C DO IPHI=IPHI1,NPHI PHID=PHI(IPHI)*DEGRAD IORI=(ITHETA-1)*NBETA*NPHI+(IBETA-1)*NPHI+IPHI C C**** Obtain propagation vectors, polarization vectors, and scattering C vectors in Target Frame for desired target rotation in Lab Frame C CALL ROTATE(A1,A2,AK1,CXE01,CXE02,BETA(IBETA), & THETA(ITHETA),PHI(IPHI),EN0R, & CXE01R,CXE02R, & MXSCA,NSCAT,ENSC,EM1,EM2, & AKSR,EM1R,EM2R) C DO J=1,3 AKR(J)=AK1*EN0R(J) ENDDO CALL NAMER(IWAV,IRAD,IORI,CFLSCA,CFLAVG) CALL GETFML(AKR,AK3,AKSR,DX, & CALPHA,CMDSOL,CMDFFT,CMDTRQ,CSHAPE, & CXADIA,CXAOFF,CXALPH,CXALOF, & CXALDS,CXALOS, & CXE,CXE01R,CXE02R,CXEPS, & CXF11,CXF12,CXF21,CXF22,CXRLOC, & CXXI,CXSC,CXSCR1,CXZC,CXZW, & EM1R,EM2R, & IBETH,IBETH1,ICOMP,ICTHM,IDVOUT,INIT, & IOCC,IORTH,IPHI,IPHI1,IPHIM,ITASK,IXYZ0, & ITNUM,MXITER, & LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0, & MXCOMP,MXCXSC,MXN3,MXNAT, & MXNX,MXNY,MXNZ,MXPHI,MXSCA,MYID, & NAT,NAT0,NAT3,NCOMP,NSCAT,NX,NY,NZ, & PHI,PHIN,PIA2, & QABS,QBKSCA,QEXT,QPHA,QSCAT,QSCAG, & QTRQAB,QTRQSC, & SCRRS1,SCRRS2,SHPAR,THETAN,TOL,TIMERS, & MXTIMERS,NTIMERS) DO J=1,2 IF(ITNUM(J).GT.ITNUMMX(J))ITNUMMX(J)=ITNUM(J) ENDDO C C Compute Mueller scattering matrix C CALL GETMUELLER(IBETA,IORTH,IPHI,ITHETA,MXBETA, & MXSCA,MXPHI,MXTHET,NSCAT, & CMDTRQ,PHIN,SM,SMORI, & S1111,S2121, & CX1121,CXE01,CXE02, & CXF11,CXF12,CXF21,CXF22, & CXS1,CXS2,CXS3,CXS4, & QABS,QABSUM,QBKSCA,QBKSUM,QEXSUM, & QEXT,QPHA,QPHSUM,QSCAG,QSCAT, & QSCGSUM,QSCSUM,QTRQAB,QTRQABSUM, & QTRQSC,QTRQSCSUM,WGTA,WGTB) C C---------------------------------------------------------------------- C**** Choose whether or not to write out scattering properties for C current orientation; conditional may be replaced if desired. C C**** If IWRKSC=1, write scattering properties of current orientation: IF(IWRKSC.EQ.1)THEN CALL WRITESCA(MXNX,MXNY,MXNZ,NX,NY,NZ, & DX,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD, & ITHETA,IBETA,IPHI,THETA,MXTHET,NTHETA, & BETA,MXBETA,NBETA,PHI,MXPHI,NPHI,NSMELTS, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & ITNUM,MXITER, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMIND1,SMIND2,SMORI, & THETAN,CX1121,CXE01,CXE01R,CXE02, & CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) ENDIF ENDDO IPHI1=1 C C*** Have completed loop over target orientation IPHI C ENDDO IBETH1=1 C C*** Have completed loop over target orientation IBETH (BETA and THETA) C*** Have completed all loops over target orientations C C*** Collect the partial sums over target orientation C from the parallel processes. This routine is a dummy C if MPI is not being used. C CALL COLSUM(IORTH,MYID,MXSCA,NSCAT, & QSCSUM,QABSUM,QEXSUM,QBKSUM,QPHSUM, & QSCGSUM,QTRQABSUM,QTRQSCSUM, & S1111,S2121,CX1121,SMORI) C C*** Write dielectric information to 'mtable' FREQ=1.E4/WAVE CXEN=SQRT(CXEPS(1)) IF(MYID.EQ.0)THEN WRITE(11,FMT=9058)WAVE,FREQ,CXEN,CXEPS(1) ENDIF C Set IORI=0 prior to calling WRITESCA in order to print orientational C averages C IORI=0 C C******** Now calculate and print out orientational averages ********** C IF(MYID.EQ.0)THEN CALL WRITESCA(MXNX,MXNY,MXNZ,NX,NY,NZ, & DX,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD, & ITHETA,IBETA,IPHI,THETA,MXTHET,NTHETA, & BETA,MXBETA,NBETA,PHI,MXPHI,NPHI,NSMELTS, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & ITNUMMX,MXITER, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMIND1,SMIND2,SMORI, & THETAN,CX1121,CXE01,CXE01R, & CXE02,CXE02R,CXEPS,CXRFR,CXF11,CXF21,CXF12, & CXF22) ENDIF ENDDO IRAD1=1 C C*** Have completed loop over target size C ENDDO C C*** Have completed loop over wavelengths C C close binary file (this perhaps should be moved to writesca) IF(MYID.EQ.0)THEN CLOSE(IOBIN) IF(CNETFLAG.NE.'NOTCDF')THEN CALL WRIMSG('DDSCAT',' close NetCDF') CNETFLAG='NCCLOS' CALL WRITESCA(MXNX,MXNY,MXNZ,NX,NY,NZ, & DX,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD, & ITHETA,IBETA,IPHI,THETA,MXTHET,NTHETA, & BETA,MXBETA,NBETA,PHI,MXPHI,NPHI,NSMELTS, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & ITNUMMX,MXITER, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMIND1,SMIND2,SMORI, & THETAN,CX1121,CXE01,CXE01R,CXE02, & CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) ENDIF ENDIF C C MPI environment shutdown: CALL MPI_FINALIZE(IERR) C STOP'Program DDSCAT terminates normally.' C 9000 FORMAT(' >DDSCAT --- ',A,/, & ' compiled with MXNX=',I3,' MXNY=',I3, & ' MXNZ=',I3) 9011 FORMAT(7X,I7,' = NAT = number of dipoles in extended target',/, & 7X,3I4,' = x,y,z length of extended target (Targ. Frame)') 9020 FORMAT(' DDSCAT --- ',A,/,' TARGET ---',A,/, & ' ',A,' --- method of solution ',/, & ' ',A,' --- prescription for polarizabilities',/, & ' ',A,' --- shape ',/, & I7,' = NAT0 = number of dipoles',/, & ' ',A,' --- source of dielectric function') 9021 FORMAT(' ',A) 9031 FORMAT(' n= (',F7.4,' , ',F7.4,'), epsilon= (',F8.4,' , ',F7.4, & ') for material',I2) 9032 FORMAT(' AEFF=',F10.5,' =',' effective radius (physical units)', & /,' WAVE=',F10.5,' = wavelength (physical units)',/, & ' K*A=',F10.5,' = 2*pi*aeff/lambda') 9043 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & 3X,'aeff',8X,'wave',7X,'Q_ext',6X,'Q_abs',6X,'Q_sca',4X, & 'g(1)=',3X,'Q_bk') 9044 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff',8X,'wave',7X,'Q_ext',7X,'Q_abs',6X,'Q_sca',5X, & 'g(1)=',3X,'Q_bk') 9045 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff wave Q_pha') 9046 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff',7X,'wave',7X,'Q_pha',7X,'Q_pol',7X,'Q_cpol') 9048 FORMAT(' wave(um) f(cm-1) Re(m) Im(m) Re(eps) ', & 'Im(eps)') 9058 FORMAT(1PE10.4,E11.4,4E10.3) C END