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

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

! #ifdef TROPNLL2DP
! 
! #define IF_NOTOK_RETURN(notok_string) call fortranlog(notok_string,len(notok_string),2); status=1; return
! #define PRINT(string) call fortranlog(string,len(string),2)
! 
! #else

#define IF_NOTOK_RETURN(notok_string) write (*,'(a)') notok_string; status=1; return
#define PRINT(string)  write (*,'(a)') trim(string)

! #endif


module MNO2RetrievalError
  !-------------------------------------------------------------------
  ! Routine "NO2ErrEstimate" for calculating the error in the 
  ! retrieval of NO2 due to various assumptions on 
  ! clouds, albedo, and NO2 profile
  !
  ! the module contains routines for the various contributions 
  ! to the error:
  !    ErrFcl     - contribution of Fresco cloud fraction errors
  !    ErrPcl     - contribution of Fresco cloud top pressure errors
  !    ErrAlb     - contribution of albedo map uncertainties 
  !    ErrProfile - contribution related to NO2 profile variations
  !                 [JDM, deleted, we use 10% of the AMF]
  !    ErrMix     - contribution related to NO2 mixing in TM3 BL 
  !    
  ! On top of this the model accounts for:
  !    - the DOAS SCD retrieval error 
  !    - the error in estimating the stratospheric NO2 "background"
  ! 
  !                        Folkert Boersma, Henk Eskes, KNMI 2002-2003
  !
  !-------------------------------------------------------------------

  use MTAmf,             only : TAmfLut
  use Mrweight,          only : rweight
  use MAmfGet,           only : GetAmf, FitWindowCentre

  implicit none  
  
  private

  public :: NO2ErrEstimate, modelStratVcdError
  public :: forecastError  

  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! Error parameters:
  
  ! vertical column forecast error (stratospheric case)
  real, parameter :: forecastError = 0.2 ! 1e15 molecules/cm2 (new choice: observation error reduced)
  ! real, parameter :: forecastError = 0.15 ! 1e15 molecules/cm2 (old choice)

  ! Error in the Fresco cloud fractions
  !JDMERR real,parameter :: frescoError = 0.05
  real, parameter :: frescoError = 0.025     

  ! Error in the Fresco cloud top pressure estimate
  real, parameter :: cloudTopError = 50.0       ! unit hPa

  ! Error in surface albedo
  !JDMERR real,parameter :: albedoError = 0.02
  real, parameter :: albedoError = 0.015

  ! Slant column retrieval (DOAS) error
  !JDMERR real,parameter :: retrievalScdError  = 0.4   	! unit 1e15 molec/cm2
  real, parameter :: retrievalScdError  = 0.55  		! unit 1e15 molec/cm2

  ! Error in estimate of the stratospheric column (from assimilation)
  !JDMERR real,parameter :: modelStratVcdError = 0.25   ! unit 1e15 molec/cm2
  real, parameter :: modelStratVcdError = 0.2   		! unit 1e15 molec/cm2

  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~

  ! TM3 covariance lookup table
  !he_cov character(160),parameter :: &
  !he_cov        CovPath= '/net/bsgi18/nobackup/users/eskes/tm3no2a/tmcovar'

  logical :: firstcall = .true.

  contains

  
  subroutine NO2ErrEstimate( lm, &
       cldfraction,cldtoppres,rclear,rcloud, &
       amflev,amfclearlev,amfcloudlev, &
       albclear,albcloud,amfgeo, &
       ix4a,iy4a,w4a, &
       at,bt,surfpres,sza,vza,aza, &
       no2collev,ltropopause, &
       amftotal,amftroptotal, &
       satvcd,satvcdtrop, &
       amfLut, &
       errVcdTotal,errVcdTropTotal,errVcdStratTotal, &
       errVcdTotalAK,errVcdTropTotalAK, &
       status )

    ! ------------------------------------------------------------------
    ! Compute an error estimate for the tropospheric vertical column and 
    ! total vertical column
    ! input:
    !    for explanation of the variables, see "ass_Hx.f90"
    ! output:
    !    errVcdTotal       : error on the total vertical NO2 column
    !    errVcdTropTotal   : error on the tropospheric NO2 column
    !    errVcdStratTotal  : error on the stratospheric NO2 column
    !    errVcdTotalAK     : error on the total vertical NO2 column, when the
    !                        averaging kernels are used 
    !                        (sum of observation + kernel error)
    !    errVcdTropTotalAK : same, for tropospheric column
    !    status            : > 0 means error
    ! ------------------------------------------------------------------

    ! use MAmfLookupTable,   only : FitWindowCentre, GetAmf

    implicit none

    ! in/out
    integer,intent(in)              :: lm
    integer,intent(in)              :: ltropopause
    integer,dimension(4),intent(in) :: ix4a, iy4a
    real,dimension(4),intent(in)    :: w4a
    real,intent(in)                 :: cldfraction,cldtoppres
    real,intent(in)                 :: rclear,rcloud,amfgeo
    real,intent(in)                 :: albclear,albcloud,surfpres
    real,intent(in)                 :: sza,vza,aza
    real,intent(in)                 :: amftotal,amftroptotal
    real,intent(in)                 :: satvcd,satvcdtrop
    real,dimension(lm),intent(in)   :: amflev,amfclearlev,amfcloudlev
    real,dimension(lm),intent(in)   :: no2collev
    real,dimension(lm+1),intent(in) :: at,bt 
    type(TAmfLut), intent(inout)    :: amfLut
    !
    real,intent(out)                :: errVcdTotal
    real,intent(out)                :: errVcdTropTotal, errVcdStratTotal
    real,intent(out)                :: errVcdTotalAK, errVcdTropTotalAK
    integer, intent(out)            :: status

    ! local
    character(*), parameter    :: rname = 'NO2ErrEstimate'
    !he_cov  real, dimension(im,jm,lm,lm) :: cov_matrix 
    !he_cov  real, dimension(im,jm,lm)    :: no2mean
    real :: errtotal_fcl,errtotal_pcl,errtotal_alb
    real :: errtroptotal_prof,errtotal_prof
    real :: errtroptotal_fcl,errtroptotal_pcl,errtroptotal_alb
    real :: err_fcl,err_pcl,err_alb,no2vcdtotal,no2vcdtroptotal
    real :: fcl_alb_correlation
    real :: erramftroptotal,erramftroptotalAK
    real :: erramftotal,erramftotalAK
    real :: errscd,errscdstrat
    real :: errvcdobs,errvcdstrat
    real :: errvcdamf,errvcdamfAK
    real :: errtroptotal_mix, errtotal_mix
    integer :: l

    ! On start, and at beginning of new month: read TM3 NO2 covariance matrix
    !he_cov  if( firstcall .or. newmonth )then
    !he_cov     call ReadCov( CovPath,month,no2mean,cov_matrix )
    !he_cov     firstcall = .false.
    !he_cov  end if  	

    status = 0

    ! initialise total slant column errors 
    errtotal_fcl = 0.0
    errtotal_pcl = 0.0
    errtotal_alb = 0.0
    no2vcdtotal = 0.0
    errtroptotal_fcl = 0.0
    errtroptotal_pcl = 0.0
    errtroptotal_alb = 0.0
    no2vcdtroptotal = 0.0

    do l = 1, lm

       call ErrFcl(cldfraction,rclear,rcloud, &
            amfclearlev(l),amfcloudlev(l),err_fcl,status)
       if ( status > 0 ) then
          IF_NOTOK_RETURN('ERROR in NO2ErrEstimate, ErrFcl') 
       end if

       call ErrPcl(l,lm,cldfraction,rclear,surfpres,cldtoppres,sza,vza,aza, &
            albcloud,amfclearlev(l),amfgeo,at,bt,amfLut,err_pcl,status)
       if ( status > 0 ) then
          IF_NOTOK_RETURN('ERROR in NO2ErrEstimate, ErrPcl') 
       end if

       call ErrAlb(l,lm,cldfraction,rcloud,surfpres,sza,vza,aza, &
            albclear,amfcloudlev(l),amfgeo,at,bt,amfLut,err_alb,status)
       if ( status > 0 ) then
          IF_NOTOK_RETURN('ERROR in NO2ErrEstimate, ErrAlb') 
       end if

       errtotal_fcl = errtotal_fcl + err_fcl*no2collev(l)
       errtotal_pcl = errtotal_pcl + err_pcl*no2collev(l)
       errtotal_alb = errtotal_alb + err_alb*no2collev(l)
       no2vcdtotal  = no2vcdtotal  +         no2collev(l)

       if(  l <= ltropopause )then
          errtroptotal_fcl = errtroptotal_fcl + err_fcl*no2collev(l)	  
          errtroptotal_pcl = errtroptotal_pcl + err_pcl*no2collev(l)
          errtroptotal_alb = errtroptotal_alb + err_alb*no2collev(l)
          no2vcdtroptotal  = no2vcdtroptotal  +         no2collev(l)
       end if 

    end do

    if ( no2vcdtotal < 1e-7 ) then
       IF_NOTOK_RETURN('ERROR in ass_err_retr.f90; no2vcdtotal is zero.') 
    end if
    errtotal_fcl = errtotal_fcl / no2vcdtotal
    errtotal_pcl = errtotal_pcl / no2vcdtotal
    errtotal_alb = errtotal_alb / no2vcdtotal

    if ( no2vcdtroptotal < 1e-7 ) then
       IF_NOTOK_RETURN('ERROR in ass_err_retr.f90; no2vcdtroptotal is zero')
    end if
    errtroptotal_fcl = errtroptotal_fcl / no2vcdtroptotal
    errtroptotal_pcl = errtroptotal_pcl / no2vcdtroptotal
    errtroptotal_alb = errtroptotal_alb / no2vcdtroptotal

    ! profile variation contribution, total column
    !he_cov  call ErrProfile(amflev,amftotal,lm, &
    !he_cov       no2mean,cov_matrix,ix4a,iy4a,w4a,errtotal_prof)
    !
    !he_cov   on average the profile related errors are about 10%
    if ( amftotal > 1.0e-7 ) then
       errtotal_prof = 0.1*amftotal
    else
       errtotal_prof = 999.0
    end if
    !
    ! profile variation contribution, troposphere only
    !he_cov  call ErrProfile(amflev,amftroptotal,ltropopause, &
    !he_cov       no2mean,cov_matrix,ix4a,iy4a,w4a,errtroptotal_prof)
    !
    !he_cov   on average the profile related errors are about 10%
    if ( amftroptotal > 0.1 ) then
       errtroptotal_prof = 0.1*amftroptotal
    else
       errtroptotal_prof = 999.0
    end if
    
    ! boundary layer mixing error contribution, troposphere
    !JDMERR call ErrMixing(lm,no2collev,at,bt,surfpres,amflev,&
    !JDMERR                 ltropopause,errtroptotal_mix)
    errtroptotal_mix = 0.0
        
    ! boundary layer mixing error contribution, total 
    !JDMERR call ErrMixing(lm,no2collev,at,bt,surfpres,amflev,&
    !JDMERR                 lm,errtotal_mix)
    errtotal_mix = 0.0

    ! Correlation between cloud fraction and albedo errors
    ! from Koelemeijer, JGR 106, 3475, 2001 (Fresco)    
    if ( abs(albcloud-albclear) < 1.0e-5 ) then          
       ! print*,'WARNING: cloud albedo and clear albedo have identical values'
       ! print*,'in ass_err_retr.f90.'
       ! print*,'Cloud albedo: ',albcloud,'Clear albedo: ',albclear
       fcl_alb_correlation =  1.0
    else	  
       fcl_alb_correlation = - ((1.0-cldfraction)*albedoError)/ &
                               ((albcloud-albclear)*frescoError)
       if ( fcl_alb_correlation >  1.0 ) fcl_alb_correlation =  1.0
       if ( fcl_alb_correlation < -1.0 ) fcl_alb_correlation = -1.0
    endif	      


    !print*,'Check 1',errtotal_fcl
    !print*,'Check 2',errtotal_pcl
    !print*,'Check 3',errtotal_alb
    !print*,'Check 4',errtotal_prof
    !print*,'Check 5',errtotal_mix
    !print*,'Check 6',fcl_alb_correlation

    ! compute total air-mass factor error 
    if (errtotal_fcl**2 + errtotal_pcl**2 + &
                        errtotal_alb**2 + errtotal_prof**2 + &
                        errtotal_mix**2 + 2.0 * fcl_alb_correlation * &
                        errtotal_fcl*errtotal_alb < 0.0 ) then

       PRINT('NO2ErrEstimate: Small error in total air-mass factor')
       erramftotal = sqrt( errtotal_fcl**2 + errtotal_pcl**2 + &
                           errtotal_alb**2 + errtotal_prof**2 + &
                           errtotal_mix**2 )			    
    else 
       erramftotal = sqrt( errtotal_fcl**2 + errtotal_pcl**2 + &
                        errtotal_alb**2 + errtotal_prof**2 + &
                        errtotal_mix**2 + 2.0 * fcl_alb_correlation * &
                        errtotal_fcl*errtotal_alb )
    endif			 

    ! compute total air-mass factor error 
    if (errtotal_fcl**2 + errtotal_pcl**2 + &
                          errtotal_alb**2 +  &
                          2.0 * fcl_alb_correlation * &
                          errtotal_fcl*errtotal_alb < 0.) then
       PRINT('NO2ErrEstimate: Small error in total AK error')
       erramftotalAK = sqrt( errtotal_fcl**2 + errtotal_pcl**2 + &
                             errtotal_alb**2)	  
    else			  
       erramftotalAK = sqrt( errtotal_fcl**2 + errtotal_pcl**2 + &
                          errtotal_alb**2 +  &
                          2.0 * fcl_alb_correlation * &
                          errtotal_fcl*errtotal_alb )
    endif			                

    ! compute total tropospheric air-mass factor error 
    if (errtroptotal_fcl**2 + errtroptotal_pcl**2 + &
                            errtroptotal_alb**2 + errtroptotal_prof**2 + &
                            errtroptotal_mix**2 + 2.0 * fcl_alb_correlation * &
                            errtroptotal_fcl*errtroptotal_alb < 0.) then
       PRINT('NO2ErrEstimate: Small error in tropospheric air-mass factor')
       erramftroptotal = sqrt( errtroptotal_fcl**2 + errtroptotal_pcl**2 + &
                               errtroptotal_alb**2 + errtroptotal_prof**2 + &
                              errtroptotal_mix**2)	    
    else			    		         
       erramftroptotal = sqrt( errtroptotal_fcl**2 + errtroptotal_pcl**2 + &
                            errtroptotal_alb**2 + errtroptotal_prof**2 + &
                            errtroptotal_mix**2 + 2.0 * fcl_alb_correlation * &
                            errtroptotal_fcl*errtroptotal_alb )
    endif			                  

    ! compute total tropospheric air-mass factor error, without profile term 
    if (errtroptotal_fcl**2 + errtroptotal_pcl**2 + &
                              errtroptotal_alb**2 + &
                              2.0 * fcl_alb_correlation * &
                              errtroptotal_fcl*errtroptotal_alb < 0.) then
       PRINT('NO2ErrEstimate: Small error in tropospheric air-mass factor')
       erramftroptotalAK = sqrt( errtroptotal_fcl**2 + errtroptotal_pcl**2 + &
                                 errtroptotal_alb**2)		      
    else		      
       erramftroptotalAK = sqrt( errtroptotal_fcl**2 + errtroptotal_pcl**2 + &
                              errtroptotal_alb**2 + &
                              2.0 * fcl_alb_correlation * &
                              errtroptotal_fcl*errtroptotal_alb )
    endif			                    

    ! error tropospheric vertical column as sum of air-mass factor error,
    ! slant column error and stratospheric slant column error

    ! measured slant column error, units 1e15 molec/cm2
    errscd = retrievalScdError  

    ! measured stratospheric slant column error, units 1e15 molec/cm2
    errscdstrat = modelStratVcdError*amfgeo

    ! error due to measurement uncertainty
    if ( amftroptotal > 0.1 ) then 
       errvcdobs = errscd/amftroptotal 
 
       ! error due to stratospheric reference
       errvcdstrat = errscdstrat/amftroptotal 

       ! error due to tropospheric AMF
       errvcdamf = erramftroptotal*satvcdtrop/amftroptotal

       ! error due to tropospheric AMF, when kernels are used
       errvcdamfAK = erramftroptotalAK*satvcdtrop/amftroptotal
    else
       errvcdobs = 999.0
       errvcdstrat = 999.0
       errvcdamf = 999.0
       errvcdamfAK = 999.0       
    end if   

    ! error in tropospheric column:

    ! tropospheric slant column error, units 1e15 molec/cm2
    if ( ( errvcdobs*errvcdobs + errvcdstrat*errvcdstrat + &
           errvcdamf*errvcdamf ) < 0.) then
       IF_NOTOK_RETURN('NO2ErrEstimate: ERROR in total tropospheric column calculation')
    else    
       errVcdTropTotal = sqrt( errvcdobs*errvcdobs + &
                               errvcdstrat*errvcdstrat + &
                               errvcdamf*errvcdamf )
    end if

    ! tropospheric slant column error, when kernels are used
    if ( ( errvcdobs*errvcdobs + errvcdstrat*errvcdstrat + &
           errvcdamfAK*errvcdamfAK) < 0. ) then
       IF_NOTOK_RETURN('NO2ErrEstimate: ERROR in total tropospheric column calculation')
    else
       errVcdTropTotalAK = sqrt( errvcdobs*errvcdobs + &
                                 errvcdstrat*errvcdstrat + &
                                 errvcdamfAK*errvcdamfAK )
    end if

    ! error in stratospheric column:

    errVcdStratTotal = modelStratVcdError

    ! error in total column:

    ! error due to the air-mass factor
    if ( amftotal > 1.0e-7 ) then
       errvcdamf = satvcd*erramftotal/amftotal

       ! error due to the air-mass factor
       errvcdamfAK = satvcd*erramftotalAK/amftotal

       ! error due to measurement uncertainty
       errvcdobs = errscd/amftotal 
    else
       errvcdamf = 999.0
       errvcdamfAK = 999.0
       errvcdobs = 999.0
    end if

    ! vertical total column error
    if ( ( errvcdobs*errvcdobs + errvcdamf*errvcdamf ) < 0. ) then
       IF_NOTOK_RETURN('NO2ErrEstimate: ERROR in total column error calculation')
    else
       errVcdTotal = sqrt( errvcdobs*errvcdobs + &
                           errvcdamf*errvcdamf )
    end if

    ! vertical total column error, when kernels are used
    if ( (errvcdobs*errvcdobs + errvcdamfAK*errvcdamfAK) < 0.) then
       IF_NOTOK_RETURN('NO2ErrEstimate: ERROR in total column error calculation')
    else
       errVcdTotalAK = sqrt( errvcdobs*errvcdobs + &
                             errvcdamfAK*errvcdamfAK )
    end if

  end subroutine NO2ErrEstimate


!*********************************************************************

  subroutine ErrFcl(cldfraction,rclear,rcloud,amfclear,amfcloud,err_fcl,status)
    !=======================================================================
    ! 
    ! Compute the uncertainty in the air-mass factor due to cloud uncertainties
    !
    ! Input:
    !    cldfraction  : cloud fraction for this pixel
    !    rclear       : radiance weight of clear part of pixel  
    !    rcloud       : radiance for cloudy part of pixel
    !    amfclear     : air-mass factor of 100% cloud-free pixel
    !    amfcloud     : air-mass factor of 100% cloudy pixel
    !        
    ! Output:
    !    err_fcl      : Uncertainty in AMF due to uncertainty in cloud fraction
    !    status       : > 0 when error occurred
    !
    !                             Folkert Boersma, KNMI, oct 2002
    !=======================================================================  

    implicit none

    ! in/out: 
    real, intent(in)     :: cldfraction, rclear, rcloud
    real, intent(in)     :: amfclear, amfcloud
    real, intent(out)    :: err_fcl  
    integer, intent(out) :: status

    ! local
    character(*), parameter    :: rname = 'ErrFcl'
    real               :: radsat1,radsat2
    ! amf for cloud fraction that is 0.05 too low
    real               :: amf_min,amf_plus
    real               :: cldfraction_min,cldfraction_plus

    cldfraction_min  = max(0.0,cldfraction-0.001)
    cldfraction_plus = min(1.0,cldfraction+0.001)
     
    radsat1=(1.0-cldfraction_min)*rclear+ &
         (cldfraction_min)*rcloud
    radsat2=(1.0-cldfraction_plus)*rclear+ &
         (cldfraction_plus)*rcloud

    ! determine new AMF with incorrect cloud fraction
    if ( radsat1 == 0. .or. radsat2 == 0.) then
       IF_NOTOK_RETURN('ERROR in err_fcl: radsat1 = 0.0')
    end if
    amf_min = (1.0-cldfraction_min)*rclear*amfclear/radsat1 &
         +cldfraction_min*rcloud*amfcloud/radsat1
    amf_plus = (1.0-cldfraction_plus)*rclear*amfclear/radsat2 &
         +cldfraction_plus*rcloud*amfcloud/radsat2
    ! determine random error in AMF due to cloud fraction uncertainty
    if ( (cldfraction_plus-cldfraction_min) < 1e-7) then
       IF_NOTOK_RETURN('ERROR in err_fcl: cldfraction_plus-cldfraction_min = 0.0')
    end if
    err_fcl = frescoerror*(amf_plus-amf_min)/ &
         (cldfraction_plus-cldfraction_min)
    !
  end subroutine ErrFcl

  
  subroutine ErrPcl(l,lm,cldfraction,rclear,surfpres,cloudpres_in,sza,vza,aza, &
       albcloud,amfclear,amfgeo,at,bt,amfLut,err_pcl,status)
    !=======================================================================
    ! 
    ! Compute the uncertainty in the air-mass factor due to cloud 
    !  pressure uncertainties
    !
    ! Input:
    !    l              : level index
    !    lm             : number of levels
    !    cloudfraction  : cloud pressure for this pixel 
    !    rclear         : radiance for cloud-free pixel
    !    surfpres       : surface pressure
    !    cloudpres_in   : cloud top pressure
    !    sza,vza,aza    : solar, viewing, azimuth angles
    !    albcloud       : cloud albedo
    !    amfclear       : air-mass factor of 100% cloud-free pixel
    !    amfgeo         : geometrical air-mass factor
    !    at,bt          : pressure level indices
    !    amfLut         : the NO2 air-mass factor LUT
    !
    ! Output:
    !    err_pcl        : Uncertainty in AMF due to cloud pressure
    !    status         : > 0 means error 
    !
    !                             Folkert Boersma, KNMI, oct 2002
    !=======================================================================  
    implicit none

    ! in/out     
    integer,intent(in) :: l, lm     
    real, intent(in)   :: cldfraction, rclear, cloudpres_in, surfpres, sza, vza, aza
    real, intent(in)   :: albcloud,amfclear,amfgeo
    real, dimension(lm+1), intent(in)  :: at,bt
    type(TAmfLut), intent(inout)    :: amfLut
    real, intent(out)  :: err_pcl    
    integer, intent(out) :: status
	  	
    ! local
    character(*), parameter    :: rname = 'ErrPcl'
    real :: cloudpres, cloudpres_plus, cloudpres_min, rcloud_plus, rcloud_min
    real :: amfcloud_plus,amfcloud_min,radsat_plus,radsat_min,amf_plus,amf_min
    real :: phybrTop,phybrBot,phybrid,amfcloudrel,abovecloudfraction

    status = 0

    cloudpres = cloudpres_in
    if ( cloudpres < 130.0 ) cloudpres = 130.0 
    
    ! Cloud pressure perturbation for finite differencing
    cloudpres_plus = min(surfpres,cloudpres + 10.0)
    
    cloudpres_min  = min(surfpres - 20.0, cloudpres - 10.0) 

    !if( cloudpres_plus > 1013 ) cloudpres_plus = 1013
    !if( cloudpres_min  > 1013 ) cloudpres_min = 1013   
    !if( cloudpres_plus < 130 )  cloudpres_plus = 130
    !if( cloudpres_min  < 130 )  cloudpres_min = 130

    ! estimate radiation intensities for cloud conditions
    call rweight( FitWindowCentre,cloudpres_plus,albcloud, &
                        sza, vza, aza, &
                        rcloud_plus )
    call rweight( FitWindowCentre,cloudpres_min,albcloud, &
                        sza, vza, aza, &
                        rcloud_min )
	  	    	  
    ! determine mid-level TM3 pressures (in hPa)
    phybrTop = 0.01*( at(l+1) + 100.0*surfpres*bt(l+1) )
    phybrBot = 0.01*( at(l)   + 100.0*surfpres*bt(l)   )
    phybrid  = 0.5*( phybrTop + phybrBot )

    ! read sensitivities from the AMF lookup table: 100% cloud case
    ! perturbed cloud top pressures
    if( cloudpres_plus >= phybrTop-1.0e-5 ) then
       if( cloudpres_plus >= phybrBot-1.0e-5 ) then
          ! model level above the cloud top
          call GetAmf( phybrid,aza,vza,sza, &
               albcloud,cloudpres_plus,amfLut,amfcloudrel,status )
          if ( status > 0 ) then
             IF_NOTOK_RETURN('ERROR returned by GetAmf, LUT range error') 
          end if
          amfcloud_plus=amfcloudrel*amfgeo
       else
          ! model level contains cloud top
          if (phybrBot-phybrTop < 1e-7 ) then
	     IF_NOTOK_RETURN('ERROR in err_pcl: phybrBot-phybrTop == 0.')
          endif
          abovecloudfraction = (cloudpres_plus-phybrTop)/(phybrBot-phybrTop)
          call GetAmf( phybrid,aza,vza,sza, &
               albcloud,phybrBot,amfLut,amfcloudrel,status )
          if ( status > 0 ) then
             IF_NOTOK_RETURN('ERROR returned by GetAmf, LUT range error') 
          end if
          !call GetAmf( phybrid,aza,vza,sza, &
          !     albcloud,phybrid+1.0e-5,amfcloudrel )
          amfcloud_plus=amfcloudrel*amfgeo*abovecloudfraction
       end if
    else
       ! model level below cloud top 
       amfcloud_plus=0.0 ! below the cloud
    end if
	  
    if( cloudpres_min >= phybrTop-1.0e-5 ) then
       if( cloudpres_min >= phybrBot-1.0e-5 ) then
          ! model level above the cloud top
          call GetAmf( phybrid,aza,vza,sza, &
               albcloud,cloudpres_min,amfLut,amfcloudrel,status )
          if ( status > 0 ) then
             IF_NOTOK_RETURN('ERROR returned by GetAmf, LUT range error') 
          end if
          amfcloud_min=amfcloudrel*amfgeo
       else
          ! model level contains cloud top
          if (phybrBot-phybrTop < 1e-7) then
	     IF_NOTOK_RETURN('ERROR in err_pcl: phybrBot-phybrTop == 0.')
	  end if
          abovecloudfraction = (cloudpres_min-phybrTop)/(phybrBot-phybrTop)
          call GetAmf( phybrid,aza,vza,sza, &
               albcloud,phybrBot,amfLut,amfcloudrel,status )
          if ( status > 0 ) then
             IF_NOTOK_RETURN('ERROR returned by GetAmf, LUT range error') 
          end if
          !call GetAmf( phybrid,aza,vza,sza, &
          !     albcloud,phybrid+1.0e-5,amfcloudrel )
          amfcloud_min=amfcloudrel*amfgeo*abovecloudfraction
       end if
    else
       ! model level below cloud top 
       amfcloud_min=0.0 ! below the cloud
    end if
	  
    ! estimated observed radiance (normalization factor)
    radsat_min =(1.0-cldfraction)*rclear+cldfraction*rcloud_min
    radsat_plus=(1.0-cldfraction)*rclear+cldfraction*rcloud_plus

    ! determine new AMF with incorrect cloud fraction
    if ( radsat_min == 0. .or. radsat_plus == 0. ) then
       IF_NOTOK_RETURN('ERROR in ErrPcl: radsat_min or radsat_plus = 0.')
    end if
    amf_min = (1.0-cldfraction)*rclear*amfclear/radsat_min &
         +cldfraction*rcloud_min*amfcloud_min/radsat_min
    amf_plus = (1.0-cldfraction)*rclear*amfclear/radsat_plus &
         +cldfraction*rcloud_plus*amfcloud_plus/radsat_plus
                    
    ! determine random error in AMF due to cloud pressure uncertainty    
    if ( cloudpres_plus-cloudpres_min < 1e-7 ) then
       err_pcl = 0.0
    else
       err_pcl = cloudtoperror*(amf_plus-amf_min)/ &
            (cloudpres_plus-cloudpres_min)
    end if	    
    	  
  end subroutine ErrPcl


  
  subroutine ErrAlb ( l,lm,cldfraction,rcloud,surfpres,sza,vza,aza, &
       albclear,amfcloud,amfgeo,at,bt,amfLut,err_alb,status )
    !=======================================================================
    ! 
    ! Compute the uncertainty in the air-mass factor due to  
    !  surface albedo uncertainties
    !
    ! Input:
    !    l              : level index
    !    lm             : number of levels
    !    cloudfraction  : cloud pressure for this pixel 
    !    rcloud         : radiance for cloudy part of pixel
    !    surfpres       : surface pressure
    !    sza,vza,aza    : solar, viewing, azimuth angles
    !    albclear       : clear albedo
    !    amfcloud       : air-mass factor of 100% cloudy pixel
    !    amfgeo         : geometrical air-mass factor
    !    at,bt          : pressure level indices
    !    amfLut         : the NO2 air-mass factor LUT
    !
    ! Output:
    !    err_alb        : Error in AMF due to albedo errors
    !    status         : > 0 means error 
    !
    !                             Folkert Boersma, KNMI, oct 2002
    !=======================================================================  
    implicit none
    ! in/out          
    integer,intent(in) :: l, lm     
    real, intent(in)   :: cldfraction,rcloud,surfpres,sza,vza,aza
    real, intent(in)   :: albclear,amfcloud,amfgeo
    real,dimension(lm+1), intent(in)  :: at, bt
    type(TAmfLut), intent(inout)    :: amfLut
    real, intent(out)  :: err_alb  
    integer, intent(out) :: status
	  	
    ! local
    character(*), parameter    :: rname = 'ErrAlb'
    real :: albedo_plus,albedo_min,rclear_plus,rclear_min
    real :: amfclear_plus,amfclear_min,radsat_plus,radsat_min,amf_plus,amf_min
    real :: phybrTop,phybrBot,phybrid

    status = 0
  	  
    ! limits check
    if ( albclear < -0.00001 .or. albclear > 1.00001 ) then
       IF_NOTOK_RETURN('ERROR n ErrAlb: albclear out of range')
    end if
    if ( l < 1 .or. l > lm ) then
       IF_NOTOK_RETURN('ERROR n ErrAlb: l out of range')
    end if

    ! Cloud pressure perturbation for finite differencing
    albedo_min  = max(0.0,albclear - 0.001)
    albedo_plus = min(1.0,albclear + 0.001)
    
    ! estimate radiation intensities for cloud-free conditions
    call rweight( FitWindowCentre,surfpres,albedo_min, &
                        sza, vza, aza, &
                        rclear_min )
    call rweight( FitWindowCentre,surfpres,albedo_plus, &
                        sza, vza, aza, &
                        rclear_plus )
	  	    	  
    ! determine mid-level TM3 pressures (in hPa)
    phybrTop = 0.01*( at(l+1) + 100.0*surfpres*bt(l+1) )
    phybrBot = 0.01*( at(l)   + 100.0*surfpres*bt(l)   )
    phybrid  = 0.5*( phybrTop + phybrBot )

    ! determine the clear amf for albedo perturbations
    call GetAmf( phybrid, aza, vza, sza, &
         albedo_min, surfpres, amfLut, amfclear_min,status )             
    if ( status > 0 ) then
       IF_NOTOK_RETURN('ERROR returned by GetAmf, LUT range error') 
    end if
    amfclear_min = amfclear_min*amfgeo	         
	  
    call GetAmf( phybrid, aza, vza, sza, &
         albedo_plus, surfpres, amfLut, amfclear_plus,status )             
    if ( status > 0 ) then
       IF_NOTOK_RETURN('ERROR returned by GetAmf, LUT range error') 
    end if
    amfclear_plus = amfclear_plus*amfgeo	         

    ! estimated observed radiance (normalization factor)
    radsat_min =(1.0-cldfraction)*rclear_min  + cldfraction*rcloud
    radsat_plus=(1.0-cldfraction)*rclear_plus + cldfraction*rcloud

    ! determine new AMF with incorrect cloud fraction
    if ( abs(radsat_min)< 0.00001 .or. abs(radsat_plus) < 0.00001 ) then
       IF_NOTOK_RETURN('ERROR in ErrAlb: radsat_min or radsat_plus = 0')
    end if
    amf_min = (1.0-cldfraction)*rclear_min*amfclear_min/radsat_min &
         +cldfraction*rcloud*amfcloud/radsat_min
    amf_plus = (1.0-cldfraction)*rclear_plus*amfclear_plus/radsat_plus &
         +cldfraction*rcloud*amfcloud/radsat_plus
                    
    ! determine random error in AMF due to cloud pressure uncertainty
    if ( albedo_plus-albedo_min < 1e-7) then
       IF_NOTOK_RETURN('ERROR in ErrAlb: albedo_plus-albedo_min = 0')
    end if
    err_alb = albedoError*(amf_plus-amf_min)/ &
         (albedo_plus-albedo_min)
    !	  
  end subroutine ErrAlb
  

  subroutine ErrMixing( lm,no2collev,at,bt,surfpres,amflev,lmax,errTotalMix,status )
    !=======================================================================
    ! 
    ! Compute an estimate of the air-mass factor error related to
    ! boundary layer mixing errors in the model. A simple estimate is 
    ! obtained by comparing the air-mass factor with constant mixing 
    ! ratio's in the lowest 4 TM3 layers ("well mixed") with a case where 
    ! all no2 of the lowest 4 layers is concentrated in layer 1 ("no mixing") 
    !
    ! Input:
    !    lm             : number of vertical levels
    !    no2collev      : original no2 profile (vector)
    !    at             : pressure parameter
    !    bt             : pressure parameter
    !    surfpres	      : surface pressure 
    !    amflev         : height-dependent air-mass factor (vector) 
    !    lmax           : top level for air-mass computation          
    ! Output:
    !    errTotalMix    : Error in AMF due to mixing errors in TM3
    !    status         : > 0 on error
    !
    !                             KFB & HJE, KNMI, March 19, 2003
    !=======================================================================
    
    implicit none

    ! in/out            
    integer,intent(in)                :: lm
    integer,intent(in)                :: lmax
    real,dimension(lm),intent(in)     :: no2collev,amflev
    real,dimension(lm+1),intent(in)   :: at,bt
    real,intent(in)                   :: surfpres
    real,intent(out)                  :: errTotalMix  
    integer, intent(out)              :: status

    ! local
    character(*), parameter    :: rname = 'ErrMixing'
    integer              :: l
    real,dimension(6)    :: p, d_p
    real                 :: d_p_layer, no2layer
    real                 :: no2slcfcPlus, no2vcdfcPlus
    real                 :: no2slcfcMin, no2vcdfcMin
    real                 :: amftotalPlus,amftotalMin

    real, dimension(:), allocatable  :: no2_mixed, no2_bottom

    ! begin code

    !if ( lm /= 19 ) then
    !   print*,'ERROR ErrMixing: routine only valid for 19-layer model'
    !   stop
    !end if

    if ( lm /= 31 .and. lm /= 38 .and. lm /= 35 .and. lm /= 34 ) then
       IF_NOTOK_RETURN('ERROR ErrMixing: routine only valid for 31, 34, 35 or 38 layer model')
    end if

    allocate ( no2_mixed(lm) )
    allocate ( no2_bottom(lm) )

    no2slcfcPlus = 0.0
    no2vcdfcPlus = 0.0
    no2slcfcMin = 0.0
    no2vcdfcMin = 0.0
    no2layer = 0.0

    no2layer = SUM(no2collev(1:5))

    do l=1,6
       ! compute pressure at bottom of layer          
       p(l) = 0.01*( at(l) + 100.0*surfpres*bt(l) )
    end do
    
    ! compute pressure difference over one layer, and 5 layers
    do l=1,5            
       d_p(l) = p(l) - p(l+1)
    end do
    d_p_layer = p(1) - p(6)

    ! CASE 1 : lowest 5 model layers well mixed

    do l = 1, lmax
       if ( l < 6 ) then 
          ! Redistribute NO2 in lowest 5 layers
          if (d_p_layer == 0.) then
	     IF_NOTOK_RETURN('ERROR in ass_err_retr.f90; d_p_layer = 0.')
	  end if
	  no2_mixed(l) = no2layer * (d_p(l)/d_p_layer)
       else
          ! Above layer 5, leave things the same
          no2_mixed(l) = no2collev(l) 
       end if

       ! compute model-predicted slant column
       ! no2_mixed = NO2 mass in grid cell (units 10^15 molecules cm^-2)
       no2slcfcPlus = no2slcfcPlus + amflev(l)*no2_mixed(l)         

       ! model vertical column (10^15 molec cm^-2)
       no2vcdfcPlus = no2vcdfcPlus + no2_mixed(l)         
    end do

    ! total air-mass factor for "mixed BL"
    if ( no2vcdfcPlus < 1.0e-7 ) then
       IF_NOTOK_RETURN('ERROR in ass_err_retr.f90; no2vcdfcPlus = 0.')
    end if
    amftotalPlus = no2slcfcPlus/no2vcdfcPlus      ! model SLC / model VCD

    ! CASE 2 : NO2 of lowest 5 model layers inserted in layer 1 ("no mixing")

    no2_bottom(1)      = no2layer
    no2_bottom(2:5)    = 0.0
    no2_bottom(6:lmax) = no2collev(6:lmax)

    ! compute model-predicted slant column
    ! no2_bottom = NO2 mass in grid cell (units 10^15 molecules cm^-2)
    do l = 1, lmax         
       ! model slant column 
       no2slcfcMin = no2slcfcMin + amflev(l)*no2_bottom(l)         

       ! model vertical column (10^15 molec cm^-2)
       no2vcdfcMin = no2vcdfcMin + no2_bottom(l)      
    end do

    ! total air-mass factor for non-mixed BL
    if ( no2vcdfcMin < 1.0e-7 ) then
       IF_NOTOK_RETURN('ERROR in ass_err_retr.f90; no2vcdfcMin = 0.')
    end if
    amftotalMin = no2slcfcMin/no2vcdfcMin     ! model SLC / model VCD

    ! 2 sigma is difference between two AMF's
    errTotalMix = 0.5*(amftotalPlus - amfTotalMin)       

    deallocate ( no2_mixed )
    deallocate ( no2_bottom )

  end subroutine ErrMixing



end module MNO2RetrievalError