      PROGRAM CON2XYZ
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER(KIND=4) (I-N)
      PARAMETER(NLONMAX=8640)
!
!...PURPOSE: CONVERTS CSR CONTOUR FILE TO GMT XYZ FORMAT FILE (LONG, LAT, VARIABLE)
!            - Extract data within user specified lat/lon bounds
!
!...NOTE :
!         (1) THE MAX. # OF DATA IN ONE ROW IS 8640
!
      CHARACTER*80 TITLE1(2), FORM1, FILE_CSRCON, FILE_GMTXYZ
!
      REAL MNLAT1,MXLAT1
      REAL MNLON1,MXLON1
      REAL MASK11,MASK12
      LOGICAL CENTER
      DATA XLATMN / -90.E0 /, XLATMX / 90.E0 /, XLONMN / 0.E0 /, XLONMX / 360.E0 /
      DATA CENTER /.TRUE./, SCALE /1.E0/
!
      DIMENSION DATA1(NLONMAX)
!
      NAMELIST /INPUT/ XLATMN,XLATMX,XLONMN,XLONMX,CENTER,SCALE

      NARGS=IARGC()
      IF(NARGS.NE.2) THEN
        PRINT*,' Usage: con2xyz CSRCON GMTXYZ < INPUT '
        PRINT*,'  Number of arguments must be 2 '
        PRINT*,' NAMELIST /INPUT/ '
        PRINT*,'    XLATMX = Maximum latitude to extract '
        PRINT*,'    XLATMN = Minimum latitude to extract '
        PRINT*,'    XLONMX = Maximum longitude to extract '
        PRINT*,'    XLONMN = Minimum longitude to extract '
        PRINT*,'    CENTER = (LOGICAL) convention - center (default) or corner '
        PRINT*,'    SCALE = Scale to apply to the data '
        STOP
      ENDIF
      CALL GETARG(1_4,FILE_CSRCON)
      CALL GETARG(2_4,FILE_GMTXYZ)

      READ(UNIT=5,NML=INPUT)

      AE=6378136.3
      PI=4.D0*DATAN(1.D0)
      DTOR=PI/180.D0
      RTOD=1.D0/DTOR

!...OPEN FILES
      OPEN (1_4, FILE=FILE_CSRCON, FORM='FORMATTED', STATUS='OLD')
      OPEN (2_4, FILE=FILE_GMTXYZ, FORM='FORMATTED' )
 
!
!...READ THE HEADERS OF THE CONTOUR FILE
      READ(1,'(A80)') TITLE1(1)
      READ(1,'(A80)') TITLE1(2)
      READ(1,100) NLAT1,NLON1
      READ(1,200) MNLAT1,MXLAT1
      READ(1,200) MNLON1,MXLON1
      READ(1,200) MASK11,MASK12
      READ(1,'(A80)') FORM1
!
      IF(NLON1.GT.NLONMAX) THEN
        PRINT*,' NLON1 GREATER THAN ',NLONMAX
        STOP
      ENDIF
!     
!...READ IN CONTOUR FILE     
      IF(CENTER) THEN
        DLAT = (MXLAT1 - MNLAT1)/FLOAT(NLAT1)
        DLON = (MXLON1 - MNLON1)/FLOAT(NLON1)
      ELSE
        DLAT = (MXLAT1 - MNLAT1)/FLOAT(NLAT1-1)
        DLON = (MXLON1 - MNLON1)/FLOAT(NLON1-1)
      ENDIF
      DATMIN=1.D99
      DATMAX=-1.D99
      DATAVE=0.
      DATRMS=0.
      DATSTD=0.
      TOTAREA = 0.D0
      WDATAVE=0.D0
      WDATRMS=0.D0
      WDATSTD=0.D0
      NPTS=0
      DO 600 I=1,NLAT1
        READ(1,FORM1) (DATA1(J), J=1,NLON1)
        DO 500 J=1,NLON1
          IF(CENTER) THEN
            RLAT = MNLAT1 + FLOAT(I-1)*DLAT + DLAT/2
            RLON = MNLON1 + FLOAT(J-1)*DLON + DLON/2
          ELSE
            RLAT = MNLAT1 + FLOAT(I-1)*DLAT 
            RLON = MNLON1 + FLOAT(J-1)*DLON 
          ENDIF
          COSPHI=COS(RLAT*DTOR)
          IF((DATA1(J).EQ.MASK11) .OR. (DATA1(J).EQ.MASK12)) GO TO 500
          IF( (RLAT.GE.XLATMN .AND. RLAT.LE.XLATMX) .AND.  &
     &        (RLON.GE.XLONMN .AND. RLON.LT.XLONMX)     )  THEN
                NPTS=NPTS+1
                DATA1J=DATA1(J)*SCALE

                DATMIN=MIN(DATMIN,DATA1J)
                DATMAX=MAX(DATMAX,DATA1J)

                DATAVE=DATAVE+DATA1J
                DATRMS=DATRMS+DATA1J**2

                AREA=(AE*DLON*DTOR)*(AE*DLAT*DTOR)*COSPHI
                TOTAREA=TOTAREA+AREA
                WDATAVE=WDATAVE+AREA*DATA1J
                WDATRMS=WDATRMS+(AREA*DATA1J)**2
                WRITE(2,300) RLON,RLAT,DATA1J
          ENDIF
 500   CONTINUE
 600   CONTINUE
       IF(NPTS.NE.0) THEN
         DATAVE=DATAVE/NPTS
         DATSTD=DATRMS-NPTS*DATAVE*DATAVE
         DATRMS=SQRT(DATRMS/NPTS)
         DATSTD=SQRT(DATSTD/NPTS)

!        WRITE(0,*) ' TOTAL AREA = ',TOTAREA
         WDATAVE=WDATAVE/TOTAREA
         WDATSTD=WDATRMS-TOTAREA*WDATAVE*WDATAVE
         WDATRMS=SQRT(WDATRMS/TOTAREA)
         WDATSTD=SQRT(WDATSTD/TOTAREA)
       ENDIF
      PRINT 501, TRIM(FILE_CSRCON),NPTS,DATAVE,DATSTD,DATMAX,DATMIN
 501  FORMAT(' N AVE STD MAX MIN : ',(A),1X,I6,4(1X,E11.5))
      PRINT 502, TRIM(FILE_CSRCON),NPTS,(WDATAVE*TOTAREA),WDATSTD,DATMAX,DATMIN
 502  FORMAT(' N WTDTOTAL WSTD MAX MIN : ',(A),1X,I6,4(1X,E11.5))
!
!...FORMAT
 100  FORMAT(16X,I5,16X,I5)
 200  FORMAT(16X,F9.4,16X,F9.4)
!300  FORMAT(3F9.3)
 300  FORMAT(2(1X,F10.4),1X,E14.6)
!...END OF GRIDIF
      STOP
      END
