ecearth3/sources/tm5mp/proj/domino/MStripeCorrection.F90

! First include the set of model-wide compiler flags 
#include "tm5.inc"

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN if (status/=0) then; TRACEBACK; return; end if
#define IF_ERROR_RETURN if (status> 0) then; TRACEBACK; return; end if
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if

! First include the set of model-wide compiler flags 
! #include "tm5.inc"

module MStripeCorrection

  ! Apply stripe correction for OMI/TROPOMI
  !
  ! Steps: 
  !  1. Read the stripe correction file, which contains:
  !     - the mean stripe correction over the past period (roughly 1 week)
  !     - the latest strip correction of the previous day
  !  2. Apply this to the slant column NO2 DOAS results
  !  3. After the vertical column has been computed:
  !     - Determine if the orbit is over the pacific
  !     - If yes, then update the stripe correction, and write to file
  !     - (Apply the new stripe correction to the SCD, VCD)
  !

  !
  !
  !                             Henk Eskes, KNMI, March-April 2016

  use GO,           only  :  gol, goErr
  use MTObsTrack,   only  :  TObsTrack
  use datetime,     only  :  date2tau, tau2date

  ! netcdf input/output makes use of the MDF module
  use MDF

  implicit none

  private

  public  ::  ReadStripeCorrection, ComputeStripeCorrection
  public  ::  stripeCorr, stripeCorrAvailable, stripeCorrDate
  public  ::  stripeCorr_averaging_time

  ! --- Storage of the stripe correction data during the run ---

  ! This (public) array stores the stripe correction = mean(scd - scd_model)
  ! which is an average over the past stripeCorr_averaging_time days
  real, dimension(:), allocatable :: stripeCorr 

  ! This (private) array stores the last stripe correction determined
  real, dimension(:), allocatable :: stripeCorr_last 

  ! Have the stripe correction parameters been read ? 
  !   should be .true. when the stripCorr is applied
  logical  ::  stripeCorrAvailable = .false.

  ! Date for which the stripe correction was determined
  integer, dimension(6)  ::  stripeCorrDate = (/ 1000, 1, 1, 0, 0, 0 /)

  ! This number defines over what period the stripe correction is averaged
  ! If set to 1.0, stripeCorr = stripeCorr_last
  ! Overwritten by value from the rc file
  real    ::  stripeCorr_averaging_time = 7.0    ! in days

  ! Do we start with a new stripe correction (missing file)
  logical  ::  stripeCorr_new = .true.

  ! Debug print statements?
  logical  ::  debug_stripe = .true.

  ! Name of the stripecorr files
  character(len=*), parameter  ::  StripeCorrFilenameBase = "NO2_StripeCorrection_"

  ! Orbit within this window are used for determining the stripe correction
  real, parameter              ::  pacific_lonmin = -163.0
  real, parameter              ::  pacific_lonmax = -137.0

  ! Minimum number of valid pixels along track (stripe correction = undefined otherwise)
  integer, parameter           ::  minNrValidPixels = 100

  ! real, parameter                    ::  undef = -1.0e30
  real, parameter              :: nf_fill_float = 9.9692099683868690e+36

  ! attributes of the files
  character(len=*), parameter  ::  dataset_author = 'Henk Eskes'
  character(len=*), parameter  ::  institution = 'KNMI'

  ! Module name
  character(len=*), parameter  ::  mname = 'MStripeCorrection'


contains


  subroutine ReadStripeCorrection ( filedir, date, nGroundPixels, status )
    !=======================================================================
    !
    ! ReadStripeCorrection: 
    !    Read the list of NO2 slant column corrections to be applied 
    !      to the viewing angles of e.g. OMI 
    !    It reads from the file of date-minus-one-day, at the start of the run
    !      Example file name: "NO2_StripeCorrection_20160131.dat"
    !        for a run starting on 2016-02-01, 00:00
    !  
    !    This routine must be called before stripe corrections are applied, 
    !      because it initialises the "stripeCorr" array
    !    The routine may be called every time step, but will only read the 
    !      stripe correction once during the first call
    !
    ! input:
    !   filedir       = directory with stripecorr files
    !   date          = the actual date-time in TM5
    !   nGroundPixels = number of across track pixels (60 for OMI)
    ! output:
    !   status        = status ( 0 means ok)
    !
    !                                           Henk Eskes, KNMI, March 2016
    !=======================================================================

    use MDF,      only : MDF_Open,      MDF_NETCDF,  MDF_READ 
    use MDF,      only : MDF_Inq_VarID, MDF_Get_Var, MDF_Close 
    use MDF,      only : MDF_Inquire_Dimension, MDF_Inq_DimID

    implicit none

    ! in/out
    character(len=*), intent(in)       :: filedir
    integer, dimension(6), intent(in)  :: date
    integer, intent(in)                :: nGroundPixels
    integer, intent(out)               :: status

    ! local
    ! character(len=*), parameter      ::  rname = mname//'/ReadStripeCorrection'
    character(len=*), parameter      ::  rname = 'ReadStripeCorrection'

    integer, dimension(6)            ::  date_yesterday
    integer                          ::  fid, dimid, id_stripecorr, id_stripecorr_last
    integer                          ::  i, ct_itau, nGP
    character(256)                   ::  filename
    character(256), dimension(4)     ::  comments 
    character(8)                     ::  datestr
    real                             ::  sc
    logical                          ::  file_exists

    real, dimension(:), allocatable  ::  stripeCorr_read

    ! start code

    ! Read only once from file at the start of the run 
    if ( .not. stripeCorrAvailable ) then

       status = 0

       ! allocate and initialise stripe correction to = 0
       allocate ( stripeCorr(nGroundPixels) )
       allocate ( stripeCorr_Last(nGroundPixels) )
       stripeCorr(:) = 0.0
       stripeCorr_last(:) = 0.0
       
       ! The stripe correction has been initialised 
       ! (even if the stripe correction file does not exist)
       stripeCorrAvailable = .true.

       ! Determine the day before the start of the run
       call date2tau ( date, ct_itau )
       ct_itau = ct_itau - 24*3600
       call tau2date ( ct_itau, date_yesterday )

       ! Construct file name, yesterday
       write ( datestr, '(i4,2i2.2)' ) date_yesterday(1:3)
       filename = trim(filedir)//'/'//StripeCorrFilenameBase//datestr//".nc"
       print '(a)', trim(rname) // ', try to read StripeCorr from "' // trim(filename) // '" '

       ! Check file availability
       inquire ( file=trim(filename), Exist=file_exists )
       if ( .not. file_exists ) then

          print '(a)', "    WARNING netcdf file not found '" // trim(filename) // "'"

          ! Check for file two days before the start of the run
          call date2tau ( date, ct_itau )
          ct_itau = ct_itau - 48*3600
          call tau2date ( ct_itau, date_yesterday )
          ! Construct file name, day before yesterday
          write ( datestr, '(i4,2i2.2)' ) date_yesterday(1:3)
          filename = trim(filedir)//'/'//StripeCorrFilenameBase//datestr//".nc"
          print '(a)', trim(rname) // ', try to read StripeCorr from "' // trim(filename) // '" '
          
          ! Check file availability, report errors 
          inquire ( file=trim(filename), Exist=file_exists )
          if ( .not. file_exists ) then
             print '(a)', "    WARNING netcdf file also not found '" // trim(filename) // "'"
             print '(a)', "    Initialising stripe correction to = 0 "
             status = -1
             return
          end if
       end if
       
       ! Open file
       call MDF_Open( trim(filename), MDF_NETCDF, MDF_READ, fid, status )
       IF_NOTOK_RETURN
       if ( debug_stripe )  print '(a,a)', '   stripeCorr file opened: ', trim(filename)

       ! get the dimension
       call MDF_Inq_DimID( fid, 'dimGroundPixel', dimid, status ) 
       IF_NOTOK_RETURN
       call MDF_Inquire_Dimension( fid, dimid, status, length=nGP ) 
       IF_NOTOK_RETURN
       if ( nGP /= nGroundPixels ) then
          print '(a)', trim(rname) // ", Error reading stripe correction, number of ground pixels incorrect"
          print '(a)', "    Initialising stripe correction to = 0 "
          return
       end if
       
       ! read the date
       call MDF_Get_Att( fid, MDF_GLOBAL, 'Date', stripeCorrDate, status )
       IF_NOTOK_RETURN
       if ( debug_stripe )  print '(a,i5,5i3)', '   stripeCorrDate = ',stripeCorrDate

       ! temp storage
       allocate ( stripeCorr_read(nGroundPixels) )

       ! Read stripe correction coefficients
       CALL MDF_Inq_VarID( fid, 'stripeCorrection', id_stripecorr, status )
       IF_ERROR_RETURN

       CALL MDF_Get_Var( fid, id_stripecorr, stripeCorr_read, status )
       IF_NOTOK_RETURN
       if ( debug_stripe )  print '(a)', '   read array stripeCorrection '

       ! if reading went well, copy
       stripeCorr(:) = stripeCorr_read(:)
       stripeCorr_read(:) = 0.0

       CALL MDF_Inq_VarID( fid, 'stripeCorrectionLast', id_stripecorr_last, status )
       IF_ERROR_RETURN

       CALL MDF_Get_Var( fid, id_stripecorr_last, stripeCorr_read, status )
       IF_NOTOK_RETURN
       if ( debug_stripe )  print '(a)', '   read array stripeCorrectionLast '

       ! if reading went well, copy
       stripeCorr_last(:) = stripeCorr_read(:)

       ! remove storage
       deallocate ( stripeCorr_read )

       ! Close file
       CALL MDF_Close( fid, status )
       IF_NOTOK_RETURN

       ! if stripeCorr file was available, and no error occurred
       stripeCorr_new = .false.

       if ( debug_stripe )  print '(a,60f7.3)', '   stripeCorr = ', stripeCorr(:)
          
    end if

  end subroutine ReadStripeCorrection



  subroutine WriteStripeCorrection ( filedir, no2Tr, lat_min, lat_max, status )
    !=======================================================================
    !
    ! WriteStripeCorrection: 
    !    Write the stripe correction data to file
    !  
    ! This routine is not public, but is called by "ComputeStripeCorrection"
    !
    ! input:
    !   filedir    = directory with stripe correction files
    !   no2Tr      = used are the file name info
    !   lat_min, lat_max  = latitude range for which the stripe correction 
    !                       was determined 
    ! output:
    !   status     = status ( 0 means the file has been created )
    !
    !                                 Henk Eskes, KNMI, March 2016, Aug 2016
    !=======================================================================

    implicit none

    ! in/out
    character(len=*), intent(in)     ::  filedir 
    type(TObsTrack), intent(in)      ::  no2Tr
    real, intent(in)                 ::  lat_min, lat_max
    integer, intent(out)             ::  status

    ! local
    ! character(len=*), parameter      ::  rname = mname//'/WriteStripeCorrection'
    character(len=*), parameter      ::  rname = 'WriteStripeCorrection'

    integer                          ::  ncid, fid
    character(256)                   ::  filename
    character(256)                   ::  commentline
    character(8)                     ::  datestr
    integer                          ::  nGroundPixels
    integer                          ::  itau_x
    integer, dimension(6)            ::  idate_x

    integer                          ::  dimid, varid_stripecorr, varid_stripecorr_last

    ! start code

    ! Determine the day to be used for the filename, e.g.
    !    stripecorr determined for orbit 24 jan, 23:52 => file date = 24 jan 
    !    stripecorr determined for orbit 25 jan, 00:12 => file date = 24 jan 
    call date2tau ( stripeCorrDate, itau_x )
    itau_x = itau_x - 6*3600
    call tau2date ( itau_x, idate_x )

    ! Construct file name
    write ( datestr, '(i4,2i2.2)' ) idate_x(1:3)
    filename = trim(filedir)//'/'//trim(StripeCorrFilenameBase)//datestr//".nc"
      
    ! The stripe correction should have been initialised 
    if ( .not. stripeCorrAvailable ) then
       print '(a)', trim(rname) // ": ERROR stripeCorr not initialised, file not created"
       status = -1
       return
    end if

    print '(a)', trim(rname) // ": now writing stripeCorr to file " // trim(filename)
       
    status = 0

    ! nr of viewing angles = nr of stripe correction values
    nGroundPixels = no2Tr%dimGroundPixel

    ! overwrite existing files (clobber)
    call MDF_Create( trim(filename), MDF_NETCDF4, MDF_REPLACE, ncid, status )
    IF_NOTOK_RETURN
       
    ! o global attributes

    call MDF_Put_Att( ncid, MDF_GLOBAL, 'Title', 'Stripe correction factors, NO2' , status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, MDF_GLOBAL, 'DatasetAuthor'    , trim(dataset_author) , status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, MDF_GLOBAL, 'Institution'      , trim(institution)    , status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, MDF_GLOBAL, 'Date'             , stripeCorrDate(1:6)  , status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, MDF_GLOBAL, 'L2_OrbitFilename' , trim(no2Tr%orbitParts(1)%filename), status)
    IF_NOTOK_MDF(fid=ncid)

    ! Write the two comment lines

    write ( commentline, '(i5,5i3,a)' ) stripeCorrDate(1:6),"  Date and time for which time-averaged slant column stripe correction was determined "
    call MDF_Put_Att( ncid, MDF_GLOBAL, 'Comment1' , trim(commentline), status)
    IF_NOTOK_MDF(fid=ncid)

    write ( commentline, '(a,f8.2,a,f8.2)' ) "Latitude range, from ",lat_min," to ", lat_max
    call MDF_Put_Att( ncid, MDF_GLOBAL, 'Comment2' , trim(commentline), status)
    IF_NOTOK_MDF(fid=ncid)

    ! o define dimensions

    call MDF_Def_Dim( ncid, 'dimGroundPixel', nGroundPixels, dimid, status )
    IF_NOTOK_MDF(fid=ncid)

    ! o define variables

    call MDF_Def_Var( ncid, 'stripeCorrection', MDF_FLOAT, (/dimid/), varid_stripecorr, status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, varid_stripecorr, 'standard_name', 'stripe correction' ,    status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, varid_stripecorr, 'long_name',     'row-dependent correction factor, mean over past orbits' ,    status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, varid_stripecorr, 'units',         '10^15 molecules cm^-2' , status)
    IF_NOTOK_MDF(fid=ncid)

    call MDF_Def_Var( ncid, 'stripeCorrectionLast', MDF_FLOAT, (/dimid/), varid_stripecorr_last, status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, varid_stripecorr_last, 'standard_name', 'stripe correction last' ,    status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, varid_stripecorr_last, 'long_name',     'row-dependent correction factor, derived from last orbit only' ,    status)
    IF_NOTOK_MDF(fid=ncid)
    call MDF_Put_Att( ncid, varid_stripecorr_last, 'units',         '10^15 molecules cm^-2' , status)
    IF_NOTOK_MDF(fid=ncid)

    !  o end definition mode

    call MDF_EndDef( ncid , status)
    IF_NOTOK_MDF(fid=ncid)

    ! o write stripe correction arrays

    call MDF_Put_Var( ncid, varid_stripecorr, stripeCorr,  status)
    IF_NOTOK_MDF(fid=ncid)

    call MDF_Put_Var( ncid, varid_stripecorr_last, stripeCorr_last,  status)
    IF_NOTOK_MDF(fid=ncid)
              
    ! close the file
    call MDF_Close( ncid, status)
    IF_NOTOK_RETURN
       

  end subroutine WriteStripeCorrection



  subroutine ComputeStripeCorrection ( filedir, date, no2Tr, nObs, modelSCD, flag, status )
    !=======================================================================
    !
    ! ComputeStripeCorrection:
    !    Check if the stripe correction should be determined for this orbit
    !      (only over the Pacific)  
    !    If so, compute NO2 slant column stripe correction
    !      and write the stripe correction data to file
    !
    ! The stripe correction is derived over the Pacific, for latitudes 
    !   between "lat_min" and "lat_max"
    !  
    ! input:
    !   filedir    = directory with stripe correction files
    !   date       = the actual date-time in TM5
    !   no2Tr      = contains the L2 slant column data product
    !   nObs       = number of observations in this orbit
    !   modelSCD   = slant column computed from the model
    !   flag       = 0 for a retrieval without errors and (snow-ice) warnings
    !
    ! output:
    !   status     = status ( 0 means the file has been created )
    !
    !                                           Henk Eskes, KNMI, March 2016
    !=======================================================================

    use physconstants,   only : pi

    implicit none

    ! in/out
    character(len=*), intent(in)          ::  filedir 
    integer, dimension(6), intent(in)     ::  date
    type(TObsTrack), intent(in)           ::  no2Tr
    integer, intent(in)                   ::  nObs
    real, intent(in), dimension(nObs)     ::  modelSCD
    integer, intent(in), dimension(nObs)  ::  flag
    integer, intent(out)                  ::  status

    ! local
    ! character(len=*), parameter        ::  rname = mname//'/ComputeStripeCorrection'
    character(len=*), parameter        ::  rname = 'ComputeStripeCorrection'

    real                               ::  dist_min, dist, lat_equator, lon_equator
    real                               ::  fit_a, fit_b 
    integer                            ::  iObs, iObs_equator, i, j
    integer                            ::  maxScanLineIndex
    integer                            ::  itau, itau_stripecorr
    integer                            ::  block_jmin, block_jmax
    real                               ::  block_lat_min, block_lat_max
    real                               ::  coeff, season_shift, year_fraction 
    integer                            ::  nGroundPixel
    integer                            ::  nValidStripeCorrs

    real, dimension(:,:), allocatable  ::  delta_scd_2d
    real, dimension(:), allocatable    ::  lat_2d
    integer, dimension(:), allocatable ::  nlat_2d
    real, dimension(:), allocatable    ::  delta_scd_aver, rownr
    integer, dimension(:), allocatable ::  n_scd_aver

    real, dimension(5) :: xtest = (/ 1.0, 2.0, 3.0, 4.0, 5.0 /), ytest = (/ 3.0, 3.0, 4.0, 4.0, nf_fill_float /)

    ! start code
    
    ! Check if this time step has already been processed
    call date2tau ( date, itau )
    call date2tau ( stripeCorrDate, itau_stripecorr )
    if ( itau == itau_stripecorr ) then
       ! A stripe correction has already been determined for this orbit
       print '(a)', trim(rname) // ": a stripe correction was already computed for this orbit" 
       return
    end if

    ! Across-track number of pixels = number of stripe correction coefficients
    nGroundPixel = no2Tr%dimGroundPixel

    ! Determine if the orbit is over the Pacific ?
    ! Nadir longitude @ equator window : -150 \pm 30
    ! OMI advances 24.72 degree longitude between orbits
    ! Criterion for accepting the orbit: When orbit longitude at the equator falls within [-163,-137] 

    ! Determine equator crossing orbit longitude 
    !   by minimising the distance to the equator and to the middle of the orbit
    !   (0.24 = roughly the width of an OMI pixel in degree)
    dist_min = 1.0e6   ! degree
    do iObs = 1, nObs
       dist = 0.24 * Abs( real(no2Tr%pixelIndex(iObs)) - (0.5*nGroundPixel+0.5) ) + Abs( no2Tr%latitude(iObs) )
       if ( dist < dist_min ) then
          iObs_equator = iObs
          lat_equator = no2Tr%latitude(iObs)
          lon_equator = no2Tr%longitude(iObs)
          dist_min = dist
       end if
    end do
    if ( debug_stripe )  print '(a,i6,2f8.2,i5)', trim(rname)//': Equator crossing coordinates: iObs, lat, lon, pixIndex = ', &
         iObs_equator, lat_equator, lon_equator, no2Tr%pixelIndex(iObs_equator)

    ! Check if this is a pacific orbit. If not, leave.
    if ( ( lon_equator < pacific_lonmin ) .or. ( lon_equator > pacific_lonmax ) ) then
       print '(a)', trim(rname) // ": this orbit is not a Pacific ocean orbit, return "
       return
    end if

    ! Now we have found an orbit over the Pacific: determine stripe correction coefficients

    print '(a)', trim(rname) // ": this orbit is a Pacific ocean orbit, now determining stripe amplitude "
    
    ! --- Re-grid the data from 1D to 2D: lon, lat, scd-scd_model ---

    ! First compute max scan line index
    maxScanLineIndex = -1
    do iObs = 1, nObs
       if ( no2Tr%scanLineIndex(iObs) > maxScanLineIndex ) then
          maxScanLineIndex = no2Tr%scanLineIndex(iObs)
       end if
    end do
    if ( debug_stripe )  print *, '   Max scan line index = ', maxScanLineIndex
    
    ! Then allocate and regrid
    allocate ( delta_scd_2d ( nGroundPixel, maxScanLineIndex ) )
    allocate ( lat_2d ( maxScanLineIndex ) )
    allocate ( nlat_2d ( maxScanLineIndex ) )
    allocate ( delta_scd_aver( nGroundPixel ) )
    allocate ( rownr( nGroundPixel ) )
    allocate ( n_scd_aver( nGroundPixel ) )
    
    ! Determine ( SCD_sat - SCD_model ) on the 2D swath grid
    delta_scd_2d(:,:) = nf_fill_float
    lat_2d(:) = 0.0
    nlat_2d(:) = 0
    do iObs = 1, nObs
       if ( iand(flag(iObs), 255) == 0 ) then
          i = no2Tr%pixelIndex(iObs)
          j = no2Tr%scanLineIndex(iObs)
          delta_scd_2d(i,j) = no2Tr%no2SLC(iObs) - modelSCD(iObs)
          lat_2d(j) = lat_2d(j) + no2Tr%latitude(iObs)
          nlat_2d(j) = nlat_2d(j) + 1
       end if
    end do
    do j = 1, maxScanLineIndex
       if ( nlat_2d(j) > 0 ) then
          lat_2d(j) = lat_2d(j) / nlat_2d(j)
       else
          lat_2d(j) = nf_fill_float
       end if
    end do

    ! print '(a,2000e10.2)', '   lat_2d = ', lat_2d(:)
    if ( debug_stripe )  print '(a,60i6)',   '    pixInd   50000+60 = ', no2Tr%pixelIndex(50000:50059)
    if ( debug_stripe )  print '(a,60f6.2)', '    no2SLC   50000+60 = ', no2Tr%no2SLC(50000:50059)
    if ( debug_stripe )  print '(a,60f6.2)', '    modelSCD 50000+60 = ', modelSCD(50000:50059)

    ! Determine block of data in latitude range
    !   winter: [30S,0],  summer: [0,30N]
    year_fraction = (real(date(2))+real(date(3))/31.0-1.0)/12.0
    season_shift = sin ( pi*year_fraction )
    block_lat_min = -30.0 + 30.0 * season_shift
    block_lat_max = 30.0 * season_shift

    if ( debug_stripe )  print *,'   date = ',date
    if ( debug_stripe )  print *,'   year fraction, season shift, latitude range = ', year_fraction, season_shift, block_lat_min, block_lat_max

    ! determine block indices
    block_jmin = 1
    block_jmax = 1
    do j = 1, maxScanLineIndex
       if ( lat_2d(j) < 0.9*nf_fill_float ) then
          if ( block_lat_min > lat_2d(j) ) block_jmin = j
          if ( block_lat_max > lat_2d(j) ) block_jmax = j
       end if
    end do
    ! print *, '   lat_2d = ', lat_2d(:)
    if ( debug_stripe )  print *,           '   block range = ',block_jmin,block_jmax
    if ( debug_stripe )  print '(a,2f8.2)', '   lat range   = ',lat_2d(block_jmin),lat_2d(block_jmax)

    ! determine mean stripe amplitude for each row over the latitude block
    nValidStripeCorrs = 0
    do i = 1, nGroundPixel
       delta_scd_aver(i) = 0.0
       n_scd_aver(i) = 0
       rownr(i) = real(i)
       do j = block_jmin, block_jmax
          if ( delta_scd_2d(i,j) < 0.9*nf_fill_float ) then
             delta_scd_aver(i) = delta_scd_aver(i) + delta_scd_2d(i,j)
             n_scd_aver(i) = n_scd_aver(i) + 1
          end if
       end do
       if ( n_scd_aver(i) >= minNrValidPixels ) then
          delta_scd_aver(i) = delta_scd_aver(i) / real(n_scd_aver(i))
          nValidStripeCorrs = nValidStripeCorrs + 1
       else
          delta_scd_aver(i) = nf_fill_float
       end if
    end do
    if ( debug_stripe )  print '(a,60f8.3)', '   delta_scd_aver(:) = ', delta_scd_aver(:)
    if ( debug_stripe )  print '(a,60i5)',   '       n_scd_aver(:) = ', n_scd_aver(:)
    if ( debug_stripe )  print '(a,i4)',     '   nr of valid stripe corrs = ', nValidStripeCorrs

    ! release storage
    deallocate ( nlat_2d )
    deallocate ( lat_2d )
    deallocate ( delta_scd_2d )

    if ( debug_stripe )  then
       ! determine linear fit
       call  linearFit ( xtest, ytest, 5, nf_fill_float, fit_a, fit_b ) 
       print '(a,2f10.4)', '   test fit_a, fit_b = ',fit_a, fit_b
    end if

    ! determine linear fit
    call  linearFit ( rownr, delta_scd_aver, nGroundPixel, nf_fill_float, fit_a, fit_b ) 
    if ( debug_stripe )  print '(a,2f10.4)', '   fit_a, fit_b = ',fit_a, fit_b
    
    ! New "stripeCorr_last":
    !   stripe correction = mean-slant-column-difference minus linear-fit 
    do i = 1, nGroundPixel
       if ( delta_scd_aver(i) < 0.9*nf_fill_float ) then
          stripeCorr_last(i) = delta_scd_aver(i) - ( fit_a * rownr(i) + fit_b )
       else
          stripeCorr_last(i) = nf_fill_float
       end if
    end do
    if ( debug_stripe )  print '(a,60f8.3)', '   stripeCorr_last(1:60) = ',stripeCorr_last(1:60)

    ! Define new stripe correction as a mix of the previous one and the new one found for this orbit
    ! with a mean averaging time
    if ( stripeCorr_averaging_time > 1.99 ) then
       if ( stripeCorr_new ) then
          ! No information from the day before, so simply copy the new stripe correction found
          stripeCorr(:) = stripeCorr_last(:)
       else
          ! mix of new stripecorr with stripecorr obtained on previous days 
          coeff = 1.0 / stripeCorr_averaging_time
          do i = 1, nGroundPixel
             if ( stripeCorr_last(i) < 0.9*nf_fill_float ) then
                if ( stripeCorr(i) < 0.9*nf_fill_float ) then
                   ! well defined new stripe corr value, and existing old one
                   stripeCorr(i) = (1.0-coeff)*stripeCorr(i) + coeff*stripeCorr_last(i)
                else
                   ! replace stripecorr when not defined
                   stripeCorr(i) = stripeCorr_last(i)
                end if
             else
                ! do not change existing stripeCorr(i)
             end if
          end do
       end if
    else
       ! The averaging time is shorter than 2 days: use instantaneous value
       stripeCorr(:) = stripeCorr_last(:)
    end if
    if ( debug_stripe )  print '(a,60f8.3)', '   stripeCorr(1:60) = ',stripeCorr(1:60)
       
    ! Date for which the stripe correction was determined
    stripeCorrDate(:) = date(:)

    ! And, finally, write the two new stripe coefficients to file
    call WriteStripeCorrection ( filedir, no2Tr, block_lat_min, block_lat_max, status )

    ! stripeCorr file is available
    stripeCorr_new = .false.

    deallocate ( delta_scd_aver )
    deallocate ( rownr )
    deallocate ( n_scd_aver )

    if ( status /= 0 ) then
       print '(a)', trim(rname) // ": ERROR writing the new stripe correction "
    end if

  end subroutine ComputeStripeCorrection



  subroutine linearFit ( x, y, n0, nf_fill_float, a, b ) 

    ! Compute simple linear fit
    !    y = a x + b
    ! For a set of n points (xi,yi) where yi <> nf_fill_float
    !
    !    nf_fill_float : large negative value indicating value is not defined
    !
    ! Solution: a = ( ave(xy) -ave(x) ave(y) ) / ( ave(x^2) - ave(x)^2 )
    !           b = ave(y) - a ave(x)
    !

    implicit none

    ! in/out
    real, dimension(n0), intent(in)  :: x, y
    integer, intent(in)              :: n0
    real, intent(in)                 :: nf_fill_float
    real, intent(out)                :: a, b

    ! local
    integer           ::  n, i 
    real              ::  sx, sy, sxx, sxy

    ! begin code
    n = 0
    sx = 0.0
    sy = 0.0
    sxx = 0.0
    sxy = 0.0
    do i = 1, n0
       if ( y(i) < 0.9 * nf_fill_float ) then
          n = n + 1
          sx = sx + x(i)
          sy = sy + y(i)
          sxx = sxx + x(i)*x(i)
          sxy = sxy + x(i)*y(i)
       end if
    end do
    if ( n > 1 ) then  ! need at least two valid points
       sx = sx/real(n)
       sy = sy/real(n)
       sxx = sxx/real(n)
       sxy = sxy/real(n)
       a = (sxy-sx*sy)/(sxx-sx*sx)
       b = sy - a*sx
    else
       a = nf_fill_float
       b = nf_fill_float
    end if

  end subroutine linearFit
    


end module MStripeCorrection