c=======================================================================
c PROGRAM
c   p1.get_report.selected_dates.f
c
c LANGUAGE
c   Fortran77
c
c-----------------------------------------------------------------------
c PURPOSE
c   to generate report of stations used to calculate precip values on
c   selected 0.25deg grid boxes for selected dates
c
c USAGE
c   p1.get_report.selected_dates.x <yyyymm>
c
c-----------------------------------------------------------------------
c INPUT FILES
c   boxes_to_check                       :
c   point_report_<yyyymmdd>              :
c   globe_uni_dly_prcp_stn_rt.<yyyymmdd> :
c
c OUTPUT FILES
c   report_<yyyymm> :
c
c INCLUDED FILES
c   oi.chk_boxes.h : NCH0
c
c-----------------------------------------------------------------------
c SUBROUTINES USED
c   distant() : 
c
c-----------------------------------------------------------------------
c COMMAND-LINE VARIABLES
c   ym : target month (optional, used in output file name)
c
c CONSTANTS
c   NB  : number of 0.25deg grid boxes to check
c   NPB : number of 0.125deg point-values used to get 0.25deg box-value [=9]
c   NP  : max number of 0.125deg grid points to check (not used here)
c   ----:------------------------------------------------------------
c   NS  : max number of stations
c   NSP : max number of station-values used to get 0.125deg point-value
c   NSB : max number of station-values used to get 0.25deg box-value
c   DLL : lat/lon resolution of 0.125deg grid
c
c LOCAL VARIABLES
c   f_box   :
c   f_fnl   :
c   f_pnt   :
c   d_stn   :
c   f_stn   :
c   ym      :
c   ln      :
c   stn()   :
c   name    :
c   sid     :
c   item()  :
c   bid     :
c   ymd     :
c   yyyy    :
c   bb      :
c   pp      :
c   sp      :
c   sb      :
c   cs      :
c   cp      :
c   csp     :
c   csb     :
c   flag    :
c   id_p()  :
c   id_b()  :
c   blat    :
c   blon    :
c   lat     :
c   lon     :
c   plat    :
c   plon    :
c   plats() :
c   plons() :
c   vs      :
c   vs_p()  :
c   vs_b()  :
c   vp      :
c   vps()   :
c   va      :
c   vb      :
c   vc      :
c   vd      :
c   vbox    :
c   w1      :
c   w2      :
c   dis     :
c   f_there :
c 
c-----------------------------------------------------------------------
c NOTE
c   1. The precip interpolation is performed on 0.125deg grid points.
c
c   2. Precip box-value of 0.25deg grid box centered at point 5 is weighted
c      average of box-values of four 0.125deg grid boxes centered at points
c      A, B, C, and D.
c
c      Precip box-value of 0.125deg grid box centered at point A is averge
c      of point-values of four 0.125deg grid points 1, 2, 4, and 5.
c
c        7---8---9
c        | C | D |
c        4---5---6
c        | A | B |
c        1---2---3
c
c      vb(A)=(vp(1)+vp(2)+vp(4)+vp(5))/4
c      vb(B)=(vp(2)+vp(3)+vp(5)+vp(6))/4
c      vb(C)=(vp(4)+vp(5)+vp(7)+vp(8))/4
c      vb(D)=(vp(5)+vp(6)+vp(8)+vp(9))/4
c
c      vb(5)=(vb(A)*wa+vb(B)*wb+vb(C)*wc+vb(D)*wd)/(wa+wb+wc+wd)
c
c      where wa=wb=cos(lat_ab), wc=wd=cos(lat_cd)
c 
c=====|==|=========|=========|=========|=========|=========|=========|==
c
      program p1_get_report_selected_dates
      implicit none

c-----------------------------------------------------------------------
c NB, NPB & NP are defined in 'oi.chk_boxes.h'
c-----------------------------------------------------------------------
      include 'oi.chk_boxes.h'

      integer NS, NSP, NSB
      real DLL
      parameter( NS=40000, NSP=10, NSB=NSP*NPB, DLL=0.125 )

      character ym*6, f_box*29, f_rpt*20, f_pnt*28
      character info_stn(NS)*51, sid*9, item(6)*10
      integer bid, yyyymmdd, yyyy, bb, pp, sp, sb, cs, cp, csp, csb
      integer flag, id_p(NSP), id_b(NSB)
      real blat, blon_w, blon
      real lat, lon, plat, plon, plats(NPB), plons(NPB)
      real vs, vs_p(NSP), vs_b(NSB), vp, vps(NPB)
      real va, vb, vc, vd, vbox, w1, w2, dis
      logical f_there, qc_run

      character*9 stnid_stn(NS)
      real rlat0_stn(NS), rlon0_stn(NS), rain0_stn(NS)

      if( iargc().lt.1 ) then
        write(*, *)
        write(*, *) '##########################################'
        write(*, *) '#                                        #'
        write(*, *) '#  parameter <yyyymm> [<qc_run>] needed  #'
        write(*, *) '#                                        #'
        write(*, *) '##########################################'
        write(*, *)
        stop
      endif
      call getarg(1, ym)

      qc_run=(iargc().ge.2)

c------------->123456789|<--
      item(1)='      Stn#'
      item(2)='  Prcp(mm)'
      item(3)=' StnID    '
      item(4)='    StnLat'
      item(5)='    StnLon'
      item(6)='  Dist(km)'

c-----------------------------------------------------------------------
c construct file names
c-----------------------------------------------------------------------
c----------->123456789|123456789|123456789<--
      f_box='boxes_to_check/boxes_to_check'
      f_rpt='report/report_yyyymm'
      f_pnt='output/point_report_yyyymmdd'

      f_rpt(15:20)=ym

      open(10, file=f_box, status='old', form='formatted')
      read(10, *)
      read(10, *)
      read(10, *)

      open(20, file=f_rpt, form='formatted')

c-----------------------------------------------------------------------
c repeat for each selected box
c-----------------------------------------------------------------------
      do bb=1, NB

        read(10, '(10x,i10,2f10.3,i10)') bid, blat, blon_w, yyyymmdd
        blon=360.0+blon_w
        yyyy=yyyymmdd/10000

        write(f_pnt(21:28), '(i8.8)') yyyymmdd

        write(20, 900) bid, blat, blon_w, blon
        write(*,  900) bid, blat, blon_w, blon
 900    format(4hBox-,i5.5,6h, Lat=,f6.3,6h, Lon=,f8.3,2h (,f7.3,1h))

        plats(1)=blat-DLL
        plats(2)=blat-DLL
        plats(3)=blat-DLL
        plats(4)=blat
        plats(5)=blat
        plats(6)=blat
        plats(7)=blat+DLL
        plats(8)=blat+DLL
        plats(9)=blat+DLL

        plons(1)=blon-DLL
        plons(2)=blon
        plons(3)=blon+DLL
        plons(4)=blon-DLL
        plons(5)=blon
        plons(6)=blon+DLL
        plons(7)=blon-DLL
        plons(8)=blon
        plons(9)=blon+DLL

c-----------------------------------------------------------------------
c repeat for each day in target month
c-----------------------------------------------------------------------
        inquire(file=f_pnt, exist=f_there)
        if( f_there ) then

          call read_stn(yyyymmdd, qc_run, stnid_stn, rlat0_stn,
     1         rlon0_stn, rain0_stn, cs, info_stn)

          open(11, file=f_pnt, status='old', form='formatted')

          cp=0
          csb=0
c-----------------------------------------------------------------------
c read in each point report until end of file
c-----------------------------------------------------------------------
 200      read(11, *, end=290)
          read(11, '(2f10.3,i10)') plon, plat, csp
          do sp=1, csp
            read(11, '(10x,i10,10x,f10.3)') id_p(sp), vs_p(sp)
          enddo
          read(11, '(f10.3)') vp

c-----------------------------------------------------------------------
c check if working point is one of the 9 points
c-----------------------------------------------------------------------
          do pp=1, NPB
            if( plon.eq.plons(pp).and.plat.eq.plats(pp) ) then
              cp=cp+1
              vps(pp)=vp
              do sp=1, csp
                do sb=1, csb
                  if( id_b(sb).eq.id_p(sp) ) goto 210
                enddo
                csb=csb+1
                id_b(csb)=id_p(sp)
                vs_b(csb)=vs_p(sp)
 210            continue
              enddo
              exit
            endif
          enddo

          goto 200

c-----------------------------------------------------------------------
c check if all 9 points have been found
c-----------------------------------------------------------------------
 290      if( cp.ne.NPB ) then
            write(*, *) '#### point count wrong: ', cp
            stop
          else

c-----------------------------------------------------------------------
c define weights
c-----------------------------------------------------------------------
            w1=blat-DLL/2.0
            w1=w1*3.1416/180.0
            w1=cos(w1)
            w2=blat+DLL/2.0
            w2=w2*3.1416/180.0
            w2=cos(w2)

c-----------------------------------------------------------------------
c define 0.125deg box-values using 0.125deg point-values
c-----------------------------------------------------------------------
            va=(vps(1)+vps(2)+vps(4)+vps(5))/4.0
            vb=(vps(2)+vps(3)+vps(5)+vps(6))/4.0
            vc=(vps(4)+vps(5)+vps(7)+vps(8))/4.0
            vd=(vps(5)+vps(6)+vps(8)+vps(9))/4.0

c-----------------------------------------------------------------------
c define 0.25deg box-value using 0.125deg box-values
c-----------------------------------------------------------------------
            vbox=(va*w1+vb*w1+vc*w2+vd*w2)/(w1+w1+w2+w2)

            write(20, 920) yyyymmdd, vbox*0.1, csb
            write(*,  920) yyyymmdd, vbox*0.1, csb
 920        format(5hDate=,i8,13h, Precip(mm)=,f6.2,11h, #StnUsed=,i2)
            write(20, 930) ('-', pp=1, 60)
            write(*,  930) ('-', pp=1, 60)
 930        format(60a)
            write(20, 940) (item(pp), pp=1, 6)
            write(*,  940) (item(pp), pp=1, 6)
 940        format(6a10)
            write(20, 930) ('-', pp=1, 60)
            write(*,  930) ('-', pp=1, 60)

            do sb=1, csb
c              read(stn(id_b(sb)), 950) sid, lat, lon, vs
 950          format(6x, a10, 2x, 2f8.2, 7x, f8.2)
              sid=stnid_stn(id_b(sb))
              lat=rlat0_stn(id_b(sb))
              lon=rlon0_stn(id_b(sb))
              vs=rain0_stn(id_b(sb))
              call distant(dis, blat, blon, lat, lon)
              write(20, 960) sb, vs*0.1, sid, lat, lon, dis
              write(*,  960) sb, vs*0.1, sid, lat, lon, dis
 960          format(i10, f10.2, 1x,a9, 3f10.2)
            enddo
          endif

          close(11)

        endif

        write(20, *)
        write(*,  *)

      enddo

      close(10)
      close(20)
      write(*, *)
      write(*, *) '  ==> ', f_rpt
      write(*, *)

      stop
      end
c
c=====|==|=========|=========|=========|=========|=========|=========|==