C*********************************************************
C***  Program to o retrieve data from MEPTV data-base  ***
C***          By CHO & PEK                             ***
C***          Modified by MORIWAKI (2001.6.19)         ***
C*********************************************************
C      MMAX: Number of the METPV observation site.  
C      JMAX: Data species.
c            J=1  --> global solar radiation
c            J=2  --> direct solar radiation
c            J=3  --> sky solar radiation
c            J=4  --> duration of sunshine
c            J=5  --> temperature
c            J=6  --> wind direction (16 compass points )
c            J=7  --> wind velocity
c            J=8  --> precipitation
c            J=9  --> amount of snow cover
C********************************************************
       IMPLICIT REAL(A-H,O-Z)
       PARAMETER(MMAX=150,JMAX=9) 
       DIMENSION  MDAT(MMAX,12,31,JMAX,24) ,RDAT(MMAX,12,31,JMAX,24)
       DIMENSION  IBNK(MMAX),INUB(MMAX),IH(MMAX)
     1           ,TLAT(MMAX),TLON(MMAX)
       DIMENSION VECU(MMAX,12,31,JMAX,1),VECV(MMAX,12,31,JMAX,1)
       DIMENSION IDAY(12) 
       DATA IDAY/31,28,31,30,31,30,31,31,30,31,30,31/
       CHARACTER FNAME*10,FILEN*11,ISTA(MMAX)*3
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---------------------------------------------------------
C**  OPEN Station list &  
c
C---------------------------------------------------------
            write(*,*) 'OPEN Station list'
          OPEN(30,FILE='METPV_station_list',STATUS='UNKNOWN')
          DO M=1,MMAX
             READ(30,999) ISTA(M),INUB(M),TLAT(M),TLON(M),IH(M)
          ENDDO
C
 999     FORMAT(A3,I5,2F10.4,I5)
          CLOSE(30)
C
C---------------------------------------------------------
C**  OPEN FILE  &  READ DATA  **
C---------------------------------------------------------
C
            write(*,*) 'OPEN DATA file'
      DO M=1,MMAX
            FILEN='MEANHOR.'//ISTA(M)
c            write(*,*) FILEN
            OPEN(10,FILE=FILEN,STATUS='OLD')
        READ(10,722)
C
        DO IM = 1,12
          DO ID=1,IDAY(IM)
            DO J=1,JMAX
              READ(10,723) MPARA,MMON,MDAY
     1       ,MHEIGHT,(MDAT(M,IM,ID,J,K),K=1,24)
     1       ,MMAXI,MMINI
     1       ,MACUM,MMEAN,MDATE
            ENDDO
          ENDDO
        ENDDO
       CLOSE(10)
      ENDDO
 722  FORMAT(51X)
 723  FORMAT(I5,I3,I2,I5,29I5)
c
C********************************************************************
c    Convert METPV value into Physical value
C********************************************************************
            write(*,*) 'Convert DATA'
            write(*,*) '  Radiation'
C Radiations
C J = 1..3
      DO M=1,MMAX
       DO IM=1,12
        DO ID=1,IDAY(IM)
         DO K=1,24
          DO J=1,3
             IRMK=INT(MOD(MDAT(M,IM,ID,J,K),10))
             MDAT(M,IM,ID,J,K)=(MDAT(M,IM,ID,J,K)-IRMK)/10
             RDAT(M,IM,ID,J,K)=MDAT(M,IM,ID,J,K)*0.01*1000000/3600
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
            write(*,*) '  duretaion and Temperature'
C Duration of sunshine & Temperature
C J = 4..5
      DO M=1,MMAX
       DO IM=1,12
        DO ID=1,IDAY(IM)
         DO K=1,24
          DO J=4,5
             IRMK=INT(MOD(MDAT(M,IM,ID,J,K),10))
             MDAT(M,IM,ID,J,K)=(MDAT(M,IM,ID,J,K)-IRMK)/10
             RDAT(M,IM,ID,J,K)=MDAT(M,IM,ID,J,K)*0.1
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
            write(*,*) '  Wind direction'
C Wind direction (16 compass points ) 
C J = 6
      DO M=1,MMAX
       DO IM=1,12
        DO ID=1,IDAY(IM)
         DO K=1,24
          DO J=6,6
             IRMK=INT(MOD(MDAT(M,IM,ID,J,K),10))
             MDAT(M,IM,ID,J,K)=(MDAT(M,IM,ID,J,K)-IRMK)/10
             RDAT(M,IM,ID,J,K)=MDAT(M,IM,ID,J,K)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
            write(*,*) '  Wind velocity'
c Wind velocity
C J = 7
      DO M=1,MMAX
       DO IM=1,12
        DO ID=1,IDAY(IM)
         DO K=1,24
          DO J=7,7
             IRMK=INT(MOD(MDAT(M,IM,ID,J,K),10))
             MDAT(M,IM,ID,J,K)=(MDAT(M,IM,ID,J,K)-IRMK)/10
             RDAT(M,IM,ID,J,K)=MDAT(M,IM,ID,J,K)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C
            write(*,*) '  Calculate vector'
c Vector of wind
      DO M=1,MMAX
       DO IM=1,12
        DO ID=1,IDAY(IM)
         DO K=1,24
            THETA = FLOAT(360/16)*RDAT(M,IM,ID,6,K)
            THETA = 90.0-THETA
            RAD   = 6.82/360. *THETA
            VECU(M,IM,ID,1,K) = -1.* COS(RAD) * RDAT(M,IM,ID,7,K)
            VECV(M,IM,ID,1,K) = -1.* SIN(RAD) * RDAT(M,IM,ID,7,K)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C
            write(*,*) '  Precipitation  &  Snow cover'
c Precipitation  &  Snow cover
C J = 8..9
      DO M=1,MMAX
       DO IM=1,12
        DO ID=1,IDAY(IM)
         DO K=1,24
          DO J=8,9
             IRMK=INT(MOD(MDAT(M,IM,ID,J,K),10))
             MDAT(M,IM,ID,J,K)=(MDAT(M,IM,ID,J,K)-IRMK)/10
             RDAT(M,IM,ID,J,K)=MDAT(M,IM,ID,J,K)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C
c--------------------------------------------
c     Judge whether the point is inside of the selected area
c    If the point is inside,  set flag IBNK =1
c    If the point is outside, set flag IBNK =0
c--------------------------------------------
          DO M=1,MMAX
                   IBNK(M) = 0
           IF((TLAT(M).GE.CLATSW).AND.(TLAT(M).LE.CLATNE).AND.
     +      (TLON(M).GE.CLONSW).AND.(TLON(M).LE.CLONNE)) THEN
                    IBNK(M) = 1
           ENDIF
         ENDDO
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C    Data output
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            write(*,*) 'Data output'
C
      OPEN(UNIT=20,FILE='METPV-OUT',STATUS='UNKNOWN')
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
                     IM = MONTH
      DO 850 M=1,MMAX
       IF (IBNK(M).EQ.1) THEN
            write(*,*) '      Data output,',INUB(M)
         WRITE(20,710) INUB(M),TLON(M),TLAT(M)
         DO 800  ID = IDYSTA,IDYEND
                            KK1 = 1
                            KK2 = 24
                            IF (ID.EQ.IDYSTA) KK1 = MHRSTA
                            IF (ID.EQ.IDYEND) KK2 = MHREND
            DO 800 K = KK1,KK2
              WRITE(20,730) 
     1        RDAT(M,IM,ID,1,K),RDAT(M,IM,ID,5,K),
     1        RDAT(M,IM,ID,6,K),RDAT(M,IM,ID,7,K),
     1        RDAT(M,IM,ID,8,K),
     1        VECU(M,IM,ID,1,K),VECV(M,IM,ID,1,K)
 800    CONTINUE
       ENDIF
 850  CONTINUE
  710   FORMAT(I5,2F12.4)
  730  FORMAT(7F12.6)
C
       CLOSE(30)
       STOP
       END
c
C
c     J=1  --> global solar radiation   (W/m2)
c     J=2  --> direct solar radiation   (W/m2)
c     J=3  --> sky solar radiation      (W/m2)
c     J=4  --> duration of sunshine     (hour)
c     J=5  --> temperature              (degC)
c     J=6  --> wind direction (16 compass points ) 
c     J=7  --> wind velocity            (m/s)
c     J=8  --> precipitation            (mm)
c     J=9  --> amount of snow cover     (cm)
