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

! 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: <o-f> = ',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: <o-a> = ',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: <o-f> = ',f8.4,'  std = ',f8.4,'  (superobs)')
    else
       write(*,1041) omfm1, sqrt(omfm2-omfm1*omfm1)
1041   format(' GetOMFHxSuper: <o-a> = ',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