Module SSMIS_Spatial_Average_Mod
!
!
! abstract:  This routine reads BUFR format SSMIS radiance 
!            (brightness temperature) files, spatially 
!            averages the data using the AAPP averaging routines
!            and writes the data back into BUFR format
!
!
! Program history log:
!    2011-11-18   collard   - Original version (for ATMS)
!    2011-12-20   eliu      - Modify to apply for SSMIS 
! 

  use kinds, only: r_kind,r_double,i_kind
  use m_uniq
  use m_distance

  implicit none     

! Declare module level parameters
  real(r_kind), parameter    :: Missing_Value=1.e11_r_double

CONTAINS 

  SUBROUTINE SSMIS_Spatial_Average(Mype, Mype_Sub, BufrSat, Method, Num_Obs, NChanl,  &
                                   FOV, Scan , Node_InOut, Time, Lat, Lon, BT_InOut, Error_status)

    IMPLICIT NONE
    
    ! Declare passed variables
    integer(i_kind) ,intent(in   ) :: Mype,Mype_Sub,BufrSat  
    integer(i_kind) ,intent(in   ) :: Method           ! 1=simple(1)  2=simple(2) 3=AAPP 
    integer(i_kind) ,intent(in   ) :: Num_Obs, NChanl
    integer(i_kind) ,intent(in   ) :: Fov(num_obs)
    integer(i_kind) ,intent(in   ) :: Scan(num_obs)
    integer(i_kind) ,intent(inout) :: Node_InOut(num_obs)
    real(r_kind)    ,intent(in   ) :: Time(Num_Obs)
    real(r_kind)    ,intent(in   ) :: Lat(Num_Obs)
    real(r_kind)    ,intent(in   ) :: Lon(Num_Obs)
    real(r_kind)    ,intent(inout) :: BT_InOut(NChanl,Num_Obs)
    integer(i_kind) ,intent(  out) :: Error_Status

    ! Declare local parameters
!   integer(i_kind), parameter :: atms1c_h_wmosatid=224
    integer(i_kind), parameter :: lninfile=15
!   integer(i_kind), parameter :: max_fov=96
    integer(i_kind), parameter :: max_fov=60
!   integer(i_kind), parameter :: max_obs=20000000
    integer(i_kind), parameter :: as_node= 1_i_kind
    integer(i_kind), parameter :: ds_node=-1_i_kind
    real(r_kind),    parameter :: scan_interval=1.90_r_kind
    real(r_kind),    parameter :: btmin=70.0_r_kind
    real(r_kind),    parameter :: btmax=320.0_r_kind
    real(r_kind),    parameter :: btbad=-9.99e11_r_kind
    real(r_kind),    parameter :: sigma = 25.0_r_kind    
!   real(r_kind),    parameter :: sigma = 50.0_r_kind  

    ! Declare local variables
    character(100) :: infile 
    character(30)  :: Cline

    integer(i_kind) :: i,iscan,ifov,ichan,nchannels,wmosatid,version
    integer(i_kind) :: is,ip,ic 
    integer(i_kind) :: iobs,ii,jj   
    integer(i_kind) :: ios,max_scan,min_scan,mintime
    integer(i_kind) :: ns1,ns2,np1,np2 
    integer(i_kind) :: nscan       
    integer(i_kind) :: ntime_scan   
    integer(i_kind) :: delta_scan   
    integer(i_kind) :: scanline_new 
    integer(i_kind) :: ncount
    integer(i_kind) :: nxaverage(nchanl),nyaverage(nchanl)
    integer(i_kind) :: channelnumber(nchanl),qc_dist(nchanl)
    integer(i_kind) :: qc_distx(nchanl),qc_disty(nchanl)
    integer(i_kind), allocatable ::  nodeinfo(:,:)
    integer(i_kind), allocatable ::  scanline(:),scanline_back(:,:)
    integer(i_kind), allocatable ::  counter(:)
    real(r_kind), allocatable ::  dt(:)
    real(r_kind), allocatable ::  latitude(:,:), longitude(:,:)
    real(r_kind), allocatable ::  time_scan(:)
    real(r_kind), allocatable ::  dist_scan(:)
    real(r_kind), allocatable ::  dist_scan_avg(:)
    real(r_kind), allocatable ::  alat(:),alon(:),blat(:),blon(:)
    real(r_kind) :: dist_scan_sum, dlat
    real(r_kind) :: time1,time2

    real(r_kind) :: sampling_distx,sampling_disty,beamwidth(nchanl) 
    real(r_kind) :: newwidth(nchanl),cutoff(nchanl)
    real(r_kind), allocatable, target :: bt_image(:,:,:)
    real(r_kind), allocatable :: bt_image_orig(:,:,:)
    real(r_kind), pointer :: bt_image1(:,:)
    real(r_kind) :: bt_mean(nchanl)
    real(r_kind) :: num  
    real(r_kind) :: t1,t2,tdiff   
    real(r_kind) :: s1,s2   
    real(r_kind) :: lat1,lat2,lon1,lon2,dist,wgt 
    real(r_kind) :: xnum,mta   
    logical      :: pix1st,pixel_avail,gaussian_wgt,node_unknown
  
    Error_Status=0

    if (Method == 1) then  ! simple averaging 1

       gaussian_wgt = .false.
       write(*,*) 'SSMIS_Spatial_Average: using method from Banghua'
       write(*,*) 'SSMIS_Spatial_Average: bufrsat = ', BufrSat
       write(*,*) 'SSMIS_Spatial_Average: Gaussian Weighted Averaging = ', gaussian_wgt 

       ! Determine scanline from time
       !==============================
       allocate(scanline(num_obs))
       t1          = time(1)  ! time for first scanline
       nscan       = 1        ! first scanline
       scanline(1) = nscan
       do iobs = 2, num_obs
          t2    = time(iobs) 
          tdiff = t2-t1
          if (tdiff >= 0.00001_r_kind) then
             nscan = nscan+1 
             t1    = t2
          endif
          scanline(iobs) = nscan
       enddo
       max_scan = maxval(scanline)
       write(*,*) 'SSMIS_Spatial_Average: max_scan = ', max_scan
       write(*,*) 'SSMIS_Spatial_Average: max_fov  = ', max_fov
       write(*,*) 'SSMIS_Spatial_Average: nchanl   = ', nchanl 

!      Allocate and initialize variables
       allocate(bt_image(max_fov,max_scan,nchanl))
       allocate(bt_image_orig(max_fov,max_scan,nchanl))
       allocate(latitude(max_fov,max_scan))
       allocate(longitude(max_fov,max_scan))
       allocate(nodeinfo(max_fov,max_scan))
       allocate(scanline_back(max_fov,max_scan))
       bt_image(:,:,:)     = 1000.0_r_kind
       bt_image_orig(:,:,:)= 1000.0_r_kind
       latitude(:,:)       = 1000.0_r_kind 
       longitude(:,:)      = 1000.0_r_kind 
       scanline_back(:,:)  = -1
       nodeinfo(:,:)       = 1000_i_kind 

!      Put data into 2D (fov vs. scanline) array
       write(*,*) 'SSMIS_Spatial_Average:  put data into array '
       do iobs = 1, num_obs
          latitude(fov(iobs),scanline(iobs))       = lat(iobs) 
          longitude(fov(iobs),scanline(iobs))      = lon(iobs) 
          bt_image(fov(iobs),scanline(iobs),:)     = bt_inout(:,iobs)
          bt_image_orig(fov(iobs),scanline(iobs),:)= bt_inout(:,iobs)
          scanline_back(fov(iobs),scanline(iobs))  = iobs
       enddo

!      Determine AS/DS node information for each scanline
       loop1: do iscan = 1, max_scan-1
          loop2: do ifov = 1, max_fov
             if (scanline_back(ifov,iscan) > 0 .and. scanline_back(ifov,iscan+1) > 0) then
                dlat = latitude(ifov,iscan+1)-latitude(ifov,iscan)
                if (dlat < 0.0_r_kind) then
                   nodeinfo(:,iscan) = ds_node 
                else
                   nodeinfo(:,iscan) = as_node 
                endif
                cycle loop1
             endif
          enddo loop2
       enddo loop1
       nodeinfo(:,max_scan) = nodeinfo(:,max_scan-1)

!      Do spatial averaging in the box centered on each fov for each channel
       write(*,*) 'SSMIS_Spatial_Average: do spatial averaging for noise reduction '
       scan_loop: do iscan = 1, max_scan 
          fov_loop: do ifov = 1, max_fov 

!            Define grid box (3 (scan direction) x 7 (satellite track dir))
             ns1 = iscan-3          
             ns2 = iscan+3          
             if (ns1 < 1) ns1=1
             if (ns2 > max_scan) ns2=max_scan
             np1 = ifov-1          
             np2 = ifov+1          
             if (np1 < 1) np1=1
             if (np2 > max_fov) np2=max_fov

             channel_loop: do ic = 1, nchanl  
                ncount = 0_i_kind
                xnum   = 0.0_r_kind
                mta    = 0.0_r_kind
                if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. bt_image_orig(np1:np1,ns1:ns2,ic) > btmax)) then 
                   bt_image(ifov,iscan,ic) = 1000.0_r_kind 
                else
                  ! Calculate distance of each fov to the center fov 
                   box_y1: do is = ns1, ns2 
                   box_x1: do ip = np1, np2 
                      lat1 = latitude(ifov,iscan)    ! lat of the center fov
                      lon1 = longitude(ifov,iscan)   ! lon of the center fov
                      lat2 = latitude(ip,is)          
                      lon2 = longitude(ip,is)
                      dist = distance(lat1,lon1,lat2,lon2) 
                      if (dist > 100.0_r_kind) cycle box_x1  ! outside the box 
                      if (gaussian_wgt) then
                         wgt = exp(-0.5_r_kind*(dist/sigma)*(dist/sigma))
                      else
                         wgt = 1.0
                      endif
                      ncount = ncount+1 
                      xnum   = xnum+wgt
                      mta    = mta +wgt*bt_image_orig(ip,is,ic)
                   enddo box_x1
                   enddo box_y1
                   bt_image(ifov,iscan,ic) = mta/xnum 
                endif
             enddo channel_loop 

          enddo fov_loop
       enddo scan_loop

       do iscan = 1, max_scan
          do ifov = 1, max_fov
             if (scanline_back(ifov,iscan) >0) then 
                bt_inout(:,scanline_back(ifov,iscan)) = bt_image(ifov,iscan,:)
                node_inout(scanline_back(ifov,iscan)) = nodeinfo(ifov,iscan)
             endif
          enddo
       enddo

!      Deallocate arrays
       deallocate(nodeinfo,scanline,scanline_back)
       deallocate(latitude,longitude)
       deallocate(bt_image_orig,bt_image)
       nullify(bt_image1)
    endif ! Method=1

!============================================================================================================

!   Simple method 2 
    if (Method == 2) then  ! simple averaging 2 

       gaussian_wgt = .false.
       write(*,*) 'SSMIS_Spatial_Average: using method from Emily'
       write(*,*) 'SSMIS_Spatial_Average: bufrsat = ', BufrSat
       write(*,*) 'SSMIS_Spatial_Avearge: Gaussian Weighted Averaging = ', gaussian_wgt 

!      Determine scanline index from time
!      Each uniq time is one uniq scanline
!      call cpu_time(time1)
       ntime_scan = nuniq(time(1:num_obs))
       allocate(scanline(num_obs))
       allocate(time_scan(ntime_scan))
       time_scan(1:ntime_scan) = time(uniq(time(1:num_obs),ntime_scan))
       do iscan = 1, ntime_scan
          where((abs(time(1:num_obs)-time_scan(iscan)))<=0.0001_r_kind) scanline = iscan
       enddo
       max_scan = maxval(scanline)
!      call cpu_time(time2)
!      write(*,*)'CPU time for determining scanline index = ', time2-time1
!      write(*,*)'determine scanline index from time'
!      write(*,*)'ntime_scan        = ', ntime_scan
!      write(*,*)'min/max scanline  = ', minval(scanline), maxval(scanline)
!      write(*,*)'min/max time      = ', minval(time), maxval(time)
!      write(*,*)'min/max time_scan = ', minval(time_scan), maxval(time_scan)

!      Put lat/lon in 2D (fov,scanline) grid to calculate distance between scanlines
!      call cpu_time(time1)
       allocate(latitude(max_fov,max_scan))
       allocate(longitude(max_fov,max_scan))
       allocate(scanline_back(max_fov,max_scan))
       latitude(:,:)      = 9999.0_r_kind   ! filled value
       longitude(:,:)     = 9999.0_r_kind   ! filled value
       scanline_back(:,:) = -1              ! filled value
       do iobs = 1, num_obs
          latitude(fov(iobs),scanline(iobs))     = lat(iobs)
          longitude(fov(iobs),scanline(iobs))    = lon(iobs)
          scanline_back(fov(iobs),scanline(iobs))= iobs
       enddo

!      Calculate distance between scanlines
!      Reassign scanline according to distance between scanlines
       allocate(dist_scan(max_fov))
       allocate(dist_scan_avg(max_scan))
       dist_scan(:)     = 0.0_r_kind
       dist_scan_avg(:) = 0.0_r_kind
       allocate(alat(max_fov),blat(max_fov),alon(max_fov),blon(max_fov))
       scanline_new = 1
       do iscan = 2, max_scan   ! loop through scanlines
          dist_scan(:) = 0.0_r_kind
          alat(:) = latitude(:,iscan-1)
          alon(:) = longitude(:,iscan-1)
          blat(:) = latitude(:,iscan)
          blon(:) = longitude(:,iscan)
!         where(alat>=1000.0_r_kind .or. alon>=1000.0_r_kind .or. blat>=1000.0_r_kind .or. blon>=1000.0_r_kind) ! missing FOVs
          where(scanline_back(:,iscan) <= 0 .or. scanline_back(:,iscan-1) <= 0) ! missing FOVs
             alat(:)=0.0_r_kind
             blat(:)=0.0_r_kind
             alon(:)=0.0_r_kind
             blon(:)=0.0_r_kind
          end where
          dist_scan            = distance(alat,alon,blat,blon)
          dist_scan_avg(iscan) = sum(dist_scan)/count(dist_scan>0.0000001_r_kind)
          delta_scan           = nint(dist_scan_avg(iscan)/12.5_r_kind)
          scanline_new         = scanline_new+delta_scan
          do ifov = 1, max_fov
             iobs = scanline_back(ifov,iscan)
             if (iobs > 0)  scanline(iobs) = scanline_new
          enddo
       enddo
       max_scan = maxval(scanline)
       min_scan = minval(scanline)
!      call cpu_time(time2)
!      write(*,*) 'CPU time for determining dist and reassign scanlines = ', time2-time1
!      write(*,*) 'SSMIS_Spatial_Average: reassigned max_scan = ', max_scan
!      write(*,*) 'SSMIS_Spatial_Average: reassigned max_fov  = ', max_fov
!      write(*,*) 'SSMIS_Spatial_Average: reassigned nchanl   = ', nchanl

       deallocate(time_scan,dist_scan,dist_scan_avg)
       deallocate(latitude,longitude,scanline_back)
       deallocate(alat,alon,blat,blon)

!      Allocate and initialize variables
       allocate(latitude(max_fov,max_scan))
       allocate(longitude(max_fov,max_scan))
       allocate(bt_image(max_fov,max_scan,nchanl))
       allocate(bt_image_orig(max_fov,max_scan,nchanl)) 
       allocate(nodeinfo(max_fov,max_scan))
       allocate(scanline_back(max_fov,max_scan))
       bt_image(:,:,:)     = 1000.0_r_kind
       bt_image_orig(:,:,:)= 1000.0_r_kind
       latitude(:,:)       = 1000.0_r_kind
       longitude(:,:)      = 1000.0_r_kind
       scanline_back(:,:)  = -1
       nodeinfo(:,:)       = 1000_i_kind

!      Put data in 2D (fov,scanline) grid to perform noise reduction
       do i=1,num_obs
          latitude(fov(i),scanline(i))       = lat(i)
          longitude(fov(i),scanline(i))      = lon(i)
          bt_image(fov(i),scanline(i),:)     = bt_inout(:,i)
          bt_image_orig(fov(i),scanline(i),:)= bt_inout(:,i) 
          scanline_back(fov(i),scanline(i))  = i
       enddo 

!      Determine AS/DS node information for each scanline
       loop3: do iscan = 1, max_scan-1
          loop4: do ifov = 1, max_fov
             if (scanline_back(ifov,iscan) > 0 .and. scanline_back(ifov,iscan+1) > 0) then
                dlat = latitude(ifov,iscan+1)-latitude(ifov,iscan)
                if (dlat < 0.0_r_kind) then
                   nodeinfo(:,iscan) = ds_node
                else
                   nodeinfo(:,iscan) = as_node
                endif
                cycle loop3
             endif
          enddo loop4
       enddo loop3
       nodeinfo(:,max_scan) = nodeinfo(:,max_scan-1)

!      write(*,*) 'ssmis_spatial_average'
!      write(*,*) 'min/max fov           = ', minval(fov), maxval(fov)
!      write(*,*) 'min/max time          = ', minval(time), maxval(time)
!      write(*,*) 'min/max scanline_back = ', minval(scanline_back), maxval(scanline_back)
!      write(*,*) 'min/max scanline      = ', minval(scanline), maxval(scanline)
!      write(*,*) 'min/max scan          = ', minval(scan), maxval(scan)

!      Do spatial averaging in the box centered on each fov for each channel
       write(*,*) 'SSMIS_Spatial_Average: do spatial averaging for noise reduction '
       scan_loop2: do iscan = 1, max_scan
          fov_loop2: do ifov = 1, max_fov

!            Define grid box (3 (scan direction) x 7 (satellite track dir))
             ns1 = iscan-3
             ns2 = iscan+3
             if (ns1 < 1) ns1=1
             if (ns2 > max_scan) ns2=max_scan
             np1 = ifov-1
             np2 = ifov+1
             if (np1 < 1) np1=1
             if (np2 > max_fov) np2=max_fov

             channel_loop2: do ic = 1, nchanl
                ncount = 0_i_kind
                xnum   = 0.0_r_kind
                mta    = 0.0_r_kind
                if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. bt_image_orig(np1:np1,ns1:ns2,ic) > btmax)) then
                   bt_image(ifov,iscan,ic) = 1000.0_r_kind
                else
                  ! Calculate distance of each fov to the center fov
                   box_y2: do is = ns1, ns2
                   box_x2: do ip = np1, np2
                      lat1 = latitude(ifov,iscan)    ! lat of the center fov
                      lon1 = longitude(ifov,iscan)   ! lon of the center fov
                      lat2 = latitude(ip,is)
                      lon2 = longitude(ip,is)
                      dist = distance(lat1,lon1,lat2,lon2)
                      if (dist > 100.0_r_kind) cycle box_x2  ! outside the box
                      if (gaussian_wgt) then
                         wgt = exp(-0.5_r_kind*(dist/sigma)*(dist/sigma))
                      else
                         wgt = 1.0
                      endif
                      ncount = ncount+1
                      xnum   = xnum+wgt
                      mta    = mta +wgt*bt_image_orig(ip,is,ic)
                   enddo box_x2
                   enddo box_y2
                   bt_image(ifov,iscan,ic) = mta/xnum
                endif
             enddo channel_loop2

          enddo fov_loop2
       enddo scan_loop2

       do iscan = 1, max_scan
          do ifov = 1, max_fov
             if (scanline_back(ifov,iscan) >0) then 
                bt_inout(:,scanline_back(ifov,iscan)) = bt_image(ifov,iscan,:)
                node_inout(scanline_back(ifov,iscan)) = nodeinfo(ifov,iscan)
             endif
          enddo
       enddo
       deallocate(nodeinfo,scanline,scanline_back)
       deallocate(latitude,longitude)
       deallocate(bt_image_orig,bt_image)
       nullify(bt_image1)
    endif


!============================================================================================================

    if (Method == 3) then  !  AAPP method

       write(*,*) 'SSMIS_Spatial_Average: using AAPP method'
       if (bufrsat == 249) infile= 'ssmis_f16_beamwidth.txt'
       if (bufrsat == 285) infile= 'ssmis_f17_beamwidth.txt'
       if (bufrsat == 286) infile= 'ssmis_f18_beamwidth.txt'

       ! Read the beamwidth requirements
       OPEN(lninfile,file=infile,form='formatted',status='old', &
            iostat=ios)
       IF (ios /= 0) THEN
          WRITE(*,*) 'Unable to open ',trim(infile)
          Error_Status=1
          RETURN
       ENDIF
       wmosatid=999
       read(lninfile,'(a30)',iostat=ios) Cline
       DO WHILE (wmosatid .NE. BufrSat .AND. ios == 0)
          DO WHILE (Cline(1:1) == '#')
             read(lninfile,'(a30)') Cline
          ENDDO
          READ(Cline,*) wmosatid
       
          read(lninfile,'(a30)') Cline
          DO WHILE (Cline(1:1) == '#')
             read(lninfile,'(a30)') Cline
          ENDDO
          READ(Cline,*) version
       
          read(lninfile,'(a30)') Cline
          DO WHILE (Cline(1:1) == '#')
             read(lninfile,'(a30)') Cline
          ENDDO
          READ(Cline,*) sampling_distx, sampling_disty
       
          read(lninfile,'(a30)') Cline
          DO WHILE (Cline(1:1) == '#')
             read(lninfile,'(a30)') Cline
          ENDDO
          READ(Cline,*) nchannels
      
          read(lninfile,'(a30)') Cline
          if (nchannels > 0) then 
             DO ichan=1,nchannels
                read(lninfile,'(a30)') Cline
                DO WHILE (Cline(1:1) == '#')
                  read(lninfile,'(a30)') Cline
                ENDDO
                READ(Cline,*) channelnumber(ichan),beamwidth(ichan), &
                      &       newwidth(ichan),cutoff(ichan),nxaverage(ichan), &
                      &       nyaverage(ichan), qc_distx(ichan), qc_disty(ichan) 
             ENDDO
          end if
          read(lninfile,'(a30)',iostat=ios) Cline
       ENDDO
       IF (wmosatid /= BufrSat) THEN
          WRITE(*,*) 'SSMIS_Spatial_Averaging: sat id not matched in ', trim(infile)
          Error_Status=1
          RETURN
       ENDIF
       CLOSE(lninfile)
 
!orig for ATMS
!      ! Determine scanline from time
!      MinTime = MINVAL(Time)
!      ALLOCATE(Scanline(Num_Obs))
!      Scanline(:)   = NINT((Time(1:Num_Obs)-MinTime)/Scan_Interval)+1
!      Max_Scan=MAXVAL(Scanline)
!      Min_Scan=MINVAL(Scanline)

!      Determine scanline index from time
!      Each uniq time is one uniq scanline 
!      call cpu_time(time1)
       ntime_scan = nuniq(time(1:num_obs))
       allocate(scanline(num_obs))
       allocate(time_scan(ntime_scan))
       time_scan(1:ntime_scan) = time(uniq(time(1:num_obs),ntime_scan))
       do iscan = 1, ntime_scan
          where((abs(time(1:num_obs)-time_scan(iscan)))<=0.0001_r_kind) scanline = iscan
       enddo
       max_scan = maxval(scanline)
       call cpu_time(time2)
!      write(*,*)'CPU time for determining scanline index = ', time2-time1 
!      write(*,*)'determine scanline index from time' 
!      write(*,*)'ntime_scan        = ', ntime_scan
!      write(*,*)'min/max scanline  = ', minval(scanline), maxval(scanline)
!      write(*,*)'min/max time      = ', minval(time), maxval(time)
!      write(*,*)'min/max time_scan = ', minval(time_scan), maxval(time_scan)

!      Put lat/lon in 2D (fov,scanline) grid to calculate distance between scanlines

!      call cpu_time(time1)
       allocate(latitude(max_fov,max_scan))
       allocate(longitude(max_fov,max_scan))
       allocate(scanline_back(max_fov,max_scan))
       latitude(:,:)      = 9999.0_r_kind   ! filled value
       longitude(:,:)     = 9999.0_r_kind   ! filled value
       scanline_back(:,:) = -1              ! filled value
       do iobs = 1, num_obs
          latitude(fov(iobs),scanline(iobs))     = lat(iobs)
          longitude(fov(iobs),scanline(iobs))    = lon(iobs)
          scanline_back(fov(iobs),scanline(iobs))= iobs
       enddo

!      Calculate distance between scanlines
!      Reassign scanline according to distance between scanlines

       allocate(dist_scan(max_fov))
       allocate(dist_scan_avg(max_scan))
       dist_scan(:)     = 0.0_r_kind
       dist_scan_avg(:) = 0.0_r_kind
       allocate(alat(max_fov),blat(max_fov),alon(max_fov),blon(max_fov))
       scanline_new = 1
       do iscan = 2, max_scan   ! loop through scanlines
          dist_scan(:) = 0.0_r_kind
          alat(:) = latitude(:,iscan-1)
          alon(:) = longitude(:,iscan-1)
          blat(:) = latitude(:,iscan)
          blon(:) = longitude(:,iscan)
!         where(alat>=1000.0_r_kind .or. alon>=1000.0_r_kind .or. blat>=1000.0_r_kind .or. blon>=1000.0_r_kind) ! missing FOVs
          where(scanline_back(:,iscan) <= 0 .or. scanline_back(:,iscan-1) <= 0) ! missing FOVs
             alat(:)=0.0_r_kind
             blat(:)=0.0_r_kind 
             alon(:)=0.0_r_kind 
             blon(:)=0.0_r_kind
          end where
          dist_scan            = distance(alat,alon,blat,blon)
          dist_scan_avg(iscan) = sum(dist_scan)/count(dist_scan>0.0000001_r_kind)
          delta_scan           = nint(dist_scan_avg(iscan)/12.5_r_kind)
          scanline_new         = scanline_new+delta_scan
          do ifov = 1, max_fov
             iobs = scanline_back(ifov,iscan)
             if (iobs > 0)  scanline(iobs) = scanline_new 
          enddo
       enddo
       max_scan = maxval(scanline)
       min_scan = minval(scanline)
!      call cpu_time(time2)
!      write(*,*)'CPU time for determining dist and reassign scanlines = ', time2-time1 
!      write(*,*) 'SSMIS_Spatial_Average: reassigned max_scan = ', max_scan
!      write(*,*) 'SSMIS_Spatial_Average: reassigned max_fov  = ', max_fov
!      write(*,*) 'SSMIS_Spatial_Average: reassigned nchanl   = ', nchanl

       deallocate(time_scan,dist_scan,dist_scan_avg)
       deallocate(latitude,longitude,scanline_back)

!      Put lat/lon in 2D (fov,scanline) grid to perform noise reduction 
!      call cpu_time(time1)
       ALLOCATE(BT_Image(Max_FOV,Max_Scan,nchanl))
       ALLOCATE(BT_Image_Orig(Max_FOV,Max_Scan,nchanl)) 
       ALLOCATE(Scanline_Back(Max_FOV,Max_Scan))
       BT_Image(:,:,:) = 1000.0_r_kind
    
       ScanLine_Back(:,:) = -1
       DO I=1,Num_Obs
          bt_image(FOV(I),Scanline(I),:)=bt_inout(:,I)
          bt_image_orig(FOV(I),Scanline(I),:)=bt_inout(:,I) 
          Scanline_Back(FOV(I),Scanline(I))=I
       END DO

!      write(*,*) 'ssmis_spatial_average' 
!      write(*,*) 'min/max fov           = ', minval(fov), maxval(fov) 
!      write(*,*) 'min/max time          = ', minval(time), maxval(time) 
!      write(*,*) 'min/max scanline_back = ', minval(scanline_back), maxval(scanline_back) 
!      write(*,*) 'min/max scanline      = ', minval(scanline), maxval(scanline) 
!      write(*,*) 'min/max scan          = ', minval(scan), maxval(scan) 
!      write(*,*) 'min/max bt_inout        ', minval(bt_inout), maxval(bt_inout) 

!      Noise reduction using FFT

       DO IChan=1,nchanl
       
          bt_image1 => bt_image(:,:,ichan)

          ! If the channel number is present in the channelnumber array we should process it 
          ! (otherwise bt_inout just keeps the same value):
          IF (ANY(channelnumber(1:nchannels) == ichan)) THEN

             CALL MODIFY_BEAMWIDTH ( MYPE, max_fov, max_scan, bt_image1, &
                  sampling_distx, sampling_disty, beamwidth(ichan), newwidth(ichan), &
                  cutoff(ichan), nxaverage(ichan), nyaverage(ichan), &
                  qc_distx(ichan), qc_disty(ichan), IOS)
          
             IF (IOS == 0) THEN
                do iscan=1,max_scan
                   do ifov=1,max_fov
                      if (Scanline_Back(IFov, IScan) > 0) &
                           bt_inout(channelnumber(ichan),Scanline_Back(IFov, IScan)) = &
                           BT_Image1(ifov,iscan)
                   end do
                end do
             ELSE
                Error_Status=1
                RETURN
             END IF
          END IF
       END DO
       DEALLOCATE(BT_Image, Scanline, Scanline_Back)
       NULLIFY(BT_Image1)
!      call cpu_time(time2)
!      write(*,*)'CPU time for noise reduction = ', time2-time1 
    endif ! method=3
    
END Subroutine SSMIS_Spatial_Average


SUBROUTINE MODIFY_BEAMWIDTH ( MYPE, nx, ny, image, sampling_distx, sampling_disty, & 
     beamwidth, newwidth, mtfcutoff, nxaverage, nyaverage, qc_distx, qc_disty, &
     Error)
     
!-----------------------------------------
! Name: $Id: modify_beamwidth.F 222 2010-08-11 14:39:09Z frna $
!
! Purpose:
!   Manipulate the effective beam width of an image. For example, convert ATMS
!   to AMSU-A-like resolution while reducing the noise.
!
! Method:
!   1) Pad the image to a power of 2 in each dimension.
! If FFT technique is to be used then: 
!   2) Assuming Gaussian beam shapes, calcluate the input and output Modulation
!      Transfer Functions (MTF).
!   3) FFT image to frequency domain (2-D).
!   4) Multiply by output MTF divided by input MTF. If a cut-off is specified
!      (when attempting to make the beam width narrower), attenuate further
!      by an exponential function - factor of 2 at the cutoff. 
!   5) FFT back to image domain 
! Finally,
!   6) Over-write the input image, with averaging if requested.
!
! COPYRIGHT
!    This software was developed within the context of the EUMETSAT Satellite
!    Application Facility on Numerical Weather Prediction (NWP SAF), under the
!    Cooperation Agreement dated 1 December 2006, between EUMETSAT and the
!    Met Office, UK, by one or more partners within the NWP SAF. The partners
!    in the NWP SAF are the Met Office, ECMWF, KNMI and MeteoFrance.
!
!    Copyright 2010, EUMETSAT, All Rights Reserved.
!
! History:
! Version    Date     Comment
!
!  1.0   22/07/2010   N.C.Atkinson
!  1.1   21/11/2011   Convert to f90. A. Collard
!
! Code Description:
!   FORTRAN 77, following AAPP standards
!
! Declarations:


      IMPLICIT NONE
! Parameters
!     INTEGER(I_KIND), PARAMETER :: nxmax=128   !Max number of spots per scan line
!     INTEGER(I_KIND), PARAMETER :: nymax=8192  !Max number of lines. Allows 6hrs of ATMS.
      INTEGER(I_KIND), PARAMETER :: nxmax=64    !Max number of spots per scan line
      INTEGER(I_KIND), PARAMETER :: nymax=16384 !Max number of lines. Allows 6hrs of ATMS.
      REAL(R_KIND) minval, maxval
      PARAMETER (minval=0.0)   !Values less than this are treated as missing
      PARAMETER (maxval=400.0) !Values greater than this are treated as missing

! Arguments
      INTEGER(I_KIND), INTENT(IN)  :: MYPE         
      INTEGER(I_KIND), INTENT(IN)  :: nx, ny         !Size of image
      REAL(R_KIND), INTENT(INOUT)  :: image(nx,ny)   !BT or radiance image
      REAL(R_KIND), INTENT(IN)     :: sampling_distx !typically degrees
      REAL(R_KIND), INTENT(IN)     :: sampling_disty !typically degrees
      REAL(R_KIND), INTENT(IN)     :: beamwidth      !ditto
      REAL(R_KIND), INTENT(IN)     :: newwidth       !ditto
      REAL(R_KIND), INTENT(IN)     :: mtfcutoff      !0.0 to 1.0
      INTEGER(I_KIND), INTENT(IN)  :: nxaverage      !Number of samples to average (or zero)
      INTEGER(I_KIND), INTENT(IN)  :: nyaverage      !Number of samples to average (or zero)
      INTEGER(I_KIND), INTENT(IN)  :: qc_distx       !Number of samples around missing data to set to 
      INTEGER(I_KIND), INTENT(IN)  :: qc_disty       !Number of samples around missing data to set to 
      INTEGER(I_KIND), INTENT(OUT) :: Error          !Error Status
       
! Local variables
      INTEGER(I_KIND) :: nxpad, nypad, dx, dy
      INTEGER(I_KIND) :: i,j,k,ix,iy, ii, jj
      INTEGER(I_KIND) :: ifirst
      INTEGER(I_KIND) :: xpow2, ypow2
      INTEGER(I_KIND) :: nxav2, nyav2, naverage
      INTEGER(I_KIND) :: deltax, minii, maxii, minjj, maxjj
      REAL(R_KIND), ALLOCATABLE :: mtfxin(:),mtfxout(:)
      REAL(R_KIND), ALLOCATABLE :: mtfyin(:),mtfyout(:)
      REAL(R_KIND) :: mtfin,mtfout,mtf_constant     
      REAL(R_KIND) :: mtfx_constant, mtfy_constant   !test
      REAL(R_KIND), ALLOCATABLE :: mtfpad(:,:)
      REAL(R_KIND), ALLOCATABLE :: imagepad(:,:)
      REAL(R_KIND), ALLOCATABLE :: work(:)
      REAL(R_KIND) :: f,df,factor
      REAL(R_KIND) :: PI, LN2, LNcsquared
      LOGICAL :: missing
      LOGICAL, ALLOCATABLE :: gooddata_map(:,:)


! End of declarations
!-----------------------------------------
      
      PI = 4.0*atan(1.0)
      LN2 = LOG(2.0)
!     MTF_Constant=-(PI/(2*sampling_dist))**2/LN2
      MTFX_Constant=-(PI/(2*sampling_distx))**2/LN2   !test
      MTFY_Constant=-(PI/(2*sampling_disty))**2/LN2   !test 
      IF (mtfcutoff .GT. 0.0) LNcsquared = LOG(mtfcutoff)**2
      nxav2 = nxaverage/2
      nyav2 = nyaverage/2
      naverage = nxaverage*nyaverage
      Error = 0

!1) Pad the image up to the nearest power of 2 in each dimension, by reversing
!the points near the edge.

      xpow2 = INT(LOG(nx*1.0)/LN2 + 1.0)
      ypow2 = INT(LOG(ny*1.0)/LN2 + 1.0)
      nxpad = 2**xpow2
      nypad = 2**ypow2
      dx = (nxpad - nx)/2
      dy = (nypad - ny)/2

      IF (nxpad .GT. nxmax) THEN
         write(*,*) 'SSMIS_Spatial_Average: nx too large, maximum allowed value is ',nxmax-1
         Error = 1
         RETURN
      END IF
      
      IF (nypad .GT. nymax) THEN
         write(*,*) 'SSMIS_Spatial_Average: ny value is ', nypad   !test 
         write(*,*) 'SSMIS_Spatial_Average: ny too large, maximum allowed value is ',nymax-1
         Error = 1
         RETURN
      END IF

      ALLOCATE(mtfxin(nxpad),mtfxout(nxpad))
      ALLOCATE(mtfyin(nypad),mtfyout(nypad))
      ALLOCATE(mtfpad(nxpad,nypad))
      ALLOCATE(imagepad(nxpad,nypad))
      ALLOCATE(work(nypad))
      ALLOCATE(gooddata_map(nxmax,nymax))

!Loop over scan positions
      DO j=dy+1,dy+ny
        DO i=dx+1,dx+nx
          if (image(i-dx,j-dy) > maxval .OR. image(i-dx,j-dy) < minval) &
               image(i-dx,j-dy) = minval - 1.0_r_kind
          imagepad(i,j) = image(i-dx,j-dy)   !Take a copy of the input data
          gooddata_map(i,j) = .TRUE.   ! Initialised for step 6)
!111      format(4i12,2f20.5,l6)
        ENDDO

!Interpolate missing points in the along-track direction

        ifirst = -1
        missing = .false.
        
        DO i=dx+1,dx+nx
          IF (.not.missing) THEN
            IF (imagepad(i,j) .GE. minval) THEN   !imagepad = -1 indicates missing data 
              ifirst = i
            ELSE
              missing = .true.
            ENDIF
          ELSE
            IF (imagepad(i,j) .GE. minval) THEN  !First good point after missing
               missing = .false.
               IF (ifirst .eq. -1) THEN
                  DO k=dx+1,i-1
                     imagepad(k,j) = imagepad(i,j)      !Constant
                  ENDDO
               ELSE
                  DO k=ifirst+1,i-1
                     factor = (i-k)*1.0/(i-ifirst)      !Interpolate
                     imagepad(k,j) = imagepad(ifirst,j)*factor + &
                          imagepad(i,j)*(1.0-factor)
                  ENDDO
               ENDIF
            ENDIF
          ENDIF
        ENDDO
        IF (missing) THEN         !Last scan is missing
          IF (ifirst .GE. 1) then
            DO k=ifirst+1,dx+nx
              imagepad(k,j) = imagepad(ifirst,j)     !Constant
            ENDDO
          ENDIF
        ENDIF          

!Continue padding the edges

        DO i=1,dx
          imagepad(i,j) = imagepad(dx+dx+2-i,j)
        ENDDO
        DO i=nx+dx+1,nxpad
          imagepad(i,j) = imagepad(nx+dx+nx+dx-i,j)
        ENDDO
     ENDDO

     DO j=1,dy
        DO i=1,nxpad
           imagepad(i,j) = imagepad(i,dy+dy+2-j)
        ENDDO
     ENDDO
     DO j=ny+dy+1,nypad
        DO i=1,nxpad
           imagepad(i,j) = imagepad(i,ny+dy+ny+dy-j)
        ENDDO
     ENDDO

!2) Compute the MTF modifications. Assume beams are Gaussian.

      IF (newwidth .GT. 0) THEN
        df = 1.0/nxpad
        DO i=1,nxpad/2+1
          f = df*(i-1)      !DC to Nyquist
!         mtfxin(i) = exp(MTF_Constant*(f*beamwidth)**2)
!         mtfxout(i) = exp(MTF_Constant*(f*newwidth)**2)
          mtfxin(i) = exp(MTFX_Constant*(f*beamwidth)**2)    !test 
          mtfxout(i) = exp(MTFX_Constant*(f*newwidth)**2)    !test
          IF (i.GT.1.AND.i.LT.nxpad/2+1) THEN
            mtfxin(nxpad-i+2) = mtfxin(i)
            mtfxout(nxpad-i+2) = mtfxout(i)
          ENDIF
        ENDDO
        df = 1.0/nypad
        DO i=1,nypad/2+1
          f = df*(i-1)      !DC to Nyquist
!         mtfyin(i) = exp(MTF_Constant*(f*beamwidth)**2)
!         mtfyout(i) = exp(MTF_Constant*(f*newwidth)**2)
          mtfyin(i) = exp(MTFY_Constant*(f*beamwidth)**2)  !test
          mtfyout(i) = exp(MTFY_Constant*(f*newwidth)**2)  !test
          IF (i.GT.1.AND.i.LT.nypad/2+1) THEN
            mtfyin(nypad-i+2) = mtfyin(i)
            mtfyout(nypad-i+2) = mtfyout(i)
          ENDIF
        ENDDO
        DO i=1,nxpad
          DO j=1,nypad
            mtfin = mtfxin(i)*mtfyin(j)
            mtfout = mtfxout(i)*mtfyout(j)
            if (mtfcutoff .GT. 0.0) THEN
              mtfpad(i,j) = (mtfout * &
                exp(-LN2/LNcsquared*(LOG(mtfout))**2))/mtfin
            else
              mtfpad(i,j) = mtfout/mtfin
            endif
          ENDDO
        ENDDO

!3) Fourier transform, line by line then column by column.
!After each FFT, points 1 to nxpad/2+1 contain the real part of the spectrum,
!the rest contain the imaginary part in reverse order.

        DO j=1,nypad
           CALL SFFTCF(imagepad(:,j),nxpad,xpow2)
        ENDDO

        DO i=1,nxpad
           DO j=1,nypad
              work(j) = imagepad(i,j)
           ENDDO
           CALL SFFTCF(work,nypad,ypow2)
           DO j=1,nypad
              imagepad(i,j) = work(j)
           ENDDO
        ENDDO

!4) Multiply the spectrum by the MTF factor

        DO j=1,nypad
           DO i=1,nxpad
            imagepad(i,j) = imagepad(i,j)*mtfpad(i,j)
          ENDDO
        ENDDO

!5) Inverse Fourier transform, column by column then line by line 

        DO i=1,nxpad
          DO j=1,nypad
            work(j) = imagepad(i,j)
          ENDDO
          CALL SFFTCB(work,nypad,ypow2)
          DO j=1,nypad
            imagepad(i,j) = work(j)
          ENDDO
        ENDDO

        DO j=1,nypad
          CALL SFFTCB(imagepad(:,j),nxpad,xpow2)
        ENDDO
     ENDIF   !New width is specified

!6) Reset missing values in gooddata_map, based on qc_dist and the values 
!   in the input image array

     ! Set the ends of the image to missing in the along track direction
     ! (doing the same across track will remove too much data)
!    gooddata_map(:,1:qc_dist)=.FALSE.
!    gooddata_map(:,ny-qc_dist+1:ny)=.FALSE.
     gooddata_map(:,1:qc_disty)=.FALSE.
     gooddata_map(:,ny-qc_disty+1:ny)=.FALSE.
     
     DO j=1,ny
        DO i=1,nx
           IF (image(i,j) <= minval) THEN
!             minjj=max(j+dy-qc_dist,0)
!             maxjj=min(j+dy+qc_dist,nymax)
              minjj=max(j+dy-qc_disty,0)
              maxjj=min(j+dy+qc_disty,nymax)
              DO jj=minjj,maxjj
!                deltax=INT(SQRT(REAL(qc_dist**2 - (jj-j-dy)**2 )))
!                minii=max(i+dx-deltax,0)
!                maxii=min(i+dx+deltax,nxmax)
                 minii=max(i+dx-qc_distx,0)
                 maxii=min(i+dx+qc_distx,nxmax)
                 DO ii=minii,maxii
                    gooddata_map(ii,jj)=.FALSE.
!345                 format(10i10)    
                 END DO
              END DO
           END IF
        END DO
     END DO

!7) Over-write the input image (points that are not missing)

     DO j=1,ny
        DO i=1,nx
           IF (gooddata_map(i+dx,j+dy)) THEN
              IF (nxav2 == 0. .AND. nyav2 == 0) THEN
                 image(i,j) = imagepad(i+dx,j+dy)
              ELSE
                 image(i,j) = 0.0             !Do averaging
                 DO ix = -nxav2,nxav2
                    DO iy = -nyav2,nyav2
                       image(i,j) = image(i,j) + imagepad(i+dx+ix,j+dy+iy)
                    ENDDO
                 ENDDO
                 image(i,j) = image(i,j)/naverage
              ENDIF
           ELSE
              image(i,j) = missing_value
           END IF
        ENDDO
     ENDDO

!8) Deallocate arrays

     DEALLOCATE(mtfxin,mtfxout)
     DEALLOCATE(mtfyin,mtfyout)
     DEALLOCATE(mtfpad)
     DEALLOCATE(imagepad)
     DEALLOCATE(work)
     DEALLOCATE(gooddata_map)

     RETURN
   END SUBROUTINE MODIFY_BEAMWIDTH

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  A real-valued, in place, split-radix FFT program
!  Real input and output in data array X
!  Length is N = 2 ** M
!  Decimation-in-time, cos/sin in second loop
!  Output in order:
!         [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
!
!  This FFT computes
!     X(k) = sum_{j=0}^{N-1} x(j)*exp(-2ijk*pi/N)
!
!
!  H.V. Sorensen, Rice University, Oct. 1985
!
!  Reference:  H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus;
!              Real-Valued Fast Fourier Transform Algorithms; IEEE
!              Trans. Acoust., Speech, Signal Process.; Vol ASSP-35,
!              June 1987, pp. 849-863.
!
!  This code was originally named RVFFT.
!
!  History:
!   21/11/2011 Converted to something resembling f90.   A.Collard
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      SUBROUTINE SFFTCF( X, N, M )

      IMPLICIT NONE

! ... Parameters ...
      REAL(R_KIND), PARAMETER :: SQRT2 = 1.4142135623730950488
      REAL(R_KIND), PARAMETER :: TWOPI = 6.2831853071795864769 

! ... Scalar arguments ...
      INTEGER(I_KIND), INTENT(IN) :: N, M
! ... Array arguments ...
      REAL(R_KIND), INTENT(INOUT) ::  X(N)
! ... Local scalars ...
      INTEGER(I_KIND)  J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8
      INTEGER(I_KIND)  N1, N2, N4, N8
      REAL(R_KIND)  XT, R1, T1, T2, T3, T4, T5, T6
      REAL(R_KIND)  A, A3, E, CC1, SS1, CC3, SS3
!
! ... Exe. statements ...
!
      IF ( N .EQ. 1 ) RETURN
!
 100  J = 1
      N1 = N - 1
      DO 104, I = 1, N1
         IF ( I .GE. J ) GOTO 101
         XT = X(J)
         X(J) = X(I)
         X(I) = XT
 101     K = N / 2
 102     IF ( K .GE. J ) GOTO 103
            J = J - K
            K = K / 2
            GOTO 102
 103     J = J + K
 104  CONTINUE
! 
      IS = 1
      ID = 4
 70   DO 60, I0 = IS, N, ID
         I1 = I0 + 1
         R1 = X(I0)
         X(I0) = R1 + X(I1)
         X(I1) = R1 - X(I1)
 60   CONTINUE
      IS = 2 * ID - 1
      ID = 4 * ID
      IF ( IS .LT. N ) GOTO 70
!
      N2 = 2
      DO 10, K = 2, M
         N2 = N2 * 2
         N4 = N2 / 4
         N8 = N2 / 8
         E = TWOPI / N2
         IS = 0
         ID = N2 * 2
 40      DO 38, I = IS, N-1, ID
            I1 = I + 1
            I2 = I1 + N4
            I3 = I2 + N4
            I4 = I3 + N4
            T1 = X(I4) + X(I3)
            X(I4) = X(I4) - X(I3)
            X(I3) = X(I1) - T1
            X(I1) = X(I1) + T1
            IF ( N4 .EQ. 1 ) GOTO 38
            I1 = I1 + N8
            I2 = I2 + N8
            I3 = I3 + N8
            I4 = I4 + N8
            T1 = ( X(I3) + X(I4) ) / SQRT2
            T2 = ( X(I3) - X(I4) ) / SQRT2
            X(I4) = X(I2) - T1
            X(I3) = - X(I2) - T1
            X(I2) = X(I1) - T2
            X(I1) = X(I1) + T2
 38      CONTINUE
         IS = 2 * ID - N2
         ID = 4 * ID
         IF ( IS .LT. N ) GOTO 40
         A = E
         DO 32, J = 2, N8
            A3 = 3 * A
            CC1 = COS(A)
            SS1 = SIN(A)
            CC3 = COS(A3)
            SS3 = SIN(A3)
            A = J * E
            IS = 0
            ID = 2 * N2
 36         DO 30, I = IS, N-1, ID
               I1 = I + J
               I2 = I1 + N4
               I3 = I2 + N4
               I4 = I3 + N4
               I5 = I + N4 - J + 2
               I6 = I5 + N4
               I7 = I6 + N4
               I8 = I7 + N4
               T1 = X(I3) * CC1 + X(I7) * SS1
               T2 = X(I7) * CC1 - X(I3) * SS1
               T3 = X(I4) * CC3 + X(I8) * SS3
               T4 = X(I8) * CC3 - X(I4) * SS3
               T5 = T1 + T3
               T6 = T2 + T4
               T3 = T1 - T3
               T4 = T2 - T4
               T2 = X(I6) + T6
               X(I3) = T6 - X(I6)
               X(I8) = T2
               T2 = X(I2) - T3
               X(I7) = - X(I2) - T3
               X(I4) = T2
               T1 = X(I1) + T5
               X(I6) = X(I1) - T5
               X(I1) = T1
               T1 = X(I5) + T4
               X(I5) = X(I5) - T4
               X(I2) = T1
 30         CONTINUE
            IS = 2 * ID - N2
            ID = 4 * ID
            IF ( IS .LT. N ) GOTO 36
 32      CONTINUE
 10   CONTINUE
      RETURN
!
! ... End of subroutine SFFTCF ...
!
   END SUBROUTINE SFFTCF

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  A real-valued, in place, split-radix IFFT program
!  Hermitian symmetric input and real output in array X
!  Length is N = 2 ** M
!  Decimation-in-frequency, cos/sin in second loop
!  Input order:
!         [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
!
!  This FFT computes
!     x(j) = (1/N) * sum_{k=0}^{N-1} X(k)*exp(2ijk*pi/N)
!
!
!  H.V. Sorensen, Rice University, Nov. 1985
!
!  Reference:  H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus;
!              Real-Valued Fast Fourier Transform Algorithms; IEEE
!              Trans. Acoust., Speech, Signal Process.; Vol ASSP-35,
!              June 1987, pp. 849-863.
!
!  This code was originally named IRVFFT.
!
!  History:
!   21/11/2011 Converted to something resembling f90.   A.Collard
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      SUBROUTINE SFFTCB( X, N, M )

      use kinds, only: r_kind,r_double,i_kind

      IMPLICIT NONE

! ... Parameters ...
      REAL(R_KIND), PARAMETER :: SQRT2 = 1.4142135623730950488
      REAL(R_KIND), PARAMETER :: TWOPI = 6.2831853071795864769 

! ... Scalar arguments ...
      INTEGER(I_KIND), INTENT(IN) :: N, M
! ... Array arguments ...
      REAL(R_KIND), INTENT(INOUT) ::  X(N)
! ... Local scalars ...
      INTEGER(I_KIND)  J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8
      INTEGER(I_KIND)  N1, N2, N4, N8
      REAL(R_KIND)  XT, R1, T1, T2, T3, T4, T5
      REAL(R_KIND)  A, A3, E, CC1, SS1, CC3, SS3
!
! ... Exe. statements ...
!
      IF ( N .EQ. 1 ) RETURN
!
      N2 = 2 * N
      DO 10, K = 1, M-1
         IS = 0
         ID = N2
         N2 = N2 / 2
         N4 = N2 / 4
         N8 = N4 / 2
         E = TWOPI / N2
 17      DO 15, I = IS, N-1, ID
            I1 = I + 1
            I2 = I1 + N4
            I3 = I2 + N4
            I4 = I3 + N4
            T1 = X(I1) - X(I3)
            X(I1) = X(I1) + X(I3)
            X(I2) = 2 * X(I2)
            X(I3) = T1 - 2 * X(I4)
            X(I4) = T1 + 2 * X(I4)
            IF ( N4 .EQ. 1 ) GOTO 15
            I1 = I1 + N8
            I2 = I2 + N8
            I3 = I3 + N8
            I4 = I4 + N8
            T1 = ( X(I2) - X(I1) ) / SQRT2
            T2 = ( X(I4) + X(I3) ) / SQRT2
            X(I1) = X(I1) + X(I2)
            X(I2) = X(I4) - X(I3)
            X(I3) = 2 * ( - T2 - T1 )
            X(I4) = 2 * ( -T2 + T1 )
 15      CONTINUE
         IS = 2 * ID - N2
         ID = 4 * ID
         IF ( IS .LT. N-1 ) GOTO 17
         A = E
         DO 20, J = 2, N8
            A3 = 3 * A
            CC1 = COS(A)
            SS1 = SIN(A)
            CC3 = COS(A3)
            SS3 = SIN(A3)
            A = J * E
            IS = 0
            ID = 2 * N2
 40         DO 30, I = IS, N-1, ID
               I1 = I + J
               I2 = I1 + N4
               I3 = I2 + N4
               I4 = I3 + N4
               I5 = I + N4 - J + 2
               I6 = I5 + N4
               I7 = I6 + N4
               I8 = I7 + N4
               T1 = X(I1) - X(I6)
               X(I1) = X(I1) + X(I6)
               T2 = X(I5) - X(I2)
               X(I5) = X(I2) + X(I5)
               T3 = X(I8) + X(I3)
               X(I6) = X(I8) - X(I3)
               T4 = X(I4) + X(I7)
               X(I2) = X(I4) - X(I7)
               T5 = T1 - T4
               T1 = T1 + T4
               T4 = T2 - T3
               T2 = T2 + T3
               X(I3) = T5 * CC1 + T4 * SS1
               X(I7) = - T4 * CC1 + T5 * SS1
               X(I4) = T1 * CC3 - T2 * SS3
               X(I8) = T2 * CC3 + T1 * SS3
 30         CONTINUE
            IS = 2 * ID - N2
            ID = 4 * ID
            IF ( IS .LT. N-1 ) GOTO 40
 20      CONTINUE
 10   CONTINUE
!
      IS = 1
      ID = 4
 70   DO 60, I0 = IS, N, ID
         I1 = I0 + 1
         R1 = X(I0)
         X(I0) = R1 + X(I1)
         X(I1) = R1 - X(I1)
 60   CONTINUE
      IS = 2 * ID - 1
      ID = 4 * ID
      IF ( IS .LT. N ) GOTO 70
!
 100  J = 1
      N1 = N - 1
      DO 104, I = 1, N1
         IF ( I .GE. J ) GOTO 101
         XT = X(J)
         X(J) = X(I)
         X(I) = XT
 101     K = N / 2
 102     IF ( K .GE. J ) GOTO 103
            J = J - K
            K = K / 2
            GOTO 102
 103     J = J + K
 104  CONTINUE
      XT = 1.0 / FLOAT( N )
      DO 99, I = 1, N
         X(I) = XT * X(I)
 99   CONTINUE
      RETURN
!
! ... End of subroutine SFFTCB ...
! 
      END SUBROUTINE SFFTCB

  SUBROUTINE realdistance(latin1,lonin1,latin2,lonin2,dist)
  !========================================================================================================
  !
  !  Purpose:
  !    To calculate geophysical distance of two points
  !
  !  Record of revisions:
  !     YYYY/MM/DD      Programmer                       Description of change
  !    ============   ==============   ===========================================================
  !     2007/03/01     Banghua Yan       Create orginal subroutine
  !                    (NOAA/NESDIS)
  !     2009/10/22     Banghua Yan       Implement to GSI package
  !
  !========================================================================================================

    use kinds, only: r_kind,r_double,i_kind
    IMPLICIT NONE

    ! Declare subroutine arguments
    REAL(r_kind) :: latin1,lonin1,latin2,lonin2
    REAL(r_kind) :: dist
    ! Declare local variables
    REAL(r_kind) :: lat1,lon1,lat2,lon2
    REAL(r_kind) :: PI,earth_radius,temp

    PI=3.14159_r_kind
    earth_radius=6378.0_r_kind
    lat1=latin1*PI/180.0_r_kind
    lon1=lonin1*PI/180.0_r_kind
    lat2=latin2*PI/180.0_r_kind
    lon2=lonin2*PI/180.0_r_kind
    temp=SIN(lat1)*SIN(lat2) + COS(lat1)*COS(lat2)*COS(lon2-lon1)
    temp=ACOS(temp)
    temp=temp*earth_radius
    IF (temp < 0.) temp=0.0_r_kind
    dist = temp
    !B.YAN (COR. 06/24/2008)
    IF (ABS(lat1-lat2) <= 0.001_r_kind .AND. ABS(lon1-lon2) <= 0.001_r_kind) temp = 0.0_r_kind

  END SUBROUTINE realdistance




END MODULE SSMIS_Spatial_Average_Mod
