C SAMPLE.F

C CONVERT YOUR DATA TO A McIDAS GR3D FILE.



C THE MAIN PROGRAM OF YOUR CONVERSION PROGRAM MUST
C BE NAMED SUBROUTINE MAIN0
C
      SUBROUTINE MAIN0
C
C THE NEXT TWO COMMENTS ARE PRINTED BY THE 'help sample' COMMAND
C ? SAMPLE program to convert data to 3D grid files
C ? sample gridf#
C
C DIMENSIONS OF 3D GRID
C NOTE NLATS, NLONS AND NHGTS MUST ALL BE AT LEAST 2
      PARAMETER (NLATS=31,NLONS=51,NHGTS=16)
C
C NUMBER OF PHYSICAL VARIABLES AND NUMBER OF TIME STEPS
C NOTE EITHER OR BOTH MAY BE EQUAL TO 1.  THAT IS, VIS-5D DOES
C NOT FORCE YOU TO HAVE MULTIPLE VARIABLES OR TIME DYNAMICS.
      PARAMETER (NVARS=5,NTIMES=100)
C
C ARRAY FOR 3D GRID DATA
      REAL*4 G(NLATS, NLONS, NHGTS)
C ARRAYS FOR GRID FILE ID AND GRID DIRECTORY
      INTEGER ID(8), IDIR(64)
C ARRAY FOR VARIABLE NAMES
      CHARACTER*4 CNAME(NVARS)
C
C LATITUDE, LONGITUDE AND HEIGHT BOUNDS FOR SPATIAL GRID
      DATA XLATS/20.0/,XLATN/50.0/
      DATA XLONE/70.0/,XLONW/120.0/
      DATA XHGTB/0.0/,XHGTT/15.0/
C
C STARTING DATE IN YYDDD AND TIME IN HHMMSS
      DATA JDAY/88035/,JTIME/020000/
C TIME STEP IN HHMMSS
      DATA JSTEP/000100/
C
C NAMES OF THE FIVE PHYSICAL VARIABLES
      DATA CNAME/'U   ', 'V   ', 'W   ', 'T   ', 'P   '/
C INITIALIZE GRID DIRECTORY TO ZEROS
      DATA IDIR/64*0/
C
C READ GRID FILE NUMBER FROM COMMAND LINE.  IPP WILL
C CONVERT THE PARAMETER # 1 TO AN INTEGER, WITH A DEFAULT
C VALUE OF 0.
      IGRIDF=IPP(1,0)
C IF ILLEGAL GRID FILE NUMBER, PRINT ERROR MESSAGE AND RETURN
      IF(IGRIDF .LT. 1 .OR. IGRIDF .GT. 9999) THEN
        CALL EDEST('BAD GRID FILE NUMBER ',IGRIDF)
        CALL EDEST('MUST BE BETWEEN 1 AND 9999 ',0)
        RETURN
      ENDIF
C
C CALCULATE GRID INTERVALS
      XLATIN=(XLATN-XLATS)/(NLATS-1)
      XLONIN=(XLONW-XLONE)/(NLONS-1)
      XHGTIN=(XHGTT-XHGTB)/(NHGTS-1)
C
C DATE AND TIME FOR FIRST TIME STEP
C IDAYS CONVERTS YYDDD FORMAT TO DAYS SINCE JAN. 1, 1900
      IDAY=IDAYS(JDAY)
C ISECS CONVERTS HHMMSS FORMAT TO SECONDS SINCE MIDNIGHT
      ISEC=ISECS(JTIME)
C
C INITIALIZE GRID IDENTIFIER TEXT TO BLANKS
C NOTE LIT CONVERTS A CHARACTER*4 TO AN INTEGER*4
      DO 10 I=1,8
10    ID(I)=LIT('    ')
C
C SET UP DIRECTORY ENTRY
C
C DIMENSIONS OF GRID
      IDIR(1)=NLATS*NLONS*NHGTS
      IDIR(2)=NLATS
      IDIR(3)=NLONS
      IDIR(4)=NHGTS
C
C LATITUDES AND LONGITUDES IN DEGREES * 10000
      IDIR(22)=4
      IDIR(23)=NINT(XLATN*10000.)
      IDIR(24)=NINT(XLONW*10000.)
      IDIR(25)=NINT(XLATIN*10000.0)
      IDIR(26)=NINT(XLONIN*10000.0)
C
C HEIGHTS IN METERS
      IDIR(31)=1
      IDIR(32)=NINT(XHGTT*1000.)
      IDIR(33)=NINT(XHGTIN*1000.)
C
C     CREATE THE GRID FILE
      CALL IGMK3D(IGRIDF, ID, NLATS*NLONS*NHGTS)
C
C     LOOP FOR TIME STEPS
      DO 200 IT=1,NTIMES
C
C        SET DATE AND TIME IN DIRECTORY ENTRY
C        IYYDDD CONVERTS DAYS SINCE JAN. 1, 1900 TO OUR YYDDD FORMAT
         IDIR(6)=IYYDDD(IDAY)
C        IHMS CONVERTS SECONDS SINCE MIDNIGHT TO OUR HHMMSS FORMAT
         IDIR(7)=IHMS(ISEC)
C
C        LOOP FOR PHYSICAL VARIABLES
         DO 190 IV=1,NVARS
C
C           SET VARIABLE NAME IN DIRECTORY ENTRY
            IDIR(9)=LIT(CNAME(IV))
C
C           *************************************************************
C           READ YOUR DATA FOR TIME STEP NUMBER IT AND VARIABLE NUMBER IV
C           INTO THE ARRAY G HERE.
C           NOTE THAT G(1,1,1) IS THE NORTH WEST BOTTOM CORNER AND
C           G(NLATS,NLONS,NHGTS) IS THE SOUTH EAST TOP CORNER.
C           MARK A GRID POINT AS 'MISSING DATA' BY SETTING IT = 1.0E35
C
C
C
C
C
C
C           *************************************************************
C
C           CALCULATE 3D GRID NUMBER
            IGRID=IV+NVARS*(IT-1)
C           WRITE DATA IN G AND DIRECTORY IN IDIR TO 3D GRID
C           NOTE WE PASS THE NEGATIVE OF THE GRID NUMBER (I.E. -IGRID)
            CALL IGPT3D(IGRIDF,-IGRID,G,NLATS,NLONS,NHGTS,IDIR,IGNO)
C     
C           END OF PHYSICAL VARIABLE LOOP
 190     CONTINUE
C     
C        INCREMENT DATE AND TIME, CONVERT JSTEP FROM HHMMSS TO SECONDS
         ISEC=ISEC+ISECS(JSTEP)
C        IF SECONDS CARRY PAST ONE DAY, ADJUST SECONDS AND DAYS
         IDAY=IDAY+ISEC/(24*3600)
         ISEC=MOD(ISEC,24*3600)
C
C     END OF TIME STEP LOOP
 200  CONTINUE
C
      RETURN
      END