PROGRAM CALLTARGET 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/1/1 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 end history C C Copyright (C) 1993,1994,1995,1996,1997 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),SHPAR(6) DOUBLE PRECISION DX,DY,DZ,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' READ(*,*)CSHAPE IF(CSHAPE.EQ.'TETRAH'.OR.CSHAPE.EQ.'tetrah')THEN CSHAPE='TETRAH' DO 1100 J=1,100 WRITE(IDVOUT,*)'Enter tetrahedron side length ', & ' (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(*,6000)SHPAR(1),NAT 1100 CONTINUE STOP ELSEIF(CSHAPE.EQ.'SPHERE'.OR.CSHAPE.EQ.'sphere')THEN CSHAPE='SPHERE' DO 1200 J=1,100 WRITE(IDVOUT,*)'Enter sphere diameter (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) RGYR2=0. XCM=0. YCM=0. ZCM=0. DO 1150 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 1150 CONTINUE 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 1200 CONTINUE STOP ELSEIF(CSHAPE.EQ.'ELLIPS'.OR.CSHAPE.EQ.'ellips')THEN CSHAPE='ELLIPS' DO 1300 J=1,100 WRITE(IDVOUT,*)'Enter ellipsoid diameters xv,yv,zv', & ' (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) XCM=0. YCM=0. ZCM=0. DO 1250 JJ=1,NAT X=IXYZ(JJ,1) Y=IXYZ(JJ,2) Z=IXYZ(JJ,3) XCM=XCM+X YCM=YCM+Y ZCM=ZCM+Z 1250 CONTINUE XCM=XCM/REAL(NAT) YCM=YCM/REAL(NAT) ZCM=ZCM/REAL(NAT) X2=0. Y2=0. Z2=0. DO 1260 JJ=1,NAT DX=IXYZ(JJ,1)-XCM DY=IXYZ(JJ,2)-YCM DZ=IXYZ(JJ,3)-ZCM X2=X2+DX*DX Y2=Y2+DY*DY Z2=Z2+DZ*DZ 1260 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 1300 CONTINUE STOP ELSEIF(CSHAPE.EQ.'CYLNDR'.OR.CSHAPE.EQ.'cylndr')THEN CSHAPE='CYLNDR' DO 1400 J=1,100 WRITE(IDVOUT,*)'Enter cylinder length (0 to stop)' READ(*,*)SHPAR(1) IF(SHPAR(1).LE.0.)STOP WRITE(IDVOUT,*)'Enter cylinder diameter' 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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 1400 CONTINUE STOP ELSEIF(CSHAPE.EQ.'RCTNGL'.OR.CSHAPE.EQ.'rctngl')THEN CSHAPE='RCTNGL' DO 1500 J=1,100 WRITE(IDVOUT,*)'Enter x-length (0 to stop)' READ(*,*)SHPAR(1) IF(SHPAR(1).LE.0.)STOP WRITE(IDVOUT,*)'Enter y-length' READ(*,*)SHPAR(2) WRITE(IDVOUT,*)'Enter z-length' READ(*,*)SHPAR(3) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 1500 CONTINUE ELSEIF(CSHAPE.EQ.'HEXGON'.OR.CSHAPE.EQ.'hexgon')THEN CSHAPE='HEXGON' DO 1600 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' READ(*,*)SHPAR(2) CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 1600 CONTINUE ELSEIF(CSHAPE.EQ.'TWOELL'.OR.CSHAPE.EQ.'twoell')THEN CSHAPE='TWOELL' DO 1700 J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 1700 CONTINUE ELSEIF(CSHAPE.EQ.'TWOAEL'.OR.CSHAPE.EQ.'twoael')THEN CSHAPE='TWOAEL' DO 1800 J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 1800 CONTINUE ELSEIF(CSHAPE.EQ.'THRELL'.OR.CSHAPE.EQ.'threll')THEN CSHAPE='THRELL' DO 1900 J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 1900 CONTINUE ELSEIF(CSHAPE.EQ.'THRAEL'.OR.CSHAPE.EQ.'thrael')THEN CSHAPE='THRAEL' DO 2000 J=1,100 WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z', & ' directions (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,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3) CALL SIZER(MXNAT,NAT,IXYZ,LXYZ) WRITE(IDVOUT,*)'NAT=',NAT 2000 CONTINUE ELSEIF(CSHAPE.EQ.'CONELL'.OR.CSHAPE.EQ.'conell')THEN CSHAPE='CONELL' DO 2100 J=1,100 WRITE(IDVOUT,*)'Enter length of outer ellipsoid in x,y,z', & 'directions (lattice units) (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 (lattice units)' 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,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 2100 CONTINUE ELSEIF(CSHAPE.EQ.'BLOCKS'.OR.CSHAPE.EQ.'blocks')THEN CSHAPE='BLOCKS' CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR, & MXNAT,SHPAR,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,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,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 and diameter b_1 of 1st spheroid' READ(*,*)SHPAR(1),SHPAR(2) WRITE(0,*)'Enter length a_2 and diameter b_2 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,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 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*(*) 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_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_null C C This version of timeit is a dummy which does not provide any C timing information. C 94/06/20 (PJF) add dtime to the formal parameters C 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 Arguments: CHARACTER CMSGTM*(*) C Local variables: CHARACTER CSTA*3,CMSGNM*70 C*********************************************************************** DATA CSTA/'ONE'/ IF(CSTA.EQ.'ONE')THEN CSTA='TWO' ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' WRITE(CMSGNM,FMT='(A,A)')'Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A)') & 'Dummy timing routine: no timings available' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END PROGRAM TSTFFT C Parameters: INTEGER MX,MY,MZ PARAMETER (MX=100, MY=100, MZ=100) C Local variables: INTEGER IN,N,NX,NY,NZ INTEGER NF235(45) COMPLEX 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 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 Copyright (C) 1994,1995 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/ PRINT *,"Calculation time for different 3-D FFT's methods" PRINT *,"================================================ " PRINT *," " 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=_______',/, & ' parameter LVR =_____',/, & ' LVR="length of vector registers" in gpfa2f,gpfa3f,gpfa5f',/, & '==========================================================',/, & ' Brenner Temperton ', & ' Temperton',/, & ' NX NY NZ (FOURX) (Old) ', & ' (GPFA)') rewind(20) DO 10 IN=1,31 N=NF235(IN) NX=N NY=N NZ=N CALL CASE(C,CW,NX,NY,NZ) 10 CONTINUE STOP END SUBROUTINE CASE(C,CW,NX,NY,NZ) C Arguments: INTEGER NX,NY,NZ COMPLEX C(NX,NY,NZ),CW(NX,NY,NZ) C Local variables: 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 C Copyright (C) 1994, 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(c,nx,ny,nz) CALL TIMEIT('BREN',T0) CALL FOURX(C,NN,3,1) CALL TIMEIT('BREN',T1) c1=c(1,1,1) c2=c(2,1,1) C--------- TEMP-OLD call gen(c,nx,ny,nz) CALL TIMEIT('OLD',T0) CALL CXFFT3(C,NX,NY,NZ,1,CW) CALL TIMEIT('OLD',T2) c3=c(1,1,1) c4=c(2,1,1) C--------- TEMP-NEW call gen(c,nx,ny,nz) CALL TIMEIT('NEW',T0) CALL CXFFT3n (C,NX,NY,NZ,1) CALL TIMEIT('NEW',T3) c5=c(1,1,1) c6=c(2,1,1) C--------- Print WRITE(21, '(1x,3i4,3F15.6)') 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,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 (a) set "include" statement in this routine to C proper location of "netcdf.inc" (b) set "ddscat.par" C variable "cnetflag" to 'ALLCDF' or 'ORICDF', (c) modify Makefile to link C with NetCDF library, e.g. (-L/netcdf/location -lnetcdf) 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 end history C Copyright (C) 1996 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),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,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,IDTHETA,IDTHETAN, + IDTHTMID,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) 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,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