ecearth3/sources/tm5mp/proj/output/user_output.F90

!
#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_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
#include "output.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      USER_OUTPUT
!
! !DESCRIPTION:  Contains calls to user-specific output routines, e.g.
!                instantaneous mix files, station output, output of flight 
!                tracks etc.
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT

  use GO,                  only : gol, goErr, goPr
  use GO,             ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start
#ifdef with_hdf4
  use user_output_noaa,    only : noaa_data
#endif
#ifdef with_m7  
  use user_output_aerocom, only : aerocom_freq, aerocom_exper, aerocom_dhour
  use user_output_aerchemmip, only : experiment,experiment_id, realm,source_type,aerchemmip_dhour, crescendo_out
  use user_output_general, only : gen_freq, gen_exper, all_chemistry
  use user_output_general, only : nCCNdiag, nSat, SuperSat
#endif
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: user_output_init
  public :: user_output_step
  public :: user_output_mean
  public :: user_output_done
  !
  ! !PRIVATE DATA MEMBERS:
  !
  logical           :: settings_data = .false.  ! signal for model settings output
  logical           :: flight_data   = .false.  ! signal for flight output
  logical           :: station_data  = .true.   ! signal for station output
  logical           :: mix_data      = .false.  ! signal for mix output
  integer           :: mix_data_dhour           ! ... and its output period 
  logical           :: mmix_data     = .false.  ! signal for mean mix output
  logical           :: aerchemmip    = .false.  ! signal for AERCHEMMIP diagnostics
  logical           :: aerchemmip_1h    = .true.  ! signal for AERCHEMMIP diagnostics
  logical           :: aerocom_data  = .false.  ! signal for AEROCOM diagnostics
  !integer           :: aerocom_dhour            ! ... and its calling period 

  logical           :: aerocom_stat  = .false.  ! signal for AEROCOM station diagnostics
  logical           :: output_pdump  = .false.  ! signal for chemistry time series output
  integer           :: output_pdump_dsec        ! ... and its lowest output period 
  logical           :: column_data              ! signal for column output
  integer           :: column_data_dtsec        ! ... and its output period (in seconds)
  logical           :: flask_data = .false.     ! signal for flask output
  logical           :: gen_data  = .false.  ! signal for GEN diagnostics
  integer           :: gen_dhour            ! ... and its calling period 
  logical           :: gen_2d  = .false.    ! signal for GEN surface and budget
  logical           :: gen_3d  = .false.    ! signal for GEN 3D tracer output
 
#ifdef with_cf_output
  logical           ::  output_cf = .false.
  integer           ::  output_cf_dhour         ! every dhour hour
#endif
  
  character(len=*), parameter  ::  mname = 'user_output'  

  !Timers:
  integer ::  itim_aero_init, itim_aero_step, itim_aero_write 
  integer ::  itim_gen_init, itim_gen_step, itim_gen_write 
!, itim_nh3, itim_sox, itim_dms, itim_ch4, itim_isop, itim_rn222

  !
  ! !REVISION HISTORY: 
  !    6 Feb 2011 - Achim Strunk - Added aerocom-2 diagnostics output. 
  !    9 Jul 2012 - Ph. Le Sager - merged (1) with version in proj/chem/base (to
  !                                account for optics data for
  !                                user_output_station) and (2) with version in
  !                                the base/trunk (to account for output_common)
  !                                and (3) with version in Climaqs proj to
  !                                incorporate pdump as replacement for retro output
  !    3 Jul 2013 - Andy Jacobson
  !        - Updated modules to use dynamic time-step weighting
  !        - Added flask, column,mmix output to trunk
  !        - Revised mix to used netCDF4 output
  !
  !   25 Aug 2015 - T. van Noije - Added switch for monthly AeroCom output 
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    USER_OUTPUT_INIT
  !
  ! !DESCRIPTION: Initialise user-specified model output (all regions)
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine user_output_init( status )
    !
    ! !USES:
    !    
    use GO,                   only : TrcFile, Init, Done, ReadRc
    use dims,	              only : ndyn_max, nregions
    use MeteoData,            only : Set, temper_dat
    use global_data,          only : rcfile
    use User_Output_Common,   only : User_Output_Common_Init
    use User_Output_Settings, only : User_Output_Settings_Init
    use user_output_mmix,     only : mmix_init
#ifdef with_hdf4    
    use user_output_station,  only : read_stationlist, init_station_output
    use user_output_noaa,     only : init_noaa_events
#endif
    use user_output_mix,      only : user_output_mix_init, mix_netcdf_attributes
    use user_output_column,   only : user_output_column_init
    use user_output_pdump,    only : output_pdump_init
    use user_output_flask,    only : user_output_flask_init
#ifdef with_cf_output
    use user_output_cf,       only : output_cf_init
#endif
! #ifdef with_optics
!     Use Optics             , only : Optics_Init, Optics_RcInit
! #endif
! #ifdef with_pm
!     Use PM                 , only : PM_Init
! #endif
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)      ::  status
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk - modified for aerocom-2 diagnostics.
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/user_output_init'
    
    ! --- local -------------------------------

    type(TrcFile)           ::  rcF
    Character (len=200)     ::  wavelengthsstring
    Character (len=200)     ::  pmsizelimits
    Character (len=300)     ::  alphastring
    integer :: n
    
    ! --- begin -------------------------------

    ! init common stuff:
    call User_Output_Common_Init( status )
    IF_NOTOK_RETURN(status=1)
    
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    
    ! ------------
    !  S E T T I N G S
    ! ------------
    ! put out settings ?
    call ReadRc( rcF, 'settings.output', settings_data, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    ! apply ?
    if ( settings_data ) then
      write (gol,'(a,":   user output settings ...")') rname; call goPr
      call User_Output_Settings_Init( rcF, status )
      IF_NOTOK_RETURN(status=1)
    end if 
    
    ! ------------
    !  S T A T I O N S
    ! ------------
    call ReadRc( rcF, 'output.station', station_data, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    if ( station_data ) then
#ifdef with_hdf4          
       call read_stationlist(status)
       IF_NOTOK_RETURN(status=1)
       call init_station_output(status)
       IF_NOTOK_RETURN(status=1)
#else
       status=1
       write (gol,'(a,":   USER_OUTPUT_STATION not available without HDF4")') rname; call goErr
       IF_NOTOK_RETURN(status=1)
#endif
    end if
    
    ! ------------
    !  N O A A
    ! ------------
#ifdef with_hdf4          
    call ReadRc( rcF, 'output.noaa', noaa_data, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    if ( noaa_data ) then
       call init_noaa_events(status)
       IF_NOTOK_RETURN(status=1)
    end if
#endif
    
    ! ------------
    !  F L I G H T
    ! ------------
    call ReadRc( rcF, 'output.flight', flight_data, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    
    ! ------------
    !  F L A S K
    ! ------------
    call ReadRc( rcF, 'output.flask', flask_data, status, default=.false.)
    IF_ERROR_RETURN(status=1)
    if(flask_data) then
       call user_output_flask_init(rcF,status)
       IF_NOTOK_RETURN(status=1)
    endif
    
    ! ------------
    !  M I X
    ! ------------
    call ReadRc( rcF, 'output.mix', mix_data, status )
    IF_NOTOK_RETURN(status=1)
    if ( mix_data ) then

      call ReadRc( rcF, 'output.mix.dhour', mix_data_dhour, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.notes', mix_netcdf_attributes%notes, status, default="")
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.disclaimer', mix_netcdf_attributes%disclaimer, status, default="")
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.email', mix_netcdf_attributes%email, status, default="")
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.url', mix_netcdf_attributes%url, status, default="")
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.institution', mix_netcdf_attributes%institution, status, default="")
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.conventions', mix_netcdf_attributes%conventions, status, default="CF-1.5")
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.mix.attributes.source', mix_netcdf_attributes%source, status, default="")
      IF_ERROR_RETURN(status=1)

      call user_output_mix_init(status)
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! ------------
    !  M M I X
    ! ------------
    ! initialise accumulation of the mean mixing ratio fields
    call ReadRc( rcF, 'output.mmix', mmix_data, status )
    IF_NOTOK_RETURN(status=1)
    if ( mmix_data ) then
       call mmix_init(status)
       IF_NOTOK_RETURN(status=1)
       
       ! require temperature then
       do n = 1, nregions
          call Set(  temper_dat(n), status, used=.true. )
       end do
    end if
    
    ! ------------
    !  A E R O C O M
    ! ------------
#ifdef with_m7
    ! initialise AEROCOM Diagnostics (hourly, daily or monthly)
    call ReadRc( rcF, 'output.aerocom', aerocom_data, status )
    IF_NOTOK_RETURN(status=1)
    
    if( aerocom_data ) then 
       call ReadRc( rcF, 'output.aerocom.freq', aerocom_freq, status, default='monthly' )
       IF_NOTOK_RETURN(status=1)
       
       call ReadRc( rcF, 'output.aerocom.exper', aerocom_exper, status, default='AP3' )
       IF_NOTOK_RETURN(status=1) 
       
       call ReadRc( rcF, 'output.aerocom.dhour', aerocom_dhour, status )
       IF_NOTOK_RETURN(status=1)

       call ReadRc( rcF, 'output.aerocom.stations', aerocom_stat, status )
       IF_NOTOK_RETURN(status=1)
    end if

    call ReadRc( rcF, 'output.general', gen_data, status )
    IF_NOTOK_RETURN(status=1)
    
    if( gen_data ) then 
       call ReadRc( rcF, 'output.general.freq', gen_freq, status, default='hourly' )
       IF_NOTOK_RETURN(status=1)
       call ReadRc( rcF, 'output.general.exper', gen_exper, status, default='AP3' )

       IF_NOTOK_RETURN(status=1) 
       !output frequency
       call ReadRc( rcF, 'output.general.dhour', gen_dhour, status )
       IF_NOTOK_RETURN(status=1)
       ! output all chemistry species or predefined smaller set (see user_output_general: variable gas_output
       call ReadRc( rcF, 'output.general.all_chemistry', all_chemistry, status, default=.false.)
       IF_NOTOK_RETURN(status=1)
       ! output 2D-fields: emi,load,surfconc,wetdep,drydep,sed,optics
       call ReadRc( rcF, 'output.general.2d', gen_2d, status )
       IF_NOTOK_RETURN(status=1)
       ! output 3D-fields: particle number/mass concentrations and diameters, chemistry species, 
       call ReadRc( rcF, 'output.general.3d', gen_3d, status )
       IF_NOTOK_RETURN(status=1)
       !call ReadRc( rcF, 'output.general.author', gen_author, status )
       !IF_NOTOK_RETURN(status=1)
       !call ReadRc( rcF, 'output.general.institute', gen_institute, status )
       !IF_NOTOK_RETURN(status=1)
       !call ReadRc( rcF, 'output.general.', gen_, status )
       !IF_NOTOK_RETURN(status=1)
       call ReadRc( rcF, 'input.diagnostic.CCN', nCCNdiag, status)
       IF_NOTOK_RETURN(status=1) 
       IF (nCCNdiag) THEN
          call ReadRc( rcF, 'input.supersat', nSat, status)
          IF_NOTOK_RETURN(status=1) 
          ALLOCATE(SuperSat(nSat))
          call ReadRc( rcF, 'supersat.values', SuperSat, status)
          SuperSat(:) = SuperSat(:) / 100.e0 ! convert from %RH
          IF_NOTOK_RETURN(status=1) 
       ENDIF
    end if
#endif    
    ! ------------
    !  A E R C H E M M I P
    ! ------------
#ifdef with_m7
    ! initialise AERCHEMMIP Diagnostics (hourly, daily or monthly)
    call ReadRc( rcF, 'output.aerchemmip', aerchemmip, status )
    IF_NOTOK_RETURN(status=1)
    if( aerchemmip ) then 
       call ReadRc( rcF, 'output.aerchemmip.1h', aerchemmip_1h, status )
       IF_NOTOK_RETURN(status=1)

       ! NOT IMPLEMENTED - 6hr- ouput is harcoded in aerchem_output.F90
       !        call ReadRc( rcF, 'output.aerchemmip.dhour', aerchemmip_dhour, status )
       !       IF_NOTOK_RETURN(status=1)
       aerchemmip_dhour=1
       
       call ReadRc( rcF, 'output.aerchemmip.experiment', experiment, status, default='AerChemMIP' )
       IF_ERROR_RETURN(status=1) 
       call ReadRc( rcF, 'output.aerchemmip.realm', realm, status, default='atmosChem' )
       IF_ERROR_RETURN(status=1) 
       call ReadRc( rcF, 'output.aerchemmip.sourcetype', source_type, status, default='AP3' )
       IF_ERROR_RETURN(status=1) 
       call ReadRc( rcF, 'output.aerchemmip.experimentid', experiment_id, status )
       IF_NOTOK_RETURN(status=1) 
        call ReadRc( rcF, 'output.crescendo', crescendo_out, status )
       IF_NOTOK_RETURN(status=1)

    end if
#endif    

    ! ------------
    !  P D U M P
    ! ------------
    ! put out in pdump format ?
    call ReadRc( rcF, 'output.pdump', output_pdump, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    ! init if necessary:
    if ( output_pdump ) then
       ! init output; return ouptut time step
       call Output_PDUMP_Init( rcF, output_pdump_dsec, status )
       IF_NOTOK_RETURN(status=1)
    end if

#ifdef with_cf_output
    ! ------------
    !  C F
    ! ------------
    ! put out in retro format ?
    call ReadRc( rcF, 'output.cf', output_cf, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    ! init if necessary:
    if ( output_cf ) then
      ! init output; return ouptut time step
      call Output_CF_Init( rcF, output_cf_dhour, status )
      IF_NOTOK_RETURN(status=1)
    end if
#endif


    ! ------------
    !  O P T I C S
    ! ------------    
    ! Optics. I just need one thing of the rc-file before closing it.
    
    ! We could remove the compiler flags with_optics and with_pm if we want to
    ! We just do with_optics if you include at least one wavelength
    ! We just do with_pm if you include at least one sizelimit
    ! Now, we simply use the compiler flags.
    
! #ifdef with_optics
!     call ReadRc( rcF, 'optics.lookup.table', lookuptable,status) ! the only line in climaqs proj. FIXME: others needed?
!     IF_NOTOK_RETURN(status=1)
!     call ReadRc( rcF, 'optics.wavelengths', wavelengthsstring,status)
!     IF_NOTOK_RETURN(status=1)
!     call ReadRc( rcF, 'optics.alphas', alphastring,status,default="")
!     IF_ERROR_RETURN(status=1)
!     Call Optics_RCInit(wavelengthsstring,alphastring,status) ! This code is so ugly, I do not want to pollute Optics_Init with this.
!     IF_NOTOK_RETURN(status=1)
!     Call Optics_Init(lookuptable,station_data,status)
!     IF_NOTOK_RETURN(status=1)
! #endif
    ! PM output.
! #ifdef with_pm
!     Call ReadRc( rcF, 'pm.sizelimits', pmsizelimits,status)
!     IF_NOTOK_RETURN(status=1)
!     ! Do not do anything if we have no stations
!     if (station_data) then
!        Call PM_Init(pmsizelimits,status)
!        IF_NOTOK_RETURN(status=1)
!     end if
! #endif


    ! ------------
    !  C O L U M N
    ! ------------    
    ! output column data?
    call ReadRc( rcF, 'output.column', column_data, status, default = .false. )
    IF_ERROR_RETURN( status=1 )

    ! apply?
    if ( column_data ) then
      call ReadRc( rcF, 'output.column.dtsec', column_data_dtsec, status , default = ndyn_max)
      IF_ERROR_RETURN(status=1)
      write (gol,*) trim(mname)//'/column_data_dtsec:', column_data_dtsec; call goPr
      call user_output_column_init( status )
      IF_NOTOK_RETURN( status=1 )
   end if
   
    call GO_Timer_Def( itim_aero_init, 'user_output_aerocom_init', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_aero_step, 'user_output_aerocom_step', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_aero_write, 'user_output_aerocom_write', status )
    IF_NOTOK_RETURN(status=1)

    call GO_Timer_Def( itim_gen_init, 'user_output_general_init', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_gen_step, 'user_output_general_step', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_gen_write, 'user_output_general_write', status )
    IF_NOTOK_RETURN(status=1)

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

    ! ok
    status = 0

  END SUBROUTINE USER_OUTPUT_INIT
  !EOC


  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    USER_OUTPUT_DONE
  !
  ! !DESCRIPTION: Finalise user-specified model output for the region given
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine user_output_done( status )
    !
    ! !USES:
    !
    use dims,                 only : nregions,newhour,newday,newmonth
    use User_Output_Settings, only : User_Output_Settings_Done
    use User_Output_Common,   only : User_Output_Common_Done
    use user_output_mmix,     only : write_mmix, mmix_Done
#ifdef with_hdf4          
    use user_output_station,  only : free_stationfields
    use user_output_noaa,     only : write_noaa_events
#endif
    use user_output_mix    ,  only : user_output_mix_done
    use user_output_column ,  only : user_output_column_done
    use user_output_flask ,   only : user_output_flask_done
#ifdef with_m7

    use user_output_aerocom,  only : output_aerocom_done, output_aerocom_write,output_aerocom_step
    use user_output_aerchemmip,  only : output_aerchemmip_done, output_aerchemmip_write, output_aerchemmip_write_hourly, output_aerchemmip_write_daily, crescendo_out, accumulate_data, output_aerchemmip_write_6hourly
    use user_output_general,  only : output_general_done, output_general_write
    use user_output_general,  only : nCCNdiag, SuperSat
#endif
    use user_output_pdump,    only : output_pdump_done
#ifdef with_cf_output
    use user_output_cf,       only : output_cf_done
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk - Added aerocom-2 case.
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/user_output_done'

    integer :: region

    ! --- begin -----------------------------

    ! settings output enabled ?
    if ( settings_data ) then
      ! done with module:
      call User_Output_Settings_Done( status )
      IF_NOTOK_RETURN(status=1)
    end if 
    
    
    ! I do the done of Optics first, since it borrows the station's hdf.
    ! I will give user_output_station the honor to finalize the station's hdf.
    ! Maybe not necessary, but it looks nicer first to finalize the sds's
    ! and then finalize the hdf.
! #ifdef with_optics
!     Call Optics_Done(station_data,status)
!     IF_NOTOK_RETURN(status=1)
! #endif
! #ifdef with_pm
!     if (station_data) then
!        Call PM_Done(status)
!        IF_NOTOK_RETURN(status=1)
!     end if
! #endif

    if ( column_data ) then
      call user_output_column_done( status )
      IF_NOTOK_RETURN( status=1 )
    endif

    ! write the mean mixing ratio fields to file
    if ( mmix_data ) then
      call write_mmix(status)
      IF_NOTOK_RETURN(status=1)
      call mmix_Done(status)
      IF_NOTOK_RETURN(status=1)      
    end if
    
#ifdef with_hdf4          
    if ( station_data) then 
       call free_stationfields(status)
       IF_NOTOK_RETURN(status=1)
    endif
#endif
    
#ifdef with_m7
    if( aerocom_data ) then 

       ! start timing:
       call GO_Timer_Start( itim_aero_write, status )
       IF_NOTOK_RETURN(status=1)

       ! write last data set before done
       do region = 1, nregions 
          call output_aerocom_write(region, aerocom_stat, status)
          IF_NOTOK_RETURN(status=1)
       end do
       call GO_Timer_End( itim_aero_write, status )
       IF_NOTOK_RETURN(status=1)

       call output_aerocom_done(aerocom_stat, status)
       IF_NOTOK_RETURN(status=1)
    end if
    
    if( gen_data ) then 
       ! start timing:
       call GO_Timer_Start( itim_gen_write, status )
       IF_NOTOK_RETURN(status=1)

       ! write last data set before done
       do region = 1, nregions 
          call output_general_write(region, status)
          IF_NOTOK_RETURN(status=1)
       end do
       call GO_Timer_End( itim_gen_write, status )
       IF_NOTOK_RETURN(status=1)

       call output_general_done( status)
       IF_NOTOK_RETURN(status=1)
    endif

    if ( nCCNdiag ) then
       deallocate(SuperSat)
    endif
    if( aerchemmip ) then 
       
       region=1
       if (newhour(1))then
          call accumulate_data(aerchemmip_dhour,.false.,status)
       end if
       !if (.not. newday)then
          call output_aerchemmip_write(region, .true.,status)
          ! write last data set before done
          IF_NOTOK_RETURN(status=1)
       !end if
       !if (.not. newday)then
          call output_aerchemmip_write_daily(region, .true.,status)
          ! write last data set before done
          IF_NOTOK_RETURN(status=1)
       !end if

       ! There should be no unfinished hours to write out
       ! so the run will end at the end of 1h time step, when
       ! output is already written out
          if (.not. newhour(1))then
             call output_aerchemmip_write_6hourly(region, .true.,status)
             IF_NOTOK_RETURN(status=1)
          end if
       if (newhour(1))then
          call output_aerchemmip_write_hourly(region, .true.,status)
          IF_NOTOK_RETURN(status=1)
       end if
       call output_aerchemmip_done(status)
       IF_NOTOK_RETURN(status=1)
       
    end if

#endif

#ifdef with_hdf4
    if ( noaa_data) then 
       call write_noaa_events(status)
       IF_NOTOK_RETURN(status=1)
    endif
#endif
    if(flask_data) then
       call user_output_flask_done(status)
      IF_NOTOK_RETURN(status=1)
    endif

    if ( mix_data ) then 
       call user_output_mix_done( status )
       IF_NOTOK_RETURN(status=1)
    endif

    if ( output_pdump ) then
       call Output_Pdump_Done( status )
       IF_NOTOK_RETURN(status=1)
    end if
    
#ifdef with_cf_output
    if ( output_cf ) then
       call Output_cf_Done( status )
       IF_NOTOK_RETURN(status=1)
    end if
#endif
    
    ! done with common stuff:
    call User_Output_Common_Done( status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  END SUBROUTINE USER_OUTPUT_DONE
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    USER_OUTPUT_STEP
  !
  ! !DESCRIPTION: Define user-specified model output for the region given
  !               Called every time step
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine user_output_step( region, status )
    !
    ! !USES:
    ! 
    use dims,                 only : itaur, newsrun, itaui, newday
#ifdef with_m7
    use dims,                 only : newmonth, newhour,newday
#endif
    use datetime,             only : tau2date
#ifdef with_hdf4          
    use user_output_station,  only : output_stationconc  
    use user_output_noaa,     only : get_noaa  
#endif
    use user_output_flight,   only : get_flightdata
    use user_output_column,   only : user_output_column_accum,user_output_column_evaluate,user_output_column_write,user_output_column_reset
    use user_output_mix,      only : user_output_mix_accum,user_output_mix_write
    use user_output_mmix,     only : accumulate_mmix
    use user_output_flask,    only : user_output_flask_sample
    use user_output_pdump   , only : Output_Pdump_Step
#ifdef with_m7
    use user_output_aerocom,  only : output_aerocom_step, output_aerocom_init
    use user_output_aerocom,  only : output_aerocom_write
    use user_output_aerchemmip,  only : output_aerchemmip_init,output_aerchemmip_write,accumulate_data, output_aerchemmip_done,output_aerchemmip_write_hourly,output_aerchemmip_write_daily, output_aerchemmip_write_6hourly
    use user_output_general,  only : output_general_step, output_general_init
    use user_output_general,  only : output_general_write
#endif
! #ifdef with_optics
!     Use Optics,               only : Optics_Step
! #endif
! #ifdef with_pm
!     Use PM,                   only : PM_Step
! #endif
#ifdef with_cf_output
    use user_output_cf      , only : Output_cf_Step
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)       ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk - modified for aerocom-2 diagnostics.
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/user_output_step'
    
    ! --- local ------------------------------
    
    integer               :: total_sec
    integer, dimension(6) :: idate_f
#ifdef with_m7
    logical               :: newaer
#endif

    ! --- begin ------------------------------
    
    call tau2date(itaur(region),idate_f)

    if ( mmix_data ) then
       call accumulate_mmix( region,status )
       IF_NOTOK_RETURN(status=1)
    endif
    

! #ifdef with_optics
!     Call Optics_Step(region,status)
!     IF_NOTOK_RETURN(status=1)
! #endif
    
#ifdef with_hdf4          
    if ( station_data ) then
       call output_stationconc(region, status)
       IF_NOTOK_RETURN(status=1)
! #ifdef with_pm
!        Call PM_Step(region,status)
!        IF_NOTOK_RETURN(status=1)
! #endif
    endif

    ! The third keyword here ("force") is for instantaneous sampling.  Default
    ! value of .false. yields a four-hour window centered on sample time.
    !
    if ( noaa_data ) call get_noaa(region,itaur(region),.false.)
#endif
    
    if ( flight_data ) then
      call get_flightdata(region, idate_f, status)
      IF_NOTOK_RETURN(status=1)
    endif
    
    if(flask_data) then
       call user_output_flask_sample(region,itaur(region),status)
      IF_NOTOK_RETURN(status=1)
    endif

    if ( mix_data ) then
      if ( modulo(itaur(region)-itaui,mix_data_dhour*3600) == 0 ) then
         call user_output_mix_write(region, status )
         IF_NOTOK_RETURN(status=1)
      endif
      call user_output_mix_accum(region, status)
      IF_NOTOK_RETURN(status=1)
    end if

    if (column_data) then
      call user_output_column_accum( region, status )
      IF_NOTOK_RETURN( status=1 )
       if ( (itaur(region)-itaui .gt. 0) .and. (modulo(itaur(region)-itaui,column_data_dtsec) == 0 )) then
          call user_output_column_evaluate( status )
          IF_NOTOK_RETURN(status=1)
          call user_output_column_write( region, status )
          IF_NOTOK_RETURN(status=1)
          call user_output_column_reset( status )
          IF_NOTOK_RETURN(status=1)
       endif
    endif


#ifdef with_m7
    ! -------------------
    ! AEROCOM DIAGNOSTICS
    ! -------------------
    IF( aerocom_data ) then 

       ! write (for previous hour, day or month) in case we reach a new hour, day or month
       select case (trim(aerocom_freq))
       case ('hourly')
          if (modulo(idate_f(4),aerocom_dhour)==0)then
             newaer = newhour(region)
          else
             newaer = .false.
          end if
       case ('daily')
          newaer = newday
       case ('monthly')
          newaer = newmonth
       end select
       
       ! start timing:
       call GO_Timer_Start( itim_aero_write, status )
       IF_NOTOK_RETURN(status=1)

       if ( newaer .and. .NOT. newsrun ) then
          call output_aerocom_write(region, aerocom_stat, status)
          IF_NOTOK_RETURN(status=1)
       end if

       call GO_Timer_End( itim_aero_write, status )
       IF_NOTOK_RETURN(status=1)

       ! start timing:
       call GO_Timer_Start( itim_aero_init, status )
       IF_NOTOK_RETURN(status=1)

       ! initialise for each new hour, day or month
       if ( newaer ) then
          call output_aerocom_init( aerocom_stat, status )
          IF_NOTOK_RETURN(status=1)
       end if
       ! stop timing
       call GO_Timer_End( itim_aero_init, status )
       IF_NOTOK_RETURN(status=1)

       ! start timing:
       call GO_Timer_Start( itim_aero_step, status )
       IF_NOTOK_RETURN(status=1)

       ! do an accumulation step 
       if ( (modulo(idate_f(4),aerocom_dhour)==0) .and. all(idate_f(5:6)==0) ) then
          
          call output_aerocom_step( region, aerocom_dhour, aerocom_stat, status )
          IF_NOTOK_RETURN(status=1)
       end if
       !stop timing
       call GO_Timer_End( itim_aero_step, status )
       IF_NOTOK_RETURN(status=1)
    end IF

    if (gen_data) then
       ! write (for previous hour, day or month) in case we reach a new hour, day or month
       select case (trim(gen_freq))
         case ('hourly')
            newaer = newhour(region)
         case ('daily')
           newaer = newday
         case ('monthly')
           newaer = newmonth
       end select
       ! start timing:
       call GO_Timer_Start( itim_gen_write, status )
       IF_NOTOK_RETURN(status=1)

      if ( newaer .and. .NOT. newsrun ) then
          call output_general_write(region, status)
          IF_NOTOK_RETURN(status=1)
       end if

       call GO_Timer_End( itim_gen_write, status )
       IF_NOTOK_RETURN(status=1)

       call GO_Timer_Start( itim_gen_init, status )
       IF_NOTOK_RETURN(status=1)

       ! initialise for each new hour, day or month
       if ( newaer ) then
          call output_general_init( status )
          IF_NOTOK_RETURN(status=1)
       end if
       ! stop timing
       call GO_Timer_End( itim_gen_init, status )
       IF_NOTOK_RETURN(status=1)

       call GO_Timer_Start( itim_gen_step, status )
       IF_NOTOK_RETURN(status=1)

       ! do an accumulation step 
       if ( (modulo(idate_f(4),gen_dhour)==0) .and. all(idate_f(5:6)==0) ) then
          call output_general_step( region, gen_dhour, status )
          IF_NOTOK_RETURN(status=1)
       end if
       !stop timing
       call GO_Timer_End( itim_gen_step, status )
       IF_NOTOK_RETURN(status=1)

    end IF
    ! -------------------
!
    if (aerchemmip) then
    
       if (newsrun) then
          call output_aerchemmip_init( status )
          IF_NOTOK_RETURN(status=1)
       end if

       !accumulate every hour
       !if ( newhour(1) .and. .not. newsrun ) then
       if ( newhour(1)  ) then
          !(modulo(idate_f(4),aerchemmip_dhour)==0) .and. all(idate_f(5:6)==0) ) then
          
          if (newsrun) then
             ! First step of the year, do instantaneous output from restarts
             ! so do only ec550aer -> .true.
             call accumulate_data(aerchemmip_dhour,.true.,status)
          else
             call accumulate_data(aerchemmip_dhour,.false.,status)
          end if
          IF_NOTOK_RETURN(status=1)
       end if

       !monthly
       if (newmonth .and. .NOT. newsrun ) then
          !if (newhour .and. .NOT. newsrun ) then
          call output_aerchemmip_write(region,  newhour(1), status)
          IF_NOTOK_RETURN(status=1)
       end if

       !hourly
       !newhour(1) means newhour at region 1 (global, now zooming expected in AERCHEMMIP)
       ! aerchemmip hourly output is a mean over 1 hour, so no writing out 
       ! initial values 
       if (newhour(1) .and. .NOT. newsrun .and. aerchemmip_1h) then
          call output_aerchemmip_write_hourly(region,  newhour(1),status)
          IF_NOTOK_RETURN(status=1)
       end if
       !6hourly
       !newhour(1) means newhour at region 1 (global, now zooming expected in AERCHEMMIP)

       ! Since 6 hourly output is instantaneous and we want to have output for hour 0 of the new year
       ! we output even on the first hour, for the new year this will be the last 
       ! value of the previous year which is read from a restart file but seems to check out
       ! see #xxx-note-x in ec-earth portal
       if (newhour(1)  .and. modulo(idate_f(4),6)==0) then
          call output_aerchemmip_write_6hourly(region,  newhour(1),status)
          IF_NOTOK_RETURN(status=1)
       end if

       !daily
       if (newday .and. .NOT. newsrun ) then
          call output_aerchemmip_write_daily(region,  newhour(1),status)
          IF_NOTOK_RETURN(status=1)
       end if

    end if
#endif

!    
    if ( output_pdump ) then
       total_sec = idate_f(4)*3600 + idate_f(5)*60 + idate_f(6)
       if ( modulo(total_sec,output_pdump_dsec) == 0 ) then
          call Output_Pdump_Step( region, idate_f, status )
          IF_NOTOK_RETURN(status=1)
       end if
    end if

#ifdef with_cf_output
    if ( output_cf ) then
      if ( (modulo(idate_f(4),output_cf_dhour)==0) .and. all(idate_f(5:6)==0) ) then
        call Output_CF_Step( region, idate_f, status )
        IF_NOTOK_RETURN(status=1)
      end if
    end if
#endif

    ! ok
    status = 0

  END SUBROUTINE USER_OUTPUT_STEP
  !EOC
    
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    USER_OUTPUT_MEAN
  !
  ! !DESCRIPTION: Calculate and write means at stations. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine user_output_mean(status)
    !
    ! !USES:
    !
    use dims, only : itau, ndyn_max
#ifdef with_hdf4          
    use user_output_station, only : evaluate_stationconc, reset_stationconc_accumulator, write_stationconc
#endif    
! #ifdef with_optics
!     Use Optics , only : Optics_Write
! #endif
! #ifdef with_pm
!     Use PM, only : PM_Write
! #endif
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)    :: status
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/user_output_mean'
  
#ifdef with_hdf4          
    IF(station_data)THEN
      IF(mod(itau, ndyn_max) == 0) THEN
         CALL evaluate_stationconc(status)
         IF_NOTOK_RETURN(status=1)
         CALL write_stationconc(status)
         IF_NOTOK_RETURN(status=1)
         CALL reset_stationconc_accumulator
      ENDIF
    ENDIF
#endif
  
! #ifdef with_optics
!     ! Write Optics as if it is a station file, so at the end of one whole time step.
!     ! We might have a parallellization problem, since we are parallel in levels or tracers,
!     ! and I do not know what is what. However, the AOD is already integrated over levels
!     ! and tracers, so it may not make a difference after all.
!     Call Optics_Write (station_data,status)
!     IF_NOTOK_RETURN(status=1)
! #endif

! #ifdef with_pm
!     ! Write particulate matter like the optics data into the station files
!     If (station_data) Then
!        Call PM_Write (status)
!        IF_NOTOK_RETURN(status=1)
!     End If
! #endif

    ! ok
    status = 0

  end subroutine user_output_mean
  
  
END MODULE USER_OUTPUT