PROGRAM CALLTARGET IMPLICIT NONE C C Program CALLTARGET is used to generate target arrays using target C generation routines employed by DDSCAT. C C B.T.Draine, Princeton Univ. Obs., 91.01.01 C History: C 91.05.24 (BTD): Modified to allow entry target names in lower case C 91.07.16 (BTD): Replaced DO...ENDDO with DO #...# CONTINUE C 92.12.16 (BTD): Added IDVSHP to argument list of TARGET C 93.01.07 (BTD): Added options TWOELL and TWOAEL C 94.05.15 (BTD): Added option CONELL C 95.12.11 (BTD): Added option BLOCKS C 96.01.04 (BTD): Added option DW1996 C 96.01.25 (BTD): Added option TWOSPH C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ here and in TARGET C Added calls to SIZER C 96.01.29 (BTD): Support selection of principal axis calculation C in TAR2SP via SHPAR(6) C 97.04.25 (BTD): Added option BLK_AN C 97.04.30 (BTD): Removed code offering choice on target axis for C HEXGON option (target axis now hardwired within C TARHEX) C 97.12.26 (BTD): Introduced lattice spacing vector DX(1-3). C Added code to input DX interactively. C Added DX to argument list in each of 15 separate C calls to TARGET. C 02.02.12 (BTD): Added options NSPHER and PRISM3 C end history C C Copyright (C) 1993,1994,1995,1996,1997,2002 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 MXNAT=max.number of dipoles in target INTEGER MXNAT PARAMETER(MXNAT=200000) C Derived parameters: INTEGER MXN3 PARAMETER(MXN3=3*MXNAT) C Variables: CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6 INTEGER IDVSHP,IDVOUT,IOSHP,J,JJ,NAT,NAT3 INTEGER LXYZ(3) INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3) REAL ERR,F,PI,RGYR2,RXY,RXZ,RYX,RYZ,RZX,RZY, & X,XCM,Y,YCM,Z,ZCM REAL A1(3),A2(3),DX(3),SHPAR(6) DOUBLE PRECISION DXL,DYL,DZL,X2,Y2,Z2 C C*********************************************************************** DATA IDVOUT/0/,IOSHP/15/ PI=4.*ATAN(1.) C Define variables CFLSHP and IDVSHP to please flint: CFLSHP='none' IDVSHP=999 1000 WRITE(IDVOUT,*)'What shape? (Enter choice in quotes) Choices:' WRITE(IDVOUT,*)' BLK_AN = construct from anisotropic cubic ', & 'blocks, input from "blocks.par"' WRITE(IDVOUT,*)' BLOCKS = construct from cubic blocks, input ', & 'from file "blocks.par"' WRITE(IDVOUT,*)' CONELL = two concentric ellipsoids' WRITE(IDVOUT,*)' CYLNDR = cylinder' WRITE(IDVOUT,*)' DW1996 = 13 cube target used by Draine & ', & 'Weingartner (1996)' WRITE(IDVOUT,*)' ELLIPS = ellipsoid' WRITE(IDVOUT,*)' HEXGON = hexagonal prism' WRITE(IDVOUT,*)' RCTNGL = rectangular solid' WRITE(IDVOUT,*)' SPHERE = sphere' WRITE(IDVOUT,*)' TETRAH = regular tetrahedron' WRITE(IDVOUT,*)' THRAEL = three collinear touching anisotropic ', & 'ellipsoids' WRITE(IDVOUT,*)' THRELL = three collinear touching ellipsoids' WRITE(IDVOUT,*)' TWOAEL = two touching anisotropic ellipsoids' WRITE(IDVOUT,*)' TWOELL = two touching isotropic ellipsoids' WRITE(IDVOUT,*)' TWOSPH = two touching spheroids (with angle ', & 'phi between symm. axes' WRITE(IDVOUT,*)' NSPHER = multisphere (union of spheres)' WRITE(IDVOUT,*)' PRISM3 = triangular prism' READ(*,*)CSHAPE WRITE(IDVOUT,*)'Specify lattice spacing, e.g.,' WRITE(IDVOUT,*)' 1 1 1 (for uniform spacing)' WRITE(IDVOUT,*)' .5 1 1 (for double resolution in x direction)' READ(*,*)DX(1),DX(2),DX(3) C C Normalize DX so that DX(1)*DX(2)*DX(3)=1. C X=(DX(1)*DX(2)*DX(3))**(1./3.) DX(1)=DX(1)/X DX(2)=DX(2)/X DX(3)=DX(3)/X WRITE(IDVOUT,*)'Renormalized lattice spacings:' WRITE(IDVOUT,*)' dx/d =',DX(1) WRITE(IDVOUT,*)' dy/d =',DX(2) WRITE(IDVOUT,*)' dz/d =',DX(3) IF(CSHAPE.EQ.'TETRAH'.OR.CSHAPE.EQ.'tetrah')THEN CSHAPE='TETRAH' DO J=1,100 WRITE(IDVOUT,*)'Enter tetrahedron side length/d ', & ' (0 to stop)' READ(*,*)SHPAR(1) SHPAR(2)=SHPAR(1) SHPAR(3)=SHPAR(1) IF(SHPAR(1).LE.0.)STOP CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(*,6000)SHPAR(1),NAT ENDDO STOP ELSEIF(CSHAPE.EQ.'SPHERE'.OR.CSHAPE.EQ.'sphere')THEN CSHAPE='SPHERE' DO J=1,100 WRITE(IDVOUT,*)'Enter sphere diameter/d (0 to stop)' READ(*,*)SHPAR(1) IF(SHPAR(1).LE.0.)stop SHPAR(2)=SHPAR(1) SHPAR(3)=SHPAR(1) CSHAPE='ELLIPS' CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) RGYR2=0. XCM=0. YCM=0. ZCM=0. DO JJ=1,NAT X=IXYZ(JJ,1) Y=IXYZ(JJ,2) Z=IXYZ(JJ,3) XCM=XCM+X YCM=YCM+Y ZCM=ZCM+Z RGYR2=RGYR2+X*X+Y*Y+Z*Z+0.25 ENDDO XCM=XCM/REAL(NAT) YCM=YCM/REAL(NAT) ZCM=ZCM/REAL(NAT) RGYR2=RGYR2/REAL(NAT)-XCM**2-YCM**2-ZCM**2 F=(5./3.)*(4.*PI/REAL(3.*NAT))**(2./3.)*RGYR2 WRITE(*,6100)SHPAR(1),SHPAR(2),SHPAR(3),NAT,F ENDDO STOP ELSEIF(CSHAPE.EQ.'ELLIPS'.OR.CSHAPE.EQ.'ellips')THEN CSHAPE='ELLIPS' DO J=1,100 WRITE(IDVOUT,*)'Enter ellipsoid diameters Dx/d,Dy/d,Dz/d', & ' (0,0,0 to stop)' READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3) IF(SHPAR(1).LE.0.)stop CSHAPE='ELLIPS' CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) XCM=0. YCM=0. ZCM=0. DO JJ=1,NAT X=IXYZ(JJ,1) Y=IXYZ(JJ,2) Z=IXYZ(JJ,3) XCM=XCM+X YCM=YCM+Y ZCM=ZCM+Z ENDDO XCM=XCM/REAL(NAT) YCM=YCM/REAL(NAT) ZCM=ZCM/REAL(NAT) X2=0. Y2=0. Z2=0. DO JJ=1,NAT DXL=IXYZ(JJ,1)-XCM DYL=IXYZ(JJ,2)-YCM DZL=IXYZ(JJ,3)-ZCM X2=X2+DXL*DXL Y2=Y2+DYL*DYL Z2=Z2+DZL*DZL ENDDO X2=X2/NAT+1./12. Y2=Y2/NAT+1./12. Z2=Z2/NAT+1./12. ERR=20.*PI*SQRT(5.*X2*Y2*Z2)/(3.*NAT)-1. RXY=SQRT(X2/Y2) RYZ=SQRT(Y2/Z2) RZX=SQRT(Z2/X2) RYX=1./RXY RZY=1./RYZ RXZ=1./RZX X2=2.*SQRT(5.*X2) Y2=2.*SQRT(5.*Y2) Z2=2.*SQRT(5.*Z2) WRITE(*,6200)SHPAR(1),SHPAR(2),SHPAR(3),NAT,X2,Y2,Z2,ERR, & RYX,RZX,RYZ,RXY,RXZ,RZY ENDDO STOP ELSEIF(CSHAPE.EQ.'CYLNDR'.OR.CSHAPE.EQ.'cylndr')THEN CSHAPE='CYLNDR' DO J=1,100 WRITE(IDVOUT,*)'Enter cylinder length/d (0 to stop)' READ(*,*)SHPAR(1) IF(SHPAR(1).LE.0.)STOP WRITE(IDVOUT,*)'Enter cylinder diameter/d' READ(*,*)SHPAR(2) WRITE(IDVOUT,*)'Enter 1,2,3 for axis along x,y,z' READ(*,*)SHPAR(3) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO STOP ELSEIF(CSHAPE.EQ.'RCTNGL'.OR.CSHAPE.EQ.'rctngl')THEN CSHAPE='RCTNGL' DO J=1,100 WRITE(IDVOUT,*)'Enter x-length/d (0 to stop)' READ(*,*)SHPAR(1) IF(SHPAR(1).LE.0.)STOP WRITE(IDVOUT,*)'Enter y-length/d' READ(*,*)SHPAR(2) WRITE(IDVOUT,*)'Enter z-length/d' READ(*,*)SHPAR(3) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO ELSEIF(CSHAPE.EQ.'HEXGON'.OR.CSHAPE.EQ.'hexgon')THEN CSHAPE='HEXGON' DO J=1,100 WRITE(IDVOUT,*)'Enter length (along symmetry axis)', & ' (0 to stop)' READ(*,*)SHPAR(1) IF(SHPAR(1).LE.0.)STOP WRITE(IDVOUT,*)'Enter (hexagon diameter=2*side width)/d' READ(*,*)SHPAR(2) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO ELSEIF(CSHAPE.EQ.'TWOELL'.OR.CSHAPE.EQ.'twoell')THEN CSHAPE='TWOELL' DO J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions/d (0 0 0 to stop)' READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3) IF(SHPAR(1).LE.0)STOP CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO ELSEIF(CSHAPE.EQ.'TWOAEL'.OR.CSHAPE.EQ.'twoael')THEN CSHAPE='TWOAEL' DO J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions/d (0 0 0 to stop)' READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3) IF(SHPAR(1).LE.0)STOP CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO ELSEIF(CSHAPE.EQ.'THRELL'.OR.CSHAPE.EQ.'threll')THEN CSHAPE='THRELL' DO J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions/d (0 0 0 to stop)' READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3) IF(SHPAR(1).LE.0)STOP CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO ELSEIF(CSHAPE.EQ.'THRAEL'.OR.CSHAPE.EQ.'thrael')THEN CSHAPE='THRAEL' DO J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions/d (0 0 0 to stop)' READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3) IF(SHPAR(1).LE.0)STOP CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT ENDDO ELSEIF(CSHAPE.EQ.'CONELL'.OR.CSHAPE.EQ.'conell')THEN CSHAPE='CONELL' DO J=1,100 WRITE(IDVOUT,*)'Enter length of outer ellipsoid in x,y,z', & 'directions/d (0 0 0 to stop)' READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3) IF(SHPAR(1).LE.0)STOP WRITE(IDVOUT,*)'Enter length of inner ellipsoid in x,y,z', & 'directions/d' READ(*,*)SHPAR(4),SHPAR(5),SHPAR(6) IF(SHPAR(4).LE.SHPAR(1).AND. & SHPAR(5).LE.SHPAR(2).AND. & SHPAR(6).LE.SHPAR(3))THEN CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSE WRITE(IDVOUT,*)'Input error: inner ellipsoid is not ', & 'contained within outer ellipsoid' ENDIF ENDDO ELSEIF(CSHAPE.EQ.'BLOCKS'.OR.CSHAPE.EQ.'blocks')THEN CSHAPE='BLOCKS' CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSEIF(CSHAPE.EQ.'BLK_AN'.OR.CSHAPE.EQ.'blk_an')THEN CSHAPE='BLK_AN' CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSEIF(CSHAPE.EQ.'DW1996'.OR.CSHAPE.EQ.'dw1996')THEN CSHAPE='DW1996' WRITE(0,*)'Enter block size (lattice units)' READ(*,*)SHPAR(1) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSEIF(CSHAPE.EQ.'TWOSPH'.OR.CSHAPE.EQ.'twosph')THEN CSHAPE='TWOSPH' WRITE(0,*)'Enter length a_1/d and diam. b_1/d of 1st spheroid' READ(*,*)SHPAR(1),SHPAR(2) WRITE(0,*)'Enter length a_2/d and diam. b_2/d of 2nd spheroid' READ(*,*)SHPAR(3),SHPAR(4) WRITE(0,*)'Enter angle phi (deg) between axes a_1 and a_2' READ(*,*)SHPAR(5) WRITE(0,*)'Enter PRINAX=0 for a_1=(1,0,0),a_2=(0,1,0) in TF' WRITE(0,*)' =1 for a_1,a_2=principal axes' READ(*,*)SHPAR(6) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSEIF(CSHAPE.EQ.'NSPHER'.OR.CSHAPE.EQ.'nspher')THEN CSHAPE='NSPHER' WRITE(0,*)'Enter DIAMX = target extent in x-direction' READ(*,*)SHPAR(1) WRITE(0,*)'Enter PRINAX=0 for a_1=(1,0,0),a_2=(0,1,0) in TF' WRITE(0,*)' =1 for a_1,a_2 = principal axes' READ(*,*)SHPAR(2) WRITE(0,*)'Enter name of sphere parameter file (e.g., ', & 'tarnsp.par [in quotes])' READ(*,*)CFLSHP CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSEIF(CSHAPE.EQ.'PRISM3'.OR.CSHAPE.EQ.'prism3')THEN CSHAPE='PRISM3' WRITE(0,*)'Enter a/d = first triangle side/d' READ(*,*)SHPAR(1) WRITE(0,*)'Enter b/a , c/a, L/a' READ(*,*)SHPAR(2),SHPAR(3),SHPAR(4) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) ELSE WRITE(IDVOUT,*)'Do not know how to make shape = ',CSHAPE GOTO 1000 ENDIF 6000 FORMAT(' AX=',F7.4,' NAT=',I7) 6100 FORMAT(' AX=',F7.4,' AY=',F7.4,' AZ=',F7.4,' NAT=',I7,' F=',F10.6) 6200 FORMAT(' XV =',F7.4,' YV =',F7.4,' ZV =',F7.4, & ' NAT=',I8,/, & ' 2a_eff=',F7.4,' 2b_eff=',F7.4,' 2c_eff=',F7.4, & ' err=',F8.7,/, & ' y:x =',F7.4,' z:x =',F7.4,' y:z =',F7.4,/, & ' x:y =',F7.4,' x:z =',F7.4,' z:y =',F7.4) END SUBROUTINE CXFFTW(CX,MX,MY,MZ,ISIGN) IMPLICIT NONE C C Arguments: C INTEGER ISIGN,MX,MY,MZ COMPLEX CX(MX,MY,MZ) C*********************************************************************** C Given: C CX = 3d array (complex) C C Returns: C CX = 3d complex fourier transform of input array CX C C Purpose: C Interface to the FFTW (fastest FFT code in the West ?) C Required routines: C FFTW package compiled in SINGLE PRECISION. C FFTW is available from http://www.fftw.org C It was written in C by Matteo Frigo and Steven G. Johnson of MIT. C To link you would typically do the following C f77 TSTFFT.o -L/usr/local/lib -lsfftw (thus the library is called libsfftw.a) C C C Advantages. (a) The code works for arbitrary N and is comparable in speed C at least on an SGI to the GPFA code of Temperton C (more tests needed) C (b) The code has several parallel processing hooks including MPI C Disadvantages C Written in C (for those of us who are FORTRAN programers) C Not always faster than the GPFA of Temperton. Sometimes slower. C C History: C 00.06.21 (PJF): Original version C 00.06.22 (BTD): minor cosmetic changes C end history C Copyright (C) 2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C WARNING: NOT TESTED IN DDSCAT YET C NOT tested on vector machines like CRAY C90 C*********************************************************************** C C Local variables: C INTEGER MXOLD,MYOLD,MZOLD DATA MXOLD,MYOLD,MZOLD/0,0,0/ SAVE MXOLD,MYOLD,MZOLD C include file (this comes with FFTW) INTEGER FFTW_FORWARD,FFTW_BACKWARD PARAMETER (FFTW_FORWARD=-1,FFTW_BACKWARD=1) INTEGER FFTW_REAL_TO_COMPLEX,FFTW_COMPLEX_TO_REAL PARAMETER (FFTW_REAL_TO_COMPLEX=-1,FFTW_COMPLEX_TO_REAL=1) INTEGER FFTW_ESTIMATE,FFTW_MEASURE PARAMETER (FFTW_ESTIMATE=0,FFTW_MEASURE=1) INTEGER FFTW_OUT_OF_PLACE,FFTW_IN_PLACE,FFTW_USE_WISDOM PARAMETER (FFTW_OUT_OF_PLACE=0) PARAMETER (FFTW_IN_PLACE=8,FFTW_USE_WISDOM=16) INTEGER FFTW_THREADSAFE PARAMETER (FFTW_THREADSAFE=128) C Constants for the MPI wrappers (not used here yet INTEGER FFTW_TRANSPOSED_ORDER, FFTW_NORMAL_ORDER INTEGER FFTW_SCRAMBLED_INPUT, FFTW_SCRAMBLED_OUTPUT PARAMETER(FFTW_TRANSPOSED_ORDER=1, FFTW_NORMAL_ORDER=0) PARAMETER(FFTW_SCRAMBLED_INPUT=8192) PARAMETER(FFTW_SCRAMBLED_OUTPUT=16384) C define pointers (whatever it is it is needed by FFTW) INTEGER P_FWD, P_RVS INTEGER N_DIM(3) SAVE P_FWD, P_RVS C C once per each mx, my, mz call C IF(MX.NE.MXOLD .OR. MY.NE.MYOLD .OR. MZ.NE.MZOLD) THEN MXOLD=MX MYOLD=MY MZOLD=MZ N_DIM(1)=MX N_DIM(2)=MY N_DIM(3)=MZ C destroy may not work for the first call, may neeed fix C (there is nothing to destroy) CALL FFTWND_F77_DESTROY_PLAN(P_FWD) CALL FFTWND_F77_DESTROY_PLAN(P_RVS) CALL FFTWND_F77_CREATE_PLAN(P_FWD,3,N_DIM,FFTW_FORWARD, $ FFTW_MEASURE+FFTW_IN_PLACE) CALL FFTWND_F77_CREATE_PLAN(P_RVS,3,N_DIM,FFTW_BACKWARD, $ FFTW_MEASURE+FFTW_IN_PLACE) ENDIF C it seems that ISIGN=1 of Temperton is equivalent to FFTW_BACKWARD in FFTW C I did not test normalization for ISIGN=-1 (can be different from DDSCAT convention) PJF IF(ISIGN.EQ.1) THEN CALL FFTWND_F77_ONE(P_RVS,CX,0) ELSE CALL FFTWND_F77_ONE(P_FWD,CX,0) ENDIF C RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C*********************************************************************** C This version of timeit uses the system call "etime" available under C SunOS and some other bsd-like Unix systems (e.g., Convex unix). C C Copyright (C) 1993,1994,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C 94.06.20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing in CMSGTM='noprint' C Arguments: CHARACTER CMSGTM*(*) REAL DTIME C C External system calls: REAL ETIME EXTERNAL ETIME C C Local variables: REAL OLDTIM CHARACTER CSTA*3,CMSGNM*70 REAL TARRAY(2) C C External subroutines: EXTERNAL WRIMSG C C Local variables to be saved: SAVE CSTA,OLDTIM C C Data statements: DATA CSTA/'ONE'/ C IF(CSTA.EQ.'ONE')THEN CSTA='TWO' OLDTIM=ETIME(TARRAY) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' DTIME=ETIME(TARRAY)-OLDTIM IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) IF(DTIME.GT.1.E4)THEN WRITE(CMSGNM,FMT='(F9.0, A)')DTIME,' CPU time (in secs)' ELSE WRITE(CMSGNM,FMT='(F9.3, A)')DTIME,' CPU time (in secs)' ENDIF CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_convex C C This version of timeit uses the system call "etime" available under C Convex unix and other bsd-like Unix systems. C C 94.06.20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments: CHARACTER CMSGTM* (*) C Local variables: REAL DTIME,OLDTIM CHARACTER CSTA*3,CMSGNM*70 REAL TARRAY(2) C External Subroutines: EXTERNAL WRIMSG C External Functions: REAL ETIME EXTERNAL ETIME C*********************************************************************** DATA CSTA/'ONE'/ IF(CSTA.EQ.'ONE')THEN CSTA='TWO' OLDTIM=ETIME(TARRAY) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' DTIME=ETIME(TARRAY)-OLDTIM IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) ' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_cray C C CRAY version. Call to second is executed. C 94.06.20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C .. Scalar Arguments .. CHARACTER CMSGTM* (*) C .. C .. Local Scalars .. REAL DTIME, OLDTIM CHARACTER CSTA*3, CMSGNM*70 C C .. C save csta, time1 C .. Data statements .. DATA CSTA/'ONE'/ C .. IF(CSTA.EQ.'ONE')THEN CSTA='TWO' call second(time1) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' CALL SECOND(TIME2) DTIME=TIME2-TIME1 IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) ' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_hp C C This version uses the HP-AUX system calls SYSCONF and TIMES C 94.06.20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C External system calls: $ALIAS sysconf = 'sysconf' (%val) INTEGER SYSCONF,TIMES EXTERNAL SYSCONF,TIMES C C Arguments: CHARACTER CMSGTM*(*) C C Local variables: INTEGER*2 NAME INTEGER CLK_TCK,I INTEGER TMS(4) REAL DTIME, OLDTIM CHARACTER CSTA*3,CMSGNM*70 C C External subroutines: EXTERNAL WRIMSG C Data statements: DATA CSTA/'ONE'/ C .. IF(CSTA.EQ.'ONE')THEN NAME=2 CLK_TCK=SYSCONF(NAME) OLDTIM=REAL(TMS(1)+TMS(2))/REAL(CLK_TCK) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' I=TIMES(TMS) DTIME=REAL(TMS(1)+TMS(2))/REAL(CLK_TCK)-OLDTIM IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) ' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_ibm6000 C 94/06/20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C C This version of timeit uses the system call "mclock" available under C RISC6000/320 IBM AIX3.1 C C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments: CHARACTER CMSGTM* (*) C C Local variables: REAL DTIME, OLDTIM CHARACTER CSTA*3, CMSGNM*70 REAL TARRAY(2) C C External Subroutines .. EXTERNAL WRIMSG C C External Functions .. REAL ETIME EXTERNAL ETIME C DATA CSTA/'ONE'/ C IF(CSTA.EQ.'ONE')THEN CSTA='TWO' I1=MCLOCK() ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' I2=MCLOCK() IALL=I2-I1 DTIME=FLOAT(IALL)/100. IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) ' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C*********************************************************************** C This version of timeit uses the system call "etime" available under C SunOS and some other bsd-like Unix systems (e.g., DEC OSF). C C Copyright (C) 1993,1994,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C 94.06.20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing in CMSGTM='noprint' C Arguments: CHARACTER CMSGTM*(*) C C External system calls: REAL ETIME EXTERNAL ETIME C C Local variables: REAL DTIME,OLDTIM CHARACTER CSTA*3,CMSGNM*70 REAL TARRAY(2) C C External subroutines: EXTERNAL WRIMSG C C Local variables to be saved: SAVE CSTA,OLDTIM C C Data statements: DATA CSTA/'ONE'/ C IF(CSTA.EQ.'ONE')THEN CSTA='TWO' OLDTIM=ETIME(TARRAY) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' DTIME=ETIME(TARRAY)-OLDTIM IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) IF(DTIME.GT.1.E4)THEN WRITE(CMSGNM,FMT='(F9.0, A)')DTIME,' CPU time (in secs)' ELSE WRITE(CMSGNM,FMT='(F9.3, A)')DTIME,' CPU time (in secs)' ENDIF CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C*********************************************************************** C This version of timeit uses the SGI IRIX system call "etime" C C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C 94.06.20 (PJF) add dtime to the formal parameters C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C Arguments: CHARACTER CMSGTM*(*) C C External system calls: REAL ETIME EXTERNAL ETIME C C Local variables: REAL DTIME,OLDTIM CHARACTER CSTA*3,CMSGNM*70 REAL TARRAY(2) C C External subroutines: EXTERNAL WRIMSG C C Local variables to be saved: SAVE CSTA,OLDTIM C C Data statements: DATA CSTA/'ONE'/ C IF(CSTA.EQ.'ONE')THEN CSTA='TWO' OLDTIM=ETIME(TARRAY) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' DTIME=ETIME(TARRAY)-OLDTIM IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) IF(DTIME.GT.1.E4)THEN WRITE(CMSGNM,FMT='(F9.0, A)')DTIME,' CPU time (in secs)' ELSE WRITE(CMSGNM,FMT='(F9.3, A)')DTIME,' CPU time (in secs)' ENDIF CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_titan C C This version of TIMEIT is for Titan computers. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C 94/06/20 (PJF) add dtime to the formal parameters C*********************************************************************** C .. Scalar Arguments .. CHARACTER CMSGTM* (*) C .. C .. Local Scalars .. REAL DTIME, OLDTIM, TIME INTEGER IAD CHARACTER CSTA*3, CMSGNM*70 C .. C .. External Subroutines .. EXTERNAL WRIMSG C .. C .. External Functions .. REAL CPUTIM EXTERNAL CPUTIM C .. C .. Data statements .. DATA CSTA/'ONE'/ C .. IF(CSTA.EQ.'ONE')THEN CSTA='TWO' OLDTIM=CPUTIM(0.) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' TIME=CPUTIM(0.) IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) ' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_vms C 94.06.20 (PJF) add dtime to the formal parameters C 96.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C This version of TIMEIT is for VMS systems. C C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. c*********************************************************************** C .. Scalar Arguments .. CHARACTER CMSGTM* (*) C .. C .. Local Scalars .. REAL DTIME, OLDTIM, TIME INTEGER IAD CHARACTER CSTA*3, CMSGNM*70 C .. C .. External Subroutines .. EXTERNAL LIB$INIT_TIMER, LIB$SHOW_TIMER, WRIMSG C .. C .. Data statements .. DATA CSTA/'ONE'/ C .. IF(CSTA.EQ.'ONE')THEN CSTA='TWO' IAD=0 CALL LIB$INIT_TIMER(IAD) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN CALL LIB$SHOW_TIMER(IAD,2) WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END PROGRAM TSTFFT IMPLICIT NONE C Parameters: INTEGER MX,MY,MZ,NF235MX PARAMETER (MX=200, MY=200, MZ=200, NF235MX=45) C Local variables: INTEGER IN,N,NX,NY,NZ INTEGER NF235(NF235MX) COMPLEX A(MX,MY,MZ),C(MX,MY,MZ),CW(MX,MY,MZ) C*********************************************************************** C Purpose: C Tests different 3D FFT implementations (CPU time). C Required routines: C cfft99f.o cxfft3.o cxfft3n.o cxfftw.o C errmsg.o fourx.o C gpfa.o timeit_sgi.o tst.o wrimsg.o gen.o C C History: C 94.06.20 (PJF): Original version C 94.12.18 (BTD): Minor changes in output formatting. C Explicit declaration of all variables. C 95.07.12 (BTD): Minor changes in formatting. C 96.11.14 (PJF): getset not needed C 00.07.05 (PJF): modified to compare Brenner, GPFA, and FFTW C 03.07.15 (BTD): modified to also do non-prime-factorable cases C Copyright (C) 1994,1995,2000,2003 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C*********************************************************************** C Elements of NF235 are numbers.le.4096 which can be factored by 2,3,5: DATA NF235/2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40, & 45,48,50,54,60,64,72,75,80,81,90,96,100,108,120,125,128,135, & 144,150,160,162,180,192,200/ WRITE(0,*)'Calculation time for different 3-D FFT methods' WRITE(0,*)'=========================', & '==========================' WRITE(0,*)' ' open(unit=20,file='err.dat',status='unknown') open(unit=21,file='res.dat',status='unknown') WRITE(21,6021) 6021 FORMAT(' CPU time (sec) for different 3-D FFT methods',/, & ' Machine=2.0 GHz Xeon, 256kB L2 cache, pgf77 compiler ',/, & ' parameter LVR = 64 ',/, & ' LVR="length of vector registers" in gpfa2f,gpfa3f,gpfa5f',/, & '============================================================', & /, & ' Brenner Temperton ', & 'Frigo & Johnson',/, & ' NX NY NZ (FOURX) (GPFA) ', & ' (FFTW)') rewind(20) C C IN=33 corresponds to 100x100x100 FFT. C On 2 GHz Xeon, this requires approx 1.6 sec for Brenner, C 0.23 sec for FFTW, C 0.18 sec for GPFA C IN=40 corresponds to 150x150x150 FFT. C On 2 GHz Xeon, this requires approx 6.9 sec for Brenner, C 0.87 sec for FFTW, C 0.68 sec for GPFA DO IN=18,45 IF(NF235(IN-1).LT.NF235(IN)-1)THEN NX=NINT(.5*(NF235(IN-1)+NF235(IN))) NY=NX NZ=NX CALL CASE(A,C,CW,NX,NY,NZ,NF235,NF235MX) ENDIF Nx=NF235(IN) Ny=NX Nz=NX CALL CASE(A,C,CW,NX,NY,NZ,NF235,NF235MX) ENDDO STOP END SUBROUTINE CASE(A,C,CW,NX,NY,NZ,NF235,NF235MX) IMPLICIT NONE C Arguments: INTEGER NX,NY,NZ,NF235MX INTEGER NF235(NF235MX) COMPLEX A(NX,NY,NZ),C(NX,NY,NZ),CW(NX,NY,NZ) C Local variables: INTEGER I,J,JX,JY,JZ,K,L,NRPT INTEGER NN(3) REAL T0,T1,T2,T3 COMPLEX C1, C2, C3, C4, C5, C6 C*********************************************************************** C Purpose: interface to 3D FFT codes (for use in FFT tests) C History: C 94.06.20 (PJF): Original version C 94.12.19 (BTD): Explicit declaration of all variables. C 03.07.15 (BTD): Modified to all use of NX,NY,NZ not factorizable C as 2^a3^b5^c C 03.07.15 (BTD): Modified to repeat for more accurate timings C Copyright (C) 1994,2003, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. ***********************************************************************C C--------- BRENNR NN(1)=NX NN(2)=NY NN(3)=NZ CALL GEN(A,NX,NY,NZ) NRPT=1+INT(1E5/(NX*NY*NZ)) CALL TIMEIT('BREN',T0) DO L=1,NRPT DO K=1,NZ DO J=1,NY DO I=1,NX C(I,J,K)=A(I,J,K) ENDDO ENDDO ENDDO CALL FOURX(C,NN,3,1) ENDDO CALL TIMEIT('BREN',T1) T1=T1/REAL(NRPT) C1=C(1,1,1) C2=C(2,1,1) C--------- GPFA C initialize c check whether NX,NY,NZ are of form 2^i3^j5^k JX=0 JY=0 JZ=0 DO J=1,NF235MX IF(NX.LE.NF235(J).AND.JX.EQ.0)JX=NF235(J) IF(NY.LE.NF235(J).AND.JY.EQ.0)JY=NF235(J) IF(NZ.LE.NF235(J).AND.JZ.EQ.0)JZ=NF235(J) ENDDO CALL CXFFT3n(C,JX,JY,JZ,1) CALL GEN(A,NX,NY,NZ) NRPT=1+INT(1E6/(JX*JY*JZ)) CALL TIMEIT('GPFA',T0) DO L=1,NRPT DO K=1,JZ DO J=1,JY DO I=1,JX C(I,J,K)=A(I,J,K) ENDDO ENDDO ENDDO CALL CXFFT3n(C,JX,JY,JZ,1) ENDDO CALL TIMEIT('GPFA',T2) T2=T2/REAL(NRPT) C5=C(1,1,1) C6=C(2,1,1) C--------- FFTW C initialize CALL CXFFTW(C,NX,NY,NZ,1) CALL GEN(A,NX,NY,NZ) NRPT=1+INT(1E6/(NX*NY*NZ)) CALL TIMEIT('FFTW',T0) DO L=1,NRPT DO K=1,NZ DO J=1,NY DO I=1,NX C(I,J,K)=A(I,J,K) ENDDO ENDDO ENDDO CALL CXFFTW(C,NX,NY,NZ,1) ENDDO CALL TIMEIT('FFTW',T3) T3=T3/REAL(NRPT) C3=C(1,1,1) C4=C(2,1,1) C--------- Print WRITE(21, '(1X,3I4,3F12.3)') NX,NY,NZ,T1,T2,T3 WRITE(20,100) NX, C1,C2 WRITE(20,100) NX, C3,C4 WRITE(20,100) NX, C5,C6 100 FORMAT(1X,I5,4G15.4) RETURN END subroutine gen(c,nx,ny,nz) C Arguments: INTEGER NX,NY,NZ complex c(nx,ny,nz) C Local variables: INTEGER I1,I2,I3 REAL F C Purpose: C generates a 3D matrix for use in 3D FFT CPU timmings do 103 i1=1, nx do 102 i2=1, ny do 101 i3=1, nz c f=i1*i2*i3 f=sin(float(i1))+sin(float(i2))+sin(float(i3)) c(i1,i2,i3)=cmplx(f*1.,f*0.5) 101 CONTINUE 102 CONTINUE 103 continue return end SUBROUTINE WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,DX,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & ITHETA,IBETA,IPHI, & IORI,AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR, & CXE01R,CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G, & QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,SM,SMORI, & PHIN,THETAN,TIMERS) C*********************************************************************** C Subroutine "writenet" writes NetCDF portable binary file C to file cnetfile ('dd.cdf') C C To use netCDF option c (a) set "include" statement in this routine to proper location of c the file "netcdf.inc" c (b) set "ddscat.par" variable "cnetflag" to 'ALLCDF' or 'ORICDF', c (c) modify Makefile to link with NetCDF library, by appropriately c defining the strings LIBNETCDF and LINKNETCDF c (see comments in the Makefile) C C Original version created by P.J.Flatau, University of California, SD C History: C 96.11.15 (PJF): First version of NetCDF (Wolff's request) C 5a6 ---> 5a7 Converted to subroutine. C 96.11.21 (PJF): Add IF(CNEFLAG.EQ.'ALLCDF')THEN and ENDIF C 96.11.21 (BTD): Removed second occurrence of MXSCA in argument list C here and call from WRITESCA C 96.11.21 (BTD): added call to NCCLOS if called with CNETFLAG='NCCLOS' C in order to get call to NCCLOS out of DDSCAT.f C 98.01.11 (BTD): added DX(3) to argument list C added variable IDDX C 03.10.30 (BTD): corrected typo (missing & in col 6) C end history C Copyright (C) 1996,1998 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C IMPORTANT warning: C NetCDF "include" statement. "include" is NOT standard FORTRAN 77. C For location of "netcdf.inc" ask your system administrator. C Flatau's system: C include '/usr/local/netcdf-2.4.2/include/netcdf.inc' C Draine's system: include '/u/flatau/netcdf-2.4.3/include/netcdf.inc' C*********************************************************************** C .. Scalar Arguments .. REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD, & THTMID,THTMXD,TOL,WAVE,XX INTEGER IBETA,ICTHM,IDNC,IORI,IORTH,IPHI,IPHIM,IRAD,ITHETA,IWAV, & MXBETA,MXCOMP,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA, & MXTHET,MXTIMERS,MXWAV,NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI, & NRAD,NSCAT,NTHETA,NTIMERS,NWAV,NX,NY,NZ CHARACTER CALPHA*6,CMDFFT*6,CMDTRQ*6,CNETFLAG*6,CSHAPE*6, & CNETFILE*14,CSTAMP*22,CDESCR*67 C .. Array Arguments .. COMPLEX CXE01R(3),CXE02R(3),CXEPS(MXCOMP),CXRFR(MXCOMP) REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),DX(3),G(2), & PHI(MXPHI),PHIN(MXSCA), & QABS(2),QAV(17),QBKSCA(2),QEXT(2),QPHA(2), & QSCAT(2),QTRQAB(3,2),QTRQSC(3,2),SM(MXSCA,4,4), & SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA), & TIMERS(MXTIMERS),WAVEA(MXWAV) INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3) C .. Local Scalars .. INTEGER ID17,ID2,ID3,IDA1,IDA2,IDAEFF,IDAEFFA,IDAK1,IDBETA, & IDBETMID,IDBETMXD,IDC22,IDC6,IDC67,IDCALPHA,IDCDESCR, & IDCMDFFT,IDCMDTRQ,IDCSHAPE,IDCSTAMP,IDCXRFR,IDDX,IDG, & IDICOMP,IDICTHM,IDIORTH,IDIPHIM,IDIX0,IDIY0,IDIZ0,IDNAT, & IDNAT0,IDNAT3,IDNBETA,IDNCOMP,IDNORI,IDNPHI,IDNRAD, & IDNSCAT,IDNTHETA,IDNTIMERS,IDNWAV,IDNX,IDNY,IDNZ,IDONE, & IDPHI,IDPHIMID,IDPHIMXD,IDPHIN,IDQABS,IDQAV,IDQBKSCA, & IDQEXT,IDQPHA,IDQSCAT,IDQTRQAB,IDQTRQSC,IDSM11,IDSM12, & IDSM13,IDSM14,IDSM21,IDSM22,IDSM23,IDSM24,IDSM31,IDSM32, & IDSM33,IDSM34,IDSM41,IDSM42,IDSM43,IDSM44,IDSMORI11, & IDSMORI12,IDSMORI13,IDSMORI14,IDSMORI21,IDSMORI22, & IDSMORI23,IDSMORI24,IDSMORI31,IDSMORI32,IDSMORI33, & IDSMORI34,IDSMORI41,IDSMORI42,IDSMORI43,IDSMORI44, & I,IDTHETA,IDTHETAN,DTHTMID,IDTHTMXD,IDTOL,IDWAVE,IDWAVEA, & IDXX,IENTRY,IERR C (defined in netcdf.inc) & NCCHAR,NCCLOB,NCFLOAT,NCLONG,NCSHORT C .. Local Arrays .. INTEGER ICOUNT(7),IDDIM(7),ISTART(7) C .. External Functions .. C (defined in netcdf.inc) INTEGER NCCRE,NCDDEF,NCVDEF C (defined in netcdf.inc) EXTERNAL NCCRE,NCDDEF,NCVDEF C .. External Subroutines .. EXTERNAL NCAPTC,NCENDF,NCVPT,NCVPT1,NCVPTC,WRIMSG C .. Save statement .. SAVE IENTRY,IDQEXT,IDQABS,IDQSCAT,IDG,IDQBKSCA,IDQPHA,IDQAV, & IDQTRQAB,IDQTRQSC,IDSM11,IDSM12,IDSM13,IDSM14,IDSM21,IDSM22, & IDSM23,IDSM24,IDSM31,IDSM32,IDSM33,IDSM34,IDSM41,IDSM42, & IDSM43,IDSM44,IDSMORI11,IDSMORI12,IDSMORI13,IDSMORI14, & IDSMORI21,IDSMORI22,IDSMORI23,IDSMORI24,IDSMORI31,IDSMORI32, & IDSMORI33,IDSMORI34,IDSMORI41,IDSMORI42,IDSMORI43,IDSMORI44, & IDWAVE,IDAEFF,IDXX,IDAK1,IDCXRFR C .. Data statements .. DATA IENTRY/0/ C 1) IENTRY=0 defines first entry to this routine. All NetCDF C dimensions will be defined in this step. This part of the C code is never executed again.NetCDF "id" variables are saved. C 2) iori=0. Routine is called for iwav, irad, itheta, ibeta, iphi C These 5 variables can be split to two groups: C ("iwav", "irad") and ("itheta", "iphi", "ibeta"). C The second group defines orientations of target and is also C defined by variable "iori". If "iori=0" the code C finished all specified ntheta, nbeta, nphi orientations C and it is time to write orientational averages. C 3) If iori .ne. 0 we cycle over orientations and iwav, irad. C Many-dimensional arrays are stored on disk. IF(CNETFLAG.EQ.'NCCLOS')THEN CALL NCCLOS(IDNC,IERR) RETURN ENDIF IF (IENTRY.EQ.0) THEN CALL WRIMSG('WRITENET',' First entry to NetCDF write') IENTRY = 1 C this dd.cdf should be "cnetfile" IDNC = NCCRE('dd.cdf',NCCLOB,IERR) C NetCDF dimensions C character variables IDC6 = NCDDEF(IDNC,'c6',6,IERR) IDC67 = NCDDEF(IDNC,'c67',67,IERR) IDC22 = NCDDEF(IDNC,'c22',22,IERR) IDCALPHA = NCVDEF(IDNC,'calpha',NCCHAR,1,IDC6,IERR) IDCDESCR = NCVDEF(IDNC,'cdescr',NCCHAR,1,IDC67,IERR) IDCMDFFT = NCVDEF(IDNC,'cmdfft',NCCHAR,1,IDC6,IERR) IDCMDTRQ = NCVDEF(IDNC,'cmdtrq',NCCHAR,1,IDC6,IERR) IDCSHAPE = NCVDEF(IDNC,'cshape',NCCHAR,1,IDC6,IERR) IDCSTAMP = NCVDEF(IDNC,'cstamp',NCCHAR,1,IDC22,IERR) C single variables IDONE = NCDDEF(IDNC,'one',1,IERR) IDIORTH = NCVDEF(IDNC,'iorth',NCLONG,1,IDONE,IERR) IDBETMID = NCVDEF(IDNC,'betmid',NCFLOAT,1,IDONE,IERR) IDBETMXD = NCVDEF(IDNC,'betmxd',NCFLOAT,1,IDONE,IERR) IDTHTMID = NCVDEF(IDNC,'thtmid',NCFLOAT,1,IDONE,IERR) IDTHTMXD = NCVDEF(IDNC,'thtmxd',NCFLOAT,1,IDONE,IERR) IDPHIMID = NCVDEF(IDNC,'phimid',NCFLOAT,1,IDONE,IERR) IDPHIMXD = NCVDEF(IDNC,'phimxd',NCFLOAT,1,IDONE,IERR) IDICTHM = NCVDEF(IDNC,'icthm',NCLONG,1,IDONE,IERR) IDIPHIM = NCVDEF(IDNC,'iphim',NCLONG,1,IDONE,IERR) IDTOL = NCVDEF(IDNC,'tol',NCFLOAT,1,IDONE,IERR) C dimensions arrays IDNX = NCDDEF(IDNC,'nx',NX,IERR) IDNY = NCDDEF(IDNC,'ny',NY,IERR) IDNZ = NCDDEF(IDNC,'nz',NZ,IERR) IDNAT0 = NCDDEF(IDNC,'nat0',NAT0,IERR) IDNAT = NCDDEF(IDNC,'nat',NAT,IERR) IDNAT3 = NCDDEF(IDNC,'nat3',NAT3,IERR) IDNSCAT = NCDDEF(IDNC,'nscat',NSCAT,IERR) IDNCOMP = NCDDEF(IDNC,'ncomp',NCOMP,IERR) IDNORI = NCDDEF(IDNC,'nori',NORI,IERR) IDNWAV = NCDDEF(IDNC,'nwav',NWAV,IERR) IDNRAD = NCDDEF(IDNC,'nrad',NRAD,IERR) IDNTHETA = NCDDEF(IDNC,'ntheta',NTHETA,IERR) IDNPHI = NCDDEF(IDNC,'nphi',NPHI,IERR) IDNBETA = NCDDEF(IDNC,'nbeta',NBETA,IERR) IDNTIMERS = NCDDEF(IDNC,'ntimers',NTIMERS,IERR) ID2 = NCDDEF(IDNC,'two',2,IERR) ID3 = NCDDEF(IDNC,'three',3,IERR) ID17 = NCDDEF(IDNC,'seventeen',17,IERR) C arrays IDWAVEA = NCVDEF(IDNC,'wavea',NCFLOAT,1,IDNWAV,IERR) CALL NCAPTC(IDNC,IDWAVEA,'long_name',NCCHAR,16, & 'wavelength array',IERR) IDAEFFA = NCVDEF(IDNC,'aeffa',NCFLOAT,1,IDNRAD,IERR) CALL NCAPTC(IDNC,IDAEFFA,'long_name',NCCHAR,22, & 'effective radius array',IERR) IDTHETA = NCVDEF(IDNC,'theta',NCFLOAT,1,IDNTHETA,IERR) CALL NCAPTC(IDNC,IDTHETA,'long_name',NCCHAR,38, & 'theta (target orientation angle array)',IERR) IDBETA = NCVDEF(IDNC,'beta',NCFLOAT,1,IDNBETA,IERR) CALL NCAPTC(IDNC,IDBETA,'long_name',NCCHAR,37, & 'beta (target orientation angle array)',IERR) IDPHI = NCVDEF(IDNC,'phi',NCFLOAT,1,IDNPHI,IERR) CALL NCAPTC(IDNC,IDPHI,'long_name',NCCHAR,36, & 'phi (target orientation angle array)',IERR) IDA1 = NCVDEF(IDNC,'a1',NCFLOAT,1,ID3,IERR) CALL NCAPTC(IDNC,IDA1,'long_name',NCCHAR,24, & 'target definition vector',IERR) IDA2 = NCVDEF(IDNC,'a2',NCFLOAT,1,ID3,IERR) CALL NCAPTC(IDNC,IDA2,'long_name',NCCHAR,24, & 'target definition vector',IERR) IDDX = NCVDEF(IDNC,'dx',NCFLOAT,1,ID3,IERR) CALL NCAPTC(IDNC,IDDX,'long_name',NCCHAR,22, & 'lattice spacing vector',IERR) IDICOMP = NCVDEF(IDNC,'icomp',NCSHORT,1,IDNAT3,IERR) CALL NCAPTC(IDNC,IDICOMP,'long_name',NCCHAR,39, & 'dipole composition definition (0-empty)',IERR) C let as split ixyz0 IDIX0 = NCVDEF(IDNC,'ix0',NCSHORT,1,IDNAT,IERR) IDIY0 = NCVDEF(IDNC,'iy0',NCSHORT,1,IDNAT,IERR) IDIZ0 = NCVDEF(IDNC,'iz0',NCSHORT,1,IDNAT,IERR) C declare "running" variables dimensions cc WRITE(IOBIN) IWAV, IRAD, IORI, (timers(i), i=1,ntimers) cc WRITE(iobin) AEFF,WAVE,XX, AK1, BETAD,THETAD,PHID, cc $ AKR, CXE01R, CXE02R, cc $ (CXRFR(i), i=1, ncomp), (cxeps(i),i=1, ncomp) cc WRITE(iobin) QEXT, QABS, QSCAT, G, QBKSCA, QPHA, QAV, cc $ QTRQAB, QTRQSC c idone=ncddef(idnc,'one',1,ierr) c idiorth=ncvdef(idnc,'iorth',nclong,1, idone,ierr) C IDDIM(1) = IDNWAV IDWAVE = NCVDEF(IDNC,'wave',NCFLOAT,1,IDDIM,IERR) IDDIM(1) = IDNRAD IDAEFF = NCVDEF(IDNC,'aeff',NCFLOAT,1,IDDIM,IERR) IDDIM(1) = IDNWAV IDDIM(2) = IDNRAD IDXX = NCVDEF(IDNC,'xx',NCFLOAT,2,IDDIM,IERR) IDAK1 = NCVDEF(IDNC,'ak1',NCFLOAT,2,IDDIM,IERR) C cxrfr is complex but NetCDF doesn't support complex: IDDIM(1) = IDNWAV IDDIM(2) = ID2 IDDIM(3) = IDNCOMP IDCXRFR = NCVDEF(IDNC,'cxrfr',NCFLOAT,3,IDDIM,IERR) IDDIM(1) = IDNWAV IDDIM(2) = IDNRAD IDDIM(3) = IDNSCAT IDSMORI11 = NCVDEF(IDNC,'smori11',NCFLOAT,3,IDDIM,IERR) IDSMORI12 = NCVDEF(IDNC,'smori12',NCFLOAT,3,IDDIM,IERR) IDSMORI13 = NCVDEF(IDNC,'smori13',NCFLOAT,3,IDDIM,IERR) IDSMORI14 = NCVDEF(IDNC,'smori14',NCFLOAT,3,IDDIM,IERR) IDSMORI21 = NCVDEF(IDNC,'smori21',NCFLOAT,3,IDDIM,IERR) IDSMORI22 = NCVDEF(IDNC,'smori22',NCFLOAT,3,IDDIM,IERR) IDSMORI23 = NCVDEF(IDNC,'smori23',NCFLOAT,3,IDDIM,IERR) IDSMORI24 = NCVDEF(IDNC,'smori24',NCFLOAT,3,IDDIM,IERR) IDSMORI31 = NCVDEF(IDNC,'smori31',NCFLOAT,3,IDDIM,IERR) IDSMORI32 = NCVDEF(IDNC,'smori32',NCFLOAT,3,IDDIM,IERR) IDSMORI33 = NCVDEF(IDNC,'smori33',NCFLOAT,3,IDDIM,IERR) IDSMORI34 = NCVDEF(IDNC,'smori34',NCFLOAT,3,IDDIM,IERR) IDSMORI41 = NCVDEF(IDNC,'smori41',NCFLOAT,3,IDDIM,IERR) IDSMORI42 = NCVDEF(IDNC,'smori42',NCFLOAT,3,IDDIM,IERR) IDSMORI43 = NCVDEF(IDNC,'smori43',NCFLOAT,3,IDDIM,IERR) IDSMORI44 = NCVDEF(IDNC,'smori44',NCFLOAT,3,IDDIM,IERR) C real 6-7 dimensional cases IDDIM(1) = IDNWAV IDDIM(2) = IDNRAD IDDIM(3) = IDNTHETA IDDIM(4) = IDNBETA IDDIM(5) = IDNPHI IDDIM(6) = ID2 IDQEXT = NCVDEF(IDNC,'qext',NCFLOAT,6,IDDIM,IERR) IDQABS = NCVDEF(IDNC,'qabs',NCFLOAT,6,IDDIM,IERR) IDQSCAT = NCVDEF(IDNC,'qscat',NCFLOAT,6,IDDIM,IERR) IDG = NCVDEF(IDNC,'g',NCFLOAT,6,IDDIM,IERR) IDQBKSCA = NCVDEF(IDNC,'qbksca',NCFLOAT,6,IDDIM,IERR) IDQPHA = NCVDEF(IDNC,'qpha',NCFLOAT,6,IDDIM,IERR) IDDIM(6) = ID17 IDQAV = NCVDEF(IDNC,'qav',NCFLOAT,6,IDDIM,IERR) IDDIM(6) = ID3 IDDIM(7) = ID2 IDQTRQAB = NCVDEF(IDNC,'qtrqab',NCFLOAT,7,IDDIM,IERR) IDQTRQSC = NCVDEF(IDNC,'qtrqsc',NCFLOAT,7,IDDIM,IERR) IDDIM(6) = IDNSCAT IF(CNETFLAG.EQ.'ALLCDF')THEN IDSM11 = NCVDEF(IDNC,'sm11',NCFLOAT,6,IDDIM,IERR) IDSM12 = NCVDEF(IDNC,'sm12',NCFLOAT,6,IDDIM,IERR) IDSM13 = NCVDEF(IDNC,'sm13',NCFLOAT,6,IDDIM,IERR) IDSM14 = NCVDEF(IDNC,'sm14',NCFLOAT,6,IDDIM,IERR) IDSM21 = NCVDEF(IDNC,'sm21',NCFLOAT,6,IDDIM,IERR) IDSM22 = NCVDEF(IDNC,'sm22',NCFLOAT,6,IDDIM,IERR) IDSM23 = NCVDEF(IDNC,'sm23',NCFLOAT,6,IDDIM,IERR) IDSM24 = NCVDEF(IDNC,'sm24',NCFLOAT,6,IDDIM,IERR) IDSM31 = NCVDEF(IDNC,'sm31',NCFLOAT,6,IDDIM,IERR) IDSM32 = NCVDEF(IDNC,'sm32',NCFLOAT,6,IDDIM,IERR) IDSM33 = NCVDEF(IDNC,'sm33',NCFLOAT,6,IDDIM,IERR) IDSM34 = NCVDEF(IDNC,'sm34',NCFLOAT,6,IDDIM,IERR) IDSM41 = NCVDEF(IDNC,'sm41',NCFLOAT,6,IDDIM,IERR) IDSM42 = NCVDEF(IDNC,'sm42',NCFLOAT,6,IDDIM,IERR) IDSM43 = NCVDEF(IDNC,'sm43',NCFLOAT,6,IDDIM,IERR) IDSM44 = NCVDEF(IDNC,'sm44',NCFLOAT,6,IDDIM,IERR) ENDIF IDPHIN = NCVDEF(IDNC,'phin',NCFLOAT,1,IDNSCAT,IERR) CALL NCAPTC(IDNC,IDPHIN,'long_name',NCCHAR,16, & 'scattering angle',IERR) IDTHETAN = NCVDEF(IDNC,'thetan',NCFLOAT,1,IDNSCAT,IERR) CALL NCAPTC(IDNC,IDTHETAN,'long_name',NCCHAR,16, & 'scattering angle',IERR) CALL NCENDF(IDNC,IERR) C write characters CALL NCVPTC(IDNC,IDCALPHA,1,6,CALPHA,6,IERR) CALL NCVPTC(IDNC,IDCDESCR,1,67,CDESCR,67,IERR) CALL NCVPTC(IDNC,IDCMDFFT,1,6,CMDFFT,6,IERR) CALL NCVPTC(IDNC,IDCMDTRQ,1,6,CMDTRQ,6,IERR) CALL NCVPTC(IDNC,IDCSHAPE,1,6,CSHAPE,6,IERR) CALL NCVPTC(IDNC,IDCSTAMP,1,22,CSTAMP,22,IERR) C write single variables: CALL NCVPT1(IDNC,IDIORTH,1,IORTH,IERR) CALL NCVPT1(IDNC,IDBETMID,1,BETMID,IERR) CALL NCVPT1(IDNC,IDBETMXD,1,BETMXD,IERR) CALL NCVPT1(IDNC,IDTHTMID,1,THTMID,IERR) CALL NCVPT1(IDNC,IDTHTMXD,1,THTMXD,IERR) CALL NCVPT1(IDNC,IDPHIMID,1,PHIMID,IERR) CALL NCVPT1(IDNC,IDPHIMXD,1,PHIMXD,IERR) CALL NCVPT1(IDNC,IDICTHM,1,ICTHM,IERR) CALL NCVPT1(IDNC,IDIPHIM,1,IPHIM,IERR) CALL NCVPT1(IDNC,IDTOL,1,TOL,IERR) C write arrays CALL NCVPT(IDNC,IDWAVEA,1,NWAV,WAVEA,IERR) CALL NCVPT(IDNC,IDAEFFA,1,NRAD,AEFFA,IERR) CALL NCVPT(IDNC,IDTHETA,1,NTHETA,THETA,IERR) CALL NCVPT(IDNC,IDBETA,1,NBETA,BETA,IERR) CALL NCVPT(IDNC,IDPHI,1,NPHI,PHI,IERR) CALL NCVPT(IDNC,IDA1,1,3,A1,IERR) CALL NCVPT(IDNC,IDA2,1,3,A2,IERR) CALL NCVPT(IDNC,IDDX,1,3,DX,IERR) CALL NCVPT(IDNC,IDICOMP,1,NAT3,ICOMP,IERR) C print, idix0 CALL NCVPT(IDNC,IDIX0,1,NAT,IXYZ0(1,1),IERR) CALL NCVPT(IDNC,IDIY0,1,NAT,IXYZ0(1,2),IERR) CALL NCVPT(IDNC,IDIZ0,1,NAT,IXYZ0(1,3),IERR) CALL NCVPT(IDNC,IDPHIN,1,NSCAT,PHIN,IERR) CALL NCVPT(IDNC,IDTHETAN,1,NSCAT,THETAN,IERR) C endif "if(ientry .eq. 0) then" END IF IF (IORI.NE.0 .AND. CNETFLAG.EQ.'ALLCDF') THEN C NetCDF option --- running variables C first 5 dimensions are always the same: C define "lower left" corner of the (many dimensional) box C some stuff will be written over and over ISTART(1) = IWAV ICOUNT(1) = 1 CALL NCVPT(IDNC,IDWAVE,ISTART,ICOUNT,WAVE,IERR) ISTART(1) = IRAD ICOUNT(1) = 1 CALL NCVPT(IDNC,IDAEFF,ISTART,ICOUNT,AEFF,IERR) ISTART(1) = IWAV ISTART(2) = IRAD ICOUNT(1) = 1 ICOUNT(2) = 1 CALL NCVPT(IDNC,IDXX,ISTART,ICOUNT,XX,IERR) CALL NCVPT(IDNC,IDAK1,ISTART,ICOUNT,AK1,IERR) ISTART(1) = IWAV ISTART(2) = 1 ISTART(3) = 1 ICOUNT(1) = 1 ICOUNT(2) = 2 ICOUNT(3) = NCOMP CALL NCVPT(IDNC,IDCXRFR,ISTART,ICOUNT,CXRFR,IERR) C ISTART(1) = IWAV ISTART(2) = IRAD ISTART(3) = ITHETA ISTART(4) = IBETA ISTART(5) = IPHI ISTART(6) = 1 ISTART(7) = 1 C define how many elements to write ICOUNT(1) = 1 ICOUNT(2) = 1 ICOUNT(3) = 1 ICOUNT(4) = 1 ICOUNT(5) = 1 ICOUNT(6) = 2 CALL NCVPT(IDNC,IDQEXT,ISTART,ICOUNT,QEXT,IERR) CALL NCVPT(IDNC,IDQABS,ISTART,ICOUNT,QABS,IERR) CALL NCVPT(IDNC,IDQSCAT,ISTART,ICOUNT,QSCAT,IERR) CALL NCVPT(IDNC,IDG,ISTART,ICOUNT,G,IERR) CALL NCVPT(IDNC,IDQBKSCA,ISTART,ICOUNT,QBKSCA,IERR) CALL NCVPT(IDNC,IDQPHA,ISTART,ICOUNT,QPHA,IERR) ICOUNT(6) = 17 CALL NCVPT(IDNC,IDQAV,ISTART,ICOUNT,QAV,IERR) ICOUNT(6) = 3 ICOUNT(7) = 2 CALL NCVPT(IDNC,IDQTRQAB,ISTART,ICOUNT,QTRQAB,IERR) CALL NCVPT(IDNC,IDQTRQSC,ISTART,ICOUNT,QTRQSC,IERR) C Muller matrix (IDL handles only 7 dimensions) - so 16 ncvpt's ICOUNT(6) = NSCAT CALL NCVPT(IDNC,IDSM11,ISTART,ICOUNT,SM(1,1,1),IERR) CALL NCVPT(IDNC,IDSM12,ISTART,ICOUNT,SM(1,1,2),IERR) CALL NCVPT(IDNC,IDSM13,ISTART,ICOUNT,SM(1,1,3),IERR) CALL NCVPT(IDNC,IDSM14,ISTART,ICOUNT,SM(1,1,4),IERR) CALL NCVPT(IDNC,IDSM21,ISTART,ICOUNT,SM(1,2,1),IERR) CALL NCVPT(IDNC,IDSM22,ISTART,ICOUNT,SM(1,2,2),IERR) CALL NCVPT(IDNC,IDSM23,ISTART,ICOUNT,SM(1,2,3),IERR) CALL NCVPT(IDNC,IDSM24,ISTART,ICOUNT,SM(1,2,4),IERR) CALL NCVPT(IDNC,IDSM31,ISTART,ICOUNT,SM(1,3,1),IERR) CALL NCVPT(IDNC,IDSM32,ISTART,ICOUNT,SM(1,3,2),IERR) CALL NCVPT(IDNC,IDSM33,ISTART,ICOUNT,SM(1,3,3),IERR) CALL NCVPT(IDNC,IDSM34,ISTART,ICOUNT,SM(1,3,4),IERR) CALL NCVPT(IDNC,IDSM41,ISTART,ICOUNT,SM(1,4,1),IERR) CALL NCVPT(IDNC,IDSM42,ISTART,ICOUNT,SM(1,4,2),IERR) CALL NCVPT(IDNC,IDSM43,ISTART,ICOUNT,SM(1,4,3),IERR) CALL NCVPT(IDNC,IDSM44,ISTART,ICOUNT,SM(1,4,4),IERR) C endif "if(iori. ne. 0 .and. cnetflag.eq.'ALLCDF') then" END IF C Here only for iori=0. Write smori to netCDF: IF (IORI.EQ.0) THEN ISTART(1) = IWAV ISTART(2) = IRAD ISTART(3) = 1 ICOUNT(1) = 1 ICOUNT(2) = 1 ICOUNT(3) = NSCAT CALL NCVPT(IDNC,IDSMORI11,ISTART,ICOUNT,SMORI(1,1,1),IERR) CALL NCVPT(IDNC,IDSMORI12,ISTART,ICOUNT,SMORI(1,1,2),IERR) CALL NCVPT(IDNC,IDSMORI13,ISTART,ICOUNT,SMORI(1,1,3),IERR) CALL NCVPT(IDNC,IDSMORI14,ISTART,ICOUNT,SMORI(1,1,4),IERR) CALL NCVPT(IDNC,IDSMORI21,ISTART,ICOUNT,SMORI(1,2,1),IERR) CALL NCVPT(IDNC,IDSMORI22,ISTART,ICOUNT,SMORI(1,2,2),IERR) CALL NCVPT(IDNC,IDSMORI23,ISTART,ICOUNT,SMORI(1,2,3),IERR) CALL NCVPT(IDNC,IDSMORI24,ISTART,ICOUNT,SMORI(1,2,4),IERR) CALL NCVPT(IDNC,IDSMORI31,ISTART,ICOUNT,SMORI(1,3,1),IERR) CALL NCVPT(IDNC,IDSMORI32,ISTART,ICOUNT,SMORI(1,3,2),IERR) CALL NCVPT(IDNC,IDSMORI33,ISTART,ICOUNT,SMORI(1,3,3),IERR) CALL NCVPT(IDNC,IDSMORI34,ISTART,ICOUNT,SMORI(1,3,4),IERR) CALL NCVPT(IDNC,IDSMORI41,ISTART,ICOUNT,SMORI(1,4,1),IERR) CALL NCVPT(IDNC,IDSMORI42,ISTART,ICOUNT,SMORI(1,4,2),IERR) CALL NCVPT(IDNC,IDSMORI43,ISTART,ICOUNT,SMORI(1,4,3),IERR) CALL NCVPT(IDNC,IDSMORI44,ISTART,ICOUNT,SMORI(1,4,4),IERR) c endif "if(iori.eq.0) then END IF RETURN END