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

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

! #ifdef TROPNLL2DP

! #define WRITE_WARNING(notok_string)   call fortranlog(notok_string,len(notok_string),2)

! #else

! Macro
#define WRITE_WARNING(notok_string)   write (*,'(a)') notok_string
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if

! #endif


module MObservationOperator
  !
  ! Computes the model predicted NO2 observation (observation operator)
  !
  ! Use:
  !    After reading a track of NO2 (and cloud data),
  !    call "GetOMFHx" to obtain a structure "ObsFcInfo" with the model 
  !    forecast, observation operator, vertical columns and related parameters
  !
  ! Contains:
  !    subroutine GetOMFHx( amfLut, no2Tr, TM5Data, no2pcf, ObsFcInfo )
  !
  !     Henk Eskes, Folkert Boersma, Bram Maasakkers KNMI, 2000-2003-2005-2012-2015

  ! printing, error handling
  use GO,               only : gol, goPr, goErr

  use pqf_module

  implicit none

  private

  public  ::  InitOMFHx, GetOMFHx

  ! --------------------------------------------------------------
  ! Temperature of the cross section used in the DOAS fit
  ! --------------------------------------------------------------

  ! Temperature of NO2 cross section used in DOAS fit = 220 K
  !    Vandaele et al., 1998
  real, parameter  :: Tdoasfit = 220.0

  ! Coefficient in the fit of the T dependence of the cross section amplitude 
  !         real, parameter  :: Tx = 11.39      ! K   This is the old T correction result, DOMINO-2
  ! Update by Marina Zara, TN-OMIE-KNMI-982, March 2016, for DOMINO-3
  !   correction factor = 1 - 0.00316 (T-Tdoasfit) + 3.39e-6 (T-Tdoasfit)^2
  real, parameter  :: Tcorrection_Acoeff = -0.00316
  real, parameter  :: Tcorrection_Bcoeff = 3.39e-6

  ! Provide error message below this value of the tropospheric AMF
  real, parameter   ::  amfTropCutOff = 0.1

  ! Undefined value
  ! real, parameter   :: undef = -999.99
  real, parameter   :: nf_fill_float = 9.9692099683868690e+36

  ! The following parameters are taken from the .rc file

  ! Is cloud radiance fraction taken from the retrieval, or calculated in GetOMFHx ?
  logical  ::  useCloudRadianceFractionFromFile

  ! Is the terrain pressure taken from the retrieval file, or calculated in GetOMFHx ?
  logical  ::  useTerrainPressureFromFile

  ! Apply stripe corrections (in the case of OMI/TROPOMI)
  logical  ::  doApplyStripeCorrection 

  ! Flag to identify first call to the routine GetOmFHx
  logical  ::  firstcall = .true.

  ! Is the albedo taken from the OMI file 
  ! logical, parameter :: useOMIAlbedo = .true.

  character(len=*), parameter  ::  mname = 'MObservationOperator'

contains


  subroutine InitOMFHx( rcFile )

    ! -----------------------------------------------------------------------
    ! Initialise the observation operator module: 
    !      read settings from .rc file
    !                                              Henk Eskes, KNMI, Aug 2016
    ! -----------------------------------------------------------------------

    use GO,                 only : TrcFile, Init, Done, ReadRc
    use MStripeCorrection,  only : stripeCorr_averaging_time    

    implicit none

    ! in/out
    character(*), intent(in)         ::  rcFile
    ! local
    type(TrcFile)    ::  rcF
    integer          ::  status
    real             ::  aveTime
    character(len=*), parameter  ::  rname = trim(mname)//'/InitOMFHx'

    ! open rcfile:
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    
    ! Is cloud radiance fraction taken from the retrieval, or calculated in GetOMFHx ?
    call ReadRc( rcF, 'domino.useCloudRadianceFractionFromFile', useCloudRadianceFractionFromFile, status, default=.true. )
    IF_NOTOK_RETURN(status=1)
    
    ! Is the terrain pressure taken from the retrieval file, or calculated in GetOMFHx ?
    call ReadRc( rcF, 'domino.useTerrainPressureFromFile', useTerrainPressureFromFile, status, default=.false. )
    IF_NOTOK_RETURN(status=1)

    ! Apply stripe corrections (in the case of OMI/TROPOMI)
    call ReadRc( rcF, 'domino.doApplyStripeCorrection', doApplyStripeCorrection, status, default=.false. )
    IF_NOTOK_RETURN(status=1)
    
    if ( doApplyStripeCorrection ) then
       ! This number defines over what period the stripe correction is averaged 
       ! If set to 1.0, stripeCorr = stripeCorr_last 
       ! Unit: days
       call ReadRc( rcF, 'domino.StripeCorrectionAveragingTime', aveTime, status, default=7.0 )
       IF_NOTOK_RETURN(status=1)
       stripeCorr_averaging_time = aveTime
    else
       print '(a)', 'InitOMFHx: No stripe correction applied'
    end if

    ! close rcfile:
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

  end subroutine InitOMFHx



  subroutine GetOMFHx( amfLut, no2Tr, TM5Data, no2pcf, ObsFcInfo, isForecast, outputdir, debug )
    !=======================================================================
    !
    ! GetOMFHx
    !     Determine Observation, model forecast and relevant parameters
    !     Data is stored in structure "ObsFcInfo"
    !
    ! Input:  
    !    amfLut     : structure with the AMF lookup-table 
    !    no2Tr      : NO2 SLC + cloud data of type(TObsTrack)
    !    TM5Data    : data structure contains all TM fields and params needed
    !    no2pcf     : 3D field of NO2 partial columns for each TM grid cell (10^15 molecules per cm^2)
    !    isForecast : true for forecast, false for analysis
    !    outputdir  : directory for stripe correction output
    !    debug      : true = provide additional debug output
    !
    ! Output: 
    !   "ObsFcInfo" - of type "TObsFcInfo"
    !            Contains output of the vertical column retrieval:
    !            grid cell indices, forecast, air mass factors 
    !            vertical columns, error estimates and 
    !            the corresponding observation operator vectors
    !
    !             Henk Eskes, Folkert Boersma, Ruud Dirksen, Bram Maasakkers
    !                                                      KNMI, 2002 - 2015
    !=======================================================================

    ! --- TM5 model arrays and variables ---
    use MTM5Data,               only : TTM5Data 

    ! --- Structures ---
    ! structure containing observation input (slant columns)
    use MTObsTrack,          	only : TObsTrack
    ! structure containing observation, forecast and the observation operator
    use MTObsFcInfo,            only : TObsFcInfo

    ! --- Assimilation interfaces ---
    use physconstants,  	only : pi, Earth_rad, grav
    use Mrweight,           	only : rweight
    use MTAmf,                  only : TAmfLut
    use MAmfGet,                only : GetAmf, GetAmf_vectorised, GetAmfGeo, FitWindowCentre
    use MGridTools,         	only : getTMCellIndex, getTMCellIndex4
    use MNO2RetrievalError, 	only : NO2ErrEstimate, modelStratVcdError

    ! --- Error message codes ---
    use pqf_module

    ! --- OMI stripe correction codes ---
    use MStripeCorrection,      only : ReadStripeCorrection, ComputeStripeCorrection
    use MStripeCorrection,      only : stripeCorr, stripeCorrAvailable 

    implicit none

    ! --- in/out ---
    type(TAmfLut), intent(inout)                                ::  amfLut
    type(TObsTrack), intent(inout)                              ::  no2Tr			
    type(TTM5Data), intent(in)                                  ::  TM5Data 
    real,dimension(TM5Data%im,TM5Data%jm,TM5Data%lm),intent(in) ::  no2pcf
    type(TObsFcInfo), intent(inout)                             ::  ObsFcInfo		
    logical, intent(in)                                         ::  isForecast
    character(len=*), intent(in)                                ::  outputdir 
    logical, intent(in)                                         ::  debug

    ! local: parameters
    integer, parameter                     :: nrtape = 55982
    real, parameter                        :: gamma = 6.5     ! K/km lapse rate
    character(*), parameter                :: rname = 'GetOMFHx'

    ! local: allocatable arrays
    real, dimension(:), allocatable       :: no2collev, amflev, amfclearlev, amfcloudlev
    integer, dimension(:), allocatable    :: ixc, iyc, ltropopause
    real, dimension(:), allocatable       :: cell_oro
    integer, dimension(:,:), allocatable  :: ix4a, iy4a
    real, dimension(:,:), allocatable     :: w4a
    real, dimension(:,:), allocatable     :: phybrid, phybrTop, phybrBot
    real, dimension(:,:), allocatable     :: amfclearrel, amfcloudrel
    real, dimension(:,:), allocatable     :: Tlayer, sigmacorrlev
    real, dimension(:), allocatable       :: sza, vza, aza
    real, dimension(:), allocatable       :: cloudfrac, cloudpres, surfpres
    real, dimension(:), allocatable       :: albcloud, albclear, cldradfraction
    real, dimension(:), allocatable       :: rclear, rcloud, amfgeo
    logical, dimension(:), allocatable    :: isSnowIce
! only used in debug
    real, dimension(:), allocatable       :: phybrid_mean, sigmacorrlev_mean

    ! local
    real                                  :: gdxc, gdyc
    real                                  :: no2gvc, no2slcfc, no2vcdfc
    real                                  :: amftotal
    real                                  :: satvcd, satscd, rsat
    real                                  :: no2slcfctrop, no2vcdfctrop, amftroptotal, amfclear
    real                                  :: no2slcfctropclear
    real                                  :: no2slcfcstrat, no2vcdfcstrat, satvcdtrop, satslctrop
    real                                  :: no2vcdfcnoghost, amftropnoghost
    real                                  :: abovecloudfraction
    real                                  :: errVcdTotal, errVcdTropTotal, errVcdStratTotal
    real                                  :: errVcdTotalAK, errVcdTropTotalAK
    real                                  :: spres_in, dT, pow, peff
    real                                  :: cos_sza, cos_vza, dphi
    real                                  :: dTemp
    integer                               :: status
    integer                               :: cntSceneScheme, nstripe
    integer                               :: NObs, iObs, iObs_print, i, i2, l
    integer                               :: nPixelsWithErrors, nSmallAmftroptotal, nPeff1050

    ! test difference between CRF computed and in the file
    real                                  :: crf_m0, crf_m1, crf_m2, crf_dif

    ! test of amf
    real,dimension(100)                   :: amf_test

    ! for error messages
    character(256)                        :: errmessage

    ! begin code

    print '(a)', 'GetOMFHx: start computing vertical columns, AMFs, kernels'

    if ( doApplyStripeCorrection .and. firstcall ) then
       ! read the OMI slant column stripe corrections from file
       ! only once when routine is called for the first time 
       call ReadStripeCorrection ( outputdir, TM5Data%date, no2Tr%dimGroundPixel, status )
       if ( status /= 0 ) then
          print '(a)', "GetOmFHx: WARNING, failed to read stripe correction file "
       end if
    end if

    print '(a,60f7.3)', "GetOmFHx: stripeCorr = ", stripeCorr(:)
    if ( doApplyStripeCorrection .and. stripeCorrAvailable ) then
       print '(a)', "GetOmFHx: stripeCorr will be applied "
    else
       print '(a)', "GetOmFHx: stripeCorr will NOT be applied "
    end if

    ObsFcInfo%count = no2Tr%count 
    NObs = no2Tr%count   ! number of observations in this track
   
    ! allocate helper arrays
    allocate(cloudfrac(Nobs))                   ; allocate(cloudpres(Nobs))          
    allocate(sza(Nobs))                         ; allocate(vza(Nobs))                
    allocate(aza(Nobs))                         ; allocate(surfpres(Nobs))           
    allocate(albcloud(Nobs))                    ; allocate(albclear(Nobs))
    allocate(rcloud(Nobs))                      ; allocate(rclear(Nobs))
    allocate(cldradfraction(Nobs))              ; allocate(isSnowIce(NObs))
    allocate(ix4a(Nobs,4))                      ; allocate(iy4a(Nobs,4))
    allocate(w4a(Nobs,4))                       ; allocate(ixc(Nobs))
    allocate(iyc(Nobs))                         ; allocate(cell_oro(Nobs))
    allocate(ltropopause(Nobs))                 ; allocate(Tlayer(Nobs,TM5Data%lm))
    allocate(phybrid(Nobs,TM5Data%lm))  	; allocate(phybrTop(Nobs,TM5Data%lm))
    allocate(phybrBot(Nobs,TM5Data%lm)) 	; allocate(amfclearrel(Nobs,TM5Data%lm))
    allocate(amfcloudrel(Nobs,TM5Data%lm))      ; allocate(amfgeo(Nobs))
    allocate(sigmacorrlev(Nobs,TM5Data%lm)) 
    allocate ( amfclearlev(TM5Data%lm) )  	; allocate(amfcloudlev(TM5Data%lm))
    allocate ( amflev(TM5Data%lm) )  
    allocate ( no2collev(TM5Data%lm) ) 

    if ( debug ) then
       ! extra mean debug ouput of sigmacorrlev
       allocate(phybrid_mean(TM5Data%lm))  	; allocate(sigmacorrlev_mean(TM5Data%lm))
    end if

    ! Copies of no2Tr
    cloudfrac(:) = no2Tr%cloudFraction(1:NObs)
    cloudpres(:) = no2Tr%cloudTopPressure(1:NObs)
    sza(:)       = no2Tr%solarZenithAngle(1:NObs)
    vza(:)       = no2Tr%viewingZenithAngle(1:NObs)
    albcloud(:)  = no2Tr%cloudAlbedo(1:NObs)

    ! Cloud model switch in case of Ice and Snow
    ! For ice or snow: use scene pressure and scene albedo instead
    ! snowIceFlag -  0: snow free; 1-100: % sea ice cover; 101: ice; 103: snow; 255: ocean
    cntSceneScheme = 0
    do iObs = 1, NObs
       ! Swith to "scene" retrieval in case of snow or ice cover > 50 %
       isSnowIce(iObs) = ( (no2Tr%snowIceFlag(iObs) > 50) .and. (no2Tr%snowIceFlag(iObs) < 200) )
       if ( isSnowIce(iObs) ) then
          cntSceneScheme = cntSceneScheme + 1
          cloudfrac(iObs) = 1.0
          cloudpres(iObs) = no2Tr%scenePressure(iObs)
          albcloud(iObs) = min(1.0,no2Tr%sceneAlbedo(iObs))
       end if
    end do
    if ( debug )  then
       print '(a,i6,a,i6)', ' Switch to scene albedo/pressure for ',cntSceneScheme,' pixels out of ',NObs
    end if

    ! Use surface albedo map provided in L2 file
    !if ( useOMIAlbedo) albclear(iObs) = no2Tr%surfaceAlbedo(iObs)
    albclear(:)  = no2Tr%surfaceAlbedo(1:NObs)

    ! Determine relative azimuth angles from the Sun and viewing azimuth
    do iObs = 1, NObs
       dphi = no2Tr%viewingAzimuthAngle(iObs) - no2Tr%solarAzimuthAngle(iObs) - 180.0;
       if ( dphi <    0.0 ) dphi = dphi + 360.0;
       if ( dphi >= 360.0 ) dphi = dphi - 360.0;
       if ( dphi >= 180.0 ) dphi = 360.0 - dphi;
       aza(iObs) = dphi                  ! in degree, range [0,180)
    end do

    iObs_print = max ( nObs/2, 1 )   ! middle of the orbit
    if ( debug )  then
       iObs = iObs_print
       print*,'NObs, iobs = ',NObs,iObs
       print*,'  cloudfrac = ',cloudfrac(iObs)
       print*,'  cloudpres = ',cloudpres(iObs)
       print*,'  cloudalb  = ',albcloud(iObs)
       print*,'  cloud CRF = ',no2Tr%cloudRadianceFraction(iObs)
       print*,'  albsurf = ',no2Tr%surfaceAlbedo(iObs)
       print*,'  terrheight = ',no2Tr%terrainHeight(iObs)
       print*,'  terrpres = ',no2Tr%terrainPressure(iObs)
       print*,'  sza = ',sza(iObs)
       print*,'  vza = ',vza(iObs)
       print*,'  aza = ',aza(iObs)
       print*,'  lat,lon = ',no2Tr%latitude(iObs),no2Tr%longitude(iObs)
       print*,'  slc = ',no2Tr%no2slc(iObs)
       print*,'  slc err = ',no2Tr%no2slcerror(iObs)
       print*,'  pixelFlag = ',no2Tr%pixelFlag(iObs)
       print*,'  snowIceFlag = ',no2Tr%snowIceFlag(iObs)
       print*,'  pixelInd = ',no2Tr%pixelIndex(iObs)
       print*,'  scanline = ',no2Tr%scanLineIndex(iObs)
       print*,'  min-max albcloud = ',minval(albcloud),maxval(albcloud)
       print*,'  min-max cloudpres = ',minval(cloudpres),maxval(cloudpres)
       print*,'  min-max sceneAlbedo = ',minval(no2Tr%sceneAlbedo),maxval(no2Tr%sceneAlbedo)
       print*,'  min-max scenePressure = ',minval(no2Tr%scenePressure),maxval(no2Tr%scenePressure)
       print*,'file = ',no2Tr%orbitParts(1)%filename
    end if

    ! initialise the flag to the flag obtained from the SLC file
    ObsFcInfo%flag(1:nObs) = no2Tr%pixelFlag(1:nObs)

    ! input range checks
    call CheckValueRanges ( nObs, cloudFrac, no2Tr%cloudRadianceFraction, albcloud, albclear, &
         sza, vza, aza, no2Tr%no2slc, no2Tr%longitude, no2Tr%latitude, ObsFcInfo%flag )

    ! clip values
    where ( albcloud > 0.99 ) albcloud = 0.99
    where ( cloudfrac > 1.0 ) cloudfrac = 1.0 
    where ( cloudfrac < 0.0 ) cloudfrac = 0.0
    where ( albclear > 1.0 ) albclear = 1.0 
    where ( albclear < 0.0 ) albclear = 0.0

    if ( debug )  then
       print *, 'TM5data, im, jm, lm ',TM5Data%im,TM5Data%jm,TM5Data%lm
       print *, 'TM5Data '
       print *, '   min max      t = ',minval(TM5Data%t),maxval(TM5Data%t)
       print *, '   min-max     ps = ',minval(TM5Data%ps),maxval(TM5Data%ps)
       print *, '   min-max    oro = ',minval(TM5Data%oro),maxval(TM5Data%oro)
       print *, '   min-max    no2 = ',minval(TM5Data%no2),maxval(TM5Data%no2)
       print *, '   min-max ltropo = ',minval(TM5Data%ltropo),maxval(TM5Data%ltropo)
    end if

    ! Compute model quantities @ observation locations
    do iObs = 1, NObs

       ! skip to next iObs if (flag = error),  accept warnings
       if ( iand(ObsFcInfo%flag(iObs), 255) /=  0 ) cycle

       ! Find the four TM5 grid points closest to the pixel centre
       call GetTMCellIndex4(TM5Data%im,TM5Data%jm,no2Tr%longitude(iObs),no2Tr%latitude(iObs),ix4a(iObs,:),iy4a(iObs,:),w4a(iObs,:))
          
       ! Find the TM5 grid cell containing the pixel
       call GetTMCellIndex(TM5Data%im,TM5Data%jm,no2Tr%longitude(iObs),no2Tr%latitude(iObs),ixc(iObs),iyc(iObs),gdxc,gdyc)
          
       ! Store the grid cell indices
       ObsFcInfo%icell(iObs)           = ixc(iObs)
       ObsFcInfo%jcell(iObs)           = iyc(iObs)

       ! Determine surface pressure (in hPa)
       surfpres(iObs) = w4a(iObs,1) * TM5Data%ps(ix4a(iObs,1),iy4a(iObs,1)) + &
            w4a(iObs,2) * TM5Data%ps(ix4a(iObs,2),iy4a(iObs,2)) + &
            w4a(iObs,3) * TM5Data%ps(ix4a(iObs,3),iy4a(iObs,3)) + &
            w4a(iObs,4) * TM5Data%ps(ix4a(iObs,4),iy4a(iObs,4))
       
       surfpres(iObs) = 0.01 * surfpres(iObs) ! Pa -> hPa
          
       ! Calculate the terrainheight according to the model (in meter "m")
       cell_oro(iObs) = w4a(iObs,1) * TM5Data%oro(ix4a(iObs,1),iy4a(iObs,1)) + &
            w4a(iObs,2) * TM5Data%oro(ix4a(iObs,2),iy4a(iObs,2)) + &
            w4a(iObs,3) * TM5Data%oro(ix4a(iObs,3),iy4a(iObs,3)) + &
            w4a(iObs,4) * TM5Data%oro(ix4a(iObs,4),iy4a(iObs,4))
       
       ! Tropoause level
       ltropopause(iObs) = TM5Data%ltropo(ixc(iObs),iyc(iObs))
          
       ! Temperature of the layers
       !   Sgphx now works similar to the call in emission_nox.F90 !JDM
       !   KFB - 01-10-2009; correct the surface pressure by using the real surface elevation of the OMI pixel
       !   Need temperature here
       do l = 1, TM5Data%lm
          ! determine layer temperature    
          Tlayer(iObs,l) = w4a(iObs,1) * TM5Data%t(ix4a(iObs,1),iy4a(iObs,1),l) + &
               w4a(iObs,2) * TM5Data%t(ix4a(iObs,2),iy4a(iObs,2),l) + &
               w4a(iObs,3) * TM5Data%t(ix4a(iObs,3),iy4a(iObs,3),l) + &
               w4a(iObs,4) * TM5Data%t(ix4a(iObs,4),iy4a(iObs,4),l) 
       end do

       if ( debug .and. (iObs == iObs_print) ) then
          print*,'Observations debug output'
          print*,'   iObs = ',iObs
          print*,'   ix4a = ',ix4a(iObs,:)
          print*,'   iy4a = ',iy4a(iObs,:)
          print*,'   ixc = ',minval(ixc),maxval(ixc)
          print*,'   min max iyc = ',minval(iyc),maxval(iyc)
          print*,'   min max w4a = ',minval(w4a),maxval(w4a)
          print*,'   min max cell_oro = ',minval(cell_oro),maxval(cell_oro)
          print*,'   min max albclear = ',minval(albclear),maxval(albclear)
          print*,'   min max surfpres = ',minval(surfpres),maxval(surfpres)
          print*,'   min max cloudpres = ',minval(cloudpres),maxval(cloudpres)
          print*,'   size modelterrainheight = ',size(ObsFcInfo%modelterrainheight)
          print*,'   size cell_oro = ',size(cell_oro)
       end if

    end do   ! loop over observations, iObs

    if ( debug ) then
       ! test profile of mean T correction
       phybrid_mean(:)      = 0.0
       sigmacorrlev_mean(:) = 0.0
    end if

    ! Count number of times pressure rescaling exceeds threshold
    nPeff1050 = 0

    if ( .not. useCloudRadianceFractionFromFile ) then
       ! Initialise CRF counters
       print *, 'Initialise CRF counters'
       crf_m0 = 0.0
       crf_m1 = 0.0
       crf_m2 = 0.0
    end if

    ! surface pressure, cloud pressure, pressure levels, radiance levels, amf_geo
    do iObs = 1, NObs

       ! skip to next iObs if error flag is set
       if ( iand(ObsFcInfo%flag(iObs), 255) /=  0 ) cycle

       if ( useTerrainPressureFromFile ) then
          ! In case the provided terrain pressure has sufficient resolution
          surfpres(iObs) = no2Tr%terrainPressure(iObs)
       else
          ! Correct the model surface pressure to account for differences in terrain height
          !    KFB - 01-10-2009; following equation (5) in Zhou et al., AMT, 2009.
          !    JDMTH, using the average terrain height instead of the value at the center of the pixel
          dT  = Tlayer(iObs,1)/( Tlayer(iObs,1) + 0.001*gamma*(cell_oro(iObs)-no2Tr%terrainHeight(iObs)) )
          pow = (-1.0*grav)/(287.*0.001*gamma) 
          peff = surfpres(iObs) * (dT**pow)
          if ( peff > 1050.0 ) then
             nPeff1050 = nPeff1050 + 1
             peff = 1050.0
          end if
          surfpres(iObs) = peff
       end if
       
       ! clip the cloud pressures, domain [surfpres, 130.]
       cloudpres(iObs) = min(surfpres(iObs),cloudpres(iObs))
       cloudpres(iObs) = max(130.01,cloudpres(iObs))
       
       ! range checks for cloud and surface pressure
       if ( (cloudpres(iObs) < 0.0) .or. (cloudpres(iObs) > 1100.0) )  ObsFcInfo%flag(iObs) = PQF_E_CLOUD_ERROR 
       if ( (surfpres(iObs) < 0.0) .or. (surfpres(iObs) > 1100.0) )    ObsFcInfo%flag(iObs) = PQF_E_GENERIC_RANGE_ERROR 
       
       ! determine mid-level TM5 pressures (in hPa), with the corrected model surface pressure
       do l = 1, TM5Data%lm          
          ! top/bottom and midlevel pressure of the layers (in hPa)
          phybrTop(iObs,l) = 0.01*( TM5Data%hyai(l+1) + 100.0*surfpres(iObs)*TM5Data%hybi(l+1) )
          phybrBot(iObs,l) = 0.01*( TM5Data%hyai(l)   + 100.0*surfpres(iObs)*TM5Data%hybi(l)   )
          phybrid(iObs,l)  = 0.5*( phybrTop(iObs,l) + phybrBot(iObs,l) ) 
       end do
       if ( phybrid(iObs,ltropopause(iObs)) > 500. ) then   
          write ( errmessage, *) 'Warning - Tropopause is too low: ', no2Tr%longitude(iObs), no2Tr%latitude(iObs), phybrid(iObs,ltropopause(iObs)) 
          WRITE_WARNING( trim(errmessage) ) 
       end if
       
       ! compute NO2 cross section temperature correction term 
       !    sigmacorr = sigma(T)/sigma(Tdoasfit)
       !    the term Tx is the result of a least-squares fit of the fitted
       !    slant column vs. cross section temperature
       do l = 1, TM5Data%lm          
          ! old formula used in DOMINO-2: (Tx = 11.39)
          !      sigmacorrlev(iObs,l) = ( Tdoasfit - Tx ) / ( Tlayer(iObs,l) - Tx )
          dTemp = Tlayer(iObs,l) - Tdoasfit
          sigmacorrlev(iObs,l) = 1.0 + Tcorrection_Acoeff * dTemp + Tcorrection_Bcoeff * dTemp * dTemp
       end do

       if ( debug ) then
          ! determine orbit-mean of the T-correction factor
          do l = 1, TM5Data%lm 
             phybrid_mean(l) = phybrid_mean(l)  + phybrid(iObs,l)
             sigmacorrlev_mean(l) = sigmacorrlev_mean(l) + sigmacorrlev(iObs,l)
          end do
       end if
       
       ! determine theoretical clear-sky and cloud-covered radiance levels, and cloud radiance fraction 
       
       ! estimate radiation intensities for cloud and no cloud conditions, used for the error estimate
       ! rcloud, rclear are needed anyway for the error estimate, 
       !    even if cldradfraction  is taken from input
       call rweight( FitWindowCentre,cloudpres(iObs),albcloud(iObs),sza(iObs),vza(iObs),aza(iObs),rcloud(iObs) )
       call rweight( FitWindowCentre,surfpres(iObs),albclear(iObs),sza(iObs),vza(iObs),aza(iObs),rclear(iObs) )
       
       ! flagging
       if ( (rclear(iObs) <= 0.0) .or. (rcloud(iObs) <= 0.0) )   ObsFcInfo%flag(iObs) = PQF_E_REFLECTANCE_RANGE_ERROR
       
       ! cloud radiance fraction
       if ( useCloudRadianceFractionFromFile ) then
          ! Use value from the retrieval
          if ( isSnowIce(iObs) ) then
             cldradfraction(iObs) = 1.0  ! Fully covered snow-ice "scene" 
          else
             cldradfraction(iObs) = no2Tr%cloudRadianceFraction(iObs)
          end if
          ! In this case:
          !      rclear = f(1-crf)/[(1-f)crf] rcloud             
          ! clipping
          if ( cldradfraction(iObs) > 1.0 ) cldradfraction(iObs) = 1.0 
          if ( cldradfraction(iObs) < 0.0 ) cldradfraction(iObs) = 0.0 
       else
          ! Estimated satellite-observed radiance (normalization factor)
          rsat = (1.0-cloudfrac(iObs))*rclear(iObs) + cloudfrac(iObs)*rcloud(iObs)
          ! Compute radiance fraction originating from the cloudy part of the pixel
          cldradfraction(iObs) = cloudfrac(iObs)*rcloud(iObs) / rsat
          ! Determine CRF difference with file (for snow/ice free pixels)
          if ( .not. isSnowIce(iObs) ) then
             !crf_dif = cldradfraction(iObs)
             crf_dif = cldradfraction(iObs) - no2Tr%cloudRadianceFraction(iObs)
             crf_m0 = crf_m0 + 1.0
             crf_m1 = crf_m1 + crf_dif
             crf_m2 = crf_m2 + crf_dif*crf_dif
          end if
       end if
       
       ! fill "amfgeo", the geometrical air mass factor
       !    NOTE: the AMF lookup table should contain AMF/AMFgeo, where AMFgeo has exactly the same form as here 
       cos_sza = cos(sza(iObs)*pi/180.0)
       cos_vza = cos(vza(iObs)*pi/180.0)
       call GetAmfGeo ( cos_vza, cos_sza, amfGeo(iObs) ) 

       if ( debug .and. (iObs == iObs_print) ) then
          print*,'Geometry debug output, iObs = ',iObs
          print*,'   no2Tr%terrainPressure(iObs) = ',no2Tr%terrainPressure(iObs)
          print*,'   surfpres(iObs)              = ',surfpres(iObs)
          print*,'   cloudpres(iObs)             = ',cloudpres(iObs)
          print'(a,34F8.2)','    phybrid                     = ',phybrid(iObs,:)
          print'(a,34F6.3)','    sigmacorrlev                = ',sigmacorrlev(iObs,:)
          print*,'   rclear,rcloud               = ',rclear(iObs),rcloud(iObs)
          print*,'   cldradfraction              = ',cldradfraction(iObs)
          print*,'   amfGeo                      = ',amfGeo(iObs)
       end if

    end do  ! iObs

    ! Determine CRF difference with file
    if ( .not. useCloudRadianceFractionFromFile ) then
       ! print *,'crf_m0,1,2 = ',crf_m0, crf_m1, crf_m2
       print '(a,F10.8,a,F10.8)', ' mean CRF difference, bias(computed-file)  = ', crf_m1/crf_m0, &
            ', var = ', sqrt( crf_m2/crf_m0 - (crf_m1/crf_m0)**2 )
    end if

    if ( nPeff1050 > 0 ) then
       write ( errmessage, '(a,i7,a)') 'WARNING GetOMFHx: peff set to 1050 for ',nPeff1050,' observations'
       WRITE_WARNING( trim(errmessage) ) 
    end if

    if ( debug ) then
       ! Mean T-correction output 
       phybrid_mean(:)      = phybrid_mean(:) / real(NObs)
       sigmacorrlev_mean(:) = sigmacorrlev_mean(:) / real(NObs)
       print*,' '
       print '(a,34F8.2)', ' mean phybrid   = ',phybrid_mean(:)
       print '(a,34F6.3)', ' mean sigmacorr = ',sigmacorrlev_mean(:)
    end if

    ! read sensitivities from the AMF lookup table: clear case    
    call GetAmf_vectorised( nObs,TM5Data%lm,phybrid,aza,vza,sza,albclear,surfpres,amfLut,amfclearrel,obsFcInfo%flag,debug )

    if ( debug ) print*,'GetAmf_vectorised done (clear)'

    if ( debug ) then
       iObs = iObs_print
       do l = 1, TM5Data%lm
          call GetAmf(phybrid(iObs,l),aza(iObs),vza(iObs),sza(iObs),albclear(iObs),surfpres(iObs),amfLut,amf_test(l),status)
       end do
       print '(a,34f6.3)', ' amf (debug) = ',amf_test(1:TM5Data%lm)
    end if

    ! read sensitivities from the AMF lookup table: 100% cloud case    
    call GetAmf_vectorised( nObs,TM5Data%lm,phybrid,aza,vza,sza,albcloud,cloudpres,amfLut,amfcloudrel,obsFcInfo%flag,.false. )

    if ( debug ) print*,'GetAmf_vectorised done (cloudy)'

    ! count the number of flagged pixels and pixels with small AMF
    nPixelsWithErrors = 0
    nSmallAmftroptotal = 0

    ! number of pixels where SCD is corrected
    nstripe = 0

    ! main retrieval loop
    mainObsLoop: do iObs = 1, NObs 

       ! skip observation when error
       if ( iand(ObsFcInfo%flag(iObs), 255) /=  0 ) then
          nPixelsWithErrors = nPixelsWithErrors + 1
          cycle
       end if       

       no2gvc       = 0.0 ; no2slcfc     = 0.0
       no2vcdfc     = 0.0 ; no2collev    = 0.0
       no2slcfctrop = 0.0 ; no2vcdfctrop = 0.0
       no2slcfctropclear = 0.0 

       ! the height-dependent AMF for clear sky 
       do l = 1, TM5Data%lm 
          amfclearlev(l) = sigmacorrlev(iObs,l)*amfclearrel(iObs,l)*amfgeo(iObs)
       end do

       ! the height-dependent AMF for 100% cloud cover
       amfcloudlev(:) = 0.0
       do l = 1, TM5Data%lm
          if ( cloudpres(iObs) >= (phybrTop(iObs,l)-1.0e-5) ) then
             if ( cloudpres(iObs) >= (phybrBot(iObs,l)-1.0e-5) ) then
                ! model level above the cloud top       
                amfcloudlev(l) = sigmacorrlev(iObs,l)*amfcloudrel(iObs,l)*amfgeo(iObs)
             else
                ! model level contains cloud top
                abovecloudfraction = (cloudpres(iObs) - phybrTop(iObs,l))/(phybrBot(iObs,l) - phybrTop(iObs,l))
                
                !Nieuwe code die clipt & melding schrijft: RuudDirksen 19 juni 2009
                spres_in = phybrBot(iObs,l)+1.0e-5
                if ( spres_in .lt. 119.82 ) then
                   write ( errmessage, '(a,3(2x,f8.4))' ) &
                        'GetOMFHx: Clipped cloud pressure; p lvl c_f c_p:',spres_in,cloudfrac(iObs),cloudpres(iObs)
                   WRITE_WARNING( trim(errmessage) ) 
                   spres_in = 119.82
                end if
                !call GetAmf(phybrid(iObs,l),aza(iObs),vza(iObs),sza(iObs),albcloud(iObs),phybrid(iObs,l)+1.0e-5,amfcloudrel(iObs,l))
                call GetAmf(phybrid(iObs,l),aza(iObs),vza(iObs),sza(iObs),albcloud(iObs),spres_in,amfLut,amfcloudrel(iObs,l),status)
                ! set flag if status = error
                if ( status > 0 )  obsFcInfo%flag(iObs) = PQF_E_LUT_RANGE_ERROR
                ! set amf for the level containing cloud top
                amfcloudlev(l) = sigmacorrlev(iObs,l)*amfcloudrel(iObs,l)*amfgeo(iObs)*abovecloudfraction
             end if
          end if
       end do  ! l, loop over vertical levels

       ! skip observation when flagged in the loop above
       if ( iand(ObsFcInfo%flag(iObs), 255) /=  0 ) then
          nPixelsWithErrors = nPixelsWithErrors + 1
          cycle
       end if

       ! loop over vertical TM5 hybrid levels
       do l = 1, TM5Data%lm

          ! determine amount of NO2 in this layer (kg/m^2)
          
          ! compute NO2 subcolumns for the observation and for layer "l" using the "no2pcf" model field
          !     unit 10^15 molec/cm^2
          no2collev(l) = w4a(iObs,1) * no2pcf(ix4a(iObs,1),iy4a(iObs,1),l) + &
                         w4a(iObs,2) * no2pcf(ix4a(iObs,2),iy4a(iObs,2),l) + &
                         w4a(iObs,3) * no2pcf(ix4a(iObs,3),iy4a(iObs,3),l) + &
                         w4a(iObs,4) * no2pcf(ix4a(iObs,4),iy4a(iObs,4),l)

          ! determine AMF weighted with cloud fraction and radiation intensity 
          ! amflev(iObs,l) = (1.0-cloudfrac(iObs))*rclear(iObs)*amfclearlev(iObs,l)/rsat(iObs) + &
          !                       cloudfrac(iObs) *rcloud(iObs)*amfcloudlev(iObs,l)/rsat(iObs)
          amflev(l) = (1.0-cldradfraction(iObs))*amfclearlev(l) + &
                           cldradfraction(iObs) *amfcloudlev(l)

          ! debug
          if ( debug .and. (iObs == iObs_print) .and. (l == TM5Data%lm) ) then
             print '(a,i6)', 'TM5 profile debug output, iObs = ',iObs
             print '(a,2f10.4)', '    crf, amfgeo  = ', cldradfraction(iObs), amfgeo(iObs) 
             print '(a,34f8.4)', '    no2collev(:) = ',no2collev(:)
             print '(a,34f8.4)', '    amfclearlev(:) = ',amfclearlev(:)
             print '(a,34f8.4)', '    amfcloudlev(:) = ',amfcloudlev(:)
             print '(a,34f8.4)', '    amflev(:)    = ',amflev(:)
          end if
	  
          ! compute model-predicted NO2 slant & vertical column
          ! no2collev = NO2 mass in grid cell (units 10^15 molecules cm^-2)
          no2slcfc = no2slcfc + amflev(l)*no2collev(l)
          no2vcdfc = no2vcdfc + no2collev(l)

          ! compute tropospheric model slant & vertical column
          if (  l <= ltropopause(iObs) ) then
             no2slcfctrop = no2slcfctrop + amflev(l)*no2collev(l)
             no2slcfctropclear = no2slcfctropclear + amfclearlev(l)*no2collev(l)
             no2vcdfctrop = no2vcdfctrop + no2collev(l)
          end if

          ! compute ghost vertical column "GVC", the vertical column below the cloud top.
          ! Burrows et al., J. Atmos. Sci., 56, 151-175, 1999.
          if ( cloudpres(iObs) >= phybrTop(iObs,l)-1.0e-5 ) then
             if ( cloudpres(iObs) >= phybrBot(iObs,l)-1.0e-5 ) then
                ! model level above the cloud top, no need to update ghost column
             else
                ! model level contains cloud top    
                abovecloudfraction  = (cloudpres(iObs) - phybrTop(iObs,l))/(phybrBot(iObs,l) - phybrTop(iObs,l))            
                no2gvc = no2gvc + (1. - abovecloudfraction )*no2collev(l)
             end if
          else
             ! model level below cloud top
             no2gvc = no2gvc + no2collev(l)
          end if
          ! ghost column only defined for cloud fractions > 0
          if ( cloudfrac(iObs) < 0.001 ) no2gvc = 0.

       end do  ! loop over TM5 hybrid levels

       ! checks to avoid division by 0
       if ( no2vcdfc <= 1.0e-8 ) then
          write ( errmessage, '(a,i6,2x,f8.2,2x,f8.2)' ) &
               'GetOMFHx: no2vcdfc <=0 for ',iObs,no2Tr%latitude(iObs),no2Tr%longitude(iObs)
          WRITE_WARNING( trim(errmessage) ) 
          ObsFcInfo%flag(iObs) = PQF_E_GENERIC_RANGE_ERROR          
          nPixelsWithErrors = nPixelsWithErrors + 1 
          cycle
       end if       
       if ( no2vcdfctrop <= 1.0e-10 ) then
          write ( errmessage, '(a,i6,2x,f8.2,2x,f8.2)' ) &
              'GetOMFHx: no2vcdfctrop <=0 for ',iObs,no2Tr%latitude(iObs),no2Tr%longitude(iObs)
          WRITE_WARNING( trim(errmessage) ) 
          ObsFcInfo%flag(iObs) = PQF_E_GENERIC_RANGE_ERROR        
          nPixelsWithErrors = nPixelsWithErrors + 1 
          cycle
       end if
       
       ! total tropospheric air-mass factor
       amftroptotal = no2slcfctrop / no2vcdfctrop

       ! too small tropospheric AMFs not allowed, provide error
       if ( amftroptotal < amfTropCutOff ) then
          nSmallAmftroptotal = nSmallAmftroptotal + 1
          nPixelsWithErrors = nPixelsWithErrors + 1
          ObsFcInfo%flag(iObs) = PQF_E_GENERIC_RANGE_ERROR         
          cycle
       end if

       ! total air-mass factor
       amftotal = no2slcfc / no2vcdfc               ! model SLC / model VCD
       
       ! total tropospheric air-mass factor, assuming the pixel is cloud free
       amfclear = no2slcfctropclear / no2vcdfctrop

       ! Observed slant column: correct for stripes
       satscd = no2Tr%no2SLC(iObs)
       if ( doApplyStripeCorrection .and. stripeCorrAvailable ) then
          ! apply the OMI slant column stripe corrections to the DOAS slant column 
          if ( stripeCorr( no2Tr%pixelIndex(iObs) ) < 0.9*nf_fill_float ) then
             ! number of pixels where SCD is corrected
             nstripe = nstripe + 1
             satscd = satscd - stripeCorr( no2Tr%pixelIndex(iObs) )
          end if
       end if

       ! retrieved total vertical NO2 column       
       satvcd = satscd / amftotal

       ! estimated NO2 vertical tropospheric column based on the tropospheric air-mass factor
       ! Total model vertical column
       no2slcfcstrat = no2slcfc - no2slcfctrop
       no2vcdfcstrat = no2vcdfc - no2vcdfctrop
       
       ! Ghost column quantities:
       ! For stratospheric clouds, take trop. column as ghost column
       if ( no2gvc > no2vcdfctrop ) no2gvc = no2vcdfctrop	  	  
       ! Calculate the model predicted NO2 VCD without the ghostcolumn contribution              
       no2vcdfcnoghost = (1.0 -  cldradfraction(iObs) )  *  no2vcdfctrop + &
                                 cldradfraction(iObs)    * (no2vcdfctrop-no2gvc)
       ! calculate the tropospheric AMF without incorporating the ghostcolumn
       amftropnoghost = no2slcfctrop / no2vcdfcnoghost

       ! --- Compute error estimates ---

       satslctrop = satscd - no2slcfcstrat
       satvcdtrop = satslctrop / amftroptotal

       call NO2ErrEstimate( TM5Data%lm, &
            cloudfrac(iObs), cloudpres(iObs), rclear(iObs), rcloud(iObs), &
            amflev(:), amfclearlev(:), amfcloudlev(:), &
            albclear(iObs), albcloud(iObs), amfgeo(iObs), &
            ix4a(iObs,:), iy4a(iObs,:), w4a(iObs,:), &
            TM5Data%hyai, TM5Data%hybi, surfpres(iObs), sza(iObs), vza(iObs), aza(iObs), &
            no2collev(:), ltropopause(iObs), &
            amftotal, amftroptotal, &
            satvcd, satvcdtrop, &
            amfLut, &
            errVcdTotal, errVcdTropTotal, errVcdStratTotal, &
            errVcdTotalAK, errVcdTropTotalAK, status )

       if ( status > 0 ) then
          write ( errmessage, '(a)' ) 'ERROR GetOMFHx: error occured in NO2ErrEstimate'
          WRITE_WARNING( trim(errmessage) ) 
          ObsFcInfo%flag(iObs) = PQF_E_GENERIC_EXCEPTION          
          nPixelsWithErrors = nPixelsWithErrors + 1 
          cycle
       end if

       ! Retrieval was successful at this point, but warnings may still occur

       ! When the cloudradiance fraction exceeds 50%, raise troposphericcolumnflag
       ! the tropospheric column value is considered unreliable because more than half of the 
       ! light/spectrum comes from above the cloud.
       ! HenkEskes: in DOMINO-3 there is no longer a warning for CRF > 0.5
       ! where (ObsFcInfo%cloudradfraction > 50.0 )     ObsFCInfo%fltrop = -1

       ! When the pixel is covered by snow the cloud parameters are unreliable. 
       ! snow covered pixels have surface albedo of 0.6
       ! if ( albclear(iObs) > 0.59 )   ObsFCInfo%flag(iObs) = ior( ObsFCInfo%flag(iObs), PQF_W_SNOW_ICE_WARNING )
       if ( isSnowIce(iObs) )  ObsFCInfo%flag(iObs) = ior( ObsFCInfo%flag(iObs), PQF_W_SNOW_ICE_WARNING )
       
       ! show results
       if ( debug .and. ( iObs == iObs_print ) ) then
          print '(a,i6,a,8i4,4f6.3,2f8.1)', 'iObs =  ', iObs, ' i, j, w, lon, lat ', ix4a(iObs,:), iy4a(iObs,:), w4a(iObs,:), no2Tr%longitude(iObs),no2Tr%latitude(iObs)
          print '(a,3f8.2)', ' vcdfc, satvcd, amf ',no2vcdfc, satvcd, amftotal
          print '(a,3f8.2)', ' vcdtrop, sig, sig_AK  ', satvcdtrop, errVcdTropTotal, errVcdTropTotalAK
          print '(a,2f8.2)', ' vcdstrat, sig ', no2vcdfcstrat, errVcdStratTotal
          print '(a,3f8.2)', ' amftotal, amftrop, amfclear ', amftotal, amftroptotal, amfclear
          print '(a,2f8.2)', ' amfstrat, amfgeo ', no2slcfcstrat / no2vcdfcstrat, amfgeo(iObs)
          print '(a,50f6.2)', ' AK ',  amflev(:)/amftotal
          print '(a,f8.2)', ' modelterrainhgt ',  cell_oro(iObs)
          print '(a,i4)', ' l_tropo ',  ltropopause(iObs)
       end if

       ! store retrieval data in structure "ObsFcInfo"

       ! vcd, amf
       ObsFcInfo%no2vcd(iObs)          = satvcd 
       ObsFcInfo%no2vcdsum(iObs)       = satvcdtrop + no2vcdfcstrat
       obsFcInfo%no2vcdsumsig(iObs)    = sqrt ( errVcdStratTotal**2 + errVcdTropTotal**2)
       ObsFcInfo%no2vcdsig(iObs)       = errVcdTotal
       ObsFcInfo%no2vcdsigak(iObs)     = errVcdTotalAK
       ObsFcInfo%no2vcdtrop(iObs)      = satvcdtrop
       ObsFcInfo%no2vcdtropsig(iObs)   = errVcdTropTotal
       ObsFcInfo%no2vcdtropsigak(iObs) = errVcdTropTotalAK
       ObsFcInfo%no2vcdstrat(iObs)     = no2vcdfcstrat
       ObsFcInfo%no2vcdstratsig(iObs)  = errVcdStratTotal
       ObsFcInfo%amf(iObs)             = amftotal
       ObsFcInfo%amftrop(iObs)         = amftroptotal
       ObsFcInfo%amfgeo(iObs)          = amfgeo(iObs)
       ObsFcInfo%amfstrat(iObs)        = no2slcfcstrat / no2vcdfcstrat
       ObsFcInfo%amfclear(iObs)        = amfclear
       ObsFcInfo%avkernel(:,iObs)      = amflev(:)/amftotal
       ObsFcInfo%ghostcol(iObs)        = no2gvc
       ! model derived quantities
       ObsFcInfo%psurf(iObs)           = surfpres(iObs)  ! in hPa
       ObsFcInfo%no2vcdfc(iObs)        = no2vcdfc
       ObsFcInfo%no2vcdfctrop(iObs)    = no2vcdfctrop
       ObsFcInfo%no2slcfc(iObs)        = no2slcfc
       ObsFcInfo%no2slcfctrop(iObs)    = no2slcfctrop
       ObsFcInfo%no2slcstrat(iObs)     = no2slcfcstrat
       ObsFcInfo%levtropopause(iObs)   = ltropopause(iObs)
       ! model orography in "m"
       ObsFcInfo%modelTerrainHeight(iObs) = cell_oro(iObs)
       ! other
       ObsFcInfo%cloudradfraction(iObs)= cldradfraction(iObs)    ! range [0,1]

    end do mainObsLoop    ! iObs: loop over observations

    ! Has the cloud radiance fraction been computed in this module?
    ObsFcInfo%cloudRadFraction_computed = (.not. useCloudRadianceFractionFromFile)

    ! number of pixels where SCD is corrected
    print '(a,i7)',  ' GetOMFHx: nr of pixels for which SCD is stripe corrected ', nstripe 
    print '(a,2i6)', ' GetOMFHx: min-max pixelIndex = ', minval(no2Tr%pixelIndex(:)) , maxval(no2Tr%pixelIndex(:)) 

    ! report on the number of errors found (pixels removed)
    if ( nPixelsWithErrors > 0 ) then
       write ( errmessage, '(a,i6)' ) &
            'GetOMFHx: WARNING, nr of pixels with error flag =', nPixelsWithErrors
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nSmallAmftroptotal > 0 ) then
       write ( errmessage, '(a,i6)' ) &
            'GetOMFHx: WARNING, nr of pixels with too small amftroptotal =', nSmallAmftroptotal
       WRITE_WARNING( trim(errmessage) ) 
    end if

    ! deallocate temporary arrays

    if (allocated(no2collev))      deallocate(no2collev) 
    if (allocated(cloudfrac))      deallocate(cloudfrac)
    if (allocated(cloudpres))      deallocate(cloudpres)
    if (allocated(sza))            deallocate(sza)
    if (allocated(vza))            deallocate(vza)
    if (allocated(aza))            deallocate(aza)
    if (allocated(surfpres))       deallocate(surfpres)
    if (allocated(albcloud))       deallocate(albcloud)
    if (allocated(albclear))       deallocate(albclear)
    if (allocated(rclear))         deallocate(rclear)
    if (allocated(rcloud))         deallocate(rcloud)
    if (allocated(cldradfraction)) deallocate(cldradfraction)
    if (allocated(isSnowIce))      deallocate(isSnowIce)
    if (allocated(amfgeo))         deallocate(amfgeo)
    if (allocated(ix4a))           deallocate(ix4a)
    if (allocated(iy4a))           deallocate(iy4a)
    if (allocated(w4a))            deallocate(w4a)
    if (allocated(ixc))            deallocate(ixc)
    if (allocated(iyc))            deallocate(iyc)
    if (allocated(cell_oro))       deallocate(cell_oro)
    if (allocated(ltropopause))    deallocate(ltropopause)
    if (allocated(Tlayer))         deallocate(Tlayer)
    if (allocated(phybrid))        deallocate(phybrid)
    if (allocated(phybrTop))       deallocate(phybrTop)
    if (allocated(phybrBot))       deallocate(phybrBot)
    if (allocated(amfclearrel))    deallocate(amfclearrel)
    if (allocated(amfcloudrel))    deallocate(amfcloudrel)
    if (allocated(sigmacorrlev))   deallocate(sigmacorrlev)
    if (allocated(amfclearlev))    deallocate(amfclearlev)
    if (allocated(amfcloudlev))    deallocate(amfcloudlev)
    if (allocated(amflev))         deallocate(amflev)

    if ( debug ) then
       ! only for extra mean ouput of sigmacorrlev
       if ( allocated(phybrid_mean) )       deallocate(phybrid_mean)
       if ( allocated(sigmacorrlev_mean) )  deallocate(sigmacorrlev_mean)
    end if

    ! new stripe correction 
    if ( doApplyStripeCorrection .and. isForecast ) then
       ! check if this is a Pacific orbit,
       ! compute a new OMI slant column stripe correction, 
       ! and save to file when a new correction was computed
       ! this is skipped in analysis mode
       call  ComputeStripeCorrection ( outputdir, TM5Data%date, no2Tr, nObs, ObsFcInfo%no2slcfc, ObsFcInfo%flag, status )
    end if

    if ( status > 0 ) then
       write ( errmessage, '(a)' ) 'ERROR GetOMFHx: error occured in ComputeStripeCorrection'
       WRITE_WARNING( trim(errmessage) ) 
    end if

    if ( debug ) print*, 'This is the end of GetOMFHx'

    firstcall = .false.

  end subroutine GetOMFHx




  subroutine CheckValueRanges ( nObs, cloudFraction, cloudRadianceFraction, albcloud, albclear, &
                                sza, vza, aza, no2slc, longitude, latitude, flag )
    !
    ! perform range checks for a list of input variables
    !

    implicit none

    integer, intent(in)                 ::  nObs
    real, intent(in), dimension(nObs)   ::  cloudFraction, cloudRadianceFraction, albcloud, albclear
    real, intent(in), dimension(nObs)   ::  sza, vza, aza, no2slc, longitude, latitude
    integer, intent(inout), dimension(nObs)  ::  flag

    integer           :: iObs
    character(256)    :: errmessage
    integer           :: nCloudFrac, nCloudRF, nAlbCloud, nAlbClear, nSZA, nVZA, nAZA, nSLC, nLon, nLat

    ! start code

    nCloudFrac = 0
    nCloudRF = 0
    nAlbCloud = 0
    nAlbClear = 0
    nSZA = 0
    nVZA = 0
    nAZA = 0
    nSLC = 0
    nLon = 0
    nLat = 0

    pixelloop: do iObs = 1, NObs

       ! only check when flag <> error 
       if ( iand(flag(iObs), 255) == 0 ) then
          
          if ( (cloudFraction(iObs) < -0.1) .or. (cloudFraction(iObs) >  1.2) ) then
             flag(iObs) = PQF_E_CLOUD_ERROR
             nCloudFrac = nCloudFrac + 1
          end if
          if ( useCloudRadianceFractionFromFile ) then
             if ( (cloudRadianceFraction(iObs) < -0.1) .or. (cloudRadianceFraction(iObs) > 1.2) ) then
                flag(iObs) = PQF_E_CLOUD_ERROR
                nCloudRF = nCloudRF + 1
             end if
          end if
          if ( (albcloud(iObs) < 0.01) .or. (albcloud(iObs) > 1.1) ) then
             flag(iObs) = PQF_E_GENERIC_RANGE_ERROR
             nAlbCloud = nAlbCloud + 1
          end if
          if ( (albclear(iObs) < -0.1) .or. (albclear(iObs) > 1.1) ) then
             flag(iObs) = PQF_E_GENERIC_RANGE_ERROR
             nAlbClear = nAlbClear + 1
          end if
          if ( (sza(iObs) > 89.0) .or. (sza(iObs) < 0.) ) then
             flag(iObs) = PQF_E_SZA_RANGE_ERROR
             nSZA = nSZA + 1
          end if
          if ( (vza(iObs) > 89.0) .or. (vza(iObs) < 0.) ) then
             flag(iObs) = PQF_E_VZA_RANGE_ERROR
             nVZA = nVZA + 1
          end if
          if ( abs(aza(iObs)) > 360.0 ) then
             flag(iObs) = PQF_E_GENERIC_RANGE_ERROR
             nAZA = nAZA + 1
          end if
          if ( (no2slc(iObs) > 1000.0) .or. (no2slc(iObs) < -1.0) ) then
             flag(iObs) = PQF_E_GENERIC_RANGE_ERROR
             nSLC = nSLC + 1
          end if
          if ( (longitude(iObs) < -360.001) .and. (longitude(iObs) > 360.001) ) then
             flag(iObs) = PQF_E_GEOLOCATION_ERROR
             nLon = nLon + 1
          end if
          if ( (latitude(iObs) < -90.0) .and. (latitude(iObs) > 90.0) ) then
             flag(iObs) = PQF_E_GEOLOCATION_ERROR
             nLat = nLat + 1
          end if

       end if

    end do pixelloop

    if ( nCloudFrac > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - cloud fraction out of range: number of pixels =', nCloudFrac
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nCloudRF > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - cloud radiance fraction out of range: number of pixels =', nCloudRF
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nAlbCloud > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - cloud albedo out of range: number of pixels =', nAlbCloud
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nAlbClear > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - clear albedo out of range: number of pixels =', nAlbClear
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nSZA > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - SZA out of range: number of pixels =', nSZA
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nVZA > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - VZA out of range: number of pixels =', nVZA
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nAZA > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - Rel azimuth angle out of range: number of pixels =', nAZA
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nSLC > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - slant column out of range: number of pixels =', nSLC
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nLon > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - longitude out of range: number of pixels =', nLon
       WRITE_WARNING( trim(errmessage) ) 
    end if
    if ( nLat > 0 ) then
       write ( errmessage, '(a,i6)') 'CheckValueRanges: WARNING - latitude out of range: number of pixels =', nLat
       WRITE_WARNING( trim(errmessage) ) 
    end if

  end subroutine CheckValueRanges

end module MObservationOperator