!=============================================================================================
!
! Name:             FilterScenes
!
! Type:             F90 Program
!
! Description:      Reads in an ascii index file and specified scene files
!                   (from a namelist), and writes out new scene file
!                   based on specified filtering (and merges into 1 single file
!                   if more than 1 are read in).
!
! Subroutines:      ReadHdrScene
!                   InitHdrScene
!                   WriteHdrScene
!                   ReadScene
!                   WriteScene
!
! Required Modules: IO_Misc
!                   IO_Scene
!                   misc
!                   Consts
!                   ErrorHandling
!
!
! History:          Kevin Garrett - IMSG Inc @ NOAA/NESDIS/STAR
!______________________________________________________________________________________________

Program FilterScenes
  USE IO_Scene
  USE IO_Misc
  USE misc
  USE utils
  USE CONSTS
  USE ErrorHandling

  !---Only explicit variable declarations allowed
  IMPLICIT NONE
  !--------------------------------------------------------------------------------------------
  CHARACTER(LEN=250), DIMENSION(:), POINTER       :: SceneFiles,SceneFilesOut
  INTEGER                                         :: nSceneFiles,nProfiles,nProfs2write,profnum
  INTEGER                                         :: iProf,iFile,iLay,iu_list,iu_sceneTemp
  INTEGER                                         :: iu_nml,iu_scenein,iu_measin,err,iu_sceneout
  INTEGER                                         :: Scene_err,Meas_err,nqc=11
  INTEGER                                         :: nTotScanLines
  INTEGER,            PARAMETER                   :: AlgSN=1030

  TYPE(Scene_type)          :: SceneTemp,SceneIn,SceneOut

  !---Namelist data
  CHARACTER(LEN=150)        :: SceneFileList,SceneFileOut,pathOut
  INTEGER                   :: NewSceneFile,NewMeasFile,SfcTyp,GeoLimits,node
  REAL                      :: minLat,maxLat,minLon,maxLon
  INTEGER                   :: FovIncrem,ScanLineIncrem,clearsky
  !-------------------------------------------------------------------------------------------
  NAMELIST /FilterCntrl/pathOut,SceneFileList,SceneFileOut,NewSceneFile,GeoLimits,minLat,maxLat,minLon,maxLon,node,SfcTyp,&
       FovIncrem,ScanLineIncrem,clearsky
  read(*,NML=FilterCntrl)
  !---Read filenames to merge
  CALL ReadList(iu_list,SceneFileList,SceneFiles,nSceneFiles,SceneFilesOut,pathOut,'SCENE')

  !---Read in Scene and Measurement file headers
  nProfs2write=0
  DO iFile=1,nSceneFiles
     CALL ReadHdrScene(iu_scenein,SceneFiles(iFile),SceneIn,nProfiles)
     print *,'Profiles to read from file ',iFile, ':  ',nProfiles
     SceneCntLoop: DO iProf=1,nProfiles
        CALL ReadScene(iu_scenein,SceneIn,Scene_err)
        IF (Scene_err .eq. Warn_EndOfFile) EXIT SceneCntLoop
        IF (Scene_err .eq. Warn_readInvalid) CYCLE SceneCntLoop

        !---Filter based on ascending/descending/both
        IF (node .ge. 0 .and. SceneIn%node .ne. node) CYCLE SceneCntLoop
        
        !---Filter based on Scanline thinning
        IF (MOD(SceneIn%iScanLine,ScanLineIncrem) /= 0) CYCLE SceneCntLoop

        !---Filter based on FOV thinning
        IF (MOD(SceneIn%iScanPos,FovIncrem) .ne. 0) CYCLE SceneCntLoop

        !---Filter based on geospatial coordinates
        IF (geoLimits .eq. 1) THEN
           IF (SceneIn%lat .lt. minLat .or. SceneIn%lat .gt. maxLat .or. &
                SceneIn%lon .lt. minLon .or. SceneIn%lon .gt. maxLon) THEN
              CYCLE SceneCntLoop
           ENDIF
        ENDIF

        !---Filter based on surface type
        IF (SfcTyp .ge. 0 .and. SceneIn%iTypSfc .ne. SfcTyp) CYCLE SceneCntLoop

!        IF (SceneIn%Emiss(1) .le. 0 ) THEN
!            PRINT*,'SceneIn%Emiss(1)  =' ,iFile, iProf
!            PRINT*,SceneIn%Emiss(1)
!        ENDIF
        IF (SceneIn%Emiss(1) .le. 0 ) CYCLE SceneCntLoop
        !---Filter out clouds/ice/rain
        IF (clearsky .eq. 1 .and. &
             ColumIntegr(SceneIn%nLay,SceneIn%Pres_lay(1:SceneIn%nLay),SceneIn%SfcPress, &
             SceneIn%CLW(1:SceneIn%nLay)) .ge. 0.05 &
             .or. ColumIntegr(SceneIn%nLay,SceneIn%Pres_lay(1:SceneIn%nLay),SceneIn%SfcPress, &
             SceneIn%Graupel(1:SceneIn%nLay)) .ge. 0.05) &
             CYCLE SceneCntLoop

        nProfs2write=nProfs2write+1
     ENDDO SceneCntLoop
     close(iu_scenein)
  ENDDO
  close(iu_scenein)
  print *,'Number of profiles to process after filtering: ',nProfs2write
  CALL ReadHdrScene(iu_sceneTemp,SceneFiles(1),SceneTemp,nProfiles)
   
  close(iu_sceneTemp)
!  PRINT*,SceneTemp%nLev,SceneTemp%nLay,SceneTemp%nChan
!  PRINT*,SceneTemp%centrfreq
!  PRINT*,SceneTemp%Polarity
!  STOP
  !---Initialize and write headers of new output files
  CALL InitHdrScene(SceneTemp%nLev,SceneTemp%nLay,SceneTemp%nChan,SceneTemp%centrfreq,SceneTemp%Polarity,    &
       SceneTemp%Pres_Lev,SceneTemp%Pres_Lay,SceneOut,SceneTemp%nParmCLW,SceneTemp%nParmRain,      &
       SceneTemp%nParmSnow,SceneTemp%nParmIce,SceneTemp%nParmGrpl,SceneTemp%nAbsorb,SceneTemp%AbsorbID, &
       SceneTemp%nqc,0,SceneTemp%nPosScan,1,AlgSN)
  CALL WriteHdrScene(iu_sceneout,SceneFileOut,SceneTemp,nProfs2write)
  profnum=0

  !---Write out profiles from each Scene File
  ScenesLoop: DO iFile=1,nSceneFiles
     CALL ReadHdrScene(iu_scenein,SceneFiles(iFile),SceneIn,nProfiles)
     ProfileLoop: DO iProf=1,nProfiles
        CALL ReadScene(iu_scenein,SceneIn,Scene_err)
        IF (Scene_err .eq. Warn_EndOfFile) EXIT ProfileLoop
        IF (Scene_err .eq. Warn_readInvalid) CYCLE ProfileLoop

        !---Filter based on ascending/descending/both
        IF (node .ge. 0 .and. SceneIn%node .ne. node) CYCLE ProfileLoop

        !---Filter based on Scanline thinning
        IF (MOD(SceneIn%iScanLine,ScanLineIncrem) /= 0) CYCLE ProfileLoop

        !---Filter based on FOV thinning
        IF (MOD(SceneIn%iScanPos,FovIncrem) /= 0) CYCLE  ProfileLoop

        !---Filter based on geospatial coordinates
        IF (geoLimits .eq. 1) THEN
           IF (SceneIn%lat .lt. minLat .or. SceneIn%lat .gt. maxLat .or. &
                SceneIn%lon .lt. minLon .or. SceneIn%lon .gt. maxLon) THEN
              CYCLE ProfileLoop
           ENDIF
        ENDIF

        !---Filter based on surface type
        IF (SfcTyp .ge. 0 .and. SceneIn%iTypSfc .ne. SfcTyp) CYCLE ProfileLoop

        IF (SceneIn%Emiss(1) .le. 0 ) CYCLE ProfileLoop
        !---Filter out clouds/ice/rain
        IF (clearsky .eq. 1 .and. &
             (ColumIntegr(SceneIn%nLay,SceneIn%Pres_lay(1:SceneIn%nLay),SceneIn%SfcPress, &
             SceneIn%CLW(1:SceneIn%nLay)) .ge. 0.05 &
             .or. (ColumIntegr(SceneIn%nLay,SceneIn%Pres_lay(1:SceneIn%nLay),SceneIn%SfcPress, &
             SceneIn%Graupel(1:SceneIn%nLay)) .ge. 0.05))) &
             CYCLE ProfileLoop

        !---Copy input Scene structure to output Scene structure
        CALL AffectScene(SceneIn,SceneOut)

        !---Check for missing data below the surface and backfill
        !---(Incomplete profiles will impact covariance matrix generation)
        SfcPressLoop: DO iLay=SceneOut%nlay,1,-1
           IF (SceneOut%Temp_lay(ilay) .gt. 0) THEN
              SceneOut%Temp_lay(ilay+1:SceneOut%nlay) = SceneOut%Temp_lay(ilay)
              SceneOut%Absorb_lay(ilay+1:SceneOut%nlay,1) = SceneOut%Absorb_lay(ilay,1)
              EXIT SfcPressLoop
           ENDIF
        ENDDO SfcPressLoop

        profnum=profnum+1
        IF (profnum .eq. nProfs2write) print *,'Profiles Written: ',profnum,'/',nProfs2write
        CALL WriteScene(iu_sceneout,SceneOut)
        IF (profnum .ge. nProfs2write) EXIT ProfileLoop
     ENDDO ProfileLoop
     close(iu_scenein)
  ENDDO ScenesLoop
  close(iu_sceneout)

End Program FilterScenes