C########1#########2#########3#########4#########5#########6#########7
C  2017-11-17   (add one-sigma level and debug of PGmin -> PGmax
C  2016-02-18   (modified from kakmin_v1.f  of  2013-01-04)
C  unix> mv kakmin.f tem.f
C  unix> tr -s '\015' '\012' < tem.f > kakmin.f
C  unix> gfortran kakmin.f
C  unix> mv a.out kakmin.out                    
C  unix> kakmin.out
C
C   after arranging March data (into February or April at 13 March)
C  unix> pgsort.out
C
C   programing note: WRITE(*) ' *** ' does not work.
C                    We have to use WRITE(6,*) ' *** ' instead!
C
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
C  kakmin.f (main) @PG/1min-raw-data/
C     1min-raw-data/kak20yymmdddaeg.min ==> (read)
C               ==> kak_10min_20yymm.csv
C
C input : minute value (Kakioka PG)
C         March data is divide into 1-10 (mm=03) and 14-31 (mm=13)
C output: every X minute value (X=10, 20, 30, 60, 90)
C         to change X, just set "intvl = 10" or "intvl = 20" etc.
C
C Required line (for 10  min):  6 x 24 x 31 = 4463 line/month
C   preface = 19 lines, title = 1 line, data = 1440 lines (=24x60)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
C 123456789-123456789-123456789-123456789-123456789-123456789-123456789-
C YYYY-MM-DD hh:mm:00.000 DOY     88888.00  88888.00  +++++.+0  88888.00
C 2011-03-17 11:13:00.000 076     88888.00  88888.00      7.30  88888.00
C 2011-03-17 11:14:00.000 076     88888.00  88888.00  99999.90  88888.00
C
C   => PGmax(144,8),
C   => PGmin(144,8),
C   YYYYMMDD hh:mm maximum %+1min +2min ... +5min -1min ... -5 min
C   YYYYMMDD hh:mm minimum %+1min +2min ... +5min -1min ... -5 min
C   as, e.g., kak_10min_20yymm.csv
C          &  kak_10min_20yymm.csv
C # NOTE-1: max (positive) before min (negative)
C # NOTE-2: decay (+1min,,,+5min) before develope (-1min,,,-5min)
C
C   comment "c1" is to take average of +-Xmin instead of +Xmin, -Xmin
C########1#########2#########3#########4#########5#########6#########7
      INTEGER JDAY(45000), YMD(45000), jmax(4500), jmin(4500), intvl
      INTEGER YYMM, YYin, MMin, DDin, NDAY, imm10
      INTEGER i, j, k, j0,j1,j2,j3,j4,j2old
      INTEGER*8 YMD0, imm
      INTEGER IDEV, ODEV, ODEV2, ITERM, OTERM
      REAL*8 PG(45000), PGmax(12, 4500), PGmin(12, 4500)
      REAL*8 PGave(2,4500), PGdev, err, err2
c      REAL*8 SDmax(12, 4500), SDmin(12, 4500)
      REAL step
      CHARACTER*11 chYMD
      CHARACTER*19 Fname
      CHARACTER*20 Fname2
      CHARACTER*5 UT(45000)
      CHARACTER*2 DUMMY
      DATA YMD/45000*0/, JDAY/45000*0/ 
     ?  , PG/45000*0./, PGmax/54000*9.99/  , PGmin/54000*9.99/ 
     ?  , PGave/9000*0./
C
C  I/O device number:  IDEV, ODEV, TERM
C
      IDEV=21                                                    
      ODEV=31
      ITERM=5
      OTERM=6
C
C  Define file
C
      intvl=10
C
C   # loop2: only works when YYMM=0 (reset loop to  = 1 for new file)
C            because loop>=2 means adding data in the same file  
C   # loop:  scan month (YYMM), and multi-month in one file
C       => YYMM=**02: first read Feb, and then read 1-10 March (**13)
C        & YYMM=**03: first read 14-31 March, and then read April
C        & YYMM=**05: read 1-31 May
C   # iday : [j3,NDAY] , j3=1 (Feb, Apr, May) or 14 (March)
C     e.g., [1,10] for YYMM=**13, and [14,31] for YYMM=**03
C   # YMD0+iday = start yyyymmdd = file name to read
C     e.g., 20110314 for YYMM=1103, and 20110301 for YYMM=1113
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  500 CONTINUE
      YYMM=0601
c      WRITE(OTERM,*) '  YYMM = 0601 / = -1 (skip, and sort data)'
c      READ(ITERM,*) YYMM
      DO 19 loop2=1,99
       IF((YYMM .LE. -1) .OR. (YYMM .GE. 1605)) THEN
         WRITE(OTERM,*) 'YYMM, loop, loop2',YYMM, loop, loop2
         GO TO 596
       ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
      DO 191 loop=1,99
       WRITE(OTERM,*) 'YYMM?  **02 or **03 or **05 or =< -1 (skip)'
       WRITE(OTERM,*) 'after second loop, three option'
       WRITE(OTERM,*) ' =0: new file with new YYMM, -1: end  '
       WRITE(OTERM,*) ' or YYMM: write on the same file (for 2005**)'
c       WRITE(OTERM,*) ' YY03: 14-31 March, YY13: 1-10 March'
       j2old = MOD(YYMM,100)
       IF(j2old .EQ. 2) THEN
         YYMM=YYMM+11
       ELSEIF(j2old .EQ. 3) THEN
         YYMM=YYMM+1
       ELSEIF(j2old .EQ. 4) THEN
         YYMM=YYMM+99
         READ(ITERM,*) YYMM
       ELSEIF(j2old .EQ. 5) THEN
         YYMM=YYMM+100
         READ(ITERM,*) YYMM
       ELSEIF(j2old .EQ. 13) THEN
         YYMM=YYMM+89
         READ(ITERM,*) YYMM
       ELSEIF(j2old .EQ. 12) THEN
         YYMM=YYMM+100
         READ(ITERM,*) YYMM
       ELSE
         YYMM=YYMM+100
         READ(ITERM,*) YYMM
       ENDIF
       IF((YYMM .LE. 600) .OR. (YYMM .GE. 1600)) THEN
         IF(loop .NE. 1) THEN
           CLOSE (ODEV)
           CLOSE (IDEV)
         ENDIF
         WRITE(OTERM,*) 'YYMM, loop, loop2',YYMM, loop, loop2
         GO TO 595
       ENDIF
       j1 = YYMM/100
       j2 = MOD(YYMM,100)
       CALL dayofmon(j1,j2, NDAY)
       j3=1
       IF(j2 .EQ. 3) j3=14
       YMD0 = (200000+YYMM)*100
       IF(j2 .EQ. 13) YMD0 = YMD0-1000
       YMD0 = YMD0+j3
       WRITE (OTERM,'(A13,i2,",",i3,"]")') '[j3,NDAY] = [',j3,NDAY
       WRITE (OTERM,'(A6,i9)') 'YMD0 =',YMD0
C########1#########2#########3#########4#########5#########6#########7
C  Open and Read the file with the format of (except header)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
C YYYY-MM-DD hh:mm:00.000 DOY     88888.00  88888.00  +++++.+0  88888.0
C
C  For every day, Header = 20 lines, and Data = 1440 lines (=24x60)
C   #  i   = minite of the day
C   #  imm = minute of the month (i.e., line number of monthly file)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  510  CONTINUE
       imm = 0
       DO 11 iday=j3,NDAY
        WRITE (chYMD,'(a3,i8)') 'kak',YMD0
        Fname=chYMD//'daef.min'
c        WRITE (OTERM,'(A7,A19)') 'open = ',Fname
        OPEN (UNIT=IDEV, FILE=Fname, STATUS='OLD')
         DO 110 i=1,20
  110      READ(IDEV,'(A2)') DUMMY
         DO 112 i=1,1440
           imm = imm+1
           READ (IDEV,415,END=519) YYin, MMin, DDin, UT(imm)
     ?       , JDAY(imm) , PG(imm)
           YMD(imm) = (((YYin*100)+MMin)*100)+DDin
           IF(YMD(imm) .NE. YMD0) GOTO 590
  112    CONTINUE
  519   CLOSE (IDEV)
        YMD0 = YMD0+1
        WRITE(OTERM,417) YMD(imm), UT(imm), PG(imm)
  11   CONTINUE
  415  FORMAT (I4,2(1x,I2),1x,a5,8x,I3,25x,F8.2) 
  417  FORMAT (I8,', ',a5,', PG =',F8.2) 
       WRITE(OTERM,*) ' Reading completed :', YYMM
C########1#########2#########3#########4#########5#########6#########7
C  Read finished (total number of PG data = imm)
C     and
C  Obtain 10-min peak values
C   #  imm   = total minutes of the month
C   #  intvl = 10 min
C   #  imm10 = imm/intvl = numbers of 10-min bins
C   #  jmax(i10), jmin(i10) = minutes registering Max/Min values of PG
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  530  CONTINUE
       imm10 = imm/intvl
       DO 130 i10=1,imm10
        PGmax(1,i10)=-9999
        PGmin(1,i10)=9999
        jmax(i10)=0
        jmin(i10)=0
  130  CONTINUE
       DO 132 i10=1,imm10
       DO 132 j=(i10*intvl)-intvl+1, i10*intvl, 1
          IF(PG(j) .LT. 9999) THEN
            IF(PG(j) .GT. PGmax(1,i10)) THEN
              PGmax(1,i10)=PG(j)
              jmax(i10)=j
            ENDIF
            IF(PG(j) .LT. PGmin(1,i10)) THEN
              PGmin(1,i10)=PG(j)
              jmin(i10)=j
            ENDIF
          ENDIF
  132  CONTINUE
       DO 134 i10=1,imm10
          IF(PGmax(1,i10) .LE. -9000) PGmax(1,i10) = 9999
          PGmax(2,i10)=PGmax(1,i10)
          PGmin(2,i10)=PGmin(1,i10)
          IF(PGmax(2,i10) .LT. 100.) PGmax(2,i10) = 9999
          IF(PGmin(2,i10) .GE. -50.) PGmin(2,i10) = 9999
  134  CONTINUE
C########1#########2#########3#########4#########5#########6#########7
C  Obtain %values at +/- 1min, 2min, 3min, 4min, 5min from the peak
C    At PGmax(3-12, i10) or PGmin(3-12, i10) for bin # i10,
C     store % values of peak when data exists
C           at +1min,,,-5min(time=j1) and at -1min,,,-5min(time=j2),
C     while assign -9.88: for end of decay or before develope
C           assign -9.99: when maximum < 100 or Minimum > 0 ==>
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
       DO 140 i10=1,imm10
       DO 140 k=1,10
          PGmax(k+2,i10)=-9.99
          PGmin(k+2,i10)=-9.99
  140  CONTINUE
C
C    When data exists i.e., when PG < 9999 =>
C     (1) ratio = PG(j1)/PGmax(1,i10), but limit to -0.01<rario<+2
C                   (-0.01 means end of decay, 2 is or error check)
C     (2) ratio(k+1) must be smaller than "ratio(k)/err"
C         => The Error range = 0.1 - 5.0 V/m (error = 1%)
C # problem is when real peak is at next/previous intvl => next section
C #                             (at i10+1 or i10-1)
C   #  imm = total minutes of the month
C             =>   1 =< j1 (=+k min from peak = decay) =< imm,
C               &  1 =< j2 (=-k min from peak = develope) =< imm,
C             =>   PG(j1) and PG(j2) exists within the monthly file
C                  (we ignore the deveple/decay across the month)
C   #  imm10 = imm/intvl = numbers of 10-min bins
C   #  i10   = scanning number of all 10-min bns.
C   #  jmax(i10), jmin(i10) = minutes registering Max/Min values of PG
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  540  CONTINUE
       err = 0.99
       DO 14 i10=1,imm10
        IF((PGmax(1,i10).GE.100.) .AND. (PGmax(1,i10).LT.9990.)) THEN
          DO 142 k=1,5
            j1=min(jmax(i10)+k,imm)
            j2=max(jmax(i10)-k,1)
            IF(PG(j1) .LT. 9999)
     ?            PGmax(k+2,i10) = max(-0.01, PG(j1)/PGmax(1,i10))
            IF(PG(j2) .LT. 9999) 
     ?            PGmax(k+7,i10) = max(-0.01, PG(j2)/PGmax(1,i10))
            IF(PGmax(k+2,i10) .GE. 2) PGmax(k+2,i10) = 2
            IF(PGmax(k+7,i10) .GE. 2) PGmax(k+7,i10) = 2
  142     CONTINUE
C
C  Extra filter for +3min, +4min, +5min (j1), -3min, -4min, -5min (j2)
C    to remove data up and downarround zero
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
          DO 143 k=1,3
            j1=min(jmax(i10)+k,imm)
            j2=max(jmax(i10)-k,1)
            j3=min(j1+2,imm)
            j4=max(j2-2,1)
            IF((PG(j1).GE.0) .AND. (PG(j3).GT.PG(j1)/err)) THEN
                  PGmax(k+4,i10) = -9.88
            ELSEIF((PG(j1).LT.0) .AND. (PG(j3).GE.0)) THEN
                  PGmax(k+4,i10) = -9.88
            ENDIF
            IF((PG(j2).GE.0) .AND. (PG(j4).GT.PG(j2)/err)) THEN
                  PGmax(k+9,i10) = -9.88
            ELSEIF((PG(j2).LT.0) .AND. (PG(j4).GE.0)) THEN
                  PGmax(k+9,i10) = -9.88
            ENDIF
  143     CONTINUE
          DO 144 k=1,4
            j1=min(jmax(i10)+k,imm)
            j2=max(jmax(i10)-k,1)
            j3=min(j1+1,imm)
            j4=max(j2-1,1)
            IF((PG(j1).GE.0) .AND. (PG(j3).GT.PG(j1)/err)) THEN
                  PGmax(k+3,i10) = -9.88
            ELSEIF((PG(j1).LT.0) .AND. (PG(j3).GE.0)) THEN
                  PGmax(k+3,i10) = -9.88
            ENDIF
            IF((PG(j2).GE.0) .AND. (PG(j4).GT.PG(j2)/err)) THEN
                  PGmax(k+8,i10) = -9.88
            ELSEIF((PG(j2).LT.0) .AND. (PG(j4).GE.0)) THEN
                  PGmax(k+8,i10) = -9.88
            ENDIF
            IF(PGmax(k+2,i10) .LT. -9.) PGmax(k+3,i10) = -9.88
            IF(PGmax(k+7,i10) .LT. -9.) PGmax(k+8,i10) = -9.88
  144     CONTINUE
        ENDIF
        IF(jmax(i10) .EQ. 0) jmax(i10)=((i10-1)*intvl)+1
C
C  The same procedure for PGmin (<0) instead of PGmax (>+100)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        IF(PGmin(1,i10) .LT. -50) THEN
          DO 146 k=1,5
            j1=min(jmin(i10)+k,imm)
            j2=max(jmin(i10)-k,1)
            IF(PG(j1) .LT. 9999) 
     ?            PGmin(k+2,i10) = max(-0.01, PG(j1)/PGmin(1,i10))
            IF(PG(j2) .LT. 9999) 
     ?            PGmin(k+7,i10) = max(-0.01, PG(j2)/PGmin(1,i10))
            IF(PGmin(k+2,i10) .GE. 2) PGmin(k+2,i10) = 2
            IF(PGmin(k+7,i10) .GE. 2) PGmin(k+7,i10) = 2
  146     CONTINUE
C
C  Extra filter for +3min, +4min, +5min (j1), -3min, -4min, -5min (j2)
C    to remove data up and downarround zero
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
          DO 147 k=1,3
            j1=min(jmin(i10)+k,imm)
            j2=max(jmin(i10)-k,1)
            j3=min(j1+2,imm)
            j4=max(j2-2,1)
            IF((PG(j1).LE.0) .AND. (PG(j3).LT.PG(j1)/err)) THEN
                  PGmin(k+4,i10) = -9.88
            ELSEIF((PG(j1).GT.0) .AND. (PG(j3).LE.0)) THEN
                  PGmin(k+4,i10) = -9.88
            ENDIF
            IF((PG(j2).LE.0) .AND. (PG(j4).LT.PG(j2)/err)) THEN
                  PGmin(k+9,i10) = -9.88
            ELSEIF((PG(j2).GT.0) .AND. (PG(j4).LE.0)) THEN
                  PGmin(k+9,i10) = -9.88
            ENDIF
  147     CONTINUE
          DO 148 k=1,4
            j1=min(jmin(i10)+k,imm)
            j2=max(jmin(i10)-k,1)
            j3=min(j1+1,imm)
            j4=max(j2-1,1)
            IF((PG(j1).LE.0) .AND. (PG(j3).LT.PG(j1)/err)) THEN
                  PGmin(k+3,i10) = -9.88
            ELSEIF((PG(j1).GT.0) .AND. (PG(j3).LE.0)) THEN
                  PGmin(k+3,i10) = -9.88
            ENDIF
            IF((PG(j2).LE.0) .AND. (PG(j4).LT.PG(j2)/err)) THEN
                  PGmin(k+8,i10) = -9.88
            ELSEIF((PG(j2).GT.0) .AND. (PG(j4).LE.0)) THEN
                  PGmin(k+8,i10) = -9.88
            ENDIF
            IF(PGmin(k+2,i10) .LT. -9.) PGmin(k+3,i10) = -9.88
            IF(PGmin(k+7,i10) .LT. -9.) PGmin(k+8,i10) = -9.88
  148     CONTINUE
        ENDIF
        IF(jmin(i10) .EQ. 0) jmin(i10)=((i10-1)*intvl)+1
   14  CONTINUE
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
C Up to here:  
C   If PGmax or PGmix cannot be defined, all PGmax(k,i) will be 9.99
C   If trend ends, PGmax(k,i) will be -9.88
C   If completely decay, PGmax(k,i) will be -0.01
C
C########1#########2#########3#########4#########5#########6#########7
C Check if the real peak in next/previous 10 min
C In that case, check how smooth the sequence is connected 
C
C   #  imm:  total number of 1-min data,
C   #  imm10 = imm/intvl: numbers of 10-min bins per month
C   #  note: re-define err = 0.9
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  550  CONTINUE
       err = 0.9
       DO 15 i10=1,imm10
        IF((PGmax(3,i10).GT.1) .OR. (PGmax(8,i10).GT.1)) THEN
            PGmax(2,i10) = -9999
          DO 151 k=1,10
            PGmax(k+2,i10) = -9.88
  151     CONTINUE
        ELSEIF(PGmax(4,i10).GT.1) THEN
c        ELSEIF((PGmax(4,i10).GT.1) .AND. (PGmax(3,i10).GT.err)) THEN
            PGmax(2,i10) = -9999
          DO 152 k=1,10
              PGmax(k+2,i10) = -9.88
  152     CONTINUE
        ELSEIF(PGmax(9,i10).GT.1) THEN
c        ELSEIF((PGmax(9,i10).GT.1) .AND. (PGmax(8,i10).GT.err)) THEN
            PGmax(2,i10) = -9999
          DO 153 k=1,10
            PGmax(k+2,i10) = -9.88
  153     CONTINUE
        ELSEIF(PGmax(5,i10) .GT. 1) THEN
          IF((PGmax(3,i10).GT.err) .AND. (PGmax(4,i10).GT.err)) THEN
              PGmax(2,i10) = -9999
            DO 154 k=1,10
              PGmax(k+2,i10) = -9.88
  154       CONTINUE
          ENDIF
        ELSEIF(PGmax(10,i10) .GT. 1) THEN
          IF((PGmax(8,i10).GT.err) .AND. (PGmax(9,i10).GT.err)) THEN
              PGmax(2,i10) = -9999
            DO 155 k=1,10
              PGmax(k+2,i10) = -9.88
  155       CONTINUE
          ENDIF
        ELSEIF(PGmax(6,i10) .GT. 1) THEN
          IF((PGmax(3,i10).GT.err) .AND. (PGmax(4,i10).GT.err)) THEN
          IF(PGmax(5,i10) .GT. err) THEN
              PGmax(2,i10) = -9999
            DO 156 k=1,10
              PGmax(k+2,i10) = -9.88
  156       CONTINUE
          ENDIF
          ENDIF
        ELSEIF(PGmax(11,i10) .GT. 1) THEN
          IF((PGmax(8,i10).GT.err) .AND. (PGmax(9,i10).GT.err)) THEN
          IF(PGmax(10,i10) .GT. err) THEN
              PGmax(2,i10) = -9999
            DO 157 k=1,10
              PGmax(k+2,i10) = -9.88
  157       CONTINUE
          ENDIF
          ENDIF
        ELSEIF(PGmax(7,i10) .GT. 1) THEN
          IF((PGmax(3,i10).GT.err) .AND. (PGmax(4,i10).GT.err)) THEN
          IF((PGmax(5,i10).GT.err) .AND. (PGmax(6,i10).GT.err)) THEN
              PGmax(2,i10) = -9999
            DO 158 k=1,10
              PGmax(k+2,i10) = -9.88
  158       CONTINUE
          ENDIF
          ENDIF
        ELSEIF(PGmax(12,i10) .GT. 1) THEN
          IF((PGmax(8,i10).GT.err) .AND. (PGmax(9,i10).GT.err)) THEN
          IF((PGmax(10,i10).GT.err) .AND. (PGmax(11,i10).GT.err)) THEN
              PGmax(2,i10) = -9999
            DO 159 k=1,10
              PGmax(k+2,i10) = -9.88
  159       CONTINUE
          ENDIF
          ENDIF
        ENDIF
   15  CONTINUE
C
C  The same procedure for PGmin (<0) instead of PGmax (>+100)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
       DO 16 i10=1,imm10
        IF((PGmin(3,i10).GT.1) .OR. (PGmin(8,i10).GT.1)) THEN
            PGmin(2,i10) = -9999
          DO 161 k=1,10
            PGmin(k+2,i10) = -9.88
  161     CONTINUE
        ELSEIF(PGmin(4,i10).GT.1) THEN
c        ELSEIF((PGmin(4,i10).GT.1) .AND. (PGmin(3,i10).GT.err)) THEN
            PGmin(2,i10) = -9999
          DO 162 k=1,10
              PGmin(k+2,i10) = -9.88
  162     CONTINUE
        ELSEIF(PGmin(9,i10).GT.1) THEN
c        ELSEIF((PGmin(9,i10).GT.1) .AND. (PGmin(8,i10).GT.err)) THEN
            PGmin(2,i10) = -9999
          DO 163 k=1,10
            PGmin(k+2,i10) = -9.88
  163     CONTINUE
        ELSEIF(PGmin(5,i10) .GT. 1) THEN
          IF((PGmin(3,i10).GT.err) .AND. (PGmin(4,i10).GT.err)) THEN
              PGmin(2,i10) = -9999
            DO 164 k=1,10
              PGmin(k+2,i10) = -9.88
  164       CONTINUE
          ENDIF
        ELSEIF(PGmin(10,i10) .GT. 1) THEN
          IF((PGmin(8,i10).GT.err) .AND. (PGmin(9,i10).GT.err)) THEN
              PGmin(2,i10) = -9999
            DO 165 k=1,10
              PGmin(k+2,i10) = -9.88
  165       CONTINUE
          ENDIF
        ELSEIF(PGmin(6,i10) .GT. 1) THEN
          IF((PGmin(3,i10).GT.err) .AND. (PGmin(4,i10).GT.err)) THEN
          IF(PGmin(5,i10) .GT. err) THEN
              PGmin(2,i10) = -9999
            DO 166 k=1,10
              PGmin(k+2,i10) = -9.88
  166       CONTINUE
          ENDIF
          ENDIF
        ELSEIF(PGmin(11,i10) .GT. 1) THEN
          IF((PGmin(8,i10).GT.err) .AND. (PGmin(9,i10).GT.err)) THEN
          IF(PGmin(10,i10) .GT. err) THEN
              PGmin(2,i10) = -9999
            DO 167 k=1,10
              PGmin(k+2,i10) = -9.88
  167       CONTINUE
          ENDIF
          ENDIF
        ELSEIF(PGmin(7,i10) .GT. 1) THEN
          IF((PGmin(3,i10).GT.err) .AND. (PGmin(4,i10).GT.err)) THEN
          IF((PGmin(5,i10).GT.err) .AND. (PGmin(6,i10).GT.err)) THEN
              PGmin(2,i10) = -9999
            DO 168 k=1,10
              PGmin(k+2,i10) = -9.88
  168       CONTINUE
          ENDIF
          ENDIF
        ELSEIF(PGmin(12,i10) .GT. 1) THEN
          IF((PGmin(8,i10).GT.err) .AND. (PGmin(9,i10).GT.err)) THEN
          IF((PGmin(10,i10).GT.err) .AND. (PGmin(11,i10).GT.err)) THEN
              PGmin(2,i10) = -9999
            DO 169 k=1,10
              PGmin(k+2,i10) = -9.88
  169       CONTINUE
          ENDIF
          ENDIF
        ENDIF
   16  CONTINUE
C
C  just check error, and no line should be "LT. -1"
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
       DO 160 i10=1,imm10
       DO 160 k=1,5
        IF(PGmax(k+2,i10) .LT. -1) PGmax(k+2,i10) = -9.88
        IF(PGmax(k+7,i10) .LT. -1) PGmax(k+7,i10) = -9.88
        IF(PGmin(k+2,i10) .LT. -1) PGmin(k+2,i10) = -9.88
        IF(PGmin(k+7,i10) .LT. -1) PGmin(k+7,i10) = -9.88
  160  CONTINUE
C########1#########2#########3#########4#########5#########6#########7
C  Obtain 10-min (or 20-min or 60min, given by "intvl") average
C   and SD over this 10-min
C   #  imm:  total number of 1-min data,
C   #  imm10 = imm/intvl: numbers of 10-min bins
C   #  i10: scanning number of all 10-min bns.
C   #  PG(j): PG values of the given minute (i) in each 10-min bin
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  570  CONTINUE
       imm10 = imm/intvl
       DO 17 i10=1,imm10
        PGave(1,i10) = 0
        j4 = 0
        DO 170 i=1,intvl
          j=((i10-1)*intvl)+i
          IF(PG(j) .GT. 9990) THEN
            j4=99
          ELSE
            PGave(1,i10) = PGave(1,i10) + PG(j)
          ENDIF
  170   CONTINUE
        PGave(1,i10) = PGave(1,i10)/intvl
        PGave(2,i10) = 0
        IF(j4 .GT. 10) THEN
          PGave(1,i10) = 9999
          PGave(2,i10) = 9999
        ELSE
          DO 172 i=1,intvl
            j=((i10-1)*intvl)+i
            PGdev = PG(j) - PGave(1,i10) 
            PGdev = PGdev*PGdev 
            PGave(2,i10) = PGave(2,i10) + PGdev 
  172     CONTINUE
          PGave(2,i10) = PGave(2,i10)/intvl
          PGave(2,i10) = sqrt(PGave(2,i10))
        ENDIF
   17  CONTINUE
C########1#########2#########3#########4#########5#########6#########7
C  Write 10 min (or 20 min or 60 min) analyses result (each month)
C   PGmax(144,8) and PGmin(144,8),
C   YYYYMMDD hh:mm maximum %+1min +2min ... +5min -1min ... -5 min
C   YYYYMMDD hh:mm minimum %+1min +2min ... +5min -1min ... -5 min
C   as, e.g., kak_10min_20yymm.csv
C          &  kak_10min_20yymm.csv
C # NOTE-1: max (positive) before min (negative)
C # NOTE-2: decay (+1min,,,+5min) before develope (-1min,,,-5min)
C 
C # NOTE-3 for loop : manual YYMM -> header only for first month
C                   : automatic YYMM -> Feb-Mar merge, Mar-Apr merge
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  580  CONTINUE
       WRITE (Fname2,'(a4,i2,a4,i6,a4)') 
     ?         'kak_',intvl,'min_', 200000+YYMM, '.csv'
       IF(loop .EQ. 1) OPEN (UNIT=ODEV, FILE=Fname2, STATUS='NEW')
       IF(loop .GE. 2) THEN
         GO TO 585
c         IF(MOD(YYMM,100) .EQ. 2) GO TO 585
c         IF(MOD(YYMM,100) .EQ. 3) GO TO 585
c         CLOSE (ODEV)
c         OPEN (UNIT=ODEV, FILE=Fname2, STATUS='NEW')
       ENDIF
        WRITE(ODEV, 481) intvl,intvl
     ?                  , 'PG_max','PG_max',1,2,3,4,5,1,2,3,4,5
     ?                  , 'PG_min','PG_min',1,2,3,4,5,1,2,3,4,5
  481  FORMAT ('YYYYMMDD,hh:mm,UTmax,UTmin,ave_',i2,'m, SD_',i2,'m'
     ?        , 2(2(', ',a6),5(', +',i1,'min'),5(', -',i1,'min'))) 
  585  CONTINUE
       imm10 = imm/intvl
       DO 18 i10=1,imm10
        WRITE(ODEV, 482)  YMD(((i10-1)*intvl)+1), UT(((i10-1)*intvl)+1)
     ?  , UT(jmax(i10)), UT(jmin(i10)), PGave(1,i10), PGave(2,i10)
     ?  , (PGmax(J,i10), J=1,12) , (PGmin(J,i10), J=1,12)
   18  CONTINUE
  482  FORMAT (I8,3(',',a5),2(', ',F6.0)
     ?        , 2(2(', ',F6.0), 5(',',F6.2) ,  5(',',F6.2))) 
c      CLOSE (ODEV)
C########1#########2#########3#########4#########5#########6#########7
C  Finish making master file of the month, toward next month
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  590  WRITE(OTERM,*) ' YYYYMMDD = ', YMD(imm)
  191 CONTINUE
  595 WRITE(OTERM,*) ' End at YYMMDD_in =', YMD0
      CLOSE (ODEV)
      CLOSE (IDEV)
   19 CONTINUE
C########1#########2#########3#########4#########5#########6#########7
C  option: makinging master file for PGmax/PGmin/PGave distribution
C  Define parameters: thre1, thre9, stepin, & intvl
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  596 CONTINUE
      step = 2.
      WRITE(OTERM,*) ' Sorting with 20% step or 200% step'
      WRITE(OTERM,*) ' steps? 1.2 = 20% , 2 = 200% , & end=0'
      READ(ITERM,*) step
      IF(step .LT. 1) GO TO 600
      intvl = 10
      CALL PGsort(step, -intvl)
      CALL PGsort(step, intvl)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
  600 CONTINUE
      WRITE(OTERM,*) ' Stop/End'
      STOP
      END                                                   
C########1#########2#########3#########4#########5#########6#########7
C  Sub-routine:  dayofmon(IYEAR,IMONTH, IDAY)
C
C    calculate numbers of day for given month
C
C    input:  IYEAR (to adjust 1980,1984,1988,1992,1996...)
C            IMONTH (1 to 12)
C
C    output: IDAY (28 to 31)
C
C########1#########2#########3#########4#########5#########6#########7
      SUBROUTINE dayofmon(IYEAR,IMONTH, NDAY)
      INTEGER IYEAR,IMONTH, NDAY
      
      SELECT CASE (IMONTH)
        CASE (1)
          NDAY=31
        CASE (2)
          NDAY=28
        CASE (3)
          NDAY=31
        CASE (4)
          NDAY=30
        CASE (5)
          NDAY=31
        CASE (6)
          NDAY=30
        CASE (7)
          NDAY=31
        CASE (8)
          NDAY=31
        CASE (9)
          NDAY=30
        CASE(10)
          NDAY=31
        CASE(11)
          NDAY=30
        CASE(12)
          NDAY=31
        CASE(13)
          NDAY=10
      END SELECT
      
      IF ((MOD(IYEAR,4).EQ.0) .AND. (IMONTH.EQ.2)) NDAY=29
      RETURN
      END
C########1#########2#########3#########4#########5#########6#########7
C  Sub-routine:  PGsort(stepin, intvl)
C    2017-11-10   (add one-sigma level and debug of PGmin -> PGmax
C    2016-2-12
C                                                 
C #   1min-raw-data/kak_10min_20yymm.csv
C               ==> (read)
C               ==> kak_10min_20yymm_mini.csv
C                or kak_10min_20yymm_maxi.csv ,
C                &  kak_10min_20yymm_aver.csv
C
C  # before sorting, must arrange monthly master file by 13 March
C  => kakmin.f main routine
C
C input : 10 min PG data (with peaks and +-5 min change), like
C 123456789-123456789-123456789-123456789-123456789-123456789-123456789-
C YYYYMMDD, hh:mm, UTmax, UTmin, ave_10, SD_10m, PG_max, +1min, +2min,
C 20060201, 00:10, 00:15, 00:10,    5.1,    1,7,   0.87,  0.87,  0.72,
C
C # NOTE-1: max (positive) before min (negative)
C # NOTE-2: decay (+1min,,,+5min) before develope (-1min,,,-5min)
C
C   parameters:
C   #  thre1 : Minimum level (e.g., -100, +100)
C              where "-" for PGmin and "+" for PGmax
C   #  thre9 : Maximum level (e.g., -2000, +2000)
C   #  stepin: range of each category (e.g., 20%=1.2, 200%=2.0)
C   #  intvl : 10 min or 20 min or 30 min (defined for master file)
C
C   store at:
C   #  PGpeak: PG value for -5min,,,+5min within PG range (category)
C   #  PGrate: relative value for -5min,,,+5min
C   #  averg : average of relative value for each of -5min,,,+5min 
C   #  SDrate: SD of (j=1-3) average and peak (before/after filter) PG
C              and (j=4-13) relative value for each of -5min,,,+5min
C   #  rateM : matrix form of averg+SDrate for 2006-2015 
C              (2005 is summation for 2006-2010 & 2012-2015)
C
C output = sort PGmin or PGmax according to threshold (thre*) values
C
C########1#########2#########3#########4#########5#########6#########7
      SUBROUTINE PGsort(stepin, intvl0)
      INTEGER YYMM, intvl, intvl0, isgn, ithre(30), icent
      INTEGER j, j0,j1,j2,j3,j4, j1old, j2old, k, k9, loop, loop2
      INTEGER*8 imm(30), jmm(30,13), irateM(30, 16, 20, 2)
      INTEGER IDEV, ODEV, ODEV2, ITERM, OTERM
      REAL*8 PGin0(13), PGpeak(30, 5000, 3), averg(30,13)
      REAL*8 SDrate(30, 13), PGdev 
      REAL PGrate(30, 5000, 10), rateM(30, 16, 20), rateM2(30, 16, 20)
      REAL thre1, thre9, step, stepin, thresh(30)
      CHARACTER*5 chsgn
      CHARACTER*8 chsgn2
      CHARACTER*20 Fname2
      CHARACTER*28 Fname3
      CHARACTER*32 Fname4
      CHARACTER*30 Fname5
      CHARACTER*26 UT(30,5000), UT0
      CHARACTER*2 DUMMY
C                                                                
C  I/O device number:  IDEV, ODEV, TERM
C
      IDEV=21                                                    
      ODEV=31
      ODEV2=32
      ODEV3=32
      ITERM=5
      OTERM=6
C
      thre1 = -100
      thre9 = -2000
      step = abs(stepin)
      intvl = abs(intvl0)
      isgn = isign(1,intvl0)
      WRITE (OTERM,*) 'sigh(intvl0) = ', isgn
   60  CONTINUE
      YYMM = 0
      IF(step .LE. 1.) GO TO 699
      IF(intvl0 .LT. -1.) THEN
        thre1 = -100
        thre9 = -2000
        chsgn = '_mini'
      ELSE IF(intvl0 .GT. 1.) THEN
        thre1 = 100
        thre9 = 2000
        chsgn = '_maxi'
      ENDIF
      icent = step*100
      WRITE (chsgn2,'(a5,i3)') chsgn, icent
      thresh(1) = thre1
      IF(intvl0 .LT. 0) THEN
        DO 202 k=2,29
          thresh(k) = thresh(k-1)*step
          IF(thresh(k) .LT. -abs(thre9))  GO TO 602
  202   CONTINUE
  602   k9=k
        thresh(k9) = -9990.
      ELSE IF(intvl0 .GT. 0) THEN
        DO 204 k=2,29
          thresh(k) = thresh(k-1)*step
          IF(thresh(k) .GT. abs(thre9))  GO TO 604
  204   CONTINUE
  604   k9=k
        thresh(k9) = 9990.
      ENDIF
      DO 208 k=1,k9
        ithre(k) = thresh(k)
  208 CONTINUE
      WRITE (OTERM,*) 'k8 and thresh(k8) = ', k9-1, ithre(k9-1)
C########1#########2#########3#########4#########5#########6#########7
C  Start of loop over all months
C  and define file name file  to read (e.g., kak_10min_201103.csv)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
      YYMM=2
      WRITE(OTERM,*) ' Month? 2/3/5/ end=0'
      READ(ITERM,*) YYMM
      loop2 = 0
      DO 29 loop=1,99
        j1old = YYMM/100
        IF((YYMM .EQ. 2) .OR. (YYMM .EQ. 3) .OR. (YYMM .EQ. 5)) THEN
          YYMM=YYMM+1100
        ELSEIF(j1old .EQ. 11) THEN
          YYMM=YYMM-500
        ELSEIF(j1old .EQ. 10) THEN
          YYMM=YYMM+200
        ELSEIF(j1old .EQ. 15) THEN
          YYMM=YYMM-1000
        ELSEIF(j1old .EQ. 5) THEN
          GO TO 690
        ELSEIF((j1old .LT. 5) .OR. (j1old .GT. 15)) THEN
          WRITE(OTERM,*) 'YYMM, loop = ',YYMM, loop
          IF(loop .NE. 1) THEN
           CLOSE (ODEV2)
          ENDIF
          GO TO 699
        ELSE
          YYMM=YYMM+100
        ENDIF
      WRITE (Fname2,'(a4,i2,a4,i6,a4)')
     ?         'kak_',intvl,'min_', 200000+YYMM, '.csv'
      OPEN (UNIT=IDEV, FILE=Fname2, STATUS='OLD')
        WRITE (OTERM,'(A7,A20)') 'open = ', Fname2
        READ(IDEV,'(A2)') DUMMY
        DO 210 k=1,k9
          imm(k) = 0
  210   CONTINUE
C
C  Read PG from kak_10min_20yymm.csv, as PGin0
C 4081     FORMAT (a26, 1x, F7.0, 8x, 86x, 2(1x,F7.0), 10(1x,F6.2))
C 4082     FORMAT (a26, 1x, F7.0, 8x, 2(1x,F7.0), 10(1x,F6.2))
C   # UT0 : first 26 characters, i.e., all UT info.
C           YYYYMMDD,hh:mm,UTmax,UTmin, => a26+1x
C   # PGin0(1) : average PG over intvl (e.g., over 10 min).
C   # (SD is skipped)
C   # PGin0(2) : maximum/minimum PG over intvl (e.g., over 10 min).
C   # PGin0(3) : same as PGin0(*,2) if >+100 or <-50 and real peak
C   # PGin0(4-13): decay 1 min to 5 min & develope -5 min to -1 min
C   # PGpeak : average, maximum, and minimum
C   # PGrate : %PG value at -5min, -4min,,,-1min, +1min, +2min,,,+5min
C   # k      : k'th threshold
C   # imm(k) : counter of 10-min bins for k'th threshold
C   # thre1  : -100 (for minimum PG) or +100 (for maximum PG)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        IF(intvl0 .LT. 0) THEN
          DO 21 i=1,90000
            READ(IDEV, 4081,END=629)  UT0,  (PGin0(J), J=1,13)
            if(i.lt.3) write(OTERM,*) UT0, (PGin0(J), J=1,2)
            DO 214 k=1,k9-1
              IF( (PGin0(3).LE.thresh(k)) .AND.
     ?            (PGin0(3).GT.thresh(k+1))         ) THEN
                imm(k) = imm(k)+1
                UT(k,imm(k)) = UT0
                PGpeak(k,imm(k),1) = PGin0(1)
                PGpeak(k,imm(k),2) = PGin0(2)
                PGpeak(k,imm(k),3) = PGin0(3)
C
C # note: we do not use PGpeak(k,imm(k),2) or PGin0(2) here
C  change order between decay (+min)and develope (-min) when store
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
                DO 212 j=1,5
                  PGrate(k,imm(k),j) = PGin0(14-j)
                  PGrate(k,imm(k),j+5) = PGin0(j+3)
  212           CONTINUE
              ENDIF
  214       CONTINUE
  21      CONTINUE
C
C  The same procedure for PGmax.
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        ELSEIF(intvl0 .GT. 0) THEN
          DO 22 i=1,90000
            READ(IDEV, 4082,END=629)  UT0,  (PGin0(J), J=1,13)
            DO 218 k=1,k9-1
              IF( (PGin0(3).GE.thresh(k)) .AND.
     ?            (PGin0(3).LT.thresh(k+1))         ) THEN
                imm(k) = imm(k)+1
                UT(k,imm(k)) = UT0
                PGpeak(k,imm(k),1) = PGin0(1)
                PGpeak(k,imm(k),2) = PGin0(2)
                PGpeak(k,imm(k),3) = PGin0(3)
                DO 216 j=1,5
                  PGrate(k,imm(k),j) = PGin0(14-j)
                  PGrate(k,imm(k),j+5) = PGin0(j+3)
  216           CONTINUE
              ENDIF
  218       CONTINUE
 4081     FORMAT (a26, 1x, F7.0, 8x, 86x, 2(1x,F7.0), 10(1x,F6.2))
 4082     FORMAT (a26, 1x, F7.0, 8x, 2(1x,F7.0), 10(1x,F6.2))
  22    CONTINUE
        ENDIF
  629 CLOSE (IDEV)
        WRITE (OTERM,'(a13,5i5)') '# of samples:', (imm(k), k=1,5)
C########1#########2#########3#########4#########5#########6#########7
C   # k       : k'th threshold
C   # imm(k)  : total number of 10-min bins for k'th threshold
C   # j2      : minute from peak (-5, -4,,,-1, +1,,,+5)
C   # jmm(k,j2): counter for k'th threshold and valid value of PGrate
C  exclude PGrate = -9.99 (not applicabel), -9.88 (end of trend)
C  include/exclude PGrate = -0.01 (reaching zero, counted as -0.01)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        DO 25 k=1,k9-1
          DO 250 j=1,13
            averg(k,j) = 0.
            jmm(k,j) = 0
  250     CONTINUE
          DO 252 j=1,10
          j2 = j+3
          DO 252 i=1,imm(k)
            IF((PGrate(k,i,j).GE.-1) .AND. (PGrate(k,i,j).LE.1)) THEN
c-see above  IF((PGrate(k,i,j).GE.0) .AND. (PGrate(k,i,j).LE.1)) THEN
              jmm(k,j2) = jmm(k,j2) + 1
              averg(k,j2) = averg(k,j2) + PGrate(k,i,j)
            ENDIF
  252     CONTINUE
          DO 254 j=1,3
          DO 254 i=1,imm(k)
            IF( (PGpeak(k,i,j).LT.9990) 
     ?          .AND. (PGpeak(k,i,j).GT.-9990) ) THEN
              jmm(k,j) = jmm(k,j) + 1
              averg(k,j) = averg(k,j) + PGpeak(k,i,j)
            ENDIF
  254     CONTINUE
          DO 256 j=1,13
  256       IF(jmm(k,j) .GE. 1) averg(k,j) = averg(k,j)/jmm(k,j)
  25    CONTINUE
C
C # note: we do not use j=2 (j=3 is the same value after filtering)
C   derive monthly standard deviation for each category each minute.
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        DO 26 k=1,k9-1
          DO 260 j=1,13
            SDrate(k,j) = 0.
  260     CONTINUE
          DO 262 j=1,10
            j2 = j+3
          DO 262 i=1,imm(k)
            IF((PGrate(k,i,j).GE.-1) .AND. (PGrate(k,i,j).LE.1)) THEN
c-see above  IF((PGrate(k,i,j).GE.0) .AND. (PGrate(k,i,j).LE.1)) THEN
              PGdev = PGrate(k,i,j) - averg(k,j2)
              SDrate(k,j2) = SDrate(k,j2) + (PGdev*PGdev)
            ENDIF
  262     CONTINUE
          DO 264 j=1,3
          DO 264 i=1,imm(k)
            IF( (PGpeak(k,i,j).LT.9990) 
     ?          .AND. (PGpeak(k,i,j).GT.-9990) ) THEN
              PGdev = PGpeak(k,i,j) - averg(k,j)
              SDrate(k,j) = SDrate(k,j) + (PGdev*PGdev)
            ENDIF
  264     CONTINUE
          DO 266 j=1,13
            IF(jmm(k,j) .GE. 1) SDrate(k,j) = SDrate(k,j)/jmm(k,j)
            IF(jmm(k,j) .GE. 1) SDrate(k,j) = sqrt(SDrate(k,j))
  266     CONTINUE
  26    CONTINUE
C########1#########2#########3#########4#########5#########6#########7
C   Write results as, e.g., kak_10min_yyyymm_maxi.csv
C                       &   kak_10min_yyyymm_mini.csv
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
      WRITE (Fname3,'(a4,i2,a4,i6,a8,a4)')
     ?        'kak_',intvl,'min_', 200000+YYMM, chsgn2,'.csv'
      OPEN (UNIT=ODEV, FILE=Fname3, STATUS='NEW')
        DO 27 k=k9-1,1,-1
          IF(intvl0 .LT. 0) THEN
            WRITE(ODEV, 4781) intvl, chsgn, 5,4,3,2,1,1,2,3,4,5
          ELSEIF(intvl0 .GT. 0) THEN
            WRITE(ODEV, 4782) intvl, chsgn, 5,4,3,2,1,1,2,3,4,5
          ENDIF
          DO 272 i=1,imm(k)
            WRITE(ODEV, 4783)  UT(k,i), (PGpeak(k,i,J), J=1,3)
     ?                   ,   (PGrate(k,i,J), J=1,10)
  272     CONTINUE
          IF(intvl0 .LT. 0) THEN
            WRITE(ODEV, 4781) intvl, chsgn, 5,4,3,2,1,1,2,3,4,5
          ELSEIF(intvl0 .GT. 0) THEN
            WRITE(ODEV, 4782) intvl, chsgn, 5,4,3,2,1,1,2,3,4,5
          ENDIF
          WRITE(ODEV, 4785) ithre(k+1), (jmm(k,J), J=1,13)
          WRITE(ODEV, 4786) ithre(k), (averg(k,J), J=1,13)
          WRITE(ODEV, 4787) ithre(k), (SDrate(k,J), J=1,13)
          WRITE(ODEV, '(2x)') 
 4781     FORMAT ('YYYYMMDD,hh:mm,UTmin,',i2,'min_ave,simp_peak'
     ?           , ',  PG', a5, 5(', -',i1,'min'),5(', +',i1,'min') ) 
 4782     FORMAT ('YYYYMMDD,hh:mm,UTmax,',i2,'min_ave,simp_peak'
     ?           , ',  PG', a5, 5(', -',i1,'min'),5(', +',i1,'min') )
 4783     FORMAT (a20, 3(',   ', F6.0), 10(',',F6.2))
 4785     FORMAT (4x,'from,',i5,', # = ',3(',   ',I6),10(',',I6))
 4786     FORMAT (4x,'to  ,',i5,', ave=',3(',   ',F6.0)
     ?                  , 10(',',F6.2)) 
 4787     FORMAT (4x,'to  ,',i5,', sd =',3(',   ',F8.2)
     ?                  , 10(',',F6.3))
  27    CONTINUE
      CLOSE (ODEV)
C
C   Write results as, e.g., kak_10min_yyyymm_mini_aver.csv
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        DO 28 k=1,k9-1
          DO 281 j=1,16
            irateM(k, j, loop, 1) = 0
            irateM(k, j, loop, 2) = 0
            rateM(k, j, loop) = 0.
            rateM2(k, j, loop) = 0.
  281     CONTINUE
          irateM(k, 1, loop, 1) = 200000+YYMM
          irateM(k, 7, loop, 1) = jmm(k,3)
          DO 283 j=1,5
            irateM(k, j+1, loop, 1) = jmm(k,j+3)
            irateM(k, j+7, loop, 1) = jmm(k,j+8)
  283     CONTINUE
          irateM(k, 13, loop, 1) = jmm(k,3)
          irateM(k, 14, loop, 1) = jmm(k,1)
          irateM(k, 15, loop, 1) = ithre(k)
          irateM(k, 16, loop, 1) = ithre(k+1)
C
C   set 1 or 0 across 6/7 boundary =INT(MIN(AN21,7)/7)
C
          irateM(k, 1, loop, 2) = irateM(k, 1, loop, 1) 
          DO 284 j=2,14
            IF(irateM(k,j,loop,1) .GE. 7) irateM(k,j,loop,2) = 1
  284     CONTINUE
          irateM(k, 15, loop, 2) = irateM(k, 15, loop, 1) 
          irateM(k, 16, loop, 2) = irateM(k, 16, loop, 1) 
C
          rateM(k, 1, loop) = 200000+YYMM
          rateM(k, 7, loop) = 1.
          DO 286 j=1,5
            rateM(k, j+1, loop) = averg(k,j+3)
            rateM(k, j+7, loop) = averg(k,j+8)
  286     CONTINUE
          rateM(k, 13, loop) = averg(k,3)
          rateM(k, 14, loop) = averg(k,1)
          rateM(k, 15, loop) = ithre(k)
          rateM(k, 16, loop) = ithre(k+1)
C
          rateM2(k, 1, loop) = 200000+YYMM
          rateM2(k, 7, loop) = 0.
          DO 288 j=1,5
            rateM2(k, j+1, loop) = SDrate(k,j+3)
            rateM2(k, j+7, loop) = SDrate(k,j+8)
  288     CONTINUE
          rateM2(k, 13, loop) = SDrate(k,3)
          rateM2(k, 14, loop) = SDrate(k,1)
          rateM2(k, 15, loop) = ithre(k)
          rateM2(k, 16, loop) = ithre(k+1)
  28    CONTINUE
        loop2 = loop
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
        IF(loop .GE. 2) GO TO 689
        WRITE (Fname4,'(a4,i2,a4,i6,a8,a8)')
     ?        'kak_',intvl,'min_', 200000+YYMM, chsgn2,'_ave.csv'
      OPEN (UNIT=ODEV2, FILE=Fname4, STATUS='NEW')
        WRITE(ODEV2, 4881) 'yyyymm', chsgn, intvl,',     range     '
     ?        , 5,4,3,2,1,1,2,3,4,5
 4881   FORMAT (a6,', item , ','PG', a5,', ', i2,'m_ave', a16
     ?        , 5(', -',i1,' min'),',  0 min',5(', +',i1,' min') )
        DO 292 k=k9-1,1,-1
          WRITE(ODEV2, 4883) 200000+YYMM, jmm(k,3), jmm(k,1)
     ?        , ithre(k), ithre(k+1), (jmm(k,J), J=4,8), jmm(k,3)
     ?        , (jmm(k,J), J=9,13)
          WRITE(ODEV2, 4884) 200000+YYMM, averg(k,3), averg(k,1)
     ?        , ithre(k), ithre(k+1), (averg(k,J), J=4,8)
     ?        , (averg(k,J), J=9,13)
          WRITE(ODEV2, 4885) 200000+YYMM, SDrate(k,3), SDrate(k,1)
     ?        , ithre(k), ithre(k+1), (SDrate(k,J), J=4,8)
     ?        , (SDrate(k,J), J=9,13)
  292   CONTINUE
 4883   FORMAT (i6,',  # = ,', 2(I6',  '), i6,' to',i6, 11(',',I7))
 4884   FORMAT (i6,', ave =,', 2(F6.0',  '), i6,' to',i6
     ?         , 5(',',F7.2), ',   1.00', 5(',',F7.2))
 4885   FORMAT (i6,', sd = ,', 2(F6.0',  '), i6,' to',i6
     ?         , 5(',',F7.4), ',   0.00', 5(',',F7.4))
  689   CONTINUE
  29  CONTINUE
  690 CONTINUE
      CLOSE (ODEV2)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
      WRITE (Fname5,'(a4,i2,a11,i1,a8,a4)')
     ?        'kak_',intvl,'min_decay_0', MOD(YYMM,100), chsgn2,'.csv'
      IF(loop2 .GE. 1)  OPEN (UNIT=ODEV, FILE=Fname5, STATUS='NEW')
        WRITE(ODEV, 4981) ' yyyymm,', 5,4,3,2,1,1,2,3,4,5, chsgn,intvl
     ?                 , ' yyyymm,', 5,4,3,2,1,1,2,3,4,5, chsgn,intvl
 4981   FORMAT (a8, ' item ', 5(', -',i1,' min')
     ?     , ',  0 min', 5(', +',i1,' min'),', PG', a5,',', i2, 'm_avr'
     ?     , ',',9x,'range ,', a8, ' item ', 5(', -',i1,' min')
     ?     , ',  0 min', 5(', +',i1,' min'),', PG', a5,',', i2, 'm_avr'
     ?     , ',',9x,'range ')
        DO 30 k=k9-1,1,-1
        DO 30 i=1,loop2
          WRITE(ODEV, 4983) (irateM(k, j, i, 1), j=1,16)
     ?                      , (rateM(k, j, i), j=1,16)
          WRITE(ODEV, 4985) (irateM(k, j, i, 2), j=1,16)
     ?                      , (rateM2(k, j, i), j=1,16)
  30    CONTINUE
 4983   FORMAT (i7,',   # ,', 11(I7,','), 2(I8,','), i6,' to',i6,','
     ?      , F7.0,', ave ,', 11(F7.2,','), 2(F8.0,',')
     ?      , F6.0,' to',F6.0)
 4985   FORMAT (i7,',>=7? ,', 11(I7,','), 2(I6,','), i6,' to',i6,','
     ?      , F7.0,',  sd ,', 11(F7.4,','), 2(F8.2,',')
     ?      , F6.0,' to',F6.0)
  699 WRITE(OTERM,*) ' End at YYMM =', YYMM
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7
      RETURN
      END                                                   
