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

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

! -----------------------------------------------------------------------------
!
! TROPOMI NL L2 Data Processors
!
! Copyright 2013-2015 KNMI
! This software is the proprietary information of KNMI
! All rights Reserved
!
! -----------------------------------------------------------------------------

module TM5IF_module

  use MTObsTrack
  use MTObsFcInfo
  use netcdf
  use pqf_module, only : PQF_E_GENERIC_EXCEPTION

  implicit none

  private

  public :: TM5IF_select, TM5IF_obstrack_read, TM5IF_obsfcinfo_write
  public :: TM5IF_time_init, TM5IF_time_to_string
  public :: TM5IF_log_info, TM5IF_log_error, TM5IF_log_debug
  public :: TM5IF_date, FILENAME_LEN

  integer, parameter          :: FILENAME_LEN = 256
  integer, parameter          :: STR_LEN = 256

  integer                     :: ncfileid
  integer                     :: ncgroupid
  character*(FILENAME_LEN)    :: ncfilename

  ! In case of QA4ECV, the unit of the columns is 10^15 molecules per cm^2
  logical                     :: use_scd_qa4ecv_units

  ! undefined
  ! real, parameter             :: nf_fill_float = 9.9692099683868690e+36
  real(4), parameter          :: nf_fill_float = 9.96921e+36

  type TM5IF_date
     integer :: year, month, day, hour, minute, second
  end type TM5IF_date

  interface TM5IF_put_var
     module procedure TM5IF_put_var_f2, TM5IF_put_var_f3, TM5IF_put_var_f4, TM5IF_put_var_i3
  end interface

  interface TM5IF_get_var
     module procedure TM5IF_get_var_f2, TM5IF_get_var_f3, TM5IF_get_var_f4, TM5IF_get_var_i1, TM5IF_get_var_i2, TM5IF_get_var_i3
  end interface

  interface TM5IF_put_att
     module procedure TM5IF_put_att_str
  end interface

  interface TM5IF_get_att
     module procedure TM5IF_get_att_i, TM5IF_get_att_str
  end interface

contains

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

  subroutine TM5IF_select ( filename, timeStart, timeStop, obstrack, before_or_after )

    !   Select L2 NO2 products.
    !   Check if product coverage time range falls in processing time range.
    !   If so, read number of scanlines used for allocation of memory for
    !   reading during pass 2.
    !
    !  before_or_after = -1 : the orbit time is before the time interval 
    !                     0 : the orbit time is inside the time interval 
    !                     1 : the orbit time is after the time interval 

    implicit none

    character*(*),intent(in)      :: filename
    type(TM5IF_date),intent(in)   :: timeStart
    type(TM5IF_date),intent(in)   :: timeStop
    type(TObsTrack),intent(inout) :: obstrack
    integer, intent(out)          :: before_or_after

    character*(STR_LEN) :: str

    integer :: i, j, n, status
    integer :: nValidPixels
    integer :: dimScanline
    integer :: dimGroundPixel
    integer :: firstimage
    integer :: secCoverageStart
    integer :: secCoverageEnd
    integer :: secCoverageMid

    type(TM5IF_date) :: timeCoverageStart
    type(TM5IF_date) :: timeCoverageEnd

    integer, allocatable :: processingQualityFlags(:,:,:)
    integer, allocatable :: processingErrorFlag(:,:,:)

    status = TM5IF_open ( filename, nf90_nowrite )
    if (status /= 0) goto 999

    status = TM5IF_get_att("time_coverage_start", str)
    if (status /= 0) goto 998

    status = TM5IF_convert_date(str, timeCoverageStart);
    if (status /= 0) goto 998

    status = TM5IF_get_att("time_coverage_end", str)
    if (status /= 0) goto 998

    status = TM5IF_convert_date(str, timeCoverageEnd);
    if (status /= 0) goto 998

    secCoverageStart = TM5IF_time_get_seconds(timeCoverageStart)
    secCoverageEnd   = TM5IF_time_get_seconds(timeCoverageEnd)
    secCoverageMid   = (secCoverageStart + secCoverageEnd) / 2

    if (secCoverageMid < TM5IF_time_get_seconds(timeStart)) then
       ! orbit before time interval
       before_or_after = -1 
       goto 998
    endif

    if (secCoverageMid > TM5IF_time_get_seconds(timeStop)) then
       ! orbit after time interval
       before_or_after = 1 
       goto 998
    endif

    ! orbit inside time interval
    before_or_after = 0 

    status = TM5IF_set_group("/PRODUCT")
    if (status /= 0) goto 998

    status = TM5IF_get_dimension("scanline", dimScanline)
    if (status /= 0) goto 998

    status = TM5IF_get_dimension("ground_pixel", dimGroundPixel)
    if (status /= 0) goto 998

    status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS")
    if (status /= 0) goto 998

    status = TM5IF_get_var("processing_quality_flags", processingQualityFlags)
    if (status /= 0) goto 998

    status = TM5IF_set_group("/PRODUCT")
    if (status /= 0) goto 998

    nValidPixels = 0

    do i = 1, dimScanline
       do j = 1, dimGroundPixel
          if (iand(processingQualityFlags(j, i, 1), 255) == 0) nValidPixels = nValidPixels + 1
       end do
    end do

    deallocate ( processingQualityFlags )

    print '(a,i6,a,i6)', ' TM5IF_select: nValidPixels = ', nValidPixels,'   total = ',dimGroundPixel*dimScanline

    n = obstrack%norbitparts + 1

    obstrack%orbitparts(n)%filename = filename
    obstrack%orbitparts(n)%nValidPixels = nValidPixels
    obstrack%orbitparts(n)%dimScanline = dimScanline

    obstrack%orbitparts(n)%starttime(1) = timeCoverageStart%year
    obstrack%orbitparts(n)%starttime(2) = timeCoverageStart%month
    obstrack%orbitparts(n)%starttime(3) = timeCoverageStart%day
    obstrack%orbitparts(n)%starttime(4) = timeCoverageStart%hour
    obstrack%orbitparts(n)%starttime(5) = timeCoverageStart%minute
    obstrack%orbitparts(n)%starttime(6) = timeCoverageStart%second

    obstrack%orbitparts(n)%endtime(1) = timeCoverageEnd%year
    obstrack%orbitparts(n)%endtime(2) = timeCoverageEnd%month
    obstrack%orbitparts(n)%endtime(3) = timeCoverageEnd%day
    obstrack%orbitparts(n)%endtime(4) = timeCoverageEnd%hour
    obstrack%orbitparts(n)%endtime(5) = timeCoverageEnd%minute
    obstrack%orbitparts(n)%endtime(6) = timeCoverageEnd%second

    !obstrack%orbitparts(n)%meantime   = TBD

    obstrack%norbitparts = n

998 continue

    status = TM5IF_close()
    if (status /= 0) goto 999

999 continue

  end subroutine TM5IF_select

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

  subroutine TM5IF_obstrack_read(obstrack)

    !   Read L2 NO2 product

    implicit none

    type(TObsTrack) :: obstrack

    character(STR_LEN) :: str
    character(STR_LEN) :: scd_unit_string
    integer :: dim_scanline
    integer :: dim_ground_pixel
    integer :: dim_corner
    integer :: dim_time
    integer :: dim_polynomial_exponents
    integer :: dim_pressure_levels
    integer :: dim_profile_layers
    integer :: i, j, k, n
    integer :: iorbit
    integer :: status

    ! must be real(4) because of reading from netcdf
    real(4), allocatable :: latitude(:,:,:)
    real(4), allocatable :: longitude(:,:,:)
    real(4), allocatable :: latitude_bounds(:,:,:,:)
    real(4), allocatable :: longitude_bounds(:,:,:,:)
    real(4), allocatable :: solar_zenith_angle(:,:,:)
    real(4), allocatable :: solar_azimuth_angle(:,:,:)
    real(4), allocatable :: viewing_zenith_angle(:,:,:)
    real(4), allocatable :: viewing_azimuth_angle(:,:,:)
    real(4), allocatable :: surface_altitude(:,:,:)
    real(4), allocatable :: surface_pressure(:,:,:)
    real(4), allocatable :: surface_albedo_no2(:,:,:)
    real(4), allocatable :: scd_no2(:,:,:)
    real(4), allocatable :: scd_no2_precision(:,:,:)
    real(4), allocatable :: chi_squared(:,:,:)
    real(4), allocatable :: cloud_fraction_no2(:,:,:)
    real(4), allocatable :: cloud_pressure(:,:,:)
    real(4), allocatable :: cloud_albedo(:,:,:)
    real(4), allocatable :: cloud_radiance_fraction_no2(:,:,:)
    real(4), allocatable :: scene_albedo(:,:,:)
    real(4), allocatable :: scene_pressure(:,:,:)
    integer, allocatable :: snow_ice_flag(:,:,:)
    integer, allocatable :: pixel_type(:,:,:)
    integer, allocatable :: time(:)
    integer, allocatable :: delta_time(:,:)
    integer, allocatable :: processing_quality_flags(:,:,:)

    obstrack%count = 0

    call TM5IF_log_info("TM5IF_obstrack_read: now reading the track (orbit parts)")

    n = 0
    do i = 1, obstrack%norbitparts
        n = n + obstrack%orbitparts(i)%nValidPixels
    end do

    call ObsTrackAllocate(obstrack, n)

    do i = 1, obstrack%norbitparts

       call TM5IF_log_info("  Reading data from " // trim(obstrack%orbitparts(i)%filename))

       status = TM5IF_open ( trim(obstrack%orbitparts(i)%filename), nf90_nowrite)
       if (status /= 0) goto 998

       status = TM5IF_get_att("orbit", iorbit)
       if (status /= 0) goto 998

       status = TM5IF_set_group("/PRODUCT")
       if (status /= 0) goto 996

       status = TM5IF_get_dimension("time", dim_time)
       if (status /= 0) goto 996

       status = TM5IF_get_dimension("scanline", dim_scanline)
       if (status /= 0) goto 996

       status = TM5IF_get_dimension("ground_pixel", dim_ground_pixel)
       if (status /= 0) goto 996

       if (dim_time /= 1) then
          call TM5IF_log_error("  Invalid value for dimension 'time'")
          goto 996
       end if

       status = TM5IF_get_var( "time", time )
       if (status /= 0) goto 996

       status = TM5IF_get_var( "delta_time", delta_time )
       if (status /= 0) goto 996

       status = TM5IF_get_var( "latitude", latitude )
       if (status /= 0) goto 996

       status = TM5IF_get_var( "longitude", longitude )
       if (status /= 0) goto 996

       status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/GEOLOCATIONS")
       if (status /= 0) goto 996

       status = TM5IF_get_var( "pixel_type", pixel_type )
       if (status /= 0) goto 996

       status = TM5IF_get_var( "latitude_bounds", latitude_bounds )
       if (status /= 0) goto 996

       status = TM5IF_get_var( "longitude_bounds", longitude_bounds )
       if (status /= 0) goto 996

       status = TM5IF_get_var( "solar_zenith_angle", solar_zenith_angle )
       if ( status /= 0 ) goto 996

       status = TM5IF_get_var( "viewing_zenith_angle", viewing_zenith_angle )
       if ( status /= 0 ) goto 996

       status = TM5IF_get_var( "solar_azimuth_angle", solar_azimuth_angle, printErrorMessage = .false. )
       if ( status /= 0 ) then
          ! try QA4ECV convention
          status = TM5IF_get_var( "relative_azimuth_angle", solar_azimuth_angle )
          if ( status /= 0 ) goto 996
          call TM5IF_log_info("  Note: using -relative_azimuth_angle- instead of -solar_azimuth_angle- ")
       end if

       status = TM5IF_get_var( "viewing_azimuth_angle", viewing_azimuth_angle, printErrorMessage = .false. )
       if ( status /= 0 ) then
          ! try QA4ECV convention
          ! set vaa = raa, so that the computation in "MObservationOperator" yields a correct value
          status = TM5IF_get_var( "relative_azimuth_angle", viewing_azimuth_angle )
          if ( status /= 0 ) goto 996
          ! re-set saa to 180, so that the computation in "MObservationOperator" yields a correct value
          solar_azimuth_angle(:,:,:) = 180.0
          call TM5IF_log_info("  Note: using -relative_azimuth_angle_surf- instead of -viewing_azimuth_angle- ")
       end if

       status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS")
       if (status /= 0) goto 996

       status = TM5IF_get_var("processing_quality_flags", processing_quality_flags)
       if (status /= 0) goto 996

       status = TM5IF_get_var("scd_no2", scd_no2)
       if (status /= 0) goto 996
       status = TM5IF_get_var_att_str("scd_no2", "units", scd_unit_string)
       if (status /= 0) goto 996
       use_scd_qa4ecv_units = .false.
       if ( trim(scd_unit_string) == "molecules cm-2" ) use_scd_qa4ecv_units = .true.
       print *, '  scd_no2: use qa4ecv units = ',use_scd_qa4ecv_units
       print *, '  scd_no2: units = ',trim(scd_unit_string)

       status = TM5IF_get_var("scd_no2_precision", scd_no2_precision)
       if (status /= 0) goto 996

       status = TM5IF_get_var( "chi_squared", chi_squared, printErrorMessage = .false. )
       if ( status /= 0 ) then
          ! QA4ECV does not provide chi_2
          ! set to fill value
          allocate ( chi_squared(dim_ground_pixel,dim_scanline,1) )
          chi_squared(:,:,:) = nf_fill_float
          call TM5IF_log_info("  Note: -chi_squared- set to fill value ")          
       end if

       status = TM5IF_get_var( "cloud_radiance_fraction_no2", cloud_radiance_fraction_no2, printErrorMessage = .false. )
       if ( status /= 0 ) then
          ! QA4ECV does not provide cloud radiance fraction: set to fill value in this case
          allocate ( cloud_radiance_fraction_no2(dim_ground_pixel,dim_scanline,1) )
          cloud_radiance_fraction_no2(:,:,:) = nf_fill_float
          call TM5IF_log_info("  Warning: -cloud_radiance_fraction_no2- set to fill value ")          
       end if

       status = TM5IF_get_var( "cloud_fraction_no2", cloud_fraction_no2, printErrorMessage = .false. )
       if ( status /= 0 ) then
          ! QA4ECV does not provide "cloud_fraction_no2", use cloud fraction instead
          status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/INPUT_DATA")
          if ( status /= 0 ) goto 996
          status = TM5IF_get_var( "cloud_fraction", cloud_fraction_no2 )
          if ( status /= 0 ) goto 996
          call TM5IF_log_info("  Note: using -cloud_fraction- instead of -cloud_fraction_no2- ")
       end if

       status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/INPUT_DATA")
       if (status /= 0) goto 996

       status = TM5IF_get_var("surface_altitude", surface_altitude)
       if (status /= 0) goto 996

       status = TM5IF_get_var("surface_pressure", surface_pressure)
       if (status /= 0) goto 996

       status = TM5IF_get_var("surface_albedo_no2", surface_albedo_no2)
       if (status /= 0) goto 996

       status = TM5IF_get_var("cloud_pressure", cloud_pressure)
       if (status /= 0) goto 996

       status = TM5IF_get_var("cloud_albedo", cloud_albedo, printErrorMessage = .false.)
       if ( status /= 0 ) then
          ! QA4ECV does not provide cloud albedo
          ! set to fill value
          allocate ( cloud_albedo(dim_ground_pixel,dim_scanline,1) )
          cloud_albedo(:,:,:) = 0.8
          call TM5IF_log_info("  Note: -cloud_albedo- set to 0.8 ")          
       end if

       status = TM5IF_get_var("snow_ice_flag", snow_ice_flag)
       if (status /= 0) goto 996

       status = TM5IF_get_var("scene_pressure", scene_pressure)
       if (status /= 0) goto 996

       status = TM5IF_get_var("scene_albedo", scene_albedo)
       if (status /= 0) goto 996

       ! note: surface_altitude_precision is a field in the file, but is not defined
       ! note: surface_classification is a field in the file, but is also not defined
       ! note: for OMI, "omi_xtrack_flags" row anomaly flags are provided. 
       !       However, these errors are included in the "processing_quality_flags"
       ! note: the "snow_ice_flag" needs to be read, because this info is not yet included
       !       in the processing_quality_flags. In the "MObservationoperator" module 
       !       the action needs to be defined for snow/ice (use scene albedo/pressure), 
       !       and the processing_quality_flags should be updated accordingly. 

       n = obstrack%count

       print *, '  dims scanline groundpix = ', dim_scanline, dim_ground_pixel

       obstrack%dimGroundPixel = dim_ground_pixel

       ! print *, 'processing_quality_flags(:, j=500, 1)',processing_quality_flags(:, 800, 1)

       obstrack%orbitNumber = iorbit

       do j = 1, dim_scanline

          do k = 1, dim_ground_pixel

             if ( iand(processing_quality_flags(k, j, 1), 255) == 0 ) then

                n = n + 1

                obstrack%orbitPartIndex(n)        = i
                obstrack%pixelIndex(n)            = k
                obstrack%scanlineIndex(n)         = j
                obstrack%subPixelNumber(n)        = 1
                obstrack%pixelFlag(n)             = processing_quality_flags(k, j, 1)
                obsTrack%subPixelNumber(n)        = pixel_type(k, j, 1)
                obstrack%latitude(n)              = latitude(k, j, 1)
                obstrack%longitude(n)             = longitude(k, j, 1)
                obstrack%cornerLatitude(:,n)      = latitude_bounds(:, k, j, 1)
                obstrack%cornerLongitude(:,n)     = longitude_bounds(:, k, j, 1)
                if ( use_scd_qa4ecv_units ) then
                   ! QA4ECV definition
                   ! units = molecules cm-2, store as (10^15 molecules cm-2) 
                   obstrack%no2SLC(n)                = scd_no2(k, j, 1) / 1.0e+15
                   ! units = molecules cm-2, store as %
                   obstrack%no2SLCError(n)           = scd_no2_precision(k, j, 1) * 100.0 / scd_no2(k, j, 1)
                else
                   ! TROPOMI definition
                   ! units = mol m-2, store as (10^15 molecules cm-2) 
                   obstrack%no2SLC(n)                = scd_no2(k, j, 1) * 6.02214e+19 / 1.0e+15
                   ! units = mol m2, store as %
                   obstrack%no2SLCError(n)           = scd_no2_precision(k, j, 1) * 100.0 / scd_no2(k, j, 1)
                end if
                obstrack%chiSquareFit(n)          = chi_squared(k, j, 1)
                obstrack%solarZenithAngle(n)      = solar_zenith_angle(k, j, 1)
                obstrack%solarAzimuthAngle(n)     = solar_azimuth_angle(k, j, 1)
                obstrack%viewingZenithAngle(n)    = viewing_zenith_angle(k, j, 1)
                obstrack%ViewingAzimuthAngle(n)   = viewing_azimuth_angle(k, j, 1)
                obstrack%cloudFraction(n)         = cloud_fraction_no2(k, j, 1)
                obstrack%cloudTopPressure(n)      = cloud_pressure(k, j, 1)
                obstrack%cloudAlbedo(n)           = cloud_albedo(k, j, 1)
                obstrack%snowIceFlag(n)           = snow_ice_flag(k, j, 1)
                obstrack%sceneAlbedo(n)           = scene_albedo(k, j, 1)
                obstrack%scenePressure(n)         = scene_pressure(k, j, 1)
                obstrack%cloudRadianceFraction(n) = cloud_radiance_fraction_no2(k, j, 1)
                obstrack%surfaceAlbedo(n)         = surface_albedo_no2(k, j, 1)
                obstrack%terrainHeight(n)         = surface_altitude(k, j, 1)
                obstrack%terrainPressure(n)       = surface_pressure(k, j, 1)

             end if 

          end do

       end do

       obstrack%count = n

996    continue

       deallocate ( latitude, stat=status )
       deallocate ( longitude, stat=status )
       deallocate ( latitude_bounds, stat=status )
       deallocate ( longitude_bounds, stat=status )
       deallocate ( solar_zenith_angle, stat=status )
       deallocate ( solar_azimuth_angle, stat=status )
       deallocate ( viewing_zenith_angle, stat=status )
       deallocate ( viewing_azimuth_angle, stat=status )
       deallocate ( surface_altitude, stat=status )
       deallocate ( surface_pressure, stat=status )
       deallocate ( surface_albedo_no2, stat=status )
       deallocate ( scd_no2, stat=status )
       deallocate ( scd_no2_precision, stat=status )
       deallocate ( chi_squared, stat=status )
       deallocate ( cloud_fraction_no2, stat=status )
       deallocate ( cloud_pressure, stat=status )
       deallocate ( cloud_albedo, stat=status )
       deallocate ( cloud_radiance_fraction_no2, stat=status )
       deallocate ( scene_albedo )
       deallocate ( scene_pressure )
       deallocate ( snow_ice_flag )
       deallocate ( pixel_type, stat=status )
       deallocate ( time, stat=status )
       deallocate ( delta_time, stat=status )
       deallocate ( processing_quality_flags, stat=status )

997    continue

       status = TM5IF_close ( )
       if ( status /= 0 ) goto 998

998    continue

    end do

    write(str, *) obstrack%count
    call TM5IF_log_info("  TM5IF_obstrack_read, loaded " // trim(adjustl(str)) // " observations")

999 continue

  end subroutine TM5IF_obstrack_read

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

  subroutine TM5IF_obsfcinfo_write ( obstrack, obsfcinfo, TM5Data )

    ! Update L2 NO2 product

    use MTM5Data,               only : TTM5Data

    implicit none

    ! in/out
    type(TObsTrack), intent(inout)  :: obstrack
    type(TObsFcInfo), intent(inout) :: obsfcinfo
    type(TTM5Data), intent(inout)   :: TM5Data

    ! local
    integer :: i, j, k, l, n, intf
    integer :: status
    real    :: f
    logical :: no_error

    ! main product
    integer, dimension(:,:,:), allocatable   :: processing_error_flag
    real(4), dimension(:,:,:), allocatable   :: amf_total
    real(4), dimension(:,:,:), allocatable   :: amf_trop
    real(4), dimension(:,:,:,:), allocatable :: averaging_kernel
    real(4), dimension(:,:,:), allocatable   :: tm5_surface_pressure
    real(4), dimension(:,:,:), allocatable   :: tm5_tropopause_layer_index
    real(4), dimension(:,:,:), allocatable   :: tropospheric_no2_vertical_column
    real(4), dimension(:,:,:), allocatable   :: tropospheric_no2_vertical_column_precision
    real(4), dimension(:,:), allocatable     :: tm5_pressure_level_a
    real(4), dimension(:,:), allocatable     :: tm5_pressure_level_b

    ! detailed results
    real(4), dimension(:,:,:), allocatable   :: amf_clear
    real(4), dimension(:,:,:), allocatable   :: amf_strat
    real(4), dimension(:,:,:), allocatable   :: ghost_column_no2
    integer, dimension(:,:,:), allocatable   :: processing_quality_flags
    real(4), dimension(:,:,:), allocatable   :: stratospheric_no2_vertical_column
    real(4), dimension(:,:,:), allocatable   :: stratospheric_no2_vertical_column_precision
    real(4), dimension(:,:,:), allocatable   :: summed_no2_total_vertical_column
    real(4), dimension(:,:,:), allocatable   :: summed_no2_total_vertical_column_precision
    real(4), dimension(:,:,:), allocatable   :: total_no2_vertical_column
    real(4), dimension(:,:,:), allocatable   :: total_no2_vertical_column_precision
    real(4), dimension(:,:,:), allocatable   :: cloud_radiance_fraction

    ! start code

    do i = 1, obstrack%norbitparts

        call TM5IF_log_info("TM5IF_obsfcinfo_write: Writing VCD NO2 data to file: ")
        call TM5IF_log_info("    " // trim(obstrack%orbitparts(i)%filename))

        status = TM5IF_open(trim(obstrack%orbitparts(i)%filename), nf90_write)
        if ( status /= 0 ) goto 998

        status = TM5IF_put_att("processing_status", "Retrieval/assimilation step complete, vertical columns and kernels added")
        if ( status /= 0 ) goto 996

        ! === Main product ===

        status = TM5IF_set_group("/PRODUCT" )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "processing_error_flag", processing_error_flag )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "amf_total", amf_total )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "amf_trop", amf_trop )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "averaging_kernel", averaging_kernel )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "tm5_pressure_level_a", tm5_pressure_level_a )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "tm5_pressure_level_b", tm5_pressure_level_b )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "tm5_surface_pressure", tm5_surface_pressure )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "tm5_tropopause_layer_index", tm5_tropopause_layer_index )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "tropospheric_no2_vertical_column", tropospheric_no2_vertical_column )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "tropospheric_no2_vertical_column_precision", tropospheric_no2_vertical_column_precision )
        if ( status /= 0 ) goto 996

        ! === Detailed Results ===

        status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS")
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "amf_clear", amf_clear )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "amf_strat", amf_strat )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "ghost_column_no2", ghost_column_no2 )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "processing_quality_flags", processing_quality_flags )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "stratospheric_no2_vertical_column", stratospheric_no2_vertical_column )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "stratospheric_no2_vertical_column_precision", stratospheric_no2_vertical_column_precision )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "summed_no2_total_vertical_column", summed_no2_total_vertical_column )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "summed_no2_total_vertical_column_precision", summed_no2_total_vertical_column_precision )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "total_no2_vertical_column", total_no2_vertical_column )
        if ( status /= 0 ) goto 996

        status = TM5IF_get_var( "total_no2_vertical_column_precision", total_no2_vertical_column_precision )
        if ( status /= 0 ) goto 996

        if ( ObsFcInfo%cloudRadFraction_computed ) then
           status = TM5IF_get_var( "cloud_radiance_fraction_no2", cloud_radiance_fraction )
           if ( status /= 0 ) goto 996
        end if

        ! QA4ECV units or TROPOMI units (NO2 columns)
        if ( use_scd_qa4ecv_units ) then
           ! QA4ECV definition
           ! units = molecules cm-2, store as (10^15 molecules cm-2) 
           f = 1.0e+15
        else
           f = 1.0e+15 / 6.02214e+19
        end if

        do n = 1, obstrack%count

            if ( obstrack%orbitPartIndex(n) == i ) then

                j = obstrack%scanlineIndex(n)
                k = obstrack%pixelIndex(n)

                ! errors/warnings
                processing_quality_flags(k, j, 1)                    = obsfcinfo%flag(n)

                ! check if first 8 bits are all = 0
                no_error = (iand(obsfcinfo%flag(n), 255) == 0) 
                if ( no_error ) then
                   processing_error_flag(k, j, 1)                       = 0   ! success
                else
                   processing_error_flag(k, j, 1)                       = 1   ! failure
                end if

                ! replace fill values only in case there was no error detected
                if ( no_error ) then

                   ! main product
                   amf_total(k, j, 1)                                   = obsfcinfo%amf(n)
                   amf_trop(k, j, 1)                                    = obsfcinfo%amftrop(n)
                   averaging_kernel(:, k, j, 1)                         = obsfcinfo%avkernel(:, n)
                   tm5_surface_pressure(k, j, 1)                        = obsfcinfo%psurf(n)
                   tm5_tropopause_layer_index(k, j, 1)                  = obsfcinfo%levtropopause(n)
                   tropospheric_no2_vertical_column(k, j, 1)            = obsfcinfo%no2vcdtrop(n) * f 
                   tropospheric_no2_vertical_column_precision(k, j, 1)  = obsfcinfo%no2vcdtropsig(n) * f
                   
                   ! detailed results
                   amf_clear(k, j, 1)                                   = obsfcinfo%amfclear(n)
                   amf_strat(k, j, 1)                                   = obsfcinfo%amfstrat(n)
                   ghost_column_no2(k, j, 1)                            = obsfcinfo%ghostcol(n) * f
                   stratospheric_no2_vertical_column(k, j, 1)           = obsfcinfo%no2vcdstrat(n) * f
                   stratospheric_no2_vertical_column_precision(k, j, 1) = obsfcinfo%no2vcdstratsig(n) * f
                   summed_no2_total_vertical_column(k, j, 1)            = obsfcinfo%no2vcdsum(n) * f
                   summed_no2_total_vertical_column_precision(k, j, 1)  = obsfcinfo%no2vcdsumsig(n) * f
                   total_no2_vertical_column(k, j, 1)                   = obsfcinfo%no2vcd(n) * f
                   total_no2_vertical_column_precision(k, j, 1)         = obsfcinfo%no2vcdsig(n) * f
                   if ( ObsFcInfo%cloudRadFraction_computed ) then
                      cloud_radiance_fraction(k, j, 1)                     = obsfcinfo%cloudradfraction(n)
                   end if

                end if

            end if

        end do

        status = TM5IF_set_group("/PRODUCT")
        if ( status /= 0 ) goto 996

        do l = 1, TM5Data%lm
           do intf = 1, 2
              tm5_pressure_level_a(intf,l) = TM5Data%hyai(l+intf-1)
              tm5_pressure_level_b(intf,l) = TM5Data%hybi(l+intf-1)
           end do
        end do

        status = TM5IF_put_var("processing_error_flag", processing_error_flag)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("amf_total", amf_total)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("amf_trop", amf_trop)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("averaging_kernel", averaging_kernel)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("tm5_pressure_level_a", tm5_pressure_level_a)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("tm5_pressure_level_b", tm5_pressure_level_b)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("tm5_surface_pressure", tm5_surface_pressure)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("tm5_tropopause_layer_index", tm5_tropopause_layer_index)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("tropospheric_no2_vertical_column", tropospheric_no2_vertical_column)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("tropospheric_no2_vertical_column_precision", tropospheric_no2_vertical_column_precision)
        if ( status /= 0 ) goto 996

        status = TM5IF_set_group("/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS")
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("amf_clear", amf_clear)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("amf_strat", amf_strat)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("ghost_column_no2", ghost_column_no2)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("processing_quality_flags", processing_quality_flags)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("stratospheric_no2_vertical_column", stratospheric_no2_vertical_column)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("stratospheric_no2_vertical_column_precision", stratospheric_no2_vertical_column_precision)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("summed_no2_total_vertical_column", summed_no2_total_vertical_column)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("summed_no2_total_vertical_column_precision", summed_no2_total_vertical_column_precision)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("total_no2_vertical_column", total_no2_vertical_column)
        if ( status /= 0 ) goto 996

        status = TM5IF_put_var("total_no2_vertical_column_precision", total_no2_vertical_column_precision)
        if ( status /= 0 ) goto 996

        if ( ObsFcInfo%cloudRadFraction_computed ) then
           status = TM5IF_put_var( "cloud_radiance_fraction_no2", cloud_radiance_fraction )
           if ( status /= 0 ) goto 996
        end if

996     continue

        deallocate( processing_error_flag, stat=status )
        deallocate( amf_total, stat=status )
        deallocate( amf_trop, stat=status )
        deallocate( averaging_kernel, stat=status )
        deallocate( tm5_pressure_level_a, stat=status )
        deallocate( tm5_pressure_level_b, stat=status )
        deallocate( tm5_surface_pressure, stat=status )
        deallocate( tm5_tropopause_layer_index, stat=status )
        deallocate( tropospheric_no2_vertical_column, stat=status )
        deallocate( tropospheric_no2_vertical_column_precision, stat=status )

        deallocate( amf_clear, stat=status )
        deallocate( amf_strat, stat=status )
        deallocate( ghost_column_no2, stat=status )
        deallocate( processing_quality_flags, stat=status )
        deallocate( stratospheric_no2_vertical_column, stat=status )
        deallocate( stratospheric_no2_vertical_column_precision, stat=status )
        deallocate( summed_no2_total_vertical_column, stat=status )
        deallocate( summed_no2_total_vertical_column_precision, stat=status )
        deallocate( total_no2_vertical_column, stat=status )
        deallocate( total_no2_vertical_column_precision, stat=status )

997    continue

       status = TM5IF_close()
       if ( status /= 0 ) goto 998

998    continue

    end do

999 continue

  end subroutine TM5IF_obsfcinfo_write

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

  integer function TM5IF_open ( filename, mode )

    implicit none

    character(len=*), intent(in) :: filename
    integer, intent(in)          :: mode

    integer ::  status

    ncfilename = filename

    call TM5IF_log_debug("Opening '" // trim(ncfilename) // "'")

    status = nf90_open(trim(ncfilename), mode, ncfileid)
    if (status /= 0) then
       call TM5IF_log_error("Failed to open file '"  // trim(ncfilename) // "'")
       TM5IF_open = 1
       return
    end if

    ncgroupid = ncfileid

    TM5IF_open = 0

  end function TM5IF_open

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

  function TM5IF_close()

    integer TM5IF_close
    integer status

    call TM5IF_log_debug("Closing '" // trim(ncfilename) // "'")

    status = nf90_close(ncfileid)
    if (status /= 0) then
       call TM5IF_log_error("Failed to close file '"  // trim(ncfilename) // "'")
       TM5IF_close = 1
       return
    end if

    TM5IF_close = 0

  end function TM5IF_close

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

  function TM5IF_set_group(name)

    integer TM5IF_set_group
    character*(*) name
    integer status
    integer ncnewgroupid
    integer i, j, n

    i = 1
    j = 1
    n = len(trim(name))

    do j = 1, n
       if (j == n) then
          status = nf90_inq_grp_ncid(ncgroupid, name(i:j), ncnewgroupid)
          if (status /= 0) then
             call TM5IF_log_error("Failed to inquire id for group '"  // trim(name(i:j)) // "'")
             TM5IF_set_group = 1
             return
          end if
          ncgroupid = ncnewgroupid
       else if (name(j:j) == "/") then
          if (i == j) then
             ncgroupid = ncfileid
          else
             status = nf90_inq_grp_ncid(ncgroupid, name(i:j-1), ncnewgroupid)
             if (status /= 0) then
                call TM5IF_log_error("Failed to inquire id for group '"  // trim(name(i:j-1)) // "'")
                TM5IF_set_group = 1
                return
             end if
             ncgroupid = ncnewgroupid
          end if
          i = j + 1
       end if
    end do

    TM5IF_set_group = 0

  end function TM5IF_set_group

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

  function TM5IF_get_dimension(name, value)

    integer TM5IF_get_dimension
    character*(*) name
    integer value
    integer status
    integer ncdimid

    status = nf90_inq_dimid(ncgroupid, name, ncdimid)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire id for dimension '"  // name // "'")
       TM5IF_get_dimension = 1
       return
    end if

    status = nf90_inquire_dimension(ncgroupid, ncdimid, len=value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire dimension '" // name // "'")
       TM5IF_get_dimension = 1
       return
    end if

    TM5IF_get_dimension = 0

  end function TM5IF_get_dimension

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

  function TM5IF_get_att_i(name, value)

    integer :: TM5IF_get_att_i
    character*(*) :: name
    integer :: value
    integer :: status

    status = nf90_get_att(ncgroupid, NF90_GLOBAL, name, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get attribute '"  // name // "'")
       TM5IF_get_att_i = 1
       return
    end if

    TM5IF_get_att_i = 0

  end function TM5IF_get_att_i

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

  function TM5IF_get_att_str(name, value)

    integer TM5IF_get_att_str
    character*(*) :: name
    character*(*) :: value
    integer status

    status = nf90_get_att(ncgroupid, NF90_GLOBAL, name, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get attribute '"  // name // "'")
       TM5IF_get_att_str = 1
       return
    end if

    TM5IF_get_att_str = 0

  end function TM5IF_get_att_str

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

  integer function TM5IF_put_att_str ( name, value )

    implicit none

    character(len=*), intent(in) :: name
    character(len=*), intent(in) :: value
    integer :: status

    status = nf90_put_att(ncgroupid, NF90_GLOBAL, name, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to put attribute '"  // name // "'")
       TM5IF_put_att_str = 1
       return
    end if

    TM5IF_put_att_str = 0

  end function TM5IF_put_att_str

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

  function TM5IF_get_var_att_str(varname, name, value)

    integer       :: TM5IF_get_var_att_str
    character*(*) :: varname
    character*(*) :: name
    character*(*) :: value

    integer :: status
    integer :: ncvarid 

    ! first, get the varid corresponding to varname
    status = nf90_inq_varid(ncgroupid, varname, ncvarid)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_att_str = 1
       return
    end if

    ! then, get the attribute of this variable
    status = nf90_get_att(ncgroupid, ncvarid, name, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get attribute '"  // name // "'")
       TM5IF_get_var_att_str = 1
       return
    end if

    TM5IF_get_var_att_str = 0

  end function TM5IF_get_var_att_str

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

  integer function TM5IF_put_var_f2(name, value)

    implicit none

    character(len=*), intent(in) :: name
    real(4), dimension(:,:), intent(inout), allocatable :: value

    integer :: status
    integer :: ncvarid
    integer :: ndims
    integer, dimension(nf90_max_var_dims) :: dimids

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_put_var_f2 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_put_var_f2 = 1
       return
    end if

    if ( ndims /= 2 ) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_put_var_f2 = 1
       return
    end if

    status = nf90_put_var(ncgroupid, ncvarid, value)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to put variable '" // name // "'")
       TM5IF_put_var_f2 = 1
       return
    end if

    TM5IF_put_var_f2 = 0

  end function TM5IF_put_var_f2


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

  integer function TM5IF_put_var_f3(name, value)

    implicit none

    character(len=*), intent(in) :: name
    real(4), dimension(:,:,:), intent(inout),allocatable :: value

    integer :: status
    integer :: ncvarid
    integer :: ndims
    integer, dimension(nf90_max_var_dims) :: dimids

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_put_var_f3 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_put_var_f3 = 1
       return
    end if

    if ( ndims /= 3 ) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_put_var_f3 = 1
       return
    end if

    status = nf90_put_var(ncgroupid, ncvarid, value)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to put variable '" // name // "'")
       TM5IF_put_var_f3 = 1
       return
    end if

    TM5IF_put_var_f3 = 0

  end function TM5IF_put_var_f3

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

  integer function TM5IF_put_var_f4(name, value)

    implicit none

    character(len=*), intent(in) :: name
    real(4), dimension(:,:,:,:), intent(inout),allocatable :: value

    integer :: status
    integer :: ncvarid
    integer :: ndims
    integer, dimension(nf90_max_var_dims) :: dimids

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_put_var_f4 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_put_var_f4 = 1
       return
    end if

    if ( ndims /= 4 ) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_put_var_f4 = 1
       return
    end if

    status = nf90_put_var(ncgroupid, ncvarid, value)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to put variable '" // name // "'")
       TM5IF_put_var_f4 = 1
       return
    end if

    TM5IF_put_var_f4 = 0

  end function TM5IF_put_var_f4

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

  integer function TM5IF_put_var_i3(name, value)

    ! unsigned byte, 3D array

    implicit none

    character(len=*), intent(in) :: name
    integer, dimension(:,:,:), intent(inout), allocatable :: value

    integer :: status
    integer :: ncvarid
    integer :: ndims
    integer, dimension(nf90_max_var_dims) :: dimids

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_put_var_i3 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_put_var_i3 = 1
       return
    end if

    if ( ndims /= 3 ) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_put_var_i3 = 1
       return
    end if

    status = nf90_put_var(ncgroupid, ncvarid, value)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to put variable '" // name // "'")
       TM5IF_put_var_i3 = 1
       return
    end if

    TM5IF_put_var_i3 = 0

  end function TM5IF_put_var_i3

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

  integer function TM5IF_get_var_f2 ( name, value, printErrorMessage )

    implicit none

    character(len=*), intent(in)  :: name
    real(4), dimension(:,:), allocatable, intent(inout) :: value
    logical, intent(in), optional :: printErrorMessage

    integer :: status
    integer :: ncvarid
    integer :: ndims, idim
    integer, dimension(nf90_max_var_dims) :: dimids
    integer, dimension(nf90_max_var_dims) :: dims
    logical :: shout    

    ! messages ?
    shout = .true.
    if ( present(printErrorMessage) ) shout = printErrorMessage

    ! variable should not have been allocated yet
    if ( allocated(value) ) deallocate(value)

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if (status /= 0) then
       if ( shout ) call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_f2 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_get_var_f2 = 1
       return
    end if

    if ( ndims /= 2 ) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_get_var_f2 = 1
       return
    end if

    do idim = 1, ndims

       status = nf90_inquire_dimension(ncgroupid, dimids(idim), len = dims(idim))
       if ( status /= 0 ) then
          call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
          TM5IF_get_var_f2 = 1
          return
       end if

    end do

    allocate(value(dims(1), dims(2)))

    status = nf90_get_var(ncgroupid, ncvarid, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get variable '" // name // "'")
       TM5IF_get_var_f2 = 1
       return
    end if

    TM5IF_get_var_f2 = 0

  end function TM5IF_get_var_f2

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

  integer function TM5IF_get_var_f3 ( name, value, printErrorMessage )

    implicit none

    character(len=*), intent(in) :: name
    real(4), dimension(:,:,:), intent(inout), allocatable :: value
    logical, intent(in), optional :: printErrorMessage

    integer :: status
    integer :: ncvarid
    integer :: ndims, idim
    integer, dimension(nf90_max_var_dims) :: dimids
    integer, dimension(nf90_max_var_dims) :: dims
    logical :: shout    

    ! messages ?
    shout = .true.
    if ( present(printErrorMessage) ) shout = printErrorMessage

    ! variable should not have been allocated yet
    if ( allocated(value) ) deallocate(value)

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if ( status /= 0 ) then
       if ( shout ) call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_f3 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if ( status /= 0 ) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_get_var_f3 = 1
       return
    end if

    if ( ndims /= 3 ) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_get_var_f3 = 1
       return
    end if

    do idim = 1, ndims

       status = nf90_inquire_dimension(ncgroupid, dimids(idim), len = dims(idim))
       if ( status /= 0 ) then
          call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
          TM5IF_get_var_f3 = 1
          return
       end if

    end do

    allocate(value(dims(1), dims(2), dims(3)))

    status = nf90_get_var(ncgroupid, ncvarid, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get variable '" // name // "'")
       TM5IF_get_var_f3 = 1
       return
    end if

    TM5IF_get_var_f3 = 0

  end function TM5IF_get_var_f3

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

  function TM5IF_get_var_f4 ( name, value, printErrorMessage )

    integer TM5IF_get_var_f4
    character*(*) name
    real(4), allocatable :: value(:,:,:,:)
    logical, intent(in), optional :: printErrorMessage

    integer status
    integer ncvarid
    integer ndims, idim
    integer, dimension(nf90_max_var_dims) :: dimids
    integer, dimension(nf90_max_var_dims) :: dims
    logical   :: shout    

    ! messages ?
    shout = .true.
    if ( present(printErrorMessage) ) shout = printErrorMessage

    ! variable should not have been allocated yet
    if ( allocated(value) ) deallocate(value)

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if (status /= 0) then
       if ( shout ) call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_f4 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_get_var_f4 = 1
       return
    end if

    if (ndims /= 4) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_get_var_f4 = 1
       return
    end if

    do idim = 1, ndims

       status = nf90_inquire_dimension(ncgroupid, dimids(idim), len = dims(idim))
       if (status /= 0) then
          call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
          TM5IF_get_var_f4 = 1
          return
       end if

    end do

    allocate(value(dims(1), dims(2), dims(3), dims(4)))

    status = nf90_get_var(ncgroupid, ncvarid, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get variable '" // name // "'")
       TM5IF_get_var_f4 = 1
       return
    end if

    TM5IF_get_var_f4 = 0

  end function TM5IF_get_var_f4

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

  function TM5IF_get_var_i1 ( name, value, printErrorMessage )

    integer TM5IF_get_var_i1
    character*(*) name
    integer, allocatable :: value(:)
    logical, intent(in), optional :: printErrorMessage

    integer status
    integer ncvarid
    integer ndims, idim
    integer, dimension(nf90_max_var_dims) :: dimids
    integer, dimension(nf90_max_var_dims) :: dims
    logical   :: shout    

    ! messages ?
    shout = .true.
    if ( present(printErrorMessage) ) shout = printErrorMessage

    ! variable should not have been allocated yet
    if ( allocated(value) ) deallocate(value)

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if (status /= 0) then
       if ( shout ) call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_i1 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_get_var_i1 = 1
       return
    end if

    if (ndims /= 1) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_get_var_i1 = 1
       return
    end if

    do idim = 1, ndims

       status = nf90_inquire_dimension(ncgroupid, dimids(idim), len = dims(idim))
       if (status /= 0) then
          call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
          TM5IF_get_var_i1 = 1
          return
       end if

    end do

    allocate(value(dims(1)))

    status = nf90_get_var(ncgroupid, ncvarid, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get variable '" // name // "'")
       TM5IF_get_var_i1 = 1
       return
    end if

    TM5IF_get_var_i1 = 0

  end function TM5IF_get_var_i1

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

  function TM5IF_get_var_i2 ( name, value, printErrorMessage )

    integer TM5IF_get_var_i2
    character*(*) name
    integer, allocatable :: value(:,:)
    logical, intent(in), optional :: printErrorMessage

    integer status
    integer ncvarid
    integer ndims, idim
    integer, dimension(nf90_max_var_dims) :: dimids
    integer, dimension(nf90_max_var_dims) :: dims
    logical   :: shout    

    ! messages ?
    shout = .true.
    if ( present(printErrorMessage) ) shout = printErrorMessage

    ! variable should not have been allocated yet
    if ( allocated(value) ) deallocate(value)

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if (status /= 0) then
       if ( shout ) call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_i2 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_get_var_i2 = 1
       return
    end if

    if (ndims /= 2) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_get_var_i2 = 1
       return
    end if

    do idim = 1, ndims

       status = nf90_inquire_dimension(ncgroupid, dimids(idim), len = dims(idim))
       if (status /= 0) then
          call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
          TM5IF_get_var_i2 = 1
          return
       end if

    end do

    allocate(value(dims(1), dims(2)))

    status = nf90_get_var(ncgroupid, ncvarid, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get variable '" // name // "'")
       TM5IF_get_var_i2 = 1
       return
    end if

    TM5IF_get_var_i2 = 0

  end function TM5IF_get_var_i2

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

  function TM5IF_get_var_i3 ( name, value, printErrorMessage )

    integer TM5IF_get_var_i3
    character*(*) name
    integer, allocatable :: value(:,:,:)
    logical, intent(in), optional :: printErrorMessage

    integer status
    integer ncvarid
    integer ndims, idim
    integer, dimension(nf90_max_var_dims) :: dimids
    integer, dimension(nf90_max_var_dims) :: dims
    logical   :: shout    

    ! messages ?
    shout = .true.
    if ( present(printErrorMessage) ) shout = printErrorMessage

    ! variable should not have been allocated yet
    if ( allocated(value) ) deallocate(value)

    status = nf90_inq_varid(ncgroupid, name, ncvarid)
    if (status /= 0) then
       if ( shout ) call TM5IF_log_error("Failed to inquire id for variable '"  // name // "'")
       TM5IF_get_var_i3 = 1
       return
    end if

    status = nf90_inquire_variable(ncgroupid, ncvarid, dimids = dimids, ndims = ndims)
    if (status /= 0) then
       call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
       TM5IF_get_var_i3 = 1
       return
    end if

    if (ndims /= 3) then
       call TM5IF_log_error("Number of dimensions not supported for variable '" // name // "'")
       TM5IF_get_var_i3 = 1
       return
    end if

    do idim = 1, ndims

       status = nf90_inquire_dimension(ncgroupid, dimids(idim), len = dims(idim))
       if (status /= 0) then
          call TM5IF_log_error("Failed to inquire dims for variable '"  // name // "'")
          TM5IF_get_var_i3 = 1
          return
       end if

    end do

    allocate(value(dims(1), dims(2), dims(3)))

    status = nf90_get_var(ncgroupid, ncvarid, value)
    if (status /= 0) then
       call TM5IF_log_error("Failed to get variable '" // name // "'")
       TM5IF_get_var_i3 = 1
       return
    end if

    TM5IF_get_var_i3 = 0

  end function TM5IF_get_var_i3

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

  function TM5IF_convert_date(string, date)

    integer TM5IF_convert_date
    character*(*) string
    type(TM5IF_date) :: date
    integer status

    if (len(string) .lt. 19) then
       TM5IF_convert_date = 1
       return
    end if

    read(string(1:4)  ,*,iostat=status,err=999) date%year
    read(string(6:7)  ,*,iostat=status,err=999) date%month
    read(string(9:10) ,*,iostat=status,err=999) date%day
    read(string(12:13),*,iostat=status,err=999) date%hour
    read(string(15:16),*,iostat=status,err=999) date%minute
    read(string(18:19),*,iostat=status,err=999) date%second

999 continue

    if (status /= 0) then
       TM5IF_convert_date = 1
       return
    end if

    TM5IF_convert_date = 0

  end function TM5IF_convert_date

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

  subroutine TM5IF_log_info(text)

    character*(*) :: text

    write(*,*) trim(text)

  end subroutine TM5IF_log_info

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

  subroutine TM5IF_log_debug(text)

    character*(*) :: text

    !write(*,*) trim(text)

  end subroutine TM5IF_log_debug

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

  subroutine TM5IF_log_error(text)

    character*(*) :: text

    write(*,*) trim(text)

  end subroutine TM5IF_log_error

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

  function TM5IF_time_get_seconds(time)

    integer TM5IF_time_get_seconds
    type(TM5IF_date) :: time

    integer,parameter,dimension(12) :: days = [31,28,31,30,31,30,31,31,30,31,30,31]
    integer,parameter :: epoch = 2000
    integer :: d, i

    d = 0
    do i = epoch, time%year - 1
       d = d + 365
       if ((mod(i, 4) .eq. 0) .and. (mod(i, 100) .ne. 0)) d = d + 1
    end do
    do i = 1, time%month - 1
       d = d + days(i)
       if ((i .eq. 2) .and. (mod(time%year, 4) .eq. 0) .and. (mod(time%year, 100) .ne. 0)) d = d + 1
    end do
    d = d + time%day - 1

    TM5IF_time_get_seconds = d * 86400 + time%hour * 3600 + time%minute * 60 + time%second

  end function TM5IF_time_get_seconds

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

  function TM5IF_time_to_string(time)

    character*(STR_LEN) TM5IF_time_to_string
    type(TM5IF_date) :: time

    write(TM5IF_time_to_string, '(I4.4,A,I2.2,A,I2.2,A,I2.2,A,I2.2,A,I2.2)') time%year, '-', time%month, '-', time%day, 'T', time%hour, ':', time%minute, ':', time%second

  end function TM5IF_time_to_string

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

  subroutine TM5IF_time_init(time, year, month, day, hour, minute, second)

    type(TM5IF_date) :: time
    integer :: year, month, day, hour, minute, second

    time%year = year
    time%month = month
    time%day = day
    time%hour = hour
    time%minute = minute
    time%second = second

  end subroutine TM5IF_time_init

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

end module TM5IF_module