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

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

! Macro
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if

! integer, parameter :: sp = selected_real_kind(6, 37)
! integer, parameter :: dp = selected_real_kind(15, 307)

!=====================================================================
! 
! Central NO2 column data assimilation module
!
! Interface call by TM5 (in assim_interface_mod.f90):
!
!         call AssimNO2( ntimestep )
!
! The subroutine AssimNO2 performs a data assimilation time step.
!
!                                        Henk Eskes, KNMI, January 2003
!                                           Version 1.0 - Februari 2003
!
!                                      J.D. Maasakkers, KNMI, July 2013
!                                              Prototype v3 - July 2013
!
!                       Henk Eskes, Mark ter Linden, KNMI, Mar-Sep 2015
!                               Interface to TM5-mp-CB05 and to TROPOMI
!======================================================================

module MAssimilation

  ! printing, error handling
  use GO,               only : gol, goPr, goErr
  ! NO2 AMF lookup table
  use MTAmf,            only : TAmfLut

  implicit none

  private

  public :: AssimNO2

  !-------------------------------------------------------------

  ! integer to count the number of calls to the routine "AssimNO2" 
  integer :: callcount  = 0
  integer :: assimcount = 0

  ! module name
  character(len=*), parameter  :: mname = 'MAssimilation'

  ! keep filenames in memory
  character(256)               :: filename_S5P_L2_NO2_filelist
  ! character(256)               :: terrain_filename
  character(256)               :: amf_filename
  logical                      :: divideByAMFGeo
  ! output directory for OmF netcdf files
  character(256)               :: outputdir

  ! NO2 AMF lookup table (keep in memory)
  type(TAmfLut)                :: amfLut

  ! If true, call the retrieval for both forecast and analysis fields
  logical  ::  computeRetrievalForAnalysis 

  ! Switch: optional to write OmF information to separate netcdf file
  logical  ::  doWriteOmF

  ! Debug: extra retrieval output
  logical  ::  debugMode

contains


  subroutine AssimNO2( rcFile, ntimestep, TM5Data, nObs, noxscaling, status )
    !======================================================================
    !
    !   Assimilation of NO2 slant column data
    !   Main routine from where all assimilation routines are called
    !
    !   rcFile     : name of the rc (settings) file
    !   ntimestep  : timestep in seconds between calls to AssimNO2
    !   TM5Data    : All TM5 fields and grid info from the interface module
    !   nObs       : nr. of observations assimilated in this timestep
    !   noxscaling : an/fc ratio for the TM5 NOx transported field
    !			                   Henk Eskes, KNMI, 2000-2002-2015
    !======================================================================

    ! --- TM5 interfaces ---

    ! TM5 model arrays and variables:
    use MTM5Data,               only : TTM5Data 
    ! TM5 date and time routines:
    use dims,			only : itau
    ! use datetime,      		only : tau2date
    use GO,                     only : TrcFile, Init, Done, ReadRc

    ! --- TROPOMI L2 interfaces ---

    ! Observational input (slant columns) is stored in this structure
    use MTObsTrack,             only : TObsTrack, ObsTrackDeallocate
    ! Storage for retrieval results (vertical columns)
    use MTObsFcInfo,            only : TObsFcInfo, ObsFcInfoAllocate, ObsFcInfoDeallocate
    ! module to read OMI NO2 slant column data (by tracks)
    use MObsRead,               only : ObsFileGet
    ! module to write OMI NO2 vertical column retrieval data (by tracks)
    use TM5IF_module,           only : TM5IF_obsfcinfo_write

    ! --- DOMINO interfaces ---

    ! compute model forecast and observation operator
    use MObservationOperator,   only : InitOMFHx, GetOMFHx
    ! construct superobservations
    use MOmFSuper,              only : GetOMFHxSuper, TObsInfoReduced, ObsInfoReduced_deallocate
    ! Kalman Filter routines
    use oi_module,              only : AssimilateColumn
    ! Terrain height module
    ! use terrain,                only : terrain_read
    ! Error module, use the model vertical column forecast error estimate 
    use MNO2RetrievalError,     only : forecastError
    ! NO2 AMF lookup table
    use MAmfRead,               only : ReadAmfLut
    ! Write OmF data to file 
    use MOmF_output,            only : OmF_Output
    ! Cloud retrieval correction
    ! use cloudcor,		only : cloudcor_read
    
    implicit none

    ! --- in-out ---

    character(*), intent(in)         ::  rcFile
    integer, intent(in)              ::  ntimestep
    type(TTM5Data), intent(inout)    ::  TM5Data 
    integer, intent(out)             ::  nObs
    real, dimension(TM5Data%im, TM5Data%jm, TM5Data%lm), intent(inout)  ::  noxscaling
    integer, intent(out)             ::  status

    ! ---  local variables ---

    character(len=*), parameter ::  rname = mname//'/AssimNO2'

    type(TrcFile)    ::  rcF

    ! the ratio between analysis and forecast (3D field)
    real, dimension(:,:,:), allocatable :: AnFcRatio
    ! array of 3D partial NO2 columns per gridbox 
    real, dimension(:,:,:), allocatable :: no2pcf
    ! array of 2D NO2 columns per lat-lon gridbox 
    real, dimension(:,:), allocatable   :: no2col_fc, no2col_an

    ! the 2D model column forecast error distribution 
    real, dimension(:,:), allocatable   :: NO2VCForecastError 

    ! ---------------------------------------------------------	
    ! Storage of track data: NO2 and cloud info
    type(TObsTrack)       :: no2Tr
    ! observations and model forecast
    type(TObsFcInfo)      :: obsFcInfo, obsAnInfo
    ! clustered superobservations
    type(TObsInfoReduced) :: obsFcInfoReduced, obsAnInfoReduced

    ! local
    integer                          :: i, l
    real                             :: slc,slcsig,tropfraction
    logical                          :: obsFound, isForecast
    integer,parameter                :: omf_tape=83552

    ! OMIFilelist !JDM
    ! integer            :: unit, ierror, system
    ! character(len=160) :: listfile, command

    ! --- begin code ---

    callcount = callcount + 1

    write(*,1111) callcount
1111   format(' ------- assimno2 ',i4,'  -------')

    ! -------------- initialise -------------------

    if ( callcount == 1 ) then

       ! open rcfile:
       call Init( rcF, rcfile, status )
       IF_NOTOK_RETURN(status=1)

       ! Read the location of the file which contains a list of L2 NO2 files
       call ReadRc( rcF, 'domino.inputlist', filename_S5P_L2_NO2_filelist, status, default='list' )
       IF_ERROR_RETURN(status=1)

       ! Read the location of the orography file
       ! call ReadRc( rcF, 'domino.terrain_filename', terrain_filename, status, default='dem.nc' ) 

       ! Read the location of the AMF lookup table file
       call ReadRc( rcF, 'domino.amf_file', amf_filename, status, default='no2_amf_lut.nc' )

       ! Divide the AMF by AMFgeo after reading 
       call ReadRc( rcF, 'domino.amf.divide.by.amfgeo', divideByAMFGeo, status, default=.false. )

       ! Optional to call observation operator after the assimilation step 
       call ReadRc( rcF, 'domino.compute.retrieval.for.analysis', computeRetrievalForAnalysis, status, default=.false. )

       ! Write OmF statistics to netcdf file
       call ReadRc( rcF, 'domino.write.OmF', doWriteOmF, status, default=.true. )

       ! Debug mode yes/no
       call ReadRc( rcF, 'domino.debug.mode', debugMode, status, default=.false. )

       ! Output directory (for OmF netcdf output)
       call ReadRc( rcF, 'output.dir', outputdir, status, default='.' )

       ! close rcfile:
       call Done( rcF, status )
       IF_NOTOK_RETURN(status=1)

       ! Fill error specification matrices
       write(*,112) ForecastError
112    format('  InitForecastError: Initialising forecast error to ',f6.2,' 1e15 molec/cm2')

       ! Read settings for the observation operator
       call InitOMFHx( rcFile )
       
       ! Here we read the 90Arc terrain height database to be able to 
       ! calculate the average terrain height of each pixel. !JDMTH
       ! call terrain_read ( terrain_filename )

       ! Read the air-mass factor lookup table from file
       print *,' Reading NO2 AMF from file : ',trim(amf_filename)
       print *,'          divide by amfgeo : ',divideByAmfGeo
       call ReadAmfLut ( amf_filename, divideByAmfGeo, amfLut, status )

       print *,' S5P list file : ',trim(filename_S5P_L2_NO2_filelist)

    end if

    ! 2D forecast error field
    allocate ( NO2VCForecastError(TM5Data%im,TM5Data%jm) )
    NO2VCForecastError(:,:) = ForecastError

    ! ---------- Allocate storage and read the track matching this model time ----------

    call ObsFileGet( filename_S5P_L2_NO2_filelist, itau-ntimestep/2,itau+ntimestep/2, no2Tr, obsFound )

    ! At this point a track of NO2 and cloud data, "no2Tr" has 
    ! been read or the track data record "no2Tr" is empty: obsFound = .false.

    if ( obsFound ) then

       write(*,'(a,i6)') 'assimNO2: No. of pixels in track ',no2Tr%count
       nObs = no2Tr%count

       ! allocate 3D partial column field and 2D column fields
       allocate ( no2pcf(TM5Data%im,TM5Data%jm,TM5Data%lm) )
       allocate ( no2col_fc(TM5Data%im,TM5Data%jm) )
       allocate ( no2col_an(TM5Data%im,TM5Data%jm) )
       
       ! Store analysis/forecast ratio, = 1 before the analysis
       noxscaling(:,:,:) = 1.0
       
       ! get 3D field of partial NO2 columns (forecast)
       call GetNO2PartialColumnField ( TM5Data, noxscaling, no2pcf, no2col_fc )
    
       ! allocate storage for retrieval output
       call ObsFcInfoAllocate ( obsFcInfo, no2Tr%count, TM5Data%lm )

       ! Calculate model forecast and observation operator 
       ! store information for individual observations in "ObsFcInfo"
       isForecast = .true.
       call GetOMFHx( amfLut, no2Tr, TM5Data, no2pcf, obsFcInfo, isForecast, outputdir, debugMode )
       
       ! Write the retrieval data to the level-2 file, based on TM5 forecast 
       call TM5IF_obsfcinfo_write( no2Tr, obsFcInfo, TM5Data )
       
       ! cluster observations into superobservations
       call GetOMFHxSuper( TM5Data%im, TM5Data%jm, TM5Data%lm, no2Tr, obsFcInfo, obsFcInfoReduced, 'f', computeRetrievalForAnalysis )

       write(*,1244) obsFcInfo%count,obsFcInfoReduced%count,TM5Data%date
1244   format(' assimno2: ',i5,' (',i4,') observ.,  Timestep ',i4,5i3)

       ! ----- Assimilation step -----

       assimcount = assimcount + 1
          
       ! allocate the helper array
       allocate ( anFcRatio(TM5Data%im,TM5Data%jm,TM5Data%lm) )
       
       print *, 'assimno2: assimilating observations'
       call AssimilateColumn ( TM5Data%im, TM5Data%jm, TM5Data%lm, NO2VCForecastError, obsFcInfoReduced, no2pcf, anFcRatio )
       
       ! Store analysis/forecast ratio
       noxscaling(:,:,:) = anFcRatio(:,:,:)
          
       ! remove helpers
       deallocate ( anFcRatio )

       ! ----- Calculate model analysis and observation operator -----

       if ( computeRetrievalForAnalysis ) then	  

          ! allocate storage for retrieval output
          call ObsFcInfoAllocate ( obsAnInfo, no2Tr%count, TM5Data%lm )
       
          ! get 3D field of partial NO2 columns: for analysis, noxscaling <> 1.0
          call GetNO2PartialColumnField ( TM5Data, noxscaling, no2pcf, no2col_an )

          ! store information for individual observations in "ObsFcInfo"
          isForecast = .false.
          call GetOMFHx ( amfLut, no2Tr, TM5Data, no2pcf, obsAnInfo, isForecast, outputdir, debugMode )

          ! cluster observations into superobservations         
          !   note: this routine allocates the storage for the superobs
          call GetOMFHxSuper ( TM5Data%im, TM5Data%jm, TM5Data%lm, no2Tr, obsAnInfo, obsAnInfoReduced, 'a', computeRetrievalForAnalysis )

       end if
          
       ! ----- Write innovation statistics to netcdf file -----
       if ( doWriteOmF ) then
          if ( computeRetrievalForAnalysis ) then
             call OmF_Output( no2Tr, obsFcInfo, obsFcInfoReduced, NO2VCForecastError, TM5Data%im, TM5Data%jm, TM5Data%lon, TM5Data%lat, no2col_fc, TM5Data%date, outputdir, status, obsAnInfo, obsAnInfoReduced, no2col_an )
          else          
             call OmF_Output( no2Tr, obsFcInfo, obsFcInfoReduced, NO2VCForecastError, TM5Data%im, TM5Data%jm, TM5Data%lon, TM5Data%lat, no2col_fc, TM5Data%date, outputdir, status )
          end if
       end if
       
       ! empty the orbit stack
       no2Tr%norbitparts = 0
       no2Tr%count       = 0

       ! remove 2D (im,jm) forecast error field
       deallocate ( NO2VCForecastError )

       ! remove storage TROPOMI L2 NO2 observations
       call ObsTrackDeallocate ( no2Tr )

       ! remove storage retrieval results
       call ObsFcInfoDeallocate ( obsFcInfo )
       if ( computeRetrievalForAnalysis ) then	  
          call ObsFcInfoDeallocate ( obsAnInfo )
       end if

       ! remove storage superobservations
       call ObsInfoReduced_deallocate ( obsFcInfoReduced )
       if ( computeRetrievalForAnalysis ) then	  
          call ObsInfoReduced_deallocate ( obsAnInfoReduced )
       end if

       ! storage for partial column field and column fields
       deallocate ( no2pcf )
       deallocate ( no2col_fc )
       deallocate ( no2col_an )
       
    else  ! no observations found

       ! write(*,'(a,t40)') '  assimno2:',TM5Data%date
       nObs = 0

    end if
    
    write(*,1112) 
1112 format(' ------- assimno2 done  -------')
    
  end subroutine AssimNO2
  




  subroutine GetNO2PartialColumnField( TM5Data, analysisScaling, no2pcf, no2col )
    !
    ! extract a 3D partial column field from the TM model field "no2"
    ! input  :  "TM5Data" structure with TM5 fields 
    !           "analysisScaling" (=1 for forecast, <>1 for analysis),
    !               scaling factors resulting from the analysis 
    ! output :  "no2pcf" in 10^15 molecules per cm^2 for each 3D grid box
    !           "no2col" total column in 10^15 molecules per cm^2 
    !
    ! Adjusted in DOMINO-3 because the interface passes mixing ratio's 
    ! while the original code is based on the mass and "rm" arrays
    !                                                     Henk Eskes, 2015 

    ! TM5 model arrays and variables:
    use MTM5Data,               only : TTM5Data 
    use physconstants,          only : pdiff2moleccm2

    implicit none

    ! -- in/out --
    type(TTM5Data), intent(in)                                   ::  TM5Data 
    real,dimension(TM5Data%im,TM5Data%jm,TM5Data%lm),intent(in)  ::  analysisScaling
    real,dimension(TM5Data%im,TM5Data%jm,TM5Data%lm),intent(out) ::  no2pcf
    real,dimension(TM5Data%im,TM5Data%jm),intent(out)            ::  no2col

    ! -- local --
    integer :: i, j, l
    real    :: pbot, ptop
    real,dimension(:), allocatable  :: no2col_zonal
    real    :: no2col_mean

    !
    ! Convert from volume mixing ratio to molecules cm^-2 
    ! using the pressure drop over a layer (in Pa)
    !
    allocate ( no2col_zonal(TM5Data%jm) )
    no2col_zonal(:) = 0.0
    no2col(:,:) = 0.0

    do i = 1, TM5Data%im
       do j = 1, TM5Data%jm
          do l = 1, TM5Data%lm
             pbot = TM5Data%hyai(l)   + TM5Data%hybi(l) * TM5Data%ps(i,j) 
             ptop = TM5Data%hyai(l+1) + TM5Data%hybi(l+1) * TM5Data%ps(i,j) 
             no2pcf(i,j,l) = pdiff2moleccm2 * ( pbot - ptop ) * TM5Data%no2(i,j,l) * analysisScaling(i,j,l)
             no2col_zonal(j) = no2col_zonal(j) + no2pcf(i,j,l)
             no2col(i,j) = no2col(i,j) + no2pcf(i,j,l)
             ! if ( i==100 .and. j==100 ) print*,'l, pbot, ptop, no2, scaling ',l,pbot, ptop, TM5Data%no2(i,j,l), analysisScaling(i,j,l)
          end do
       end do
    end do
    no2col_zonal(:) = no2col_zonal(:) / real(TM5Data%im)
    no2col_mean     = sum(no2col_zonal) / real(TM5Data%jm)

    if ( isnan(no2col_mean) ) then 
       print '(a)', "ERROR in GetNO2PartialColumnField: NaN found in the TM5 NO2 field"
       stop "stop because of NaN"
    end if

!    print *, 'no2col_zonal -90, -60, .. 90 = ',no2col_zonal(1),no2col_zonal(30),no2col_zonal(60),no2col_zonal(90),no2col_zonal(120),no2col_zonal(150),no2col_zonal(180)
!    print *, 'no2col at (59,7) and (34,38) ',sum(no2pcf(59,7,:)),sum(no2pcf(34,38,:))
!    print *, 'no2pcf profile at 5E, 53N ',no2pcf(185,143,:)

    print '(a,f10.3)', ' GetNO2PartialColumnField, mean NO2 column = ', no2col_mean

    deallocate ( no2col_zonal )

  end subroutine GetNO2PartialColumnField



end module MAssimilation