!$$$ subprogram documentation block
!               .      .    .
! subprogram:	valid			validate the obs and penalty values
!     prgmmr:   safford			date:  2009-12
!
! abstract:	This module contains code to read a given satellite's 
!               base file and then validate new obs(count) and penalty
!               values by comparing them to the baseline values.  
!
! program history log:
!	2009-12-07 safford - initial coding
!
! contains:
!
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$

module valid

  implicit none
  
  private

! --- module routines
  public :: load_base
  public :: validate_count
  public :: validate_penalty 

! --- module parameters
  integer, parameter  :: funit = 17
  real,parameter :: rmiss      = -999.0

! --- module vars
  logical :: base_loaded = .FALSE.
  integer :: nchan, nregion, j, k, dummy

  real,allocatable,dimension(:,:)::  avg_count, sdv_count 
  real,allocatable,dimension(:,:)::  min_count, max_count 
  real,allocatable,dimension(:,:)::  avg_penalty, sdv_penalty 
  real,allocatable,dimension(:,:)::  min_penalty, max_penalty 

  integer,allocatable,dimension(:):: nu_chan
 
  contains

    !--------------------------------------------------------------------------
    !  load the base file for the given satellite
    !
    !  Note a change in validation for iasi and cris sources.
    !  Per Kristen (2/20/20):
    !   With correlated error, the contribution to penalty calculation in the
    !   radmon is no longer mathematically correct. The calculated penalty is
    !   expected to be larger than what it actually is. Both IASI instruments 
    !   will use correlated error over land and sea, and both CrIS instruments 
    !   will use correlated error over sea in V16. You should therefore 
    !   suppress/modify the warnings for IASI and CrIS. 
    !
    !  The agreed upon solution is to multiply the historical max penalty
    !  values, making iasi 5x and cris 6x.
    !---------------------------------------------------------------------------

    subroutine load_base( satname, iret )
     

      !--- interface
      character(20), intent( in )	:: satname
      integer, intent( out )		:: iret     

      !--- variables
      character(20) fname
      character(20) test_satname
      character(10) base_date

      integer fios, multiply_by
      integer chan, region

      logical fexist 


      !--- initialization
      iret   = -1
      fname  = trim(satname) // '.base'
      fexist = .FALSE.

      write(*,*) '--> valid, satname, fname = ', satname, fname

      !--- verify file exists and open the file
      inquire( file = fname, exist = fexist )
      if( fexist .eqv. .FALSE. ) then
         fios = -1 
      else 
         open( UNIT=funit, FILE=fname, IOSTAT=fios )
         write(*,*) ' fios = ', fios
      end if

      if( fios == 0 ) then
         !--- read the file header
         read(funit,*) test_satname, base_date, nchan, nregion
         write(*,*) ' test_satname, base_date, nchan, nregion = ', test_satname, base_date, nchan, nregion
        
         allocate( avg_count(nchan,nregion), sdv_count(nchan,nregion), &
              avg_penalty(nchan,nregion), sdv_penalty(nchan,nregion), &
              min_count(nchan,nregion), max_count(nchan,nregion), &
              min_penalty(nchan,nregion), max_penalty(nchan,nregion) )

         ! --- set all missing
         do k=1,nregion
            do j=1,nchan
               avg_count(j,k) = rmiss
               sdv_count(j,k) = rmiss
               min_count(j,k) = rmiss
               max_count(j,k) = rmiss
               avg_penalty(j,k) = rmiss
               sdv_penalty(j,k) = rmiss
               min_penalty(j,k) = rmiss
               max_penalty(j,k) = rmiss
            end do
         end do

         do k=1,nregion
            do j=1,nchan
               read(funit,*) chan, region, avg_count(j,k), sdv_count(j,k), min_count(j,k), max_count(j,k), &
                             avg_penalty(j,k), sdv_penalty(j,k), min_penalty(j,k), max_penalty(j,k)
            end do
         end do

         iret = 0 
         base_loaded = .TRUE.
      else
         write(*,*) 'unable to load fname for data error checking'
      end if

    end subroutine load_base


    !---------------------------------------------------------------
    !  validate a count
    !     given a count value for a channel and region, determine
    !     if the count is within +/- 2*sdv 
    !
    !     iret         0 = normal
    !                 -1 = invalid channel
    !                 -2 = invalid region
    !                  1 = base file wasn't loaded, unable to validate
    !---------------------------------------------------------------
    subroutine validate_count( channel, region, count, valid, avg_cnt, iret )

      !--- interface
      integer, intent( in )		:: channel
      integer, intent( in )		:: region
      real, intent( in )                :: count
      logical, intent( out )            :: valid
      real, intent( out )               :: avg_cnt
      integer, intent( out )		:: iret

      !--- vars
      real cnt, hi, lo, sdv2

      write(*,*) '--> validate_count, channel, region, count ', channel, region, count
      !--- initialize vars
      iret = 0
      cnt = count
      valid = .FALSE.
      avg_cnt = 0.00 

      if( base_loaded .eqv. .TRUE. ) then

         if( channel < 1 .OR. channel > nchan ) then
            iret = -1
            write(*,*) 'Warning:  In validate_count attempt to validate channel out of range', channel
            valid = .TRUE.
         else if( region < 1 .OR. region > nregion ) then
            iret = -2
            write(*,*) 'Warnig:  In validate_count attempt to validate region out of range', region
            valid = .TRUE.
         else
            !---------------------------------------------------------------------
            !  all unassimilated channels in the base files will have an rmiss
            !  value and are considered valid for verification purposes
            !
            avg_cnt = avg_count( channel, region )
         
            if( avg_count(channel,region) < 0.0 ) then
               valid = .TRUE.
            else
               !------------------------------------------------------------------
               !  Consider any count valid if:
               !    cnt is within 2 sdvs from avg 
               !
               sdv2 = 2 * sdv_count( channel, region )
!               lo = avg_count(channel,region) - sdv2
               lo = avg_count(channel,region) * 0.67

               if( cnt >= lo ) then
                  valid = .TRUE.
               end if

 	       if( valid .eqv. .FALSE. ) then
                  write(*,*) 'LOW COUNT:  cnt, lo, min_count, avg_cnt = ', & 
			 cnt, lo, min_count(channel,region), avg_cnt
               end if

            end if

         end if


      else 
         !---------------------------------------------------------
         !  base file wasn't loaded, so validation wasn't possible 
         !
         write(*,*) 'base file not loaded, unable to validate count'
         iret = 1 
      end if 

      write (*,*) '<-- valid, iret=', valid, iret
    end subroutine validate_count


    !-------------------------------------------------------------
    !  validate a penalty value
    !     given a penalty value for a channel and region, determine
    !     if the penalty is within +/- 2*sdv 
    !
    !     iret         0 = normal
    !                 -1 = invalid channel
    !                 -2 = invalid region
    !-------------------------------------------------------------
    subroutine validate_penalty( channel, region, penalty, valid, bound, iret )

      !--- interface
      integer, intent( in )		:: channel
      integer, intent( in )		:: region
      real, intent( in )                :: penalty
      logical, intent( out )            :: valid
      real, intent( out )               :: bound
      integer, intent( out )		:: iret

      !--- vars
      real hi, lo, sdv2

      write(*,*) '--> validate_penalty, channel, region, penalty ', channel, region, penalty

      !--- initialize vars
      iret = 0 
      valid = .FALSE.
      bound = rmiss

      if( base_loaded .eqv. .TRUE. .AND. nchan > 1 ) then
         if( channel < 1 .OR. channel > nchan ) then
            iret = -1
            write(*,*) 'Warning:  In validate_penalty attempt to validate channel out of range', channel
            valid = .TRUE.
         else if( region < 1 .OR. region > nregion ) then
            iret = -2
            write(*,*) 'Warning:  In validate_penalty attempt to validate region out of range', region
            valid = .TRUE.
         else
            !
            !  all unassimilated channels in the base files will have an rmiss
            !  value and are considered valid for verification purposes
            !
            bound = max_penalty( channel,region )

            if( avg_penalty(channel,region) < 0.0 ) then
               valid = .TRUE.
            else
               sdv2 = 2 * sdv_penalty( channel, region )
               hi = avg_penalty(channel,region) + sdv2
               lo = avg_penalty(channel,region) - sdv2

               !
               !  Consider any penalty value valid if:
               !    penalty is 2 sdvs from avg or
               !    penalty is greater than the established max range for 
               !      chan,region by 20% or more --> idea here is to stage
               !      two levels of warning, say a 20% exceedence and then
               !      maybe a 25% or 30% exceedence for a yellow and red
               !      level warning.
               !
               if( (penalty >= 0.0) .AND. penalty <= hi ) then
                  valid = .TRUE.
               else if( (penalty > 0.0) .AND. &
                   (penalty <= (max_penalty( channel,region )*1.2) )  ) then
                  valid = .TRUE.
               end if
  
            end if
         end if

         if ( valid .eqv. .FALSE. ) then
            write(*,*) ' BAD:  penalty, avg_penalty(channel,region), sdv2, hi, lo, bound = ', penalty, avg_penalty(channel,region), sdv2, hi, lo, bound
         end if
         write (*,*) '<-- valid, iret=', valid, iret
      else 
         !--- base file was not loaded, or nchan was 0 so return 
         !--- a warning that validation isn't possible
         write (*,*) 'Warning:  base file not loaded or nchan < 1, nchan= ', nchan
         iret = 1 
      end if 
    end subroutine validate_penalty

    
end module valid