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

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

! macro defs
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if

!=====================================================================
! 
! Central interface to link the assimilation code to the TM5 model
!
! The retrieval-assimilation will be performed on the root node
! The steps performed:
!       - Gather all arrays needed by the assimilation
!       - call AssimNO2( ntimestep )
!       - Scatter the arrays for the next TM5 time step
!
! The interface is designed such that the same fields are passed as in the 
!   NRT interface, with TM5 fields read from file
!
! The subroutine AssimNO2 performs a data assimilation time step.
!
!                                        Henk Eskes, KNMI, March 2015
!                                        Interface to TM5-mp-CB05
!======================================================================

module assim_interf_mod

  ! Definition of structure to store TM5 data
  use MTM5Data,         only : TTM5Data

  implicit none

  private

  public ::  AssimNO2_interface, AssimNO2_interface_init, AssimNO2_interface_done
  ! public ::  TM5Data

  ! This structure stores all info passed from TM5 to the assimilation 
  type(TTM5Data)  :: TM5Data

  ! scaling factor, the result of the assimilation 
  real, dimension(:,:,:), allocatable  :: noxscaling

  ! This array stores the 3D NOx mass field 
  real, dimension(:,:,:), allocatable :: store_rm_nox
  ! and these arrays the slopes
#ifdef slopes
  real, dimension(:,:,:), allocatable :: store_rxm_nox
  real, dimension(:,:,:), allocatable :: store_rym_nox
  real, dimension(:,:,:), allocatable :: store_rzm_nox  
#endif
  ! Arrays used to convert the surface pressure and orography from 3D to 2D
  real, dimension(:,:,:), allocatable :: sp_temp
  real, dimension(:,:,:), allocatable :: oro_temp
  real, dimension(:,:,:), allocatable :: gph_temp    ! to store geopotential height

  character(len=*), parameter  ::  mname = 'assim_interf_mod'

  ! ODIN HNO3 / O3 climatology related
  logical             ::  doNudgeOdinHNO3 
  character(len=255)  ::  OdinHNO3File
  ! This array stores the 3D HNO3 mass field (to be updated by the HNO3/O3 nudging) 
  real, dimension(:,:,:), allocatable :: store_rm_hno3

  ! Used for debug printing
  logical, parameter  :: showDebugPrints = .false.

contains



  subroutine AssimNO2_interface_init ( rcFile )
    !======================================================================
    !
    !   Initialise storage space for TM5 3D global arrays
    !
    !      rcFile     : name of the rc (settings) file
    !                                          Henk Eskes, KNMI, March 2015
    !======================================================================

    ! For printing information 
    use GO,               only : gol, goPr, goErr
    ! For reading keys from the .rc file
    use GO,               only : TrcFile, Init, Done, ReadRc
    ! To identify root proces
    use partools,         only : isRoot 
    ! for dimensions
    use meteodata,        only : global_lli, levi  
    ! to read ODIN HNO3/O3 climatology
    use MOdinHNO3,        only : ReadOdinHNO3

    implicit none

    ! --- in/out
    character(*), intent(in)     ::  rcFile

    ! --- local
    character(len=*), parameter  ::  rname = mname//'/AssimNO2_interface_init'

    type(TrcFile)    ::  rcF

    integer :: imr, jmr, lmr, n
    integer :: status

    ! --- start code ---

    if ( isRoot ) then 
       write(gol,'("Assim_interface: initialise interface")')  ; call goPr
    end if

    ! entire region grid size
    n = 1    ! global region
    imr = global_lli(n)%nlon
    jmr = global_lli(n)%nlat
    lmr = levi%nlev

    ! Fill fixed-dim constants
    TM5Data%im = imr
    TM5Data%jm = jmr
    TM5Data%lm = lmr

    ! Allocate arrays on root, dummy arrays elsewhere
    if ( isRoot ) then

       allocate ( TM5Data%lon(imr) )
       allocate ( TM5Data%lat(jmr) )
       allocate ( TM5Data%hyai(lmr+1) )
       allocate ( TM5Data%hybi(lmr+1) )
       allocate ( TM5Data%hyam(lmr) )
       allocate ( TM5Data%hybm(lmr) )
       allocate ( TM5Data%ps(imr,jmr) )
       allocate ( TM5Data%oro(imr,jmr) )
       allocate ( TM5Data%t(imr,jmr,lmr) )
       allocate ( TM5Data%ltropo(imr,jmr) )
       allocate ( TM5Data%no2(imr,jmr,lmr) )

       allocate ( noxscaling(imr,jmr,lmr) )

       allocate ( store_rm_nox(imr,jmr,lmr) )
#ifdef slopes
       allocate ( store_rxm_nox(imr,jmr,lmr) )
       allocate ( store_rym_nox(imr,jmr,lmr) )
       allocate ( store_rzm_nox(imr,jmr,lmr) )
#endif
       allocate ( sp_temp(imr,jmr,1) )
       allocate ( oro_temp(imr,jmr,1) )
       allocate ( gph_temp(imr,jmr,lmr+1) )

    else

       allocate ( TM5Data%lon(1) )
       allocate ( TM5Data%lat(1) )
       allocate ( TM5Data%hyai(1) )
       allocate ( TM5Data%hybi(1) )
       allocate ( TM5Data%hyam(1) )
       allocate ( TM5Data%hybm(1) )
       allocate ( TM5Data%ps(1,1) )
       allocate ( TM5Data%oro(1,1) )
       allocate ( TM5Data%t(1,1,1) )
       allocate ( TM5Data%ltropo(1,1) )
       allocate ( TM5Data%no2(1,1,1) )

       allocate ( noxscaling(1,1,1) )

       allocate ( store_rm_nox(1,1,1) )
#ifdef slopes
       allocate ( store_rxm_nox(1,1,1) )
       allocate ( store_rym_nox(1,1,1) )
       allocate ( store_rzm_nox(1,1,1) )
#endif
       allocate ( sp_temp(1,1,1) )
       allocate ( oro_temp(1,1,1) )
       allocate ( gph_temp(1,1,1) )

    end if

    if ( isRoot ) then 

       ! open rcfile:
       call Init( rcF, rcfile, status )
       IF_NOTOK_RETURN(status=1)
       
       ! Apply the alternative ODIN HNO3/O3 climatology ?
       call ReadRc( rcF, 'domino.doNudgeOdinHNO3', doNudgeOdinHNO3, status, default=.false. )
       IF_ERROR_RETURN(status=1)
       
       if ( doNudgeOdinHNO3 ) then
          ! Read the .nc file name (path/filename) which contains the ODIN HNO3/O3 climatology
          call ReadRc( rcF, 'domino.odinclimatologyfile', OdinHNO3File, status, default='list' )
          IF_ERROR_RETURN(status=1)
       end if

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

    end if

    ! Read the ODIN climatology of HNO3 and O3
    if ( doNudgeOdinHNO3 ) then

       if ( isRoot ) then 

          write(gol,'("AssimNO2_interface: initialise ODIN HNO3 climatology")')  ; call goPr
          
          ! read the ODIN climatology
          call ReadOdinHNO3 ( OdinHNO3File, jmr, lmr, status )
          IF_NOTOK_RETURN(status=1)

          ! Allocate the memory for the hno3 mass to be adjusted
          allocate ( store_rm_hno3(imr,jmr,lmr) )

       else  ! processor is not "root"

          ! Allocate the memory for the hno3 mass to be adjusted
          allocate ( store_rm_hno3(1,1,1) )

       end if

    end if

  end subroutine AssimNO2_interface_init



  subroutine AssimNO2_interface_done ( )
    !======================================================================
    !
    !   Release storage space for TM5 3D global arrays
    !                                          Henk Eskes, KNMI, March 2015
    !======================================================================

    ! For printing information 
    use GO,               only : gol, goPr, goErr
    ! To identify root proces
    use partools,         only : isRoot 
    ! Odin-HNO3 climatology
    use MOdinHNO3,        only : DoneOdinHNO3

    implicit none

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

    ! --- start code ---

    if ( isRoot ) then 
       write(gol,'("Assim_interface: deallocate arrays")') ; call goPr
    end if

    if ( allocated(TM5Data%lon)       )  deallocate ( TM5Data%lon )
    if ( allocated(TM5Data%lat)       )  deallocate ( TM5Data%lat )
    if ( allocated(TM5Data%hyai)      )  deallocate ( TM5Data%hyai )
    if ( allocated(TM5Data%hybi)      )  deallocate ( TM5Data%hybi )
    if ( allocated(TM5Data%hyam)      )  deallocate ( TM5Data%hyam )
    if ( allocated(TM5Data%hybm)      )  deallocate ( TM5Data%hybm )
    if ( allocated(TM5Data%ps)        )  deallocate ( TM5Data%ps )
    if ( allocated(TM5Data%oro)       )  deallocate ( TM5Data%oro )
    if ( allocated(TM5Data%t)         )  deallocate ( TM5Data%t )
    if ( allocated(TM5Data%ltropo)    )  deallocate ( TM5Data%ltropo )
    if ( allocated(TM5Data%no2)       )  deallocate ( TM5Data%no2 )

    if ( allocated(noxscaling)        )  deallocate ( noxscaling )

    if ( allocated(store_rm_nox)      )  deallocate ( store_rm_nox )
#ifdef slopes
    if ( allocated(store_rxm_nox)     )  deallocate ( store_rxm_nox )
    if ( allocated(store_rym_nox)     )  deallocate ( store_rym_nox )
    if ( allocated(store_rzm_nox)     )  deallocate ( store_rzm_nox )
#endif
    if ( allocated(sp_temp)           )  deallocate ( sp_temp )
    if ( allocated(oro_temp)          )  deallocate ( oro_temp )
    if ( allocated(gph_temp)          )  deallocate ( gph_temp )
       

    ! ODIN climatology of HNO3 and O3
    if ( doNudgeOdinHNO3 ) then

       ! Remove ODIN-HNO3 storage
       call DoneOdinHNO3 ( )

       ! De-allocate storage space of the gathered 3D HNO3 field
       if ( allocated(store_rm_hno3) )  deallocate ( store_rm_hno3 )

    end if

  end subroutine AssimNO2_interface_done



  subroutine AssimNO2_interface( rcFile, ntimestep, status )
    !======================================================================
    !
    !   Gather - Assimilate - Scatter
    !                                          Henk Eskes, KNMI, March 2015
    !======================================================================

    ! Observation vector dimensions:
    use tm5_distgrid,     only : dgrid, Get_DistGrid, Gather, Scatter
    ! NOx tracer indices
    use chem_param,       only : fscale, ino, ino2, ino3, in2o5, ihno4, inox
    use chem_param,       only : ihno3, io3
    ! Tracer concentrations: transported and chemistry-only
    use tracer_data,      only : mass_dat, chem_dat
    ! Meteo fields
    use meteodata,        only : temper_dat, oro_dat, sp_dat, m_dat, gph_dat 
    ! grid/level info
    use meteodata,        only : global_lli, levi
    ! contains "dxyp"
    use global_data,      only : region_dat
    ! For printing information 
    use GO,               only : gol, goPr, goErr
    ! To identify root proces
    use partools,         only : isRoot 
    ! For computing "dxyp"
    use binas,            only : ae, grav   ! earth radius (m), gravitation constant
    use dims,             only : dx, dy, gtor, idate
    ! Routine to perform an NO2 retrieval + assimilation step
    use MAssimilation,    only : assimNO2
    ! to nudge to ODIN HNO3/O3 climatology
    use MOdinHNO3,        only : NudgeOdinHNO3
    
    implicit none

    real, dimension(:,:,:), allocatable  :: tile_no2_field

    ! For ODIN climatology nudging of HNO3
    real, dimension(:,:,:), allocatable  :: tile_o3_field, tile_hno3_field
    real, dimension(:,:,:), allocatable  :: mixratio_3d_hno3, mixratio_3d_o3
    real, dimension(:,:), allocatable    :: mixratio_zonal_hno3, mixratio_zonal_o3

    character(len=*), intent(in)      ::  rcfile                 ! name of rc file
    integer, intent(out)    ::  status
    integer, intent(in)     ::  ntimestep

    ! --- local varables ---

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

    integer :: imr, jmr, lmr
    integer :: i1, i2, j1, j2, i, j, l, n
    integer :: nObs
    real    :: mix_no, mix_no2, mix_no3, mix_n2o5, mix_hno4, mix_nox_transport
    real    :: mix_hno3, mix_o3
    real    :: noxratio, chemnox
    real    :: dxx, dyy, lat 

    ! --- start code ---
       
    ! entire region grid size
    n = 1    ! global region
    imr = global_lli(n)%nlon
    jmr = global_lli(n)%nlat
    lmr = levi%nlev

    if ( isRoot ) then 

       ! Fill grid and time properties that are not distributed
       TM5Data%date(:) = idate(:)

       TM5Data%im = imr
       TM5Data%jm = jmr
       TM5Data%lm = lmr
       
       TM5Data%lon(:) = global_lli(n)%lon_deg(:)
       TM5Data%lat(:) = global_lli(n)%lat_deg(:)
       
       TM5Data%hyai(1:lmr+1) = levi%a(0:lmr)
       TM5Data%hybi(1:lmr+1) = levi%b(0:lmr)
       TM5Data%hyam(1:lmr) = levi%fa(1:lmr)
       TM5Data%hybm(1:lmr) = levi%fb(1:lmr)

       if ( showDebugPrints ) then
          print*,'date ',TM5Data%date
          print*,'im, jm, lm ',TM5Data%im,TM5Data%jm,TM5Data%lm
          print*,'lon range ',TM5Data%lon(1),TM5Data%lon(imr)
          print*,'lat range ',TM5Data%lat(1),TM5Data%lat(jmr)
          print*,'hyai range ',TM5Data%hyai(1),TM5Data%hyai(lmr+1)
          print*,'hybi range ',TM5Data%hybi(1),TM5Data%hybi(lmr+1)
          print*,'hyam range ',TM5Data%hyam(1),TM5Data%hyam(lmr)
          print*,'hybm range ',TM5Data%hybm(1),TM5Data%hybm(lmr)
       end if

    end if

    ! *** Gather all arrays needed ***

    if ( isRoot ) then 
       write(gol,'("Assim_interface: Gather NO2, NOx and meteo fields")') ; call goPr
    end if

    ! First, determine updated NO2 after NOx transport (distributed on all arrays)
    !   All NOx species are scaled with the same factor

    call Get_DistGrid( dgrid(n), I_STRT=i1,  I_STOP=i2, J_STRT=j1,  J_STOP=j2 ) 

    allocate( tile_no2_field(i1:i2, j1:j2, lmr) )

    do l = 1, lmr
       do j = j1, j2 
          do i = i1, i2
             ! Get NOx after chemistry from sum of non-transported components
             ! "fscale(ispec) = xmair/ra(ispec)"
             mix_no   = chem_dat(n)%rm(i,j,l,ino) * fscale(ino) / m_dat(n)%data(i,j,l)
             mix_no2  = chem_dat(n)%rm(i,j,l,ino2) * fscale(ino2) / m_dat(n)%data(i,j,l)
             mix_no3  = chem_dat(n)%rm(i,j,l,ino3) * fscale(ino3) / m_dat(n)%data(i,j,l)
             mix_n2o5 = chem_dat(n)%rm(i,j,l,in2o5) * fscale(in2o5) / m_dat(n)%data(i,j,l)
             mix_hno4 = chem_dat(n)%rm(i,j,l,ihno4) * fscale(ihno4) / m_dat(n)%data(i,j,l)
             ! NOx after chemistry, rescaled with molecule mass
             chemnox = mix_no + mix_no2 + mix_no3 + 2.0*mix_n2o5 + mix_hno4
             ! NOx after time stepping, rescaled with molecule mass
             mix_nox_transport = mass_dat(n)%rm(i,j,l,inox) * fscale(inox) / m_dat(n)%data(i,j,l)
             ! Get ratio of transported and non-transported NOx
             if ( chemnox > 1.0e-15 ) then
                noxratio = mix_nox_transport/chemnox
             else
                noxratio = 1.0
             end if
             ! Re-scale NO2 mixing ratio field to concentration after transport
             tile_no2_field(i,j,l) = noxratio*mix_no2 
             if ( showDebugPrints ) then
                if (l==1 .and. j==j1 .and. i==295 ) &
                     print '(a,3i4,f8.4,e13.4)',' i,j,l,nox_ratio,nox = ',i,j,l,noxratio,mix_nox_transport
             end if
          end do
       end do
    end do

    if ( showDebugPrints ) then
       if ( isRoot ) then
          print '(a,4i4)','i1, i2, j1, j2 ',i1, i2, j1, j2
          print '(a,2e10.3)','fscale no, no2 ',fscale(ino),fscale(ino2)
          print '(a,34e10.2)', 'Tile no2 (i2,j2,:) ', tile_no2_field(i2,j2,:)
       end if
    end if

    ! Store NO2 field tiles into global NO2 array (volume mixing ratio)
    call Gather( dgrid(n), tile_no2_field, TM5Data%no2, 0, status)
    IF_NOTOK_RETURN(status=1)

    deallocate( tile_no2_field )


    ! *** Nudging to ODIN climatology, gather the zonal means of HNO3 and O3 ***

    if ( doNudgeOdinHNO3 ) then

       write(gol,'("Assim_interface: compute ODIN HNO3 nudging factors ")') ; call goPr

       n = 1
       call Get_DistGrid( dgrid(n), I_STRT=i1,  I_STOP=i2, J_STRT=j1,  J_STOP=j2 ) 

       ! allocate the tiles
       allocate( tile_hno3_field(i1:i2, j1:j2, lmr) )
       allocate( tile_o3_field(i1:i2, j1:j2, lmr) )

       do l = 1, lmr
          do j = j1, j2 
             do i = i1, i2
                ! Get O3 and HNO3 mixing ratios
                ! "fscale(ispec) = xmair/ra(ispec)"
                tile_hno3_field(i,j,l) = mass_dat(n)%rm(i,j,l,ihno3) * fscale(ihno3) / m_dat(n)%data(i,j,l)
                tile_o3_field(i,j,l)   = mass_dat(n)%rm(i,j,l,io3) * fscale(io3) / m_dat(n)%data(i,j,l)
             end do
          end do
       end do

       if ( showDebugPrints ) then
          if ( isRoot ) then
             print '(a,4i4)','i1, i2, j1, j2 ',i1, i2, j1, j2
             print '(a,2e10.3)','fscale hno3, o3 ',fscale(ihno3),fscale(io3)
             print '(a,34e10.2)', 'Tile hno3 (i2,j2,:) ', tile_hno3_field(i2,j2,:)
             print '(a,34e10.2)', 'Tile o3 (i2,j2,:)   ', tile_o3_field(i2,j2,:)
          end if
       end if

       ! Allocate the memory for the hno3 / o3 3D mixing ratio fields on root
       if ( isRoot ) then
          allocate ( mixratio_3d_hno3(imr,jmr,lmr) )
          allocate ( mixratio_3d_o3(imr,jmr,lmr) )
       else  
          allocate ( mixratio_3d_hno3(1,1,1) )
          allocate ( mixratio_3d_o3(1,1,1) )
       end if
       
       ! Store O3, HNO3 field tiles into global arrays (volume mixing ratio)
       call Gather( dgrid(n), tile_hno3_field, mixratio_3d_hno3, 0, status)
       IF_NOTOK_RETURN(status=1)
       call Gather( dgrid(n), tile_o3_field, mixratio_3d_o3, 0, status)
       IF_NOTOK_RETURN(status=1)
       
       ! remove the tiles
       deallocate( tile_hno3_field )
       deallocate( tile_o3_field )

       ! Store HNO3 mass field tiles into global HNO3 array (kg)
       call Gather( dgrid(n), mass_dat(n)%rm(:,:,:,ihno3), store_rm_hno3, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if ( isRoot ) then

          ! Allocate zonal mean arrays
          allocate ( mixratio_zonal_hno3(jmr,lmr) )
          allocate ( mixratio_zonal_o3(jmr,lmr) )

          ! Determine 2D zonal mean mixing ratio fields of HNO3 and O3
          do l = 1, lmr
             do j = 1, jmr
                mixratio_zonal_hno3(j,l) = 0.0
                mixratio_zonal_o3(j,l) = 0.0
                do i = 1, imr
                   mixratio_zonal_hno3(j,l) = mixratio_zonal_hno3(j,l) + mixratio_3d_hno3(i,j,l)
                   mixratio_zonal_o3(j,l) = mixratio_zonal_o3(j,l) + mixratio_3d_o3(i,j,l)
                end do
                mixratio_zonal_hno3(j,l) = mixratio_zonal_hno3(j,l) /imr
                mixratio_zonal_o3(j,l) = mixratio_zonal_o3(j,l) / imr
             end do
          end do

          if ( showDebugPrints ) then
             print '(a,3i4,i5,5i3,i6)', 'assim_interf: imr, jmr, lmr, idate, timestep = ', imr, jmr, lmr, idate, ntimestep
             print '(a,34e10.2)', 'assim_interf: mixratio_zonal_hno3 j 120 = ',mixratio_zonal_hno3(120,:)
             print '(a,34e10.2)', 'assim_interf: mixratio_zonal_o3   j 120 = ',mixratio_zonal_o3(120,:)
             print '(a,34e10.2)', 'assim_interf: store_rm_hno3  i 1 j 120 = ',store_rm_hno3(1,120,:)
             print '(a)', 'assim_interf: calling NudgeOdinHNO3'
          end if

          ! Determine scaled HNO3 field
          call NudgeOdinHNO3 ( imr, jmr, lmr, idate, ntimestep, mixratio_zonal_hno3, mixratio_zonal_o3, store_rm_hno3 )

          ! Remove zonal mean arrays
          deallocate ( mixratio_zonal_hno3 )
          deallocate ( mixratio_zonal_o3 )

       end if   ! isRoot

       ! release storage of the 3D mixing ratio fields
       deallocate ( mixratio_3d_hno3 )
       deallocate ( mixratio_3d_o3 )

    end if     ! doNudgeOdinHNO3


    ! *** Gather other TM5 fields needed ***

    ! Store NOx mass field tiles into global NOx array (kg)
    call Gather( dgrid(n), mass_dat(n)%rm(:,:,:,inox), store_rm_nox, mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)

    ! Same for the slopes
#ifdef slopes 
    call Gather( dgrid(n), mass_dat(n)%rxm(:,:,:,inox), store_rxm_nox, mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)
    call Gather( dgrid(n), mass_dat(n)%rym(:,:,:,inox), store_rym_nox, mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)
    call Gather( dgrid(n), mass_dat(n)%rzm(:,:,:,inox), store_rzm_nox, mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)
#endif

    ! Surface pressure (Pa): sp(imr,jmr,1) 
    call Gather( dgrid(n), sp_dat(n)%data, sp_temp, sp_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)

    ! Orography (m2/s2): oro(imr,jmr,1), divide by "g" for meters 
    call Gather( dgrid(n), oro_dat(n)%data, oro_temp, oro_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)

    ! Temperature field (K): temper(imr,jmr,lmr) 
    call Gather( dgrid(n), temper_dat(n)%data, TM5Data%t, temper_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)

    ! Geopotential height (m): gph(imr,jmr,lmr+1) 
    call Gather( dgrid(n), gph_dat(n)%data, gph_temp, gph_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)

    if ( isRoot ) then

       ! Store surface pressure 
       TM5Data%ps(:,:) = sp_temp(:,:,1)

       ! Store model orography
       TM5Data%oro(:,:) = oro_temp(:,:,1)/grav

       ! Determine tropopause level and store
       do i = 1, imr
          do j = 1, jmr
             TM5Data%ltropo(i,j) = ltropo(lmr,TM5Data%hyai,TM5Data%hybi,TM5Data%t(i,j,:),gph_temp(i,j,:))
          end do
       end do

       ! Initialise assimilation scaling factor
       noxscaling(:,:,:) = 1.0

       ! *** Call the assimilation ***

       ! Assimilation is done only on the root processor

       if ( showDebugPrints ) then
          i = 350
          j = 120
          print '(a,i4,a,i4,a)', "Gathered data for cell (",i,", ",j,":, 1) " 
          print '(a,10e9.2)', "  NO2  ", TM5Data%no2(i, j:j+9, 1)
          print '(a,10e9.2)', "  NOx  ", store_rm_nox(i, j:j+9, 1)
          print '(a,10e9.2)', "  NOx slope x  ", store_rxm_nox(i, j:j+9, 1)
          print '(a,10e9.2)', "  NOx slope y  ", store_rym_nox(i, j:j+9, 1)
          print '(a,10e9.2)', "  NOx slope z  ", store_rzm_nox(i, j:j+9, 1)
          print '(a,10f9.1)', "  temp ", TM5Data%t(i, j:j+9, 1)
          print '(a,10f9.1)', "  oro  ", TM5Data%oro(i, j:j+9)
          print '(a,180i3)' , "  ltropo_lat  ", TM5Data%ltropo(i, 1:180)
          print '(a,180i3)' , "  ltropo_lon  ", TM5Data%ltropo(1:360, j)
          print '(a,180f5.0)', "  ptropo ", TM5Data%hyai(TM5Data%ltropo(i, 1:180))*0.01+TM5Data%hybi(TM5Data%ltropo(i, 1:180))*1013.0
          print '(a,35f6.0)', "  plevs ", TM5Data%hyai(:)*0.01+TM5Data%hybi(:)*1013.0
          print '(a,10f9.0)', "  sp   ", TM5Data%ps(i, j:j+9)*0.01
       end if

       write(gol,'("Assim_interface: Fields have been gathered, now calling assimNO2")') ; call goPr

       call assimNO2( rcFile, ntimestep, TM5Data, nObs, noxscaling, status )
       IF_NOTOK_RETURN(status=1)

       if ( nObs > 0 ) then

          ! *** Finally, scatter the analysis NOx field ***
          
          write(gol,'("Assim_interface: Update NOx field, slopes, and scatter")') ; call goPr
          
          ! Apply assimilation scaling to NOx (only) ...
          store_rm_nox(:,:,:) = store_rm_nox(:,:,:) * noxscaling(:,:,:)
          ! Scale also slopes 
#ifdef slopes 
          store_rxm_nox(:,:,:) = store_rxm_nox(:,:,:) * noxscaling(:,:,:)
          store_rym_nox(:,:,:) = store_rym_nox(:,:,:) * noxscaling(:,:,:)
          store_rzm_nox(:,:,:) = store_rzm_nox(:,:,:) * noxscaling(:,:,:)
#endif
       end if   ! nObs > 0

    end if  ! isRoot

    ! ... and scatter NOx back to the MPI tiles
    call Scatter( dgrid(n), mass_dat(n)%rm(:,:,:,inox), store_rm_nox, &
         mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)

#ifdef slopes 
    call Scatter( dgrid(n), mass_dat(n)%rxm(:,:,:,inox), store_rxm_nox, &
         mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)
    call Scatter( dgrid(n), mass_dat(n)%rym(:,:,:,inox), store_rym_nox, &
         mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)
    call Scatter( dgrid(n), mass_dat(n)%rzm(:,:,:,inox), store_rzm_nox, &
         mass_dat(n)%halo, status)
    IF_NOTOK_RETURN(status=1)
#endif

    ! *** Scatter result of the HNO3 nudging to ODIN climatology ***

    if ( doNudgeOdinHNO3 ) then
          
       write(gol,'("Assim_interface: now scattering updated HNO3 concentration (ODIN nudging) ")') ; call goPr
       call Scatter( dgrid(n), mass_dat(n)%rm(:,:,:,ihno3), store_rm_hno3, &
            mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)
    end if

  end subroutine AssimNO2_interface


  !-----------------------------------------------------------------------
  !                    TM5                                                !
  !-----------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  ltropo
  !
  ! !DESCRIPTION: determine tropopause level from temperature gradient
  !	          reference: WMO /Bram Bregman
  !\\
  !\\
  ! !INTERFACE:
  !
  integer function ltropo(lmx,at,bt,tx,gphx)
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)                :: lmx
    real,dimension(lmx+1), intent(in) :: at, bt
    real,dimension(lmx),intent(in)    :: tx
    real,dimension(lmx+1),intent(in)  :: gphx 
    !
    ! !REVISION HISTORY: 
    !	programmed by		 fd 01-11-1996 
    !   changed to function      fd 12-07-2002
    !   16 Jul 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    integer :: ltmin,ltmax,l
    real    :: dz,dt

    !   ltropo is highest tropospheric level
    !   above is defined as stratosphere

    ! min tropopause level 450 hPa (ca. 8 km)
    ltmin=lvlpress(lmx,at,bt,45000.,98400.)
    ! max tropopause level 70 hPa (ca. 20 km)
    ltmax=lvlpress(lmx,at,bt,7000.,98400.)
    ltropo=ltmin
    do l=ltmin,ltmax
       dz=(gphx(l)-gphx(l-2))/2. 
       dt=(tx(l)-tx(l-1))/dz
       if ( dt < -0.002) ltropo=l !wmo 85 criterium for stratosphere
       ! increase upper tropospheric level
    end do !l

  end function ltropo
  !EOC
  
  !---------------------------------------------------------------------------
  !                    TM5                                                    !
  !---------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  lvlpress
  !
  ! !DESCRIPTION: lvlpress determines the index of the level with a pressure
  !               greater than press i.e. in height below press
  !               determines level independent of vertical resolution lm 
  !               uses hybrid level coefficients to determine lvlpress
  !\\
  !\\
  ! !INTERFACE:
  !
  integer function lvlpress(lm,at,bt,press,press0)
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: lm
    real, dimension(lm+1), intent(in) :: at, bt
    real,    intent(in) :: press
    real,    intent(in) :: press0
    !
    ! !REVISION HISTORY: 
    !   programmed by           Peter van Velthoven, 6-november-2000
    !   16 Jul 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real             :: zpress0, zpress
    integer          :: l,lvl

    if ( press0 == 0.0 ) then
       zpress0 = 98400.
    else
       zpress0 = press0
    endif
    lvl = 1
    ! l increases from bottom (l=1) to top (l=lm)
    do l = 1, lm
       zpress =  (at(l)+at(l+1) + (bt(l)+bt(l+1))*zpress0)*0.5
       if ( press < zpress ) then
          lvl = l
       endif
    end do
    lvlpress = lvl

  end function lvlpress
  !EOC
  

end module assim_interf_mod