*****************************************************************
* PROGRAM:    globe_ssmi_pen_def_ssmi_pen_v0711.lnx.f
*
* LANGUAGE:   Fortran77
*
* MACHINE:    Any LINUX BOX
*----------------------------------------------------------------*
* PURPOSE:    to define GrADS formatted file of gridded pen ssmi data
*----------------------------------------------------------------*
* USAGE:
*----------------------------------------------------------------*
* INPUT FILES:    Unit 10 == availability of SSM/I daily files
*                      20 == yearly file of daily SSM/I 
*                      30 == SSM/I data mask (avoid snow contamination) 
*
* OUTPUT FILES:   Unit 50 == GrADS formatted pentad SSM/I on 2.5deg lat/lon
*----------------------------------------------------------------*
* SUBROUTINES USED:   NONE
*
* FUNCTIONS USED:     (all Standard Fortran 77 functions)
*
*----------------------------------------------------------------*
* INPUT VARIABLES:    kyy == year of the target date
*                     kmm == month of the target date
*                     kdd == day of the target date
*                     kjj == Julian Day
*
* LOCAL VARIABLES:    kstatus      == data availability
*                     rmask        == SSMI data mask 
*                     rain0        == SSMI daily precp     (2.5deg)
*                     samp0        == SSMI daily sampling  (2.5deg)
*                     rain         == SSMI pentad precp    (2.5deg)
*                     samp         == SSMI pentad sampling (2.5deg)  
*
*----------------------------------------------------------------*
* DATE              2007-11-16 
*----------------------------------------------------------------*
* MODIFICATIONS:
******************************************************************

      character yyyymmdd*8, f_msk*35, f_dly*54, f_pen*48
      integer yyyy, mm, dd, days(12), m, jjj

      dimension    rain0(144,72), samp0(144,72), 
     #             rmask(144,72), kcheck (144),  
     #             rain (144,72), samp (144,72), ndat (144,72)     

      data days / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /

c     ************************************************************************
c     *                                                                      *
c     *   1.  to get the target date                                         *
c     *                                                                      *
c     ************************************************************************

      if( iargc().lt.1 ) then
        write(*, *)
        write(*, *) '#################################'
        write(*, *) '#                               #'
        write(*, *) '#  parameter <yyyymmdd> needed  #'
        write(*, *) '#                               #'
        write(*, *) '#################################'
        write(*, *)
        stop
      endif

      call getarg(1, yyyymmdd)
      read(yyyymmdd, '(i4.4,2i2.2)') yyyy, mm, dd

      if( mod(yyyy,4).eq.0 )   days(2)=29
      if( mod(yyyy,100).eq.0 ) days(2)=28
      if( mod(yyyy,400).eq.0 ) days(2)=29

      jjj=dd
      do m=1, mm-1
        jjj=jjj+days(m)
      enddo

      if( days(2).eq.29 ) then
        kleap=1
      else
        kleap=0
      endif

c     ************************************************************************
c     *                                                                      *
c     *   2.  to determine if it is the date to define the pentad            *
c     *                                                                      *
c     ************************************************************************

      if( jjj.le.59 )  then 
        kdays = jjj
      else 
        kdays = jjj - kleap 
      end  if  
      if( mod (kdays,5).ne.0 ) then
        write(*, *) '## NOT end of pentad: ', yyyy, mm, dd, jjj
        stop
      endif
 
      kpen = kdays / 5 
      write(*, *) yyyy, mm, dd, jjj, ' pentad: ', kpen
      if  (kpen.le.11)  then 
        kday1 = (kpen-1)*5 + 1 
        kday2 = kday1 + 4 
      else  if  (kpen.eq.12)  then 
        kday1 = 56 
        kday2 = 60 + kleap 
      else 
        kday1 = (kpen-1)*5 + 1 + kleap  
        kday2 = kday1 + 4 
      end  if  
c
c
c     ************************************************************************
c     *                                                                      *
c     *   3.  to open the daily and pentad ssmi/i files                      *
c     *                                                                      *
c     ************************************************************************

c----------->123456789|123456789|123456789|123456789|123456789|1234<--
      f_dly='../output/SSMI/yyyy/ppt/globe_dly_2.5deg_ssmi.lnx.yyyy' 
      f_pen='../output/SSMI/yyyy/ppt/ssmi_pen_2.5deg.lnx.yyyy'
      f_msk='../library/nesdis_ssmi_mask_pr1.lnx'

      write(f_dly(16:19), '(i4.4)') yyyy
      write(f_dly(51:54), '(i4.4)') yyyy

      write(f_pen(16:19), '(i4.4)') yyyy
      write(f_pen(45:48), '(i4.4)') yyyy


      open(10, file=f_dly, status='old', access='direct', recl=144*72*4)
      open(11, file=f_msk, status='old', access='direct', recl=144*72*4)

      open(20, file=f_pen, access='direct', recl=144*72*4)

c     ************************************************************************
c     *                                                                      *
c     *   4.  to define the pentad gpi                                       *
c     *                                                                      *
c     ************************************************************************
c
      do 4001 jj = 1, 72 
      do 4001 ii = 1, 144 
        ndat (ii,jj) = 0 
        rain (ii,jj) = 0.0 
        samp (ii,jj) = 0.0 
 4001 continue
c
      do 4003 kday = kday1, kday2  
        kinp1 = (kday-1)*2 + 1 
        kinp2 = (kday-1)*2 + 2  
        read  (10,rec=kinp1)   rain0  
        read  (10,rec=kinp2)   samp0  
        do 4002 jj =  1, 72  
        do 4002 ii =  1, 144  
          if  (rain0(ii,jj).ge.0.0.and. 
     #         samp0(ii,jj).gt.0.0)  then 
            ndat (ii,jj) = ndat (ii,jj) + 1 
            rain (ii,jj) = rain (ii,jj) + rain0 (ii,jj) 
            samp (ii,jj) = samp (ii,jj) + samp0 (ii,jj) 
          end  if    
 4002   continue          
        write(*, *) 'pentad: ', kpen, ' day: ', kday
 4003 continue 
c
      do 4004 jj = 1, 72 
      do 4004 ii = 1, 144 
        if  (ndat(ii,jj).ge.1)  then 
          rain (ii,jj) = rain (ii,jj) / ndat (ii,jj)  
        else 
          rain (ii,jj) = -999.0 
          samp (ii,jj) = 0.0 
        end  if  
 4004 continue  
       write  (6,*)  'finished defining pentad SSMI !!' 
c
c
c     ************************************************************************
c     *                                                                      *
c     *   5.  to mask out areas with no values                               *
c     *                                                                      *
c     ************************************************************************
c
      read  (11,rec=kpen)   rmask 
      do 5001 jj = 1, 72 
      do 5001 ii = 1, 144 
        if  (rmask(ii,jj).eq.0.0)  then 
          rain (ii,jj) = -999.0  
          samp (ii,jj) = 0.0 
        end  if  
 5001 continue 
      write  (6,*)  'finished applying mask on pentad SSMI !!' 
c
c
c     ************************************************************************
c     *                                                                      *
c     *   6.  to kick out error                                              *
c     *                                                                      *
c     ************************************************************************
c
c     6.1  to kick out large values 
      do 6001 jj = 1, 16   
      do 6001 ii = 1, 144 
        if  (rain(ii,jj).gt.10.0)  then 
          rain (ii,jj) = -999.0 
          samp (ii,jj) = 0.0 
        end  if  
 6001 continue 
c
c     6.2  to kick out suspicious values over mid- and high latitudes 
      do 6002 ii = 1, 144 
          kcheck (ii) = 1 
 6002 continue  
      do 6003 ii = 1, 144 
      do 6003 jj = 71, 53, -1  
        if  (kcheck(ii).eq.1.and.
     #       rain(ii,jj).gt.10.0)  then 
          rain (ii,jj) = -999.0 
          samp (ii,jj) = 0.0
        end  if  
        i1 = ii - 1 
        if  (i1.eq.0)  then 
          i1 = 144 
        end  if 
        i2 = ii 
        i3 = ii + 1 
        if  (i3.eq.145)  then 
          i3 = 1 
        end  if  
        j1 = jj - 1 
        j2 = jj 
        j3 = jj + 1 
        if  (rain(i1,j2).ge.0.0.and.
     #       rain(i3,j2).ge.0.0.and.
     #       rain(i1,j3).ge.0.0.and.
     #       rain(i2,j3).ge.0.0.and.
     #       rain(i3,j3).ge.0.0.and. 
     #       rain(i1,j1).ge.0.0.and. 
     #       rain(i2,j1).ge.0.0.and. 
     #       rain(i3,j1).ge.0.0)  then 
          kcheck (ii) = 0 
        end  if  
 6003 continue 
      write  (6,*)  'finished kicking out suspicious values !!' 
c
c
c     ************************************************************************
c     *                                                                      *
c     *   7.  to output                                                      *
c     *                                                                      *
c     ************************************************************************

      write(20, rec=(kpen-1)*2+1)   rain
      write(20, rec=(kpen-1)*2+2)   samp

      close(10)
      close(11)
      close(20)

      stop
      end