ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c    Read SDP data and output   by Pek                                c
c                 Modified by MORIWAKI (2001.6.19)                    c
c                                                                     c
c      MMAX: Number of the observation site.                          c
c      JMAX: Number of data in a day.                                 c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c-----------------------------------------------
c    Parameter shuld be choiced by user
c-----------------------------------------------
!         PARAMETER(MMAX=155,JMAX=8)    ! if the year is before 1990
       PARAMETER(MMAX=155,JMAX=24)   ! if the year is after 1991
c
c-----------------------------------------------
c    Format shuld be choiced by user
c-----------------------------------------------
c---   before 1990  ----
!  52    FORMAT(I3,I1,10x,I4,3I2,8(I6,6x,I5,4x,I4,I3,I4,11x),152x)
c---   after 1991   ----
  52    FORMAT(I3,I1,10x,I4,3I2,24(I6,6x,I5,4x,I4,I3,I4,24x),152x)
c----------------------------------------------------------------------
               IMPLICIT REAL*4 (A-H,O-Z)
          DIMENSION  KNUB(MMAX)
     &              ,KMON(12)
     &              ,KDN(12),KDAY(12,31)
     &              ,ISTA(MMAX),XLAT(MMAX),XLON(MMAX),IZX(MMAX)
     &              ,KDATP(MMAX,12,31,JMAX)
     &              ,KDATT(MMAX,12,31,JMAX)
     &              ,KDATH(MMAX,12,31,JMAX)
     &              ,KDATD(MMAX,12,31,JMAX)
     &              ,KDATW(MMAX,12,31,JMAX)
     &              ,DATP(MMAX,12,31,JMAX)
     &              ,DATT(MMAX,12,31,JMAX)
     &              ,DATH(MMAX,12,31,JMAX)
     &              ,DATD(MMAX,12,31,JMAX)
     &              ,DATW(MMAX,12,31,JMAX)
     &              ,VECU(MMAX,12,31,JMAX)
     &              ,VECV(MMAX,12,31,JMAX)
     &              ,IBNK(MMAX)
C
          CHARACTER FILEN*11
          CHARACTER FN*4,FNB*4
c-----------------------------------------------
c    Parameter shuld be set by user
c-----------------------------------------------
      MONTH  = 8
      IDYSTA = 5
      IDYEND = 6
      MHRSTA = 21
      MHREND = 21
      CLATSW = 35.0
      CLONSW = 134.0
      CLATNE = 36.0
      CLONNE = 135.0
c-----------------------------------------------
      IF (JMAX.EQ.8) THEN
        IF((MOD(MHRSTA,3).NE.0).or.(MOD(MHREND,3).NE.0)) THEN
          WRITE(*,*) 'Error !'
          WRITE(*,*) 'Please set MHRSTA and MHREND multiple of 3.'
          STOP
        ENDIF
c
        MHRSTA = INT(MHRSTA/3)
        MHREND = INT(MHREND/3)
      ENDIF
c
ccccccccccccccccccccccccccccccccccccccccc
c Open site information file.
c      ISTA: Station number
c      XLAT: Latitude , XLON: Longitude
c      IZX : Elevation 
cccccccccccccccccccccccccccccccccccccccccc
c
          OPEN(10,FILE='sdp_station_list.txt',
     1                        STATUS='OLD',FORM='FORMATTED')
          DO M=1,MMAX
             READ(10,80) ISTA(M),XLAT(M),XLON(M),IZX(M)
          ENDDO
          CLOSE(10,STATUS='KEEP')
 80   FORMAT(I3,2F10.4,I5) 
c
cccccccccccccccccccccccccccccccccccccccccc
C   Read Data
cccccccccccccccccccccccccccccccccccccccccc
       WRITE(*,*) 'Reading SDP data'
c---------------------------------------------------------
              FILEN='SDP1991.TXT'  
c---------------------------------------------------------
              WRITE(*,*) FILEN
            OPEN(20,FILE=FILEN,STATUS='OLD')
c
       FN=FILEN(4:8)
       FNB='1991'
       IF(((FN.GE.FNB).AND.(JMAX.EQ.8)).or.
     +     ((FN.LT.FNB).AND.(JMAX.EQ.24))) THEN
         WRITE(*,*) ' Error!'
         WRITE(*,*) ' File and Parameter are not mach.'
         STOP
       ENDIF
C
cccccccccccccccccccccccccccccccccccccccccc
C  INPUT INITIAL VALUE
c      M     : Oberevation Site
c      IM    : Month
c      ID    : Day
c      J     : Time series
cccccccccccccccccccccccccccccccccccccccccc
C
           DO 500     M = 1,MMAX
           DO 500    IM = 1,12
           DO 500    ID = 1,31
           DO 500     J = 1,JMAX
C
           KDATP(M,IM,ID,J) = -999.
           KDATT(M,IM,ID,J) = -999.
           KDATH(M,IM,ID,J) = -999.
           KDATD(M,IM,ID,J) = -999.
           KDATW(M,IM,ID,J) = -999.
  500     CONTINUE
C
cccccccccccccccccccccccccccccccccccccccccc
c  KDATP : Pressure
c  KDATT : Temperature
c  KDATH : Relative Humidity
c  KDATD : Wind direction
c  KDATW : Wind velocity
cccccccccccccccccccccccccccccccccccccccccc
c
           DO 600     M = 1,MMAX
           DO 600    IM = 1,12
                     ID = 1
C
         READ(20,52,END=700) KNUB(M),X,KYEAR,
     &     KMON(IM),KDAY(IM,ID),KDN(IM),
     &     ((KDATP(M,IM,ID,J),KDATT(M,IM,ID,J),KDATH(M,IM,ID,J),
     &       KDATD(M,IM,ID,J),KDATW(M,IM,ID,J)),J=1,JMAX)
C
           DO 600    ID = 2,KDN(IM)
c
         READ(20,52,END=700) KNUB(M),X,KYEAR,
     &     KMON(IM),KDAY(IM,ID),KDN(IM),
     &     ((KDATP(M,IM,ID,J),KDATT(M,IM,ID,J),KDATH(M,IM,ID,J),
     &       KDATD(M,IM,ID,J),KDATW(M,IM,ID,J)),J=1,JMAX)
C
  600    CONTINUE
  700    CONTINUE
C        
C
C
C
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c    Transform SDP value to Physical value
c  &    Error checking
ccccccccccccccccccccccccccccccccccccccccccccccccccc
C
           DO 850     M = 1,MMAX
           DO 800    IM = 1,12
           DO 800    ID = 2,KDN(IM)
           DO 800     J = 1,JMAX
c
c-------------------------------------------------------c
c If the 1st order of data is not 8, the data is error. c
c When the data is error, input -999.                   c
c-------------------------------------------------------c
c-- Pressure --
          IRMK=INT(MOD(KDATP(M,IM,ID,J),10))
          IF ((IRMK.NE.8).or.(KDATP(M,IM,ID,J).EQ.-999))  THEN
            KDATP(M,IM,ID,J) = -999
          ELSE
            DATP(M,IM,ID,J) = FLOAT(KDATP(M,IM,ID,J)-IRMK)/10.*0.1
          ENDIF
c-- Temperature --
          IRMK=INT(MOD(KDATT(M,IM,ID,J),10))
          IF ((IRMK.NE.8).or.(KDATT(M,IM,ID,J).EQ.-999))  THEN
            KDATT(M,IM,ID,J) = -999
          ELSE
            DATT(M,IM,ID,J) = FLOAT(KDATT(M,IM,ID,J)-IRMK)/10.*0.1
          ENDIF
c-- Relative Humidity --
          IRMK=INT(MOD(KDATH(M,IM,ID,J),10))
          IF ((IRMK.NE.8).or.(KDATH(M,IM,ID,J).EQ.-999))  THEN
            KDATH(M,IM,ID,J) = -999
          ELSE
            DATH(M,IM,ID,J) = FLOAT(KDATH(M,IM,ID,J)-IRMK)/10.
          ENDIF
c-- Wind Speed --
          IRMK=INT(MOD(KDATW(M,IM,ID,J),10))
          IF ((IRMK.NE.8).or.(KDATW(M,IM,ID,J).EQ.-999))  THEN
            KDATW(M,IM,ID,J) = -999
          ELSE
            DATW(M,IM,ID,J) = FLOAT(KDATW(M,IM,ID,J)-IRMK)/10.*0.1
          ENDIF
c-- Wind Direction --
          IRMK=INT(MOD(KDATD(M,IM,ID,J),10))
          IF ((IRMK.NE.8).or.(KDATD(M,IM,ID,J).EQ.-999))  THEN
            KDATD(M,IM,ID,J) = -999
          ELSE
            DATD(M,IM,ID,J) = FLOAT(KDATD(M,IM,ID,J)-IRMK)/10.
          ENDIF
c
cccccccccccccccccccccccccccccccccccc
c  calculate the vector
cccccccccccccccccccccccccccccccccccc
c
          THETA = FLOAT(360/16)*DATD(M,IM,ID,J)
          THETA = 90.0-THETA
          RAD   = 6.82/360. *THETA
c
c-- U --
          VECU(M,IM,ID,J) = -1.* COS(RAD) * DATW(M,IM,ID,J)
c-- V --
          VECV(M,IM,ID,J) = -1.* SIN(RAD) * DATW(M,IM,ID,J)
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c  +----+---+---+---+---+---+---+---+---+---+
c  |DATD| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
c  +----+---+---+---+---+---+---+---+---+---+
c  |DIR | - |NNE| NE|ENE| E |ESE| SE|SSE| S |
c  +----+---+---+---+---+---+---+---+---+---+
c
c  +----+---+---+---+---+---+---+---+---+
c  |DATD| 9 | 10| 11| 12| 13| 14| 15| 16|
c  +----+---+---+---+---+---+---+---+---+
c  |DIR |SSW| SW|WSW| W |WNW| NW|NNW| N |
c  +----+---+---+---+---+---+---+---+---+
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c-- if wind data is error, input 0 to VECU & VECV
c
          IF ((KDATW(M,IM,ID,J).EQ.-999).or.(KDATD(M,IM,ID,J).EQ.-999))
     1                                                        THEN
             VECU(M,IM,ID,J) = -999.
             VECV(M,IM,ID,J) = -999.
          ENDIF
c
c-- if wind data is error, input 0 to VECU & VECV
c
          IF (KDATD(M,IM,ID,J).EQ.0)  THEN
             VECU(M,IM,ID,J) = -999.
             VECV(M,IM,ID,J) = -999.
          ENDIF
c
c
c
  800    CONTINUE
c---------------------------------------------
c  Judge whether observation point is within indicated area
c    If yes,  set flag IBNK =1
c    If no,   set flag IBNK =0
c---------------------------------------------
                   IBNK(M) = 0
         IF((XLAT(M).GE.CLATSW).AND.(XLAT(M).LE.CLATNE).AND.
     +      (XLON(M).GE.CLONSW).AND.(XLON(M).LE.CLONNE)) THEN
                   IBNK(M) = 1
         ENDIF
c
  850    CONTINUE
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C    Data output
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
        OPEN(30,FILE='SDP-OUT',STATUS='UNKNOWN',FORM='FORMATTED')
c       write(*,*) 'OK1'
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
                     IM = MONTH
           DO 900     M = 1,MMAX
       IF (IBNK(M).EQ.1) THEN
         WRITE(30,710) ISTA(M),XLON(M),XLAT(M)
               DO    ID = IDYSTA,IDYEND
                            JJ1 = 1
                            JJ2 = JMAX
                            IF (ID.EQ.IDYSTA) JJ1 = MHRSTA
                            IF (ID.EQ.IDYEND) JJ2 = MHREND
                  DO  J = JJ1,JJ2
         WRITE(30,711) DATW(M,IM,ID,J),DATD(M,IM,ID,J),
     &               VECU(M,IM,ID,J),VECV(M,IM,ID,J),
     &               DATT(M,IM,ID,J),DATH(M,IM,ID,J)
                  ENDDO
               ENDDO
       ENDIF
 900       CONTINUE
c--
 710   FORMAT(I5,2F12.4)
 711   FORMAT(6E14.6)
c--
      CLOSE(20)
      CLOSE(30)
c--
          STOP
          END

