;****************************************************************************
; Author  :  Josh Grant, COAPS
;
; Description : This program takes a divergence file which is in the format
;               previously specified for this project and grids it based on
;               the week of the swath.  Each weekly value is then used to
;               calculate the weekly average of each week in the NSCAT
;               mission.
;
; **Warning** : Due to poor planning on my part and lack of time the bulk of
;               this program is not commented like it should be.
;****************************************************************************

;@/usr/people/whalley/lidl/jday2date.pro
;@/usr/people/whalley/lidl/date2jday.pro

PRO WRITEWEEK, grid, totals, week, weeklyDirectory
    stringFormat = STRCOMPRESS( '(i' + STRING( 2 ) + '.' $
                   + STRING( 2 ) + ')', /REMOVE_ALL )
    IF ( week LT 15 ) THEN BEGIN
      file_name = weeklyDirectory + STRTRIM( 1996, 1 ) + $
        STRING( week + 38, format = stringFormat ) + ".week.avg"
    ENDIF ELSE BEGIN
      file_name = weeklyDirectory + STRTRIM( 1997, 1 ) + $
        STRING( week - 14, format = stringFormat ) + ".week.avg"
    ENDELSE
    FOR lon = 0, 359 DO BEGIN
      FOR lat = 0, 179 DO BEGIN
        IF ( totals( week, lon, lat ) NE 0 ) THEN BEGIN
          grid( week, lon, lat ) = grid( week, lon, lat ) / $
            totals( week, lon, lat )
        ENDIF
      ENDFOR
    ENDFOR
    OPENW, output, file_name, /GET_LUN
    FOR lat = 0, 179 DO BEGIN
      printf, output, FORMAT = '(360F10.3)', grid( week, *, lat )
    ENDFOR
    FREE_LUN, output
    print, "Finished " + file_name     
END

  base = "/Net/NSCAT/people/grant/diverg_proj/"
  swathsDirectory = base + "div2grid/data/daily_swath/"
  weeklyDirectory = base + "div2grid/data/gridded/weekly/"
  
  tempFile = "temp2.div"

  miss = 99999.0
  time = LON64ARR( 18000 )
  data = FLTARR( 3, 48, 18000 )

  grid = FLTARR( 42, 360, 180 )
  totals = FLTARR( 42, 360, 180 )

  ny = 0
  nx = 0

  grid( *, *, * ) = miss
  totals( *, *, * ) = 0.0
  index = 0
  ptweek = 0  
  start = 357
  periodstart = 259
  endday = 366
  FOR year = 1997, 1997 DO BEGIN
    IF ( year EQ 1997 ) THEN BEGIN
      start = 171
      endday = 180
    ENDIF
    FOR number = start, endday DO BEGIN
      stringFormat = STRCOMPRESS( '(i' + STRING( 3 ) + '.' $
                    + STRING( 3 ) + ')', /REMOVE_ALL )

      file_name = swathsDirectory + STRTRIM( year, 1 ) + $
        STRING( number, format = stringFormat ) + '.div.gz'

      READFILE, file_name, tempFile, data, time, ny, nx, error
      IF ( year EQ 1996 ) THEN index = number - periodstart
      IF ( year EQ 1997 ) THEN index = number + 366 - periodstart

;      print, 'index = ' + STRTRIM(index, 1)
      IF ( error EQ 0 ) THEN BEGIN
        check = 0
        IF ( index EQ 287 ) THEN ptweek = 40 ELSE ptweek = FIX( index / 7 )
        min2date, time( 0 ), ptyear, ptmonth1, ptday1, pthour, ptmin
        min2date, time( ny - 1 ), ptyear, ptmonth2, ptday2, pthour, ptmin
        IF ( ptday1 NE ptday2 ) THEN check = 1

;        print, 'ptweek = ' + STRTRIM(ptweek, 1)
        FOR col = 0, nx - 1 DO BEGIN
          FOR row = 0, ny - 1 DO BEGIN
            IF ( data( 2, col, row ) NE miss ) THEN BEGIN
              grid_lat = FIX( data( 0, col, row ) + 90.5 )
              grid_lon = FIX( data( 1, col, row ) + 0.5 )
              IF ( check EQ 1 ) THEN BEGIN
                min2date, time( row ), ptyear, ptmonth, ptday, pthour, ptmin
                jday = DATE2JDAY( ptyear + 1900, ptmonth, ptday )
                IF ( ptyear + 1900 EQ 1996 ) THEN jday = jday - periodstart $
                  ELSE jday = jday + 366 - periodstart
                ptweek = FIX( jday / 7 )
                IF ( ptweek GT 40 ) THEN ptweek = 40
              ENDIF

              IF ( grid_lon GT 359 ) THEN grid_lon = 359
              IF ( grid_lat GT 179 ) THEN grid_lat = 179
              IF ( grid( ptweek, grid_lon, grid_lat ) EQ miss ) THEN BEGIN
                grid( ptweek, grid_lon, grid_lat ) = data( 2, col, row )
                totals( ptweek, grid_lon, grid_lat ) = 1.0
              ENDIF ELSE BEGIN
                grid( ptweek, grid_lon, grid_lat ) = $
                  grid( ptweek, grid_lon, grid_lat ) + data( 2, col, row )
                totals( ptweek, grid_lon, grid_lat ) = $
                  totals( ptweek, grid_lon, grid_lat ) + 1.0
              ENDELSE
            ENDIF
          ENDFOR
        ENDFOR

        print, "Finished " + file_name
      ENDIF

      IF ( ptweek NE 0 AND ptweek NE 41 ) THEN BEGIN
        print, 'ptweek = ' + STRTRIM(ptweek, 1)
        print, 'index MOD 7 = ' + STRTRIM(index MOD 7, 1)
        IF ( index MOD 7 EQ 2 ) THEN BEGIN
            ;print, 'index = ' + index
          WRITEWEEK, grid, totals, ptweek - 1, weeklyDirectory
        ENDIF
      ENDIF

    ENDFOR
  ENDFOR

  WRITEWEEK, grid, totals, 40, weeklyDirectory

END