!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (rcode/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (rcode> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
#include "output.inc"  
!
!----------------------------------------------------------------------------
!                    TM5                                                    !
!----------------------------------------------------------------------------
!BOP
!
! !MODULE:      USER_OUTPUT_MMIX
!
! !DESCRIPTION: Handles initialisation, accumulation and output of monthly
!               mean tracers mixing ratio, temperature, and pressure.
!               For run shorter than a month, average is done over the
!                duration of the run.
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_MMIX
  !
  ! !USES:
  !
  use GO, only : gol, goPr, goErr, goBug, goLabel
  use dims !, mname_dims=>mname
  use chem_param, only: ntrace, nstd

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: mmix_Init, mmix_Done  ! object life cycle
  public :: accumulate_mmix       ! accumulate data
  public :: reset_mmix            ! reset data
  public :: write_mmix            ! write data 
  public :: mmix_dat, w_mmix
  !
  ! !PUBLIC TYPES:
  !
  type, public :: mmix_data
     real,dimension(:,:,:,:),pointer     :: rmmix
     real,dimension(:,:,:,:),pointer     :: std_mmix
     real,dimension(:,:,:),pointer       :: tempm
     real,dimension(:,:),pointer         :: presm
  end type mmix_data

  type(mmix_data), dimension(nregions), target :: mmix_dat  ! accumulated data
  real,            dimension(nregions)         :: w_mmix    ! accumulated weight

  integer, dimension(6) :: startdate
  
  character(len=*), parameter  ::  mname = 'user_output_mmix'
  !
  ! !REVISION HISTORY: 
  !   16 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MMIX_INIT
  !
  ! !DESCRIPTION: Allocate and set to 0 fields to accumulate.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE MMIX_INIT(rcode)
    !
    ! !USES:
    !
    use tm5_distgrid,  only : dgrid, Get_DistGrid
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  rcode
    !
    ! !REVISION HISTORY: 
    !   16 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/Init_mmix'

    integer :: region, i1, i2, j1, j2

    call goLabel(rname)

    do region=1,nregions

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

       allocate(mmix_dat(region)%rmmix   (i1:i2, j1:j2, lm(region), ntrace) )
       allocate(mmix_dat(region)%std_mmix(i1:i2, j1:j2, lm(region), nstd)   )
       allocate(mmix_dat(region)%tempm   (i1:i2, j1:j2, lm(region))         )
       allocate(mmix_dat(region)%presm   (i1:i2, j1:j2)                     )
       w_mmix(region) = 0.0

       mmix_dat(region)%rmmix    = 0.0
       mmix_dat(region)%std_mmix = 0.0
       mmix_dat(region)%presm    = 0.0
       mmix_dat(region)%tempm    = 0.0
    end do
    
    startdate = idatei ! start of model run
    
    call goLabel(); rcode=0

  end subroutine mmix_Init
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MMIX_DONE
  !
  ! !DESCRIPTION: deallocate data
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE MMIX_DONE(rcode)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  rcode
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/mmix_Done'

    integer :: region

    call goLabel( rname )

    do region = 1,nregions
       deallocate ( mmix_dat(region)%rmmix   )
       deallocate ( mmix_dat(region)%std_mmix)
       deallocate ( mmix_dat(region)%tempm   )
       deallocate ( mmix_dat(region)%presm   )
    end do
    
    call goLabel(); rcode=0
    
  END SUBROUTINE MMIX_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESET_MMIX
  !
  ! !DESCRIPTION: reset to zero all accumulated fields
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESET_MMIX( rcode )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  rcode
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/reset_mmix'

    integer :: region

    call goLabel( rname )

    do region=1,nregions

       w_mmix(region) = 0.0

       mmix_dat(region)%rmmix    = 0.0
       mmix_dat(region)%std_mmix = 0.0
       mmix_dat(region)%presm    = 0.0
       mmix_dat(region)%tempm    = 0.0

    enddo
    
    startdate = idate
    
    ! ok
    call goLabel(); rcode=0

  END SUBROUTINE RESET_MMIX
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   ACCUMULATE_MMIX
  !
  ! !DESCRIPTION: accumulates fields
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ACCUMULATE_MMIX( region, rcode)
    !
    ! !USES:
    !
    !use dims,          only : lm, okdebug, ndyn, ndyn_max
    use chem_param,    only : istd , nstd, ntrace, ntracet, ntrace_chem
    use global_data,   only : region_dat, mass_dat,  chem_dat
    use meteodata  ,   only : sp_dat, temper_dat, m_dat
    use partools,      only : isRoot
    use tm5_distgrid,  only : dgrid, Get_DistGrid
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)     :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  rcode
    !
    ! !REVISION HISTORY: 
    !   16 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/accumulate_mmix'

    real,dimension(:,:,:,:),pointer      :: rm, rmmix, std_mmix, rmc
    real,dimension(:,:,:),  pointer      :: t, tempm, m
    real,dimension(:,:,:),  pointer      :: p
    real,dimension(:,:),    pointer      :: presm
    integer,dimension(:,:), pointer      :: zoomed
    integer         :: n, i1, i2, j1, j2, itrace, i,j,l
    real            :: weight

    ! start
    call goLabel( rname )
    
    if (newmonth.and.(.not.newsrun)) then
       call write_mmix(rcode)
       IF_NOTOK_RETURN(rcode=1)
       call reset_mmix(rcode)
       IF_NOTOK_RETURN(rcode=1)
    endif
    
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    ! debug
    if (region /= 1) then
       print*, rname//" - unexpected region - FIXME ZOOM"; rcode=1
       IF_NOTOK_RETURN(rcode=1)
    end if


    ! input
    rm  => mass_dat(region)%rm(i1:i2,j1:j2,:,:)
    if ( ntrace_chem > 0 ) rmc => chem_dat(region)%rm(i1:i2,j1:j2,:,:)
    
    m  =>      m_dat(region)%data(i1:i2,j1:j2,:)
    t  => temper_dat(region)%data(i1:i2,j1:j2,:)
    p  =>     sp_dat(region)%data(i1:i2,j1:j2,:)

    ! output
    rmmix    => mmix_dat(region)%rmmix
    std_mmix => mmix_dat(region)%std_mmix
    tempm    => mmix_dat(region)%tempm
    presm    => mmix_dat(region)%presm

    weight = float(ndyn)/float(ndyn_max)

    ! mmix....
    do n=1, ntracet
       rmmix(:,:,:,n) = rmmix(:,:,:,n) + weight*rm(:,:,:,n)/m(:,:,:)
    end do


    ! non-transported tracers, if any
    if ( ntrace_chem > 0 ) then
       do n= ntracet+1, ntracet+ntrace_chem
          rmmix(:,:,:,n) = rmmix(:,:,:,n) + weight*rmc(:,:,:,n-ntracet)/m(:,:,:)
       end do
    end if

    ! stdt deviations
    do n=1,nstd
       itrace = istd(n)
       std_mmix(:,:,:,n) = std_mmix(:,:,:,n) + weight*(rm(:,:,:,itrace)/m(:,:,:))**2
    end do

    ! met fields
    tempm = tempm + weight*t
    presm = presm + weight*p(:,:,1)
    
    w_mmix(region) = w_mmix(region) + weight

    if ( okdebug ) print*, 'accumulate_mmix: region ',region, &
         '; w_mmix',w_mmix(region)

    nullify(m)
    nullify(rm)
    nullify(t)
    nullify(p)
    nullify(rmmix)
    nullify(std_mmix)
    nullify(presm)
    nullify(tempm)
    
    ! ok
    call goLabel(); rcode=0

  END SUBROUTINE ACCUMULATE_MMIX
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    WRITE_MMIX
  !
  ! !DESCRIPTION: Write to file mean fields defined in this module, plus all
  !               essential model parameters 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE WRITE_MMIX( rcode )
    !
    ! !USES:
    !
    use chem_param
#ifdef with_hdf4
    use io_hdf
#endif
    use datetime,      only : tstamp
    use global_data,   only : outdir
    use ParTools,      only : Par_Barrier, isRoot
    use User_Output_Common, only : User_Output_Check_Overwrite
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)   ::  rcode
    !
    ! !REVISION HISTORY: 
    !   16 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/write_mmix'

#ifdef with_hdf4
    integer           :: istat, sfsnatt, sfscatt, io, sfstart
    integer           :: sfend
#endif
    integer           :: region
    integer           :: i,j,l,n,k,ind
    !integer,dimension(6) :: 
    character(len=12) :: name
    character(len=200):: FFilename
    
    ! start
    call goLabel( rname )

    REG : do region=nregions,1,-1
       
#ifdef with_hdf4
       IOroot : if ( isRoot ) then

          write (FFilename,'(a,"/mmix_",i4.4,3i2.2,"_",i4.4,3i2.2,"_",a,".hdf")') &
!               trim(outdir), idatei(1:4), idatee(1:4),trim(region_name(region)) 
               trim(outdir), startdate(1:4), idate(1:4),trim(region_name(region)) 

          ! check existence ...
          call User_Output_Check_Overwrite( trim(FFilename), rcode )
          IF_NOTOK_RETURN(rcode=1)

          ! open
          io = sfstart(FFilename,DFACC_CREATE)

          if ( io < 0 ) then
             write (gol,'("While starting mmix file")'); call goErr
             write (gol,'("Filename:",a)') trim(Ffilename); call goErr
             write (gol,'("IIO unit:",i10)') io; call goErr
             write (gol,'("in program")'); call goErr; rcode=1; return
          end if

          !write (gol,'("write_mmix: io unit",i10)') io; call goPr

          istat = sfsnatt(io,'itau',  DFNT_INT32, 1, itau)
          istat = sfsnatt(io,'nregions',  DFNT_INT32, 1, nregions)
          istat = sfscatt(io,'region_name',  DFNT_CHAR, &
               len_trim(region_name(region)), trim(region_name(region)) )
          istat = sfsnatt(io,'im',    DFNT_INT32, 1, im(region))
          istat = sfsnatt(io,'jm',    DFNT_INT32, 1, jm(region))
          istat = sfsnatt(io,'lm',    DFNT_INT32, 1, lm(region))
          istat = sfsnatt(io,'dx',    DFNT_FLOAT64, 1, dx/xref(region))
          istat = sfsnatt(io,'dy',    DFNT_FLOAT64, 1, dy/yref(region))
          istat = sfsnatt(io,'dz',    DFNT_FLOAT64, 1, dz/zref(region))
          istat = sfsnatt(io,'xbeg',  DFNT_INT32, 1, xbeg(region))
          istat = sfsnatt(io,'xend',  DFNT_INT32, 1, xend(region))
          istat = sfsnatt(io,'ybeg',  DFNT_INT32, 1, ybeg(region))
          istat = sfsnatt(io,'yend',  DFNT_INT32, 1, yend(region))
          istat = sfsnatt(io,'zbeg',  DFNT_INT32, 1, zbeg(region))
          istat = sfsnatt(io,'zend',  DFNT_INT32, 1, zend(region))
          if(region/=1) then
             istat = sfsnatt(io,'ibeg',  DFNT_INT32, 1, ibeg(region))
             istat = sfsnatt(io,'iend',  DFNT_INT32, 1, iend(region))
             istat = sfsnatt(io,'jbeg',  DFNT_INT32, 1, jbeg(region))
             istat = sfsnatt(io,'jend',  DFNT_INT32, 1, jend(region))
             istat = sfsnatt(io,'lbeg',  DFNT_INT32, 1, lbeg(region))
             istat = sfsnatt(io,'lend',  DFNT_INT32, 1, lend(region))
          end if
          istat = sfsnatt(io,'xref',  DFNT_INT32, 1, xref(region))
          istat = sfsnatt(io,'yref',  DFNT_INT32, 1, yref(region))
          istat = sfsnatt(io,'zref',  DFNT_INT32, 1, zref(region))
          istat = sfsnatt(io,'tref',  DFNT_INT32, 1, tref(region))
          istat = sfsnatt(io,'ntrace',DFNT_INT32, 1, ntrace)
          istat = sfsnatt(io,'ntracet',DFNT_INT32, 1, ntracet)
          istat = sfsnatt(io,'nstd',DFNT_INT32, 1, nstd)
          istat = sfsnatt(io,'idate' ,DFNT_INT32, 6, idate)
          istat = sfsnatt(io,'istart',  DFNT_INT32, 1, istart)
          istat = sfsnatt(io,'ndyn_max',  DFNT_INT32, 1, ndyn_max)
          istat = sfsnatt(io,'nconv', DFNT_INT32, 1, nconv)
          istat = sfsnatt(io,'ndiag', DFNT_INT32, 1, ndiag)
          istat = sfsnatt(io,'nchem', DFNT_INT32, 1, nchem)
          istat = sfsnatt(io,'nsrce', DFNT_INT32, 1, nsrce)
          istat = sfsnatt(io,'nread', DFNT_INT32, 1, nread)
          istat = sfsnatt(io,'nwrite',DFNT_INT32, 1, nwrite)
          istat = sfsnatt(io,'ninst', DFNT_INT32, 1, ninst)
          istat = sfsnatt(io,'ncheck',DFNT_INT32, 1, ncheck)
          istat = sfsnatt(io,'ndiff', DFNT_INT32, 1, ndiff)
          istat = sfsnatt(io,'itaui',    DFNT_INT32, 1, itaui)
          istat = sfsnatt(io,'itaue',    DFNT_INT32, 1, itaue)
          istat = sfsnatt(io,'itaut',    DFNT_INT32, 1, itaut)
          istat = sfsnatt(io,'itau0',    DFNT_INT32, 1, itau0)
          istat = sfsnatt(io,'idatei' ,  DFNT_INT32, 6, idatei)
          istat = sfsnatt(io,'idatee' ,  DFNT_INT32, 6, idatee)
          istat = sfsnatt(io,'idatet' ,  DFNT_INT32, 6, idatet)
          istat = sfsnatt(io,'idate0' ,  DFNT_INT32, 6, idate0)
          istat = sfsnatt(io,'icalendo' ,DFNT_INT32, 1, icalendo)
          istat = sfsnatt(io,'iyear0' ,  DFNT_INT32, 1, iyear0)
          istat = sfsnatt(io,'julday0' , DFNT_INT32, 1, julday0)
          istat = sfsnatt(io,'ndiagp1' , DFNT_INT32, 1, ndiagp1)
          istat = sfsnatt(io,'ndiagp2' , DFNT_INT32, 1, ndiagp2)
          istat = sfsnatt(io,'nstep'   , DFNT_INT32, 1, nstep)
          istat = sfsnatt(io,'cpu0'   ,  DFNT_FLOAT64, 1, cpu0)
          istat = sfsnatt(io,'cpu1'   ,  DFNT_FLOAT64, 1, cpu1)
          istat = sfsnatt(io,'ra'     ,  DFNT_FLOAT64, ntracet, ra)
          istat = sfsnatt(io,'fscale' ,  DFNT_FLOAT64, ntrace, fscale)
          istat = sfscatt(io,'names'  ,  DFNT_CHAR, ntrace*tracer_name_len, names)
          istat = sfsnatt(io,'areag'  ,  DFNT_FLOAT64, 1, areag)
          istat = sfsnatt(io,'czeta'  ,  DFNT_FLOAT64, 1, czeta)
          istat = sfsnatt(io,'czetak'  , DFNT_FLOAT64, 1, czetak)
          istat = sfscatt(io,'xlabel'  , DFNT_CHAR, 160, xlabel)
          istat = sfsnatt(io,'istd'    , DFNT_INT32, nstd, istd)
          istat = sfsnatt(io,'newyr'   , DFNT_INT32, 1, newyr)
          istat = sfsnatt(io,'newmonth', DFNT_INT32, 1, newmonth)
          istat = sfsnatt(io,'newday'  , DFNT_INT32, 1, newday)
          istat = sfsnatt(io,'newsrun' , DFNT_INT32, 1, newsrun)
!          istat = sfsnatt(io,'cdebug'  , DFNT_INT32, 1, cdebug)
          istat = sfsnatt(io,'limits'  , DFNT_INT32, 1, limits)
          istat = sfsnatt(io,'nregions'  , DFNT_INT32, 1, nregions)
          istat = sfsnatt(io,'w_mmix'  , DFNT_FLOAT64, 1, w_mmix(region))
          istat = sfsnatt(io,'at'  , DFNT_FLOAT64,lm(1)+1, at)
          istat = sfsnatt(io,'bt'  , DFNT_FLOAT64,lm(1)+1, bt)
#ifdef slopes
#ifndef secmom
          istat = sfscatt(io,'adv_scheme'  , DFNT_CHAR, 5, 'slope')
#else
          istat = sfscatt(io,'adv_scheme'  , DFNT_CHAR, 5, '2nd_m')
#endif
#endif
          istat = sfsnatt(io,'nsplitsteps'  , DFNT_INT32, 1, nsplitsteps)
          istat = sfscatt(io,'splitorder'  , DFNT_CHAR,  nsplitsteps, splitorder)
       end if IOroot
       
       CALL WRITEMMIX_REGION( region, rcode )
       IF_NOTOK_RETURN(rcode=1)

       call par_barrier

       if ( isRoot ) then
          write(gol,'("write mmix: sfend returns",i4)') sfend(io) ; call goPr
       endif
#else
    write (gol,'("not compiled with hdf4")'); call goErr
    TRACEBACK; rcode=1; return
#endif // with_hdf4

    END DO REG

    ! ok
    call goLabel(); rcode=0


  CONTAINS

    !------------------------------------------------------------------------
    !                    TM5                                                !
    !------------------------------------------------------------------------
    !BOP
    !
    ! !IROUTINE:    WRITEMMIX_REGION
    !
    ! !DESCRIPTION: Write to file mean fields defined in this module
    !\\
    !\\
    ! !INTERFACE:
    !
    SUBROUTINE WRITEMMIX_REGION( region, rcode ) 
      !
      ! !USES:
      !
      use global_data,   only : region_dat
      use tm5_distgrid,  only : dgrid, Get_DistGrid, gather
      use ParTools,      only : isRoot
      !
      ! !INPUT PARAMETERS:
      !
      integer,intent(in)      :: region
      !
      ! !INPUT/OUTPUT PARAMETERS:
      !
      integer, intent(inout)  :: rcode
      !
      ! !REVISION HISTORY: 
      !   17 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
      !
      ! !REMARKS:
      !
      !EOP
      !---------------------------------------------------------------------
      !BOC
      character(len=*), parameter ::  rname = mname//'/write_mmix_region'

      real,dimension(:,:,:,:),pointer        :: rmmix, std_mmix
      real,dimension(:,:,:)  ,pointer        :: tempm
      real,dimension(:,:)    ,pointer        :: presm
      real,dimension(:,:,:,:),allocatable    :: rmmix_glb
      real,dimension(:,:,:,:),allocatable    :: std_mmix_glb
      integer :: imr, jmr, lmr, nsend, i1, i2, j1, j2, n
      real    :: ahelp,ahelp1

      ! start
      call goLabel( rname )

#ifdef with_hdf4

      ! global and local indices
      imr = im(region)
      jmr = jm(region)
      lmr = lm(region)
      call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
      
      ! array to gather data
      if (isRoot) then 
         allocate(   rmmix_glb(imr,jmr,lmr,ntrace) )
         allocate(std_mmix_glb(imr,jmr,lmr,nstd)   )
      else
         allocate(   rmmix_glb(1,1,1,1) )
         allocate(std_mmix_glb(1,1,1,1) )
      end if

      ! shortcut
      rmmix    => mmix_dat(region)%rmmix
      std_mmix => mmix_dat(region)%std_mmix
      tempm    => mmix_dat(region)%tempm
      presm    => mmix_dat(region)%presm

      ! 
      !write(gol,'("writemmix_region: w_mmix",f12.2)') w_mmix(region); call goPr      
      if ( w_mmix(region) > 2 ) then

         ! scale stdt deviation tracers
         do k= 1,nstd
         do l= 1,lmr
         do j=j1,j2
         do i=i1,i2
            
            n = istd(k)
            ahelp = rmmix(i,j,l,n)

            ahelp1= fscale(n) * &
                 (std_mmix(i,j,l,k) - ahelp*ahelp/w_mmix(region) ) / &
                 (w_mmix(region)-1)

            std_mmix(i,j,l,k)=max(1e-35,ahelp1)

         end do
         end do
         end do
         end do

         ! scale tracers
         do n=1,ntrace
            rmmix(:,:,:,n) = fscale(n) * rmmix(:,:,:,n)/w_mmix(region)
         end do

         ! scale temperature and pressure
         tempm = tempm/w_mmix(region)
         presm = presm/w_mmix(region)

         ! GATHER AND WRITE METEO
         call gather(dgrid(region), tempm, rmmix_glb(:,:,:,1), 0, rcode)         
         IF_NOTOK_RETURN(rcode=1)
         
         if (isRoot) then
            call io_write3d_32d(io,imr,'LON'//trim(region_name(region)), &
                 jmr,'LAT'//trim(region_name(region)),lmr, &
                 'HYBRID',rmmix_glb(1:imr,1:jmr,1:lmr,1),'tempm',idate)
         end if
         
         call gather(dgrid(region), presm, rmmix_glb(:,:,1,1), 0, rcode)
         IF_NOTOK_RETURN(rcode=1)
         
         if (isRoot) then
            call io_write2d_32d(io,imr,'LON'//trim(region_name(region)), &
                 jmr,'LAT'//trim(region_name(region)), &
                 rmmix_glb(1:imr,1:jmr,1,1),'presm',idate)
         end if
         
         ! gather & write mmix
         call gather(dgrid(region), rmmix, rmmix_glb, 0, rcode)
         IF_NOTOK_RETURN(rcode=1)

         if (isRoot) then
            do n=1,ntrace
               name=names(n)
               call io_write3d_32d(io,imr,'LON'//trim(region_name(region)), &
                    jmr,'LAT'//trim(region_name(region)),lmr,'HYBRID', &
                    rmmix_glb(1:imr,1:jmr,1:lmr,n),name,idate)
            end do
         end if

         ! gather and write stdt dev
         call gather(dgrid(region), std_mmix, std_mmix_glb, 0, rcode)
         IF_NOTOK_RETURN(rcode=1)
         
         if (isRoot) then
            do k=1,nstd
               n = istd(k)
               name = 'std_'//names(n)
               call io_write3d_32d(io,imr,'LON'//trim(region_name(region)), &
                    jmr,'LAT'//trim(region_name(region)),lmr,'HYBRID', &
                    std_mmix_glb(1:imr,1:jmr,1:lmr,k),name,idate)
            end do
         end if

      else
         write(gol,'("writemmix_region: NO WRITING! just avoided division by zero!")') ; call goPr
      end if
      
      
      deallocate(rmmix_glb)
      deallocate(std_mmix_glb)
      nullify(rmmix)
      nullify(std_mmix)
      nullify(tempm)
      nullify(presm)
#else
      write (gol,'("not compiled with hdf4")'); call goErr
      TRACEBACK; rcode=1; return
#endif // with_hdf4

      ! ok
      call goLabel(); rcode=0

    END SUBROUTINE WRITEMMIX_REGION
    !EOC

  END SUBROUTINE WRITE_MMIX
  !EOC

END MODULE USER_OUTPUT_MMIX