!### macro's #####################################################
!
#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"
!
!#################################################################
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      USER_OUTPUT_MIX
!
! !DESCRIPTION: write instantaneous output fields as total column in Dobson units.
!\\
!\\
! !INTERFACE: 
!
module user_output_mix
  use dims,  only : nregions
  use GO, only : gol, goErr

  implicit none
  ! interface

  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: output_mix, write_mix
  public :: output_mix_init, output_mix_close
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*),  parameter ::  mname = 'user_output_mix'
  !
  ! !PRIVATE TYPES:
  !
  type mixfile
     real,dimension(:,:,:,:),pointer         :: field  ! 3D fields to write
     real,dimension(:,:),pointer             :: ps     ! the surface pressure
     real,dimension(:,:),pointer             :: wda    ! number of averages
     character(len=20),dimension(:),pointer  :: namefield
     character(len=300) :: fname
     integer  :: funit = 0                   ! unit number 
     integer,dimension(:),allocatable  :: var_sds        ! dataset ID 
     real     :: lon_first                   ! x geometry field
     real     :: lon_inc                     ! x geometry field
     real     :: lat_first                   ! y geometry field
     real     :: lat_inc                     ! y geometry field
     integer  :: lon_n                       ! x dimension requested field
     integer  :: lat_n                       ! y dimension requested field
     integer,dimension(6)  :: idate_start    ! date of collection!
     integer,dimension(6)  :: idate_end      ! date of collection!
  end type mixfile

  type(mixfile),dimension(nregions),save  :: mixf

  ! generate output if true:
  logical :: okdebug = .false.
  !
  ! !REVISION HISTORY: 
  !   13 Jun 2013 - W. Peters - added 3D options taken from release/2.0 code
  !   10 Jul 2012 - Ph. Le Sager - okdebug from dims; rid of escape_tm; protex
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------
contains
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_MIX_INIT
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !



  subroutine output_mix_init(status)
    !
    !
    use dims,        only : nregions, newsrun
    use dims,        only : im, jm, lm, region_name, parent,  itau
    use dims,        only : xbeg, xend, ybeg, yend, dx, dy, xref, yref
    use dims,        only : zbeg, zend, dz, zref, idatei
    use chem_param,  only : names

    use ParTools,    only : myid, ntracetloc, ntracet_ar

    implicit none
    ! IO
    integer, intent(out) :: status
    ! local
    integer               :: region, offsetn
    integer               :: i,j,n 
    integer, dimension(6) :: idater
    character(len=200)   :: FFilename
    character(len=2)   :: mynos

    offsetn = sum(ntracet_ar(0:myid-1))
    
    regionloop: do region = 1, nregions

          !  initialize the mix output: instantaneous

          allocate(mixf(region)%field(im(region),jm(region),lm(region),ntracetloc+1))  ! +1 for gph
          allocate(mixf(region)%ps(im(region),jm(region)))
          allocate(mixf(region)%wda(im(region),jm(region)))
          allocate(mixf(region)%namefield(ntracetloc+1))
          allocate(mixf(region)%var_sds(ntracetloc+1))
          do i=1, ntracetloc
             write(mynos,'(i2.2)') i
             mixf(region)%namefield(i) =  trim(names(i+offsetn))//'_'//trim(mynos)
          enddo
          mixf(region)%namefield(ntracetloc+1) = 'gph'

          mixf(region)%lon_first = xbeg(region) + 0.5*dx/xref(region)
          mixf(region)%lat_first = ybeg(region) + 0.5*dy/yref(region)
          mixf(region)%lat_inc   = dy/yref(region)
          mixf(region)%lon_inc   = dx/xref(region)
          mixf(region)%lon_n     = im(region)
          mixf(region)%lat_n     = jm(region)
          mixf(region)%ps        = 0.0
          mixf(region)%field     = 0.0
          mixf(region)%wda       = 0.0

          mixf(region)%idate_start = idatei  ! initialize start date of interval

    end do  regionloop ! nregions

    status = 0

  end subroutine output_mix_init



  subroutine output_mix_close(status)
    !-----------------------------------------------------------
    !
    !-----------------------------------------------------------
    use dims,        only : nregions

    integer, intent(out) :: status
    ! local
    integer               :: region



    ! start
    do region=1,nregions

       call write_mix(region,status)

    enddo   !regions. Do not yet deallocate because parent mix is needed for updates!

    do region=1,nregions
       deallocate(mixf(region)%field)   ! free memory
       deallocate(mixf(region)%ps)   ! free memory
       deallocate(mixf(region)%wda)   ! free memory
       deallocate(mixf(region)%namefield)   ! free memory
       deallocate(mixf(region)%var_sds)   ! free memory
    enddo   !regions


  end subroutine output_mix_close



  subroutine output_mix(region,status)
    !-----------------------------------------------------------
    !
    !-----------------------------------------------------------
    use tracer_data, only : mass_dat
    use global_data, only : region_dat
    use meteodata  , only : m_dat
    use meteo      , only : sp_dat,gph_dat, m_dat
    use dims,        only : im,jm,lm, itaur, ndyn, tref, isr, jsr, ier, jer, at, bt, ndyn, ndyn_max
    use chem_param,  only : names, ra, fscale, ntracet, ntrace_chem
    use binas    ,   only : Avog,Dobs,pi, xmair
    use datetime,    only : tau2date, date2tau
    use toolbox ,    only : escape_tm
    use ParTools,    only : myid, ntracetloc, ntracet_ar

    implicit none

    ! in/out
    integer, intent(in) :: region
    integer, intent(out) :: status

    ! local
    integer               :: io,istat,i,j, imr, jmr, lmr, n, l
    integer               :: sfstart, sfscatt, sfsnatt, sfend
    integer, dimension(6) :: idater
    integer               :: itau_lt, dtime, itau_out12, itau_out00, offsetn, ditau, dtime_max
    real                  :: longitude, weight

    ! MPI arrays to gather fields.
    real,dimension(:,:,:,:),               pointer     :: rm
    real,dimension(:,:,:),                 pointer     :: p
    real,dimension(:,:,:),                 pointer     :: gph
    real,dimension(:,:,:),                 pointer     :: m
    real,dimension(:),                     pointer     :: dxyp
    integer,dimension(:,:),                pointer     :: zoomed

    zoomed => region_dat(region)%zoomed

    dtime = ndyn/tref(region)   ! basic time step (seconds)
    dtime_max = ndyn_max/tref(region)   ! basic time step (seconds)
    weight = float(ndyn)/float(ndyn_max)
    offsetn = sum(ntracet_ar(0:myid-1))

    ! collect output FOR LT
    imr = im(region) ; jmr = jm(region) ; lmr = lm(region)


    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
    p=> sp_dat(region)%data
    gph=> gph_dat(region)%data
    dxyp=> region_dat(region)%dxyp


!  fill fields for averaging:

    do j=jsr(region),jer(region)
       do i=isr(region), ier(region)
          if(zoomed(i,j) /= region) cycle
              ! average output over day:
              mixf(region)%wda(i,j) = mixf(region)%wda(i,j) + weight
              do l=1,lm(region)
                do n=1,ntracetloc
                   mixf(region)%field(i,j,l,n) = &
                      mixf(region)%field(i,j,l,n) + weight*rm(i,j,l,n)/m(i,j,l)*(xmair/ra(offsetn+n))
                enddo
                mixf(region)%field(i,j,l,ntracetloc+1) = &
                     mixf(region)%field(i,j,l,ntracetloc+1) + weight*gph(i,j,l)
             enddo
             mixf(region)%ps(i,j) = &
                mixf(region)%ps(i,j) + weight*p(i,j,1)
       enddo
    enddo
    nullify(m)    ! reset pointers
    nullify(rm)    ! reset pointers
    nullify(p)
    nullify(gph)
    nullify(dxyp)
    nullify(zoomed)

    status = 0

  end subroutine output_mix

subroutine update_parent_mix(region)

  use dims
  use global_data,  only: region_dat
  use toolbox,      only : escape_tm
  use chem_param,   only: ntrace

  IMPLICIT NONE

  !__IO___________________________________________________________________

  INTEGER, INTENT(IN)   :: region

  !__LOCAL_VARIABLES______________________________________________________

  REAL, DIMENSION(:,:,:,:), POINTER    :: mean_mix       ! column value m  - mean
  REAL, DIMENSION(:,:), POINTER    :: mean_ps      ! column value ps  - mean

  REAL, DIMENSION(:,:,:,:), POINTER    :: p_mean_mix       ! column value m  - mean
  REAL, DIMENSION(:,:), POINTER    :: p_mean_ps      ! column value ps  - mean

  REAL, DIMENSION(:), POINTER    :: dxyp

  INTEGER :: i,j
  INTEGER :: ip, jp
  INTEGER :: my_parent
  INTEGER :: xref_,yref_,zref_

  INTEGER :: imp,jmp
  INTEGER :: imr,jmr

  INTEGER :: ic,jc

  INTEGER :: iox,ioy1,ioy2
  REAL    :: w,wtot



  !__START_SUBROUTINE______________________________________________________



  IF (region==1) RETURN


  imr = im(region)
  jmr = jm(region)

  ! determine parent
  my_parent = parent(region)

  xref_ = xref(region)/xref(my_parent)
  yref_ = yref(region)/yref(my_parent)

  imp = im(region)/xref_
  jmp = jm(region)/yref_

  ! check calculated imp, jmp, lmp
  IF(ibeg(region) < iend(region) .and. imp /= iend(region)-ibeg(region)+1) THEN
    call escape_tm('stopped in update_parent_columns')
  ENDIF
    
  IF(jmp.ne.jend(region)-jbeg(region)+1) THEN
    call escape_tm('stopped in update_parent_columns')
  ENDIF  
 

  mean_mix  => mixf(region)%field
  mean_ps  => mixf(region)%ps

  p_mean_mix  => mixf(my_parent)%field
  p_mean_ps => mixf(my_parent)%ps

  dxyp=> region_dat(region)%dxyp

  iox = isr(region)/xref_
  ioy1 = jsr(region)/yref_
  ioy2 = (jm(region)-jer(region)+1)/yref_


  DO jp=jbeg(region)+ioy1, jend(region)-ioy2  
    jc = (jp-jbeg(region))*yref_
    DO ip=ibeg(region)+iox, iend(region)-iox
      ic = (ip-ibeg(region))*xref_

      p_mean_ps(ip,jp) = 0.0
      p_mean_mix(ip,jp,:,:) = 0.0


      wtot = 0.0
      DO j=1,yref_
        w = dxyp(jc+j)
        DO i=1,xref_

          p_mean_mix(ip,jp,:,:) = p_mean_mix(ip,jp,:,:)  + mean_mix(ic+i,jc+j,:,:)*w
          p_mean_ps(ip,jp) = p_mean_ps(ip,jp)  + mean_ps(ic+i,jc+j)*w

          wtot = wtot+w
        ENDDO
      ENDDO
      p_mean_mix(ip,jp,:,:) = p_mean_mix(ip,jp,:,:)/wtot
      p_mean_ps(ip,jp) = p_mean_ps(ip,jp)/wtot


    ENDDO
  ENDDO


  NULLIFY(mean_mix)  
  NULLIFY(mean_ps)  

  NULLIFY(p_mean_mix)  
  NULLIFY(p_mean_ps)  

END SUBROUTINE update_parent_mix

!===========================================================================================================
!===========================================================================================================

  subroutine write_mix(region,status)
    !-----------------------------------------------------------
    !
    !-----------------------------------------------------------
    use global_data, only : mass_dat, region_dat
    use meteo      , only : sp_dat
    use dims,        only : im,jm,lm, itaur, itaui, ndyn, tref, isr, jsr, ier, jer, at, bt, ndyn, ndyn_max
    use dims       , only : nregions, region_name
    use dims       , only : dx, xref, xbeg, xend, ibeg, iend
    use dims       , only : dy, yref, ybeg, yend, jbeg, jend
    use dims       , only : dz, zref, zbeg, zend, lbeg, lend
    use chem_param,  only : names, ra, fscale, ntracet, ntrace_chem
    use binas    ,   only : Avog,Dobs,pi, xmair
    use datetime,    only : tau2date, date2tau
    use toolbox ,    only : escape_tm
    use ParTools,    only : myid, ntracetloc, ntracet_ar
    use MDF      , only : MDF_Open
    use MDF      , only : MDF_Create, MDF_Close, MDF_EndDef
    use MDF      , only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ
    use MDF      , only : MDF_Put_Att
    use MDF      , only : MDF_Def_Dim
    use MDF      , only : MDF_Inq_Dimid, MDF_Inquire_Dimension
    use MDF      , only : MDF_Inq_VarID
    use MDF      , only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var

    implicit none

    ! in/out
    integer, intent(in) :: region
    integer, intent(out) :: status

    character(len=*), parameter  ::  rname = mname//'write_mix'
    ! local
    integer               :: io,istat,i,j, imr, jmr, lmr, n, l
    integer, dimension(6) :: idater
    integer               :: itau_lt, dtime, itau_out12, itau_out00, offsetn, ditau, dtime_max
    real                  :: longitude, weight

    integer,dimension(:,:),                pointer     :: zoomed

    integer :: hid
    integer :: dim_lon_id, dim_lat_id, dim_lev_id
    integer :: var_ps, var_temp, var_lon, var_lat, var_lev
    integer, dimension(6) :: idatei,idatee,idatef
    character(len=1024)   :: attstring
    character(len=200)    :: dsname,attname,comp,pp

    if(itaur(region) == itaui) return

    zoomed => region_dat(region)%zoomed

    dtime = ndyn/tref(region)   ! basic time step (seconds)
    dtime_max = ndyn_max/tref(region)   ! basic time step (seconds)

    ! go to middle of interval 
    call tau2date(itaur(region),idater)

    mixf(region)%idate_end = idater 

    ! Create file plus header

    call mix_file_init ! open mix file

    if ( okdebug ) print *, 'output_mix: Writing mix file output',region

    !WP! First divide each element by weight, then update_parents, and finally
    !write the output.

    do j=jsr(region),jer(region)
       do i=isr(region), ier(region)
          if(zoomed(i,j) /= region) cycle
          if(mixf(region)%wda(i,j)>0) &
          mixf(region)%ps(i,j)  =  mixf(region)%ps(i,j)/mixf(region)%wda(i,j)
       enddo
    enddo
    do n=1,ntracetloc+1
      do l=1,lm(region)     
         do j=jsr(region),jer(region)
            do i=isr(region),ier(region)
               if(zoomed(i,j) /= region) cycle
               if(mixf(region)%wda(i,j)>0) &
               mixf(region)%field(i,j,l,n) =  mixf(region)%field(i,j,l,n)/   &
                                                   mixf(region)%wda(i,j)
            enddo
         enddo
      enddo   
    enddo   

    !WP! Added explicit call to update_parent to ensure that inner domain of
    !parents is filled with appropriate values. This was missing before. Note
    !that the order in which to write mix files *MUST* be from smallest zoom
    !to parents.

    call update_parent_mix(region)

    call MDF_Put_Var( hid, var_lon, (/(xbeg(region)+i*dx/xref(region),i=0,im(region)-1)/), status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Var( hid, var_lat, (/(ybeg(region)+i*dy/yref(region),i=0,jm(region)-1)/), status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Var( hid, var_lev, (/(i,i=1,lm(region))/), status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Var( hid, var_ps, mixf(region)%ps , status )
    IF_NOTOK_RETURN(status=1)

    do n=1,ntracetloc+1

      var_temp = mixf(region)%var_sds(n)
      call MDF_Put_Var( hid, var_temp, mixf(region)%field(:,:,:,n),status)
      IF_ERROR_RETURN(status=1)

    enddo  !n

    call MDF_Close( hid, status )
    IF_NOTOK_RETURN(status=1)

    mixf(region)%ps        = 0.0
    mixf(region)%field     = 0.0
    mixf(region)%wda       = 0.0
    mixf(region)%var_sds   = 0
    mixf(region)%idate_start = mixf(region)%idate_end ! end of interval becomes start of next interval  

    nullify(zoomed)

   contains

    subroutine mix_file_init
    use dims,        only : nregions
    use dims,        only : im, jm, lm, region_name, parent,  itaur, itau
    use dims,        only : xbeg, xend, ybeg, yend, dx, dy, xref, yref
    use dims,        only : zbeg, zend, dz, zref, at, bt
    use global_data, only : region_dat, outdir
    use datetime,    only : tau2date
    use partools,    only : myid, ntracetloc, ntracet_ar
    use MDF      , only : MDF_Open
    use MDF      , only : MDF_Create, MDF_Close, MDF_EndDef
    use MDF      , only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ
    use MDF      , only : MDF_Put_Att
    use MDF      , only : MDF_Def_Dim
    use MDF      , only : MDF_Inq_Dimid, MDF_Inquire_Dimension
    use MDF      , only : MDF_Inq_VarID
    use MDF      , only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var

    integer               :: sfstart, sfscatt, sfsnatt, sfend
    integer, dimension(6) :: idater, idatee
    character(len=12)      :: idatestext,idateetext ! arj change, was len=10
    character(len=2)      :: myids,myno
    integer		  :: io,istat

    write(myids,'(i2.2)') myid
     
    idater=mixf(region)%idate_start  ! use start of interval in file name
    idatee=mixf(region)%idate_end  ! use end of interval

    write(idatestext, '(i4,4i2.2)') idater(1), idater(2), idater(3), idater(4), idater(5)  ! arj change, was up to idater(4)
    write(idateetext, '(i4,4i2.2)') idatee(1), idatee(2), idatee(3), idatee(4), idatee(5)  ! arj change, was up to idater(4)
    mixf(region)%fname = trim(outdir)//'/molefractions_'//region_name(region)//'_'//idatestext//'_'//idateetext//'.nc'

    call MDF_Create( trim(mixf(region)%fname), MDF_NETCDF, MDF_REPLACE, hid, status )
    IF_NOTOK_RETURN(status=1)

    ! Global Attributes 

    write(attstring,'(i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC")') idater
    call MDF_Put_Att(hid, MDF_GLOBAL,"model_start_date",values=trim(attstring),status=status)
    IF_NOTOK_RETURN(status=1)

    write(attstring,'(i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC")') idatee
    call MDF_Put_Att(hid, MDF_GLOBAL,"model_end_date",values=trim(attstring),status=status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att( hid,MDF_GLOBAL, 'ibeg', ibeg(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'iend', iend(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'jbeg', jbeg(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'jend', jend(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'lbeg', lbeg(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'lend', lend(region), status )

    call MDF_Put_Att( hid,MDF_GLOBAL, 'xref', xref(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'yref', yref(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'zref', zref(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'tref', tref(region), status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'at', at, status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'bt', bt, status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'idate_start', mixf(region)%idate_start, status )
    call MDF_Put_Att( hid,MDF_GLOBAL, 'idate_end', mixf(region)%idate_end, status )


    ! define dimensions:                                                                                                                                                                                            
    call MDF_Def_Dim( hid, 'lon', im(region) , dim_lon_id, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Def_Dim( hid, 'lat', jm(region) , dim_lat_id, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Def_Dim( hid, 'lev', lm(region) , dim_lev_id, status )
    IF_NOTOK_RETURN(status=1)

    ! dimension variables:                                                                                                                                                                                                    
    call MDF_Def_Var( hid, 'lon', MDF_FLOAT, (/dim_lon_id/), var_lon, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lon,"units",values="degrees East",status=status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lon,"comment",values="center gridbox",status=status)
    IF_NOTOK_RETURN(status=1)

    call MDF_Def_Var( hid, 'lat', MDF_FLOAT, (/dim_lat_id/), var_lat, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lat,"units",values="degrees North",status=status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lat,"comment",values="center gridbox",status=status)
    IF_NOTOK_RETURN(status=1)

    call MDF_Def_Var( hid, 'lev', MDF_FLOAT, (/dim_lev_id/), var_lev, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lat,"units",values="unitless",status=status)
    IF_NOTOK_RETURN(status=1)

    ! define other variables: 

    call MDF_Def_var( hid,'ps' ,MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_ps, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid,var_ps,'idate_start',mixf(region)%idate_start,status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid,var_ps,'idate_end',mixf(region)%idate_end,status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid,var_ps,'units','Pa',status)
    IF_ERROR_RETURN(status=1)

    do n=1,ntracetloc+1

        dsname=mixf(region)%namefield(n)

        call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id,dim_lev_id/), var_temp, status )
        IF_NOTOK_RETURN(status=1)
        mixf(region)%var_sds(n) = var_temp

        call MDF_Put_Att(hid, var_temp,'idate_start',mixf(region)%idate_start,status)
        IF_NOTOK_RETURN(status=1)
        call MDF_Put_Att(hid, var_temp,'idate_end',mixf(region)%idate_end,status)
        IF_NOTOK_RETURN(status=1)
        if(n .le. ntracetloc) then
           call MDF_Put_Att(hid, var_temp,'units','mole/mole',status)
           IF_NOTOK_RETURN(status=1)
        else
           call MDF_Put_Att(hid, var_temp,'units','m',status)
           IF_NOTOK_RETURN(status=1)
        endif

    enddo

    ! finished definition:
    call MDF_EndDef( hid, status )
    IF_NOTOK_RETURN(status=1)

    end subroutine mix_file_init

  end subroutine write_mix
  
end module user_output_mix
