! First include the set of model-wide compiler flags #include "tm5.inc" module MOmFSuper ! ! Cluster NO2 observations and observation operators into ! superobservation for TM gridcells ! ! Henk Eskes, KNMI, 2002 ! modified Nov 2004, aug 2015 implicit none private public :: GetOMFHxSuper, TObsInfoReduced public :: ObsInfoReduced_allocate, ObsInfoReduced_deallocate ! ---------------------------------------------------- ! "TObsInfoReduced" contains the ! observation, forecast and the observation operator ! for the superobservations ! ! count : number of grid cells with NO2 observations ! obsi : pixel longitude grid index ! obsj : pixel latitude grid index ! nomf : Nr of obs in grid cell ! slcobs : grid-mean slant column observation ! slcfc : grid-mean slant column forecast ! vcdfc : grid-mean column forecast ! sigobs : grid-mean observation error ! hobs : observation operator (kernel) vector for ! each observed grid cell ! ! NOTE: the slant columns are divided by the geometrical air-mass factor ! ---------------------------------------------------- type TObsInfoReduced integer :: count integer,dimension(:),allocatable :: obsi, obsj, nomf real,dimension(:),allocatable :: slcobs, slcfc, vcdfc, sigobs real,dimension(:,:),allocatable :: hobs end type TObsInfoReduced ! the "usePixel" flag is shared between the forecast and analysis superobservations ! it list the individual observations used to construct the superobservations logical, dimension(:), allocatable :: usePixel contains subroutine GetOMFHxSuper( im, jm, lm, no2Tr, ObsFcInfo, ObsInfoReduced, anFc, computeRetrievalForAnalysis ) !======================================================================= ! ! GetOMFHxReduced ! Determine a reduced set of Observation minus Forecast ! innovations and observation operators, defined on the ! model grid ! ! Data is stored in structure "ObsInfoReduced" ! ! Input: ! im, jm, lm : model grid dimensions ! no2Tr : Observation info: geometry, scd etc ! ObsFcInfo : Observation, forecast and observation operator, for ! all pixels (forecast or analysis) ! AnFc : 'a' = analysis, 'f' = forecast ! ! Output: ! "ObsInfoReduced" - of type "TObsInfoReduced" ! Contains grid cell indices, OmF innovations, ! the observation errors and ! the corresponding observation operator vectors ! the superobservations are slant columns / geometrical air mass factor ! ! ! Henk Eskes, KNMI, 2002 !======================================================================= use MTObsTrack, only : TObsTrack use MTObsFcInfo, only : TObsFcInfo implicit none ! in/out integer, intent(in) :: im, jm, lm type(TObsTrack),intent(in) :: no2Tr type(TObsFcInfo),intent(in) :: ObsFcInfo type(TObsInfoReduced),intent(inout) :: ObsInfoReduced character(1), intent(in) :: anFc logical, intent(in) :: computeRetrievalForAnalysis ! parameters to accept/reject observations real,parameter :: maxSZA = 80.0 ! maximum solar zenith angle real,parameter :: maxLatRelative = 5 ! no assimilation near Npole !JDM1132 !real,parameter :: maxLatIndex = jm-2 ! no assimilation near Npole !JDM1132 real,parameter :: minLatIndex = 5 ! no assimilation near Spole !JDM1132 ! real,parameter :: minLatIndex = 3 ! no assimilation near Spole !JDM1132 ! real,parameter :: maxCloudFrac = 0.8 ! maximum cloud fraction real,parameter :: maxCloudFrac = 1.01 ! maximum cloud fraction real,parameter :: minCloudFrac = -0.01 ! minimum cloud fraction !real,parameter :: maxOmFdiff = 4.0 ! maximum difference VCD obs-forec (1e15 molec/cm2) real,parameter :: maxOmFdiff = 10.0 ! maximum difference VCD obs-forec (1e15 molec/cm2) ! stratos slc error / amf_geo (1e15 molec/cm2) ! real,parameter :: strat_error = 0.25 ! old value, resulting in (a-f)/(o-f) about 30% real,parameter :: strat_error = 0.15 ! corresponding to a superobs SLC noise level of about 0.4 for an AMF of 2.5 ! tropos slc error / amf_geo (1e15 molec/cm2) real,parameter :: trop_error = 4.00 ! estimated average observation correlation inside a TM grid box real,parameter :: obscorfac=1.0 ! el. [0..1] ! local integer, dimension(:,:), allocatable :: grid_index integer, dimension(:,:), allocatable :: grid_count ! integer :: i, j, l, ipix, ipixred, ipixredmax, iplusj integer :: nnegative, nbadpix, nlargeslcsig, nsza, npole, ncloud real :: sig,snorm,sigcor real :: omf,omfm1,omfm2,slcgeotrop,slcgeo real :: omf_m1_pixels, omf_m2_pixels integer :: omf_count real :: maxLatIndex ! logical :: QualityFlag ! ! observation selection - remove the following: ! 1) high solar zenith angles ! 2) high cloud fractions ! 3) negative observations ! 4) Pixels with errors ! ! Consistency check: both structures should contain the same nr of observations if ( ObsFcInfo%count /= no2Tr%count ) then stop 'ERROR GetOMFHxSuper, number of observations inconsistent !' end if if ( anFc == 'f' ) then ! UsePixel is used for both the forecast and analysis sweep, and is allocated at forecast print *, 'Allocate usePixel' allocate ( usePixel(no2Tr%count) ) end if ! Determine which pixels to use for the superobs (usePixel) and provide feedback ! on how many pixels were removed if ( anFc == 'f' ) then maxLatIndex = jm + 1 - maxLatRelative nsza = 0 npole = 0 ncloud = 0 nnegative = 0 nlargeslcsig = 0 nbadpix = 0 print *, 'GetOMFHxSuper: WARNING, keeping only 1 in 2 superobs, i+j=odd' do ipix = 1, no2Tr%count OmF = (no2Tr%no2slc(ipix)-ObsFcInfo%no2slcfc(ipix)) & / ObsFcInfo%amfgeo(ipix) QualityFlag = ( Abs(OmF) <= maxOmFdiff ) i = ObsFcInfo%icell(ipix) j = ObsFcInfo%jcell(ipix) ! Speed up: Only assimilate 1 in 2 grid cells, sum(i+j) odd, checkerboard iplusj = i+j+1 usePixel(ipix) = & ( no2Tr%solarZenithAngle(ipix) < maxSZA ) .and. & ( ObsFcInfo%jcell(ipix) <= maxLatIndex ) .and. & ( ObsFcInfo%jcell(ipix) >= minLatIndex ) .and. & ( no2Tr%cloudFraction(ipix) <= maxCloudFrac ) .and. & ( no2Tr%cloudFraction(ipix) >= minCloudFrac ) .and. & ( no2Tr%no2slc(ipix) > 0.0 ) .and. & ( no2Tr%no2slcError(ipix) < 100.0 ) .and. & ( iand(ObsFcInfo%flag(ipix), 255) == 0 ) .and. & ( mod(iplusj,2) == 0 ) .and. & QualityFlag if ( no2Tr%solarZenithAngle(ipix) >= maxSZA ) nsza = nsza + 1 if ( ( ObsFcInfo%jcell(ipix) > maxLatIndex ) .or. & ( ObsFcInfo%jcell(ipix) < minLatIndex ) ) npole = npole + 1 if ( ( no2Tr%cloudFraction(ipix) > maxCloudFrac ) .or. & ( no2Tr%cloudFraction(ipix) < minCloudFrac ) ) ncloud = ncloud + 1 if ( no2Tr%no2slc(ipix) <= 0.0 ) nnegative = nnegative + 1 if ( no2Tr%no2slcError(ipix) >= 100.0 ) nlargeslcsig = nlargeslcsig + 1 if ( .not. QualityFlag ) nbadpix = nbadpix + 1 end do ! warning/info messages if ( nsza > 0 ) then write(*,'(a,i5,a,i5,a)')' GetOMFHxSuper: ', & nsza,' of ',ObsFcInfo%count,' pixels have too large SZA' end if if ( npole > 0 ) then write(*,'(a,i5,a,i5,a)')' GetOMFHxSuper: ', & npole,' of ',ObsFcInfo%count,' pixels rejected: close to the poles' end if if ( ncloud > 0 ) then write(*,'(a,i5,a,i5,a)')' GetOMFHxSuper: ', & ncloud,' of ',ObsFcInfo%count,' pixels with cloud fraction outside range' end if if ( nnegative > 0 ) then write(*,'(a,i5,a,i5,a)')' GetOMFHxSuper: ', & nnegative,' of ',ObsFcInfo%count,' pixels are negative' end if if ( nlargeslcsig > 0 ) then write(*,'(a,i5,a,i5,a)')' GetOMFHxSuper: ', & nlargeslcsig,' of ',ObsFcInfo%count,' pixels have error > 100% ' end if if ( nbadpix > 0 ) then write(*,'(a,i5,a,f6.2)')' GetOMFHxSuper: Quality control, ', & nbadpix,' pixels found with abs(OmF) > ',maxOmFdiff end if end if ! ( anFc == 'f' ) ! At this point the list of pixels to use ("usePixel") is known ! Write mean OmF statistics omf_m1_pixels = 0.0 omf_m2_pixels = 0.0 omf_count = 0 do ipix = 1, no2Tr%count if ( usePixel(ipix) ) then omf = (no2Tr%no2slc(ipix)-ObsFcInfo%no2slcfc(ipix)) & / ObsFcInfo%amfgeo(ipix) omf_m1_pixels = omf_m1_pixels + omf omf_m2_pixels = omf_m2_pixels + omf*omf omf_count = omf_count + 1 end if end do omf_m1_pixels = omf_m1_pixels / real(max(1,omf_count)) omf_m2_pixels = omf_m2_pixels / real(max(1,omf_count)) print '(a,i8)', ' GetOMFHxSuper: number of pixels used to build superobs =',omf_count if ( anFc == 'f' ) then write(*,103) omf_m1_pixels, sqrt(omf_m2_pixels-omf_m1_pixels*omf_m1_pixels) 103 format(' GetOMFHxSuper: = ',f8.4,' std = ',f8.4,' (individual pixels)' ) else write(*,104) omf_m1_pixels, sqrt(omf_m2_pixels-omf_m1_pixels*omf_m1_pixels) 104 format(' GetOMFHxSuper: = ',f8.4,' std = ',f8.4,' (individual pixels)') end if ! determine number of pixels in each cell, and total number of superobservations allocate ( grid_index(im,jm) ) allocate ( grid_count(im,jm) ) print *,'grid index and count allocated' grid_index(:,:) = 0 grid_count(:,:) = 0 ipixredmax = 0 do ipix = 1, no2Tr%count if ( usePixel(ipix) ) then if ( iand(ObsFcInfo%flag(ipix), 255) == 0 ) then i = ObsFcInfo%icell(ipix) j = ObsFcInfo%jcell(ipix) if ( grid_count(i,j) == 0 ) then ipixredmax = ipixredmax + 1 grid_index(i,j) = ipixredmax end if grid_count(i,j) = grid_count(i,j) + 1 else ! the analysis may have more flagged pixels than the forecast ! discard also those pixels usePixel(ipix) = .false. end if end if end do write (*, 105) ipixredmax 105 format(' GetOMFHxSuper: number of superobs =',i6) ! Number of superobs is now known: allocate storage call ObsInfoReduced_allocate ( ipixredmax, lm, ObsInfoReduced ) ! calculate obs, forecast and observation operators for superobservations ObsInfoReduced%count = ipixredmax ObsInfoReduced%nomf(:) = 0 ObsInfoReduced%slcobs(:) = 0.0 ObsInfoReduced%slcfc(:) = 0.0 ObsInfoReduced%vcdfc(:) = 0.0 ObsInfoReduced%sigobs(:) = 0.0 ObsInfoReduced%hobs(:,:) = 0.0 do ipix = 1, no2Tr%count if ( usePixel(ipix) ) then i = ObsFcInfo%icell(ipix) j = ObsFcInfo%jcell(ipix) ipixred = grid_index(i,j) ObsInfoReduced%obsi(ipixred) = i ObsInfoReduced%obsj(ipixred) = j ObsInfoReduced%nomf(ipixred) = & ObsInfoReduced%nomf(ipixred) + 1 ObsInfoReduced%slcobs(ipixred) = & ObsInfoReduced%slcobs(ipixred) + & no2Tr%no2slc(ipix) / ObsFcInfo%amfgeo(ipix) ObsInfoReduced%slcfc(ipixred) = & ObsInfoReduced%slcfc(ipixred) + & ObsFcInfo%no2slcfc(ipix) / ObsFcInfo%amfgeo(ipix) ObsInfoReduced%vcdfc(ipixred) = & ObsInfoReduced%vcdfc(ipixred) + & ObsFcInfo%no2vcdfc(ipix) ! attribute large error to tropospheric part of the ! observation operator slcgeotrop = ObsFcInfo%no2slcfctrop(ipix) slcgeo = ObsFcInfo%no2slcfc(ipix) if ( slcgeo <= 0. ) then print*,'ERROR in ass_Hx_super.f90; slcgeo negative.' stop end if sig = ( trop_error * slcgeotrop + & strat_error * ( slcgeo - slcgeotrop ) ) / slcgeo ObsInfoReduced%sigobs(ipixred) = & ObsInfoReduced%sigobs(ipixred) + sig do l = 1, lm ObsInfoReduced%hobs(l,ipixred) = & ObsInfoReduced%hobs(l,ipixred) + & ObsFcInfo%avkernel(l,ipix)*ObsFcInfo%amf(ipix) / & ObsFcInfo%amfgeo(ipix) end do end if end do ! normalise slc, standard deviation and observation operator ! determine O-F statistics omfm1 = 0.0 omfm2 = 0.0 do ipixred = 1, ipixredmax snorm = 1.0/real(max(1,ObsInfoReduced%nomf(ipixred))) ObsInfoReduced%slcobs(ipixred) = ObsInfoReduced%slcobs(ipixred)*snorm ObsInfoReduced%slcfc(ipixred) = ObsInfoReduced%slcfc(ipixred)*snorm ObsInfoReduced%vcdfc(ipixred) = ObsInfoReduced%vcdfc(ipixred)*snorm sigcor = sqrt( (obscorfac+(1.0-obscorfac)*snorm) ) ObsInfoReduced%sigobs(ipixred) = & ObsInfoReduced%sigobs(ipixred)*snorm*sigcor do l = 1, lm ObsInfoReduced%hobs(l,ipixred) = & ObsInfoReduced%hobs(l,ipixred)*snorm end do omf = ObsInfoReduced%slcobs(ipixred)-ObsInfoReduced%slcfc(ipixred) omfm1 = omfm1+omf omfm2 = omfm2+omf*omf end do ! omfm1 = omfm1 / real(max(1,ipixredmax)) omfm2 = omfm2 / real(max(1,ipixredmax)) if ( anFc == 'f' ) then write(*,1031) omfm1, sqrt(omfm2-omfm1*omfm1) 1031 format(' GetOMFHxSuper: = ',f8.4,' std = ',f8.4,' (superobs)') else write(*,1041) omfm1, sqrt(omfm2-omfm1*omfm1) 1041 format(' GetOMFHxSuper: = ',f8.4,' std = ',f8.4,' (superobs)') end if ! release storage deallocate ( grid_index ) deallocate ( grid_count ) if ( computeRetrievalForAnalysis ) then ! only deallocate "usePixel" after the analysis step if ( anFc /= 'f' ) then print *, 'Deallocate usePixel' deallocate ( usePixel ) end if else ! Forecast retrieval only: deallocate print *, 'Deallocate usePixel' deallocate ( usePixel ) end if end subroutine GetOMFHxSuper subroutine ObsInfoReduced_allocate( NObsSuper, lm, ObsInfoReduced ) ! ! Create storage for the superobservations in (empty) ObsInfoReduced ! Henk Eskes, KNMI, Aug 2015 ! implicit none ! in/out integer, intent(in) :: NObsSuper, lm type(TObsInfoReduced),intent(inout) :: ObsInfoReduced ! start code if ( allocated(ObsInfoReduced%hobs) ) & stop 'ERROR subroutine ObsInfoReduced_allocate: already allocated!' allocate ( ObsInfoReduced%obsi(NObsSuper) ) allocate ( ObsInfoReduced%obsj(NObsSuper) ) allocate ( ObsInfoReduced%nomf(NObsSuper) ) allocate ( ObsInfoReduced%slcobs(NObsSuper) ) allocate ( ObsInfoReduced%slcfc(NObsSuper) ) allocate ( ObsInfoReduced%vcdfc(NObsSuper) ) allocate ( ObsInfoReduced%sigobs(NObsSuper) ) allocate ( ObsInfoReduced%hobs(lm,NObsSuper) ) end subroutine ObsInfoReduced_allocate subroutine ObsInfoReduced_deallocate( ObsInfoReduced ) ! ! Remove storage for the superobservations in (empty) ObsInfoReduced ! Henk Eskes, KNMI, Aug 2015 ! implicit none ! in/out type(TObsInfoReduced),intent(inout) :: ObsInfoReduced ! start code if ( allocated ( ObsInfoReduced%obsi ) ) deallocate ( ObsInfoReduced%obsi ) if ( allocated ( ObsInfoReduced%obsj ) ) deallocate ( ObsInfoReduced%obsj ) if ( allocated ( ObsInfoReduced%nomf ) ) deallocate ( ObsInfoReduced%nomf ) if ( allocated ( ObsInfoReduced%slcobs ) ) deallocate ( ObsInfoReduced%slcobs ) if ( allocated ( ObsInfoReduced%slcfc ) ) deallocate ( ObsInfoReduced%slcfc ) if ( allocated ( ObsInfoReduced%vcdfc ) ) deallocate ( ObsInfoReduced%vcdfc ) if ( allocated ( ObsInfoReduced%sigobs ) ) deallocate ( ObsInfoReduced%sigobs ) if ( allocated ( ObsInfoReduced%hobs ) ) deallocate ( ObsInfoReduced%hobs ) end subroutine ObsInfoReduced_deallocate end module MOmFSuper