! 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