!#################################################################
!
! handles initialisation, accumulation and output of mean mixing ratio data
!
!### 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"
!
!#################################################################

module user_output_mmix

  use GO, only : gol, goPr, goErr, goBug, goLabel
  use dims, mname_dims=>mname
  use chem_param, only: ntrace, nstd

  implicit none

  private

  public :: mmix_Init
  public :: mmix_Done
  public :: accumulate_mmix 
  public :: reset_mmix
  public :: write_mmix
  public :: mmix_data, mmix_dat, w_mmix

  type 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

  real,dimension(nregions)                     :: w_mmix

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

contains



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


  subroutine mmix_Init(status)

    implicit none

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

    !__LOCAL_VARIABLES_______________________________________________________

    integer :: region

    ! --- const ------------------------------

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

    !__START_SUBROUTINE______________________________________________________

    call goLabel(rname)

    do region=1,nregions
       allocate(mmix_dat(region)%rmmix   (im(region),jm(region),lm(region),ntrace))
       allocate(mmix_dat(region)%std_mmix(im(region),jm(region),lm(region),nstd))
       allocate(mmix_dat(region)%tempm   (im(region),jm(region),lm(region)))
       allocate(mmix_dat(region)%presm   (im(region),jm(region)))
       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

    call goLabel()
    ! ok
    call goLabel(); status=0

  end subroutine mmix_Init

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


  subroutine mmix_Done(status)
    implicit none

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

    !__LOCAL_VARIABLES_______________________________________________________

    integer :: region

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

    !__START_SUBROUTINE______________________________________________________

    call goLabel( rname )


    do region = 1,nregions
       deallocate ( mmix_dat(region)%rmmix   )
       deallocate ( mmix_dat(region)%std_mmix)
       deallocate ( mmix_dat(region)%tempm   )  !WP! over global domain
       deallocate ( mmix_dat(region)%presm   )
    end do
    ! ok
    call goLabel(); status=0
  end subroutine mmix_Done


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


  subroutine reset_mmix(status)

    implicit none

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


    !__LOCAL_VARIABLES_______________________________________________________

    integer :: region

    ! --- const ------------------------------

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

    !__START_SUBROUTINE______________________________________________________

    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
    ! ok
    call goLabel(); status=0

  end subroutine reset_mmix

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

  subroutine accumulate_mmix(region,status)
    !
    ! MK JUL2003
    ! in calling this routine the data are either paralel 
    ! over tracers or paralel over levels
    !
    use dims,        only: isr,ier,jsr,jer,lm,okdebug,lm, ndyn, ndyn_max
    use chem_param,   only: fscale, istd , nstd, ntrace, ntracet
    use global_data, only: mass_dat, region_dat
    use meteo      , only : sp_dat, temper_dat
    use meteodata  , only : m_dat
#ifdef MPI  
    use mpi_const
#endif
    use partools   , only : offsetl_k => offsetl

    implicit none

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

    ! local
    real,dimension(:,:,:,:),pointer           :: rm, rmmix, std_mmix
    real,dimension(:,:,:),  pointer           :: t, tempm, m
    real,dimension(:,:,:),  pointer           :: p
    real,dimension(:,:),    pointer           :: presm
    integer,dimension(:,:),    pointer        :: zoomed
    integer         :: i,j,l,n, itrace,lmr,nt
    integer   :: communicator,offsetl,offsetn,lglob,nglob,root_id
    real            :: weight

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

    ! start
    call goLabel( rname )


    call accumulate_mmix_short(region,status)   ! gather short lived species
    IF_NOTOK_RETURN(status=1)

#ifdef MPI
    which_par=previous_par(region)

    if(which_par == 'tracer'.and.ntracetloc == 0) return 
    if(which_par == 'levels'.and.lmloc == 0) return  !WP! 
#endif

#ifdef with_zoom
    zoomed => region_dat(region)%zoomed
#endif

    m => m_dat(region)%data

#ifdef MPI
    if(which_par == 'tracer') then
#endif

       rm => mass_dat(region)%rm_t
       lmr = lm(region)
       offsetl = 0
#ifdef MPI  
       nt=ntracetloc
       communicator=com_trac  !WP! assign com_trac as communicator
       root_id=root_t
       offsetn=sum(ntracet_ar(0:myid-1) )  !offset for global value of n
#else
       nt=ntracet
       offsetn=0
#endif
#ifdef MPI
    else if ( which_par == 'levels' ) then

       rm => mass_dat(region)%rm_k
       lmr = lmloc
       offsetl = offsetl_k
       nt=ntracet
       communicator=com_lev  !WP! assign com_lev as communicator
       root_id=root_k
       offsetn=0    ! no offset for tracers

    end if
#endif

    rmmix => mmix_dat(region)%rmmix
    std_mmix => mmix_dat(region)%std_mmix
    tempm => mmix_dat(region)%tempm
    t => temper_dat(region)%data
    p => sp_dat(region)%data
    presm => mmix_dat(region)%presm

    weight = float(ndyn)/float(ndyn_max)
    do n=1,nt
       nglob=n+offsetn  !WP! offset is zero on level domain
       do l=1,lmr
          lglob=l+offsetl  !WP! offset is zero on tracer domain
          do j = jsr(region),jer(region)
             do i = isr(region), ier(region)
#ifdef with_zoom
                if(zoomed(i,j) /= region) cycle
#endif
                rmmix(i,j,lglob,nglob) = rmmix(i,j,lglob,nglob) + &
                     weight*rm(i,j,l,n)/m(i,j,lglob)
             end do
          end do
       end do
    end do

    do n=1,nstd
       itrace = istd(n)
       !CMKbug   if(which_par == 'levels') itrace=itrace-offsetn 

#ifdef MPI  
       !WP! calculate local itrace  
       if ( which_par == 'tracer' ) itrace=itrace-offsetn
       !WP! only when itrace is meaningful on this PE
       if ( itrace <= 0 .or. itrace > nt ) cycle
#endif
       do l=1,lmr
          lglob=l+offsetl
          do j = jsr(region),jer(region)
             do i = isr(region), ier(region)
#ifdef with_zoom 
                if(zoomed(i,j) /= region) cycle
#endif
                std_mmix(i,j,lglob,n) = std_mmix(i,j,lglob,n) + &
                     weight*(rm(i,j,l,itrace)/m(i,j,lglob))**2
            end do
          end do
       end do
    end do

    do l=1,lm(region)  !WP! over global domain
       do j = jsr(region),jer(region)
          do i = isr(region), ier(region)
#ifdef with_zoom
             if(zoomed(i,j) /= region) cycle
#endif
             tempm(i,j,l) = tempm(i,j,l) + weight*t(i,j,l)
          end do
       end do
    end do
    do j = jsr(region),jer(region)
       do i = isr(region), ier(region)
#ifdef with_zoom
          if(zoomed(i,j) /= region) cycle
#endif
          presm(i,j) = presm(i,j) + weight*p(i,j,1)
       end do
    end do


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

#ifdef with_zoom
    nullify(zoomed)
#endif
    nullify(m)
    nullify(rm)
    nullify(t)
    nullify(p)
    nullify(rmmix)
    nullify(std_mmix)
    nullify(presm)
    nullify(tempm)

    ! ok
    call goLabel(); status=0

  end subroutine accumulate_mmix



  subroutine accumulate_mmix_short(region,status)
    !
    ! The short lived compounds are always paralel over levels
    !
    use dims,          only : isr,ier,jsr,jer,lm,okdebug,lm, ndyn, ndyn_max
    use global_data,   only : region_dat, mass_dat
    use meteodata    , only : m_dat
    use tracer_data  , only : chem_dat
    use chem_param,     only : fscale, istd , nstd, ntrace, ntracet
#ifdef MPI  
    use mpi_const
#endif
    use partools, only : lmloc, offsetl
    use partools, only : which_par

    implicit none

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

    ! local
    real,dimension(:,:,:,:),pointer :: rmmix, std_mmix
    real,dimension(:,:,:,:),pointer :: rmc
    real,dimension(:,:,:),pointer   :: m
    integer,dimension(:,:),pointer  :: zoomed
    integer                         :: n, i, j, l, lmr, lmm, lglob
    real                            :: weight

    character(len=*), parameter ::  rname = mname//'/accumulate_mmix_short'
    ! start
    call goLabel( rname )


#ifdef with_zoom
    zoomed => region_dat(region)%zoomed
#endif
    rmc => chem_dat(region)%rm_k
    rmmix => mmix_dat(region)%rmmix
    std_mmix => mmix_dat(region)%std_mmix
    lmr = lmloc

    m => m_dat(region)%data
    lmm = offsetl
    weight = float(ndyn)/float(ndyn_max)
    do n= ntracet+1, ntrace
       do l=1,lmr
          lglob=l+offsetl  !WP! offset is zero on tracer domain
          do j = jsr(region),jer(region)
             do i = isr(region), ier(region)
#ifdef with_zoom
                if ( zoomed(i,j) /= region ) cycle
#endif
                rmmix(i,j,lglob,n) = rmmix(i,j,lglob,n) + &
                     weight*rmc(i,j,l,n)/m(i,j,l+lmm)
             end do
          end do
       end do
    end do

#ifdef with_zoom
    nullify(zoomed)
#endif
    nullify(rmc)
    nullify(rmmix)
    nullify(std_mmix)
    nullify(m)

    ! ok
    call goLabel(); status=0

  end subroutine accumulate_mmix_short


  subroutine update_mmix_parent(region, status)
    !
    ! pass the mmix information  from child to parent
    ! 

    use global_data,  only : region_dat
    use chem_param,    only : ntrace, istd, nstd, ntracet
    implicit none

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

    ! local
    real,dimension(:,:,:,:),pointer  :: rmmix, rmmixp, std_mmix, std_mmixp
    real,dimension(:,:,:),pointer    :: tempm,tempmp
    real,dimension(:,:),pointer      :: presm,presmp
    real,dimension(:),pointer        :: dxyp

    integer :: my_parent,xref_,yref_,zref_
    integer :: ip,jp,ic,jc,i,j,l,n,n1
    integer :: imp,jmp,lmp,imr,jmr,lmr
    integer :: iox,ioy1,ioy2
    real    :: w,wtot

    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/update_mmix_parent'
    
    ! start
    call goLabel( rname )


    if ( region == 1 ) return

    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    my_parent = parent(region)
    xref_ = xref(region)/xref(my_parent)
    yref_ = yref(region)/yref(my_parent)
    zref_ = zref(region)/zref(my_parent)

    if ( okdebug ) then
       print *,'update_mmix_parent: my_parent=',my_parent, &
            ' x-,y-,zref_: ',xref_,yref_,zref_
    end if

    imp = im(region)/xref_
    jmp = jm(region)/yref_
    lmp = lm(region)/zref_

    if ( ibeg(region) < iend(region) .and. &
         imp /= iend(region)-ibeg(region)+1 ) then
           write (gol,'("Grid setup error")'); call goErr
           write (gol,'("in program")'); call goErr; status=1; return
    end if

    if ( jmp /= jend(region)-jbeg(region)+1 ) then
           write (gol,'("Grid setup error")'); call goErr
           write (gol,'("in program")'); call goErr; status=1; return
    end if

    if ( lmp /= lend(region)-lbeg(region)+1 ) then
           write (gol,'("Grid setup error")'); call goErr
           write (gol,'("in program")'); call goErr; status=1; return
    end if

    rmmix => mmix_dat(region)%rmmix
    std_mmix => mmix_dat(region)%std_mmix
    presm => mmix_dat(region)%presm
    tempm => mmix_dat(region)%tempm
    dxyp => region_dat(region)%dxyp

    rmmixp => mmix_dat(my_parent)%rmmix
    std_mmixp => mmix_dat(my_parent)%std_mmix
    presmp => mmix_dat(my_parent)%presm
    tempmp => mmix_dat(my_parent)%tempm

    iox = isr(region)/xref_
    ioy1 = jsr(region)/yref_
    ioy2 = (jm(region)-jer(region)+1)/yref_
    write(gol, *) 'update_mmix_parent: iox,ioy1,ioy2 ', iox,ioy1,ioy2, &
         'parent..,child', my_parent, region; call goPr

    do l=1,lm(region)
       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_
             tempmp(ip,jp,l) = 0.0
             wtot = 0.0
             do j=1,yref_
                w = dxyp(jc+j)
                do i=1,xref_
                   tempmp(ip,jp,l) = tempmp(ip,jp,l)  + tempm(ic+i,jc+j,l)*w
                   wtot = wtot+w
                end do
             end do
             tempmp(ip,jp,l) = tempmp(ip,jp,l)/wtot
          end do
       end do
    end do

    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_
          presmp(ip,jp) = 0.0
          wtot = 0.0
          do j=1,yref_
             w = dxyp(jc+j)
             do i=1,xref_
                presmp(ip,jp) = presmp(ip,jp)  + presm(ic+i,jc+j)*w
                wtot = wtot+w
             end do
          end do
          presmp(ip,jp) = presmp(ip,jp)/wtot
       end do
    end do

    do n=1,ntrace
      do l=1,lm(region)
        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_
            rmmixp(ip,jp,l,n) = 0.0
            wtot = 0.0
            do j=1,yref_
              w = dxyp(jc+j)
              do i=1,xref_
                rmmixp(ip,jp,l,n) = rmmixp(ip,jp,l,n) + &
                       rmmix(ic+i,jc+j,l,n)*w
                wtot = wtot+w
              end do
            end do
            rmmixp(ip,jp,l,n) = rmmixp(ip,jp,l,n)/wtot
         end do
       end do
     end do
    end do

    do n=1,nstd
      do l=1,lm(region)
        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_
            std_mmixp(ip,jp,l,n) = 0.0
            wtot = 0.0
            do j=1,yref_
              w = dxyp(jc+j)
              do i=1,xref_
                std_mmixp(ip,jp,l,n) = std_mmixp(ip,jp,l,n)  + &
                      std_mmix(ic+i,jc+j,l,n)*w
                wtot = wtot+w
              end do
            end do
            std_mmixp(ip,jp,l,n) = std_mmixp(ip,jp,l,n)/wtot
          end do
        end do
      end do
    end do

    nullify(rmmix)
    nullify(std_mmix)
    nullify(presm)
    nullify(tempm)
    nullify(dxyp)

    nullify(rmmixp)
    nullify(std_mmixp)
    nullify(presmp)
    nullify(tempmp)

    ! ok
    call goLabel(); status=0

  end subroutine update_mmix_parent


  subroutine write_mmix(status)
    !------------------------------------------------------------------
    ! save all essential model parameters and fields on unit kdisk
    !
    !                                       v 8.4
    !  modified for HDF output
    !------------------------------------------------------------------

    use chem_param
    use io_hdf
    use datetime,     only: tstamp
    use global_data, only : outdir

    use ParTools, only : myid, root_t, Par_Barrier
    use User_Output_Common, only : User_Output_Check_Overwrite

    implicit none

    ! --- in/out -----------------------------
    
    integer, intent(inout)      ::  status


    !__LOCAL_VARIABLES_______________________________________________________

    integer           :: istat, sfsnatt, sfscatt, io, sfstart
    integer           :: sfend
    integer           :: region
    integer           :: i,j,l,n,k,ind
    character(len=12) :: name

    character(len=200)   :: FFilename

    ! --- const ------------------------------

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

    !__START_SUBROUTINE______________________________________________________

    call goLabel( rname )


    do region=nregions,1,-1

       if ( myid == root_t ) then
         ! filename:
         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))
         ! check existence ...
         call User_Output_Check_Overwrite( trim(FFilename), status )
         IF_NOTOK_RETURN(status=1)
         ! oepn:
         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; status=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, ntracet)
         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*8, 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 !all pe's from here
       call writemmix_region(region, status)
       IF_NOTOK_RETURN(status=1)

       call par_barrier

       if ( myid == root_t ) then
          write(gol,'("write_mmix: sfend returns",i4)') sfend(io) ; call goPr
       endif

    end do



    call tstamp(kmain,itau,'savemmix')

    if ( cdebug ) then
       call tstamp(kdebug,itau,'savemmix')
    end if
    call par_barrier

    ! ok
    call goLabel(); status=0


  contains

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

    subroutine writemmix_region(region, status)

#ifdef MPI 
      use mpi_const, only : my_real, mpi_sum, com_trac, ierr
#endif
      use global_data, only: region_dat

      implicit none

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

      ! local
      real,dimension(:,:,:,:),pointer        :: rmmix, std_mmix
      real,dimension(:,:,:)  ,pointer        :: tempm
      real,dimension(:,:)    ,pointer        :: presm
#ifdef with_zoom
      integer,dimension(:,:) ,pointer        :: zoomed
#endif
      real,dimension(:,:,:,:),allocatable    :: rmmix_sum
      real,dimension(:,:,:,:),allocatable    :: std_mmix_sum
      integer imr,jmr,lmr,nsend
      real    :: ahelp,ahelp1

     ! --- const ------------------------------
 
     character(len=*), parameter ::  rname = mname//'/write_mmix_region'
    

      ! start
    call goLabel( rname )


      imr = im(region) ; jmr = jm(region) ; lmr = lm(region)
      allocate(rmmix_sum(imr,jmr,lmr,ntrace))
      allocate(std_mmix_sum(imr,jmr,lmr,nstd))

      rmmix => mmix_dat(region)%rmmix
      std_mmix => mmix_dat(region)%std_mmix
      tempm => mmix_dat(region)%tempm
      presm => mmix_dat(region)%presm
#ifdef with_zoom
      zoomed => region_dat(region)%zoomed
#endif

      write(gol,'("writemmix_region: w_mmix",f12.2)') w_mmix(region); call goPr
      if ( w_mmix(region) > 2 ) then
#ifdef MPI   
         nsend=imr*jmr*lmr*ntrace
         call mpi_allreduce( rmmix, rmmix_sum, nsend, &
              my_real, mpi_sum, com_trac, ierr )
         nsend=imr*jmr*lmr*nstd
         call mpi_allreduce( std_mmix, std_mmix_sum, nsend, &
              my_real, mpi_sum, com_trac, ierr )
#else
         rmmix_sum = rmmix
         std_mmix_sum = std_mmix
#endif
         do k=1,nstd
            n = istd(k)
            do l=1,lmr
               do j=jsr(region),jer(region)
                  do i=isr(region),ier(region)
#ifdef with_zoom
                     if(zoomed(i,j)/=region) cycle
#endif
                     ahelp = rmmix_sum(i,j,l,n)
                     ahelp1= fscale(n)* (std_mmix_sum(i,j,l,k) - &
                          ahelp*ahelp/w_mmix(region))/(w_mmix(region)-1)
                     std_mmix_sum(i,j,l,k)=max(1e-35,ahelp1)
                  end do
               end do
            end do
         end do
         do n=1,ntrace
            do l=1,lmr
               do j=jsr(region),jer(region)
                  do i=isr(region),ier(region)
#ifdef with_zoom
                     if(zoomed(i,j)/=region) cycle
#endif
                     rmmix_sum(i,j,l,n) = fscale(n)* &
                          rmmix_sum(i,j,l,n)/w_mmix(region)
                  end do
               end do
            end do
         end do
         do l=1,lmr
            do j=jsr(region),jer(region)
               do i=isr(region),ier(region)
#ifdef with_zoom
                  if(zoomed(i,j)/=region) cycle
#endif
                  tempm(i,j,l) = tempm(i,j,l)/w_mmix(region)
               end do
            end do
         end do
         do j=jsr(region),jer(region)
            do i=isr(region),ier(region)
#ifdef with_zoom
               if(zoomed(i,j)/=region) cycle
#endif
               presm(i,j) = presm(i,j)/w_mmix(region)
            end do
         end do

         if ( myid == root_t ) then
            rmmix = rmmix_sum      ! CORRECTED OCT2003 CMK
            ! pass the information of the core zoom to the parent
            ! (which has no info yet!) only root!
            call update_mmix_parent(region, status)
            IF_NOTOK_RETURN(status=1)
            call io_write3d_32d(io,imr,'LON'//trim(region_name(region)), &
                 jmr,'LAT'//trim(region_name(region)),lmr, &
                 'HYBRID',tempm(1:imr,1:jmr,1:lmr),'tempm',idate)
            call io_write2d_32d(io,imr,'LON'//trim(region_name(region)), &
                 jmr,'LAT'//trim(region_name(region)), &
                 presm(1:imr,1:jmr),'presm',idate)
            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_sum(1:imr,1:jmr,1:lmr,n),name,idate)
            end do

            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_sum(1:imr,1:jmr,1:lmr,k),name,idate)
            end do

         end if ! all PEs from here

      end if ! nmmix>2

      deallocate(rmmix_sum)
      deallocate(std_mmix_sum)
      nullify(rmmix)
      nullify(std_mmix)
      nullify(tempm)
      nullify(presm)
#ifdef with_zoom
      nullify(zoomed)
#endif

    ! ok
    call goLabel(); status=0

    end subroutine writemmix_region

  end subroutine write_mmix

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



end module user_output_mmix
