! 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