! 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