;****************************************************************************
; Author         : Josh Grant, COAPS
; Procedure Name : CALC_GRID_DIV
;
; Description : This procedure calculates the divergence of gridded vector 
;               dataset.  It is assumed that each u and v array has the
;               dimensions 360x180 (the entire globe).  The divergence is
;               calculated using the formula
;                         du/dx + dv/dy                                      
;               Where du is the wind magnitude change in the u direction, dx
;               is the change in x direction, dv is the wind magnitude change
;               in the v directions, and dy is the change is y direction.
; Parameters :
;  u     - a 360x180 array containing the wind magnitudes in the u direction.
;  v     - a 360x180 array containing the wind magnitudes in the v direction.
;  delta - the number of degrees that will be used to find the change is x
;          and y.  The actual dx and dy value used will be 2*delta.
;  miss  - the missing value in the both u and v arrays.
;  diverg- the variable that is updated by the procedure and will contain the
;          divergence values.  This program is currently set up so that this
;          variable must be initialized as a 360x180 array before for the
;          procedure is called.
;
;****************************************************************************

PRO CALC_GRID_DIV, u, v, delta, miss, diverg

  intdelta = FIX(delta)

  du = 0
  dv = 0
  dy = 0
  dx = 0
  diverg( *, * ) = miss
  umissing = WHERE( u EQ miss, umessage )
  vmissing = WHERE( v EQ miss, vmessage )
  
  FOR lon = 0, 359 DO BEGIN
    FOR lat = 0, 179 DO BEGIN
      IF ( u[ lon, lat ] NE miss AND v[ lon, lat ] NE miss ) THEN BEGIN
        IF ( lat LT intdelta OR lat GT 179 - intdelta ) THEN BEGIN
          dv = miss
        ENDIF ELSE BEGIN
          IF ( v[ lon, lat + intdelta ] NE miss AND $
               v[ lon, lat - intdelta ] NE miss ) THEN BEGIN
            dv = v[ lon, lat + intdelta ] - v[ lon, lat - intdelta ]
          ENDIF ELSE BEGIN
            dv = miss
          ENDELSE
        ENDELSE

        upper = 0
        lower = 0
        IF ( lon LT intdelta ) THEN BEGIN
          upper = lon + intdelta
          lower = 360 + (lon - intdelta)
        ENDIF ELSE BEGIN
          IF ( lon GT 359 - intdelta ) THEN BEGIN
            upper = lon + intdelta - 360
            lower = lon - intdelta
          ENDIF ELSE BEGIN
            upper = lon + intdelta
            lower = lon - intdelta
          ENDELSE
        ENDELSE
        IF ( u[ upper, lat ] NE miss AND u[ lower, lat ] NE miss ) THEN BEGIN
          du = u[ upper, lat ] - u[ lower, lat ]
        ENDIF ELSE BEGIN
          du = miss
        ENDELSE

        rad = ( ACOS( -1 ) ) / 180.0
        dy = 2 * delta * .111000
        dx = 2 * COS( ( lat - 90.0 ) * rad ) * delta * .111000

        IF ( du EQ miss OR dv EQ miss ) THEN BEGIN
          diverg[ lon, lat ] = miss
        ENDIF ELSE BEGIN
          diverg[ lon, lat ] = du / dx + dv / dy
        ENDELSE
        
      ENDIF
    ENDFOR
  ENDFOR

;  IF ( umessage GT 0 ) THEN BEGIN
;    diverg[ umissing ] = miss
;  ENDIF
;  IF ( vmessage GT 0 ) THEN BEGIN
;    diverg[ vmissing ] = miss
;  ENDIF
END

PRO READ_NSCAT_GRID, filename, u, v, err

  OPENR, input, filename, /GET_LUN, ERROR = err
  IF ( err EQ 0 ) THEN BEGIN
    temp = FLTARR( 360 )
    FOR y = 0, 179 DO BEGIN
      READF, input, FORMAT = '(360F10.3)', temp
      u( *, y ) = temp
    ENDFOR
    FOR y = 0, 179 DO BEGIN
      READF, input, FORMAT = '(360F10.3)', temp
      v( *, y ) = temp
    ENDFOR
    FREE_LUN, input
  ENDIF
  
END

PRO WRITE_GRID_DIV, file_name, diverg
  
  OPENW, output, file_name, /GET_LUN
  FOR lat = 0, 179 DO BEGIN
    printf, output, FORMAT = '(360F10.3)', diverg( *, lat )
  ENDFOR
  FREE_LUN, output
  print, "Finished " + file_name     

END