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

!#################################################################
!
! MODULE:  user_output_mix  *************** DUMMY TM5-MP VERSION *******************
! 
! PUBLIC MEMBER FUNCTIONS:
!    user_output_mix_init
!    user_output_mix_write
!    user_output_mix_accum
!    user_output_mix_done
!
! DESCRIPTION: 
!   Write mixing ratio tracer fields. Called from user_output.F90. 
!
! REVISION HISTORY: 
!   Wouter Peters
!
!   Mike Trudeau, Oct 2011 
!   Modified to write NetCDF4 output in NOAA/ESRL CarbonTracker
!   Release format.
!
!   Mike trudeau, Feb 2012
!   Modified to output temperature, orography, gph at model level
!   boundaries, and remove "halo" grid cells.
!
!   Andy Jacobson, 13 Jun 2012
!    Add date dimension to temperature variable in output file
!    Change nc variable names: temp to temperature, press to pressure.
!
!   Andy Jacobson, 17 Jun 2012
!    Dymamic time step weighting in averages.
!
!   Andy Jacobson, 3 Jul 2013
!    NetCDF file attributes read from rc file
!
!### macros #####################################################
!
#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_mix

!   use dims,                                   only : nregions
!   use GO,                                     only : gol, goErr, goPr

  implicit none

  ! --- in/out ------------------------

  private

  public :: user_output_mix_init, user_output_mix_done, user_output_mix_write, user_output_mix_accum
  public :: mix_netcdf_attributes

  ! --- const ----------------------------------
    
  character(len=*), parameter                 :: mname = 'user_output_mix'

  ! --- var ---------------------------

!   type fields 
!     real, dimension(:,:,:,:), pointer         :: field ! 4D fields to write
!     real, dimension(:,:), pointer             :: ps ! the surface pressure
!     real, dimension(:,:), pointer             :: orography ! ground height
!     real, dimension(:,:,:), pointer           :: temperature ! air temperature
!     real, dimension(:,:,:), pointer           :: gph ! geopotential height
!     real, dimension(:,:), pointer             :: wda ! number of averages
!     character(len=20), dimension(:), pointer  :: namefield
!     character(len=1024)                       :: fname
!     integer, dimension(6)                     :: idate_start
!     integer, dimension(6)                     :: idate_end
!     integer                                   :: funit = 0 
!     real                                      :: lon_first
!     real                                      :: lon_inc
!     real                                      :: lat_first
!     real                                      :: lat_inc 
!     integer                                   :: lon_n 
!     integer                                   :: lat_n
!   end type fields 

!   type(fields), dimension(nregions), save     :: ncfile

  type netcdf_attribute_t
     character(len=256)		              :: notes 
     character(len=1024)	              :: disclaimer
     character(len=256) 	              :: email
     character(len=256) 	              :: url
     character(len=256) 	              :: institution
     character(len=256) 	              :: conventions
     character(len=256) 	              :: source 
  end type netcdf_attribute_t

   type(netcdf_attribute_t)		      :: mix_netcdf_attributes

contains

  subroutine user_output_mix_init(status)

!     use dims,        only : nregions
!     use dims,        only : im, jm, lm
!     use dims,        only : xbeg, ybeg, dx, dy, xref, yref
!     use dims,        only : idatei
!     use chem_param,  only : names

!     use partools,    only : myid, ntracetloc, ntracet_ar

!     implicit none

!     ! --- in/out ------------------------

     integer, intent(out)  :: status
!     integer               :: region, offsetn
!     integer               :: i, j, n 

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

!     offsetn = sum(ntracet_ar(0:myid-1))
    
!     ! initialize the mix output
!     regionloop: do region = 1, nregions

!       allocate(ncfile(region)%field(im(region),jm(region),lm(region),ntracetloc))
!       allocate(ncfile(region)%ps(im(region),jm(region)))
!       allocate(ncfile(region)%orography(im(region),jm(region)))
!       allocate(ncfile(region)%temperature(im(region),jm(region),lm(region)))
!       allocate(ncfile(region)%gph(im(region),jm(region),lm(region)+1))
!       allocate(ncfile(region)%wda(im(region),jm(region)))
!       allocate(ncfile(region)%namefield(ntracetloc))

!       do i=1, ntracetloc
!         ncfile(region)%namefield(i) =  names(i+offsetn)
!       enddo

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

!       ncfile(region)%idate_start = idatei

!     end do  regionloop

    ! ok
    status = 0

  end subroutine user_output_mix_init

  subroutine user_output_mix_done(status)

!     use dims,             only : nregions
 
!     implicit none

!     ! --- in/out ------------------------

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

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

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

!     ! do not yet deallocate because parent mix is needed for updates
! !    do region = 1, nregions
!     do region = nregions, 1, -1  ! ARJ, 17 May 2012
!        call user_output_mix_write( region, status )
!        IF_NOTOK_RETURN(status=1)
!     end do

!     do region = 1, nregions
!       deallocate(ncfile(region)%field) 
!       deallocate(ncfile(region)%ps)
!       deallocate(ncfile(region)%wda)
!       deallocate(ncfile(region)%namefield)
!       deallocate(ncfile(region)%orography)
!       deallocate(ncfile(region)%temperature)
!       deallocate(ncfile(region)%gph)
!     enddo 

    ! ok
    status = 0

  end subroutine user_output_mix_done

  subroutine user_output_mix_accum(region, status)

!     use global_data, only : mass_dat, region_dat
!     use meteo,       only : sp_dat, gph_dat, m_dat, oro_dat, temper_dat
!     use dims,        only : im, jm, lm
!     use dims,        only : isr, jsr, ier, jer
!     use dims,        only : tref, ndyn, ndyn_max
!     use binas,       only : xmair
!     use ParTools,    only : myid, ntracetloc, ntracet_ar
!     use chem_param,  only : ra,uscale

!     implicit none

!     ! --- in/out ------------------------

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

!     integer                           :: i ,j, l, n
!     integer                           :: imr, jmr, lmr
!     integer                           :: dtime, offsetn, dtime_max
!     real                              :: 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
!     real, dimension(:,:), pointer     :: orography 
!     real, dimension(:,:,:), pointer   :: temperature 

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

!     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
!     orography => oro_dat(region)%data(:,:,1)
!     temperature => temper_dat(region)%data

!     ! 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
!           ncfile(region)%wda(i,j) = ncfile(region)%wda(i,j) + weight

!           do l = 1, lm(region)
!             do n = 1, ntracetloc
!               ncfile(region)%field(i,j,l,n) = &
!                 ncfile(region)%field(i,j,l,n) + uscale(n)*weight*rm(i,j,l,n)/m(i,j,l)*(xmair/ra(offsetn+n))
!             enddo

!             ncfile(region)%temperature(i,j,l) = &
!               ncfile(region)%temperature(i,j,l) + weight*temperature(i,j,l)
!           enddo

!           do l = 1, lm(region) + 1
!             ncfile(region)%gph(i,j,l) = &
!               ncfile(region)%gph(i,j,l) + weight*gph(i,j,l)
!           enddo

!           ncfile(region)%ps(i,j) = &
!             ncfile(region)%ps(i,j) + weight*p(i,j,1)

!           ncfile(region)%orography(i,j) = &
!             ncfile(region)%orography(i,j) + weight*orography(i,j)
!       enddo
!     enddo

!     nullify(rm) 
!     nullify(m)
!     nullify(p)
!     nullify(gph)
!     nullify(dxyp)
!     nullify(zoomed)
!     nullify(orography)
!     nullify(temperature)
!     nullify(gph)

    ! ok
    status = 0

  end subroutine user_output_mix_accum

!   subroutine output_update_parent(region)

!     use dims
!     use global_data,  only: region_dat
!     use toolbox,      only: escape_tm

!     implicit none

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

!     ! --- var ---------------------------

!     real, dimension(:,:,:,:), pointer  :: mean_mix ! column value m - mean
!     real, dimension(:,:,:), pointer      :: mean_gph ! column value gph - mean
!     real, dimension(:,:), pointer        :: mean_orography ! column value orography - mean
!     real, dimension(:,:,:), pointer      :: mean_temperature ! column value temperature - 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_gph ! column value gph - mean
!     real, dimension(:,:), pointer        :: p_mean_orography ! column value orography - mean
!     real, dimension(:,:,:), pointer      :: p_mean_temperature ! column value temperature - 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_
!     integer                              :: imp, jmp
!     integer                              :: imr, jmr
!     integer                              :: ic, jc
!     integer                              :: iox, ioy1, ioy2
!     real                                 :: w, wtot

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

!     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 => ncfile(region)%field(:,:,:,:)
!     mean_gph => ncfile(region)%gph
!     mean_orography => ncfile(region)%orography
!     mean_temperature => ncfile(region)%temperature
!     mean_ps => ncfile(region)%ps

!     p_mean_mix  => ncfile(my_parent)%field(:,:,:,:)
!     p_mean_gph  => ncfile(my_parent)%gph
!     p_mean_orography  => ncfile(my_parent)%orography
!     p_mean_temperature  => ncfile(my_parent)%temperature
!     p_mean_ps => ncfile(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
!         p_mean_gph(ip,jp,:) = 0.0
!         p_mean_orography(ip,jp) = 0.0
!         p_mean_temperature(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_gph(ip,jp,:) = p_mean_gph(ip,jp,:)  + mean_gph(ic+i,jc+j,:)*w
!             p_mean_orography(ip,jp) = p_mean_orography(ip,jp)  + mean_orography(ic+i,jc+j)*w
!             p_mean_temperature(ip,jp,:) = p_mean_temperature(ip,jp,:)  + mean_temperature(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_gph(ip,jp,:) = p_mean_gph(ip,jp,:)/wtot
!         p_mean_orography(ip,jp) = p_mean_orography(ip,jp)/wtot
!         p_mean_temperature(ip,jp,:) = p_mean_temperature(ip,jp,:)/wtot
!         p_mean_ps(ip,jp) = p_mean_ps(ip,jp)/wtot

!       enddo
!     enddo

!     nullify(mean_mix)  
!     nullify(mean_gph)  
!     nullify(mean_orography)  
!     nullify(mean_temperature)  
!     nullify(mean_ps)  

!     nullify(p_mean_mix)  
!     nullify(p_mean_gph)  
!     nullify(p_mean_orography)  
!     nullify(p_mean_temperature)  
!     nullify(p_mean_ps)  

!   end subroutine output_update_parent

  subroutine user_output_mix_write(region, status)

!     use global_data, only : region_dat
!     use dims,        only : im, jm, lm
!     use dims,        only : itaur, tref, at, bt, ndyn, ndyn_max
!     use dims,        only : isr, jsr, ier, jer
!     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 dims,        only : region_name
!     use chem_param,  only : ntracet, names
!     use datetime,    only : tau2date, date2tau, idate2ddate
!     use ParTools,    only : myid, root, ntracetloc
!     use global_data, only : outdir

!     use MDF,         only : MDF_Create, MDF_Close, MDF_EndDef
!     use MDF,         only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_UNLIMITED
!     use MDF,         only : MDF_INT, MDF_FLOAT, MDF_DOUBLE, MDF_CHAR
!     use MDF,         only : MDF_Put_Att
!     use MDF,         only : MDF_Def_Dim
!     use MDF,         only : MDF_Def_Var, MDF_Put_Var

!     implicit none

!     ! --- in/out ------------------------

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

!     integer                          :: dimid_lon, dimid_lat, dimid_lvl, dimid_bnd, dimid_date, dimid_cal, dimid_ntracetloc, dimid_char8
!     integer                          :: varid_lon, varid_lat, varid_lvl, varid_bnd
!     integer                          :: varid_press
!     integer			     :: varid_tracer,varid_tracernm
!     integer                          :: varid_elapsed, varid_dec, varid_cal_int, varid_date_int
!     integer                          :: varid_gph
!     integer                          :: varid_orography, varid_temperature
!     integer                          :: fid

!     integer                          :: i, j, l
!     integer                          :: imr, jmr, lmr
!     integer                          :: n
!     integer, dimension(6)            :: idatee
!     integer                          :: dtime, dtime_max

!     integer                          :: itau_start
!     integer                          :: itau_end
!     integer                          :: itau_avg
!     integer                          :: itau_ref
!     integer                          :: itau_elapsed
!     real*8                           :: ddate_avg 
!     real*8                           :: days_elapsed 
!     integer, dimension(6)            :: idate_ref
!     integer, dimension(6)            :: idate_avg 
!     real*4                           :: fillval_r4 = -1e34
!     real*8                           :: fillval_r8 = -1e34

!     integer, dimension(8)            :: isysdate

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

!     integer                          :: iHalo, imHalo, jHalo, jmHalo

!     type ncType
!       real, dimension(:,:), allocatable     :: orography 
!       real, dimension(:,:,:), allocatable   :: temperature
!       real, dimension(:,:,:), allocatable   :: pressure
!       real, dimension(:,:,:), allocatable   :: gph 
!       real, dimension(:,:,:,:), allocatable :: field
!       real, dimension(:), allocatable       :: lonMean 
!       real, dimension(:), allocatable       :: latMean
!       real, dimension(:), allocatable       :: lonEdge 
!       real, dimension(:), allocatable       :: latEdge
!     end type ncType
    
!     type(ncType), dimension(nregions) :: ncout
!     character(len=256) 		              :: progstring
!     character(len=256) 		              :: sysdate
!     character(len=256) 		              :: history
    
!     ! --- const ---------------------------

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

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

!     if (myid .ne. root ) 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), idatee)

!     ncfile(region)%idate_end = idatee

!     ! 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 (ncfile(region)%wda(i,j) > 0) then 

!           ! surface pressure and orography
!           ncfile(region)%ps(i,j) = ncfile(region)%ps(i,j)/ncfile(region)%wda(i,j)
!           ncfile(region)%orography(i,j) = ncfile(region)%orography(i,j)/ncfile(region)%wda(i,j)

!           ! tracer fields and air temperature
!           do l = 1, lm(region)
!             do n = 1, ntracetloc
!               ncfile(region)%field(i,j,l,n) = ncfile(region)%field(i,j,l,n) / ncfile(region)%wda(i,j)
!             enddo

!             ncfile(region)%temperature(i,j,l) = ncfile(region)%temperature(i,j,l) / ncfile(region)%wda(i,j)
!           enddo

!           ! geopotential height
!           do l = 1, lm(region) + 1
!             ncfile(region)%gph(i,j,l) = ncfile(region)%gph(i,j,l)  / ncfile(region)%wda(i,j)
!           enddo
!         endif
!       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 output_update_parent(region)

!     call date_and_time(values = isysdate)
!     write (sysdate, '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",i2.2," UTC")') &
!            isysdate(1), isysdate(2), isysdate(3), isysdate(5), isysdate(6), isysdate(7)


!     call getarg (0, progstring, status)
!     write(history,'("Created ",a," by ",a,".")') trim(sysdate),trim(progstring)

!     ! date/time conversions
!     call date2tau(ncfile(region)%idate_start, itau_start)
!     call date2tau(ncfile(region)%idate_end, itau_end)

!     itau_avg = nint(5.0D-1 * dble(itau_start + itau_end))

!     call tau2date(itau_avg, idate_avg)
!     ddate_avg = idate2ddate(idate_avg)

!     idate_ref = (/2000, 1, 1, 0, 0, 0/)

!     call date2tau(idate_ref, itau_ref)
!     itau_elapsed = itau_avg - itau_ref
!     days_elapsed = dble(itau_elapsed) / 86400.0D+0

!     ! define index offsets to elimate "halo" cells
!     if ( region == 1 ) then
!       iHalo = 0 
!       jHalo = 0
!     else
!       iHalo = 3 
!       jHalo = 2 
!     endif

!     imHalo = im(region) - iHalo * 2
!     jmHalo = jm(region) - jHalo * 2

!     allocate(ncout(region)%orography(imHalo, jmHalo))
!     allocate(ncout(region)%temperature(imHalo, jmHalo, lm(region)))
!     allocate(ncout(region)%pressure(imHalo, jmHalo, lm(region)+1))
!     allocate(ncout(region)%gph(imHalo, jmHalo, lm(region)+1))
!     allocate(ncout(region)%field(imHalo, jmHalo, lm(region), ntracetloc))
!     allocate(ncout(region)%lonMean(imHalo))
!     allocate(ncout(region)%latMean(jmHalo))
!     allocate(ncout(region)%lonEdge(imHalo+1))
!     allocate(ncout(region)%latEdge(jmHalo+1))

!     ! trim "halo" cells
!     do i = 1, imHalo
!       do j = 1, jmHalo
!         do l = 1, lm(region)
!           do n = 1, ntracetloc
!             ncout(region)%field(i,j,l,n) = ncfile(region)%field(i+iHalo,j+jHalo,l,n)
!           enddo
!           ncout(region)%temperature(i,j,l) = ncfile(region)%temperature(i+iHalo,j+jHalo,l)
!         enddo
!         ncout(region)%orography(i,j) = ncfile(region)%orography(i+iHalo,j+jHalo)
!       enddo
!     enddo

!     ! Subtract the background. This is a CarbonTracker-CO2-specific operation.
!     !    do n = 2, ntracetloc
!     !      ncout(region)%field(:,:,:,n) = ncout(region)%field(:,:,:,n) - ncout(region)%field(:,:,:,1) 
!     !    enddo

!     ! compute 3D pressure field at model level boundaries & trim "halo" cells
!     do i = 1, imHalo
!       do j = 1, jmHalo
!         do l = 1, lm(region) + 1
!           ncout(region)%pressure(i,j,l) = at(l) + bt(l) * ncfile(region)%ps(i+iHalo,j+jHalo)
!           ncout(region)%gph(i,j,l) = ncfile(region)%gph(i+iHalo,j+jHalo,l)
!         enddo
!       enddo
!     enddo

!     ! longitude arrays
!     do i = 1, imHalo + 1 
!       ncout(region)%lonEdge(i) = xbeg(region) + (i + iHalo - 1) * dx / xref(region)
!     enddo
!     do i = 1, imHalo 
!       ncout(region)%lonMean(i) = (ncout(region)%lonEdge(i) + ncout(region)%lonEdge(i+1)) / 2.0
!     enddo

!     ! latitude arrays
!     do j = 1, jmHalo + 1
!       ncout(region)%latEdge(j) = ybeg(region) + (j + jHalo - 1) * dy / yref(region)
!     enddo
!     do j = 1, jmHalo
!       ncout(region)%latMean(j) = (ncout(region)%latEdge(j) + ncout(region)%latEdge(j+1)) / 2.0
!     enddo

!     ! create new file
!     write (ncfile(region)%fname, '(a,"/molefrac_",a,"_",i4.4,3i2.2,"_",i4.4,3i2.2,".nc")') &
!            trim(outdir), trim(region_name(region)), ncfile(region)%idate_start(1:4), ncfile(region)%idate_end(1:4)

!     call MDF_Create( trim(ncfile(region)%fname), MDF_NETCDF, MDF_REPLACE, fid, status )
!     IF_NOTOK_RETURN(status=1)

!     ! global attributes 
!     if(len_trim(mix_netcdf_attributes%notes) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "notes", values=trim(mix_netcdf_attributes%notes), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(mix_netcdf_attributes%disclaimer) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "disclaimer", values=trim(mix_netcdf_attributes%disclaimer), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(mix_netcdf_attributes%email) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "email", values=trim(mix_netcdf_attributes%email), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(mix_netcdf_attributes%url) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "url", values=trim(mix_netcdf_attributes%url), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(mix_netcdf_attributes%institution) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "institution", values=trim(mix_netcdf_attributes%institution), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(mix_netcdf_attributes%conventions) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "conventions", values=trim(mix_netcdf_attributes%conventions), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(history) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "history", values=trim(history), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
!     if(len_trim(mix_netcdf_attributes%source) .gt. 0) then
!        call MDF_Put_Att( fid, MDF_GLOBAL, "source", values=trim(mix_netcdf_attributes%source), status=status )
!        IF_NOTOK_RETURN(status=1)
!     end if
    
!     ! define dimensions
!     call MDF_Def_Dim( fid, 'date', MDF_UNLIMITED, dimid_date, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'char8', 8, dimid_char8, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'ntracers', ntracetloc, dimid_ntracetloc, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'calendar_components', 6, dimid_cal, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'lon', imHalo, dimid_lon, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'lat', jmHalo, dimid_lat, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'level', lm(region) , dimid_lvl, status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Dim( fid, 'boundary', lm(region)+1, dimid_bnd, status )
!     IF_NOTOK_RETURN(status=1)

!     ! dimension variables
!     call MDF_Def_Var( fid, 'date', MDF_DOUBLE, (/dimid_date/), varid_elapsed, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_elapsed, "units", values="days since 2000-01-01 00:00:00 UTC", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_elapsed, "long_name", values="date", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'decimal_date', MDF_DOUBLE, (/dimid_date/), varid_dec, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_dec, "units", values="years", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_dec, "_FillValue", values=fillval_r8, status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'calendar_components', MDF_INT, (/dimid_cal/), varid_cal_int, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_cal_int, "units", values="none", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_cal_int, "long_name", values="calendar_components", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'date_components', MDF_INT, (/dimid_cal, dimid_date/), varid_date_int, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_date_int, "units", values="none", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_date_int, "long_name", values="Integer value calendar components of UTC date.", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_date_int, "comment", values="year, month, day, hour, minute, second", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'lon', MDF_DOUBLE, (/dimid_lon/), varid_lon, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lon, "units", values="degrees_east", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lon, "long_name", values="lon", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lon, "actual_range", values=(/ncout(region)%lonEdge(1), ncout(region)%lonEdge(imHalo+1)/), status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'lat', MDF_DOUBLE, (/dimid_lat/), varid_lat, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lat, "units", values="degrees_north", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lat, "long_name", values="lat", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lat, "actual_range", values=(/ncout(region)%latEdge(1), ncout(region)%latEdge(jmHalo+1)/), status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'level', MDF_INT, (/dimid_lvl/), varid_lvl, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lvl, "units", values="none", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lvl, "long_name", values="level", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lvl, "positive", values="up", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'boundary', MDF_INT, (/dimid_bnd/), varid_bnd, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "units", values="none", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "long_name", values="boundary", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "positive", values="up", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, "tracer", MDF_FLOAT, (/dimid_lon, dimid_lat, dimid_lvl, dimid_ntracetloc,dimid_date/), varid_tracer, deflate_level=9, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "units", values="unspecified", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_tracer, "_FillValue", values=fillval_r4, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_tracer, "long_name", values="tracer_mixing_ratio", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, "tracer_names", MDF_CHAR, (/dimid_char8, dimid_ntracetloc/), varid_tracernm, status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'pressure', MDF_FLOAT, (/dimid_lon, dimid_lat, dimid_bnd, dimid_date/), varid_press, deflate_level=9, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_press, "units", values="Pa", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_press, "_FillValue", values=fillval_r4, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_press, "long_name", values="air_pressure", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_press, "comment", values="air pressure at level boundaries", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'gph', MDF_FLOAT, (/dimid_lon, dimid_lat, dimid_bnd, dimid_date/), varid_gph, deflate_level=9, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_gph, "units", values="m", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_gph, "_FillValue", values=fillval_r4, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_gph, "long_name", values="geopotential_height", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_gph, "comment", values="geopotential height at level boundaries", status=status )
!     IF_NOTOK_RETURN(status=1)

!     ! note: Reported as geopotential for consistency with previous CT releases. Divide by 9.80665 for surface height in meters
!     ! source: http://www.ecmwf.int/products/data/archive/data_faq.html#geopotential
!     ! 21 Feb 2012, M. Trudeau
!     call MDF_Def_Var( fid, 'oro', MDF_FLOAT, (/dimid_lon, dimid_lat/), varid_orography, deflate_level=9, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_orography, "units", values="m^2/s^2", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_orography, "long_name", values="orography", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_orography, "standard_name", values="surface geopotential", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'temperature', MDF_FLOAT, (/dimid_lon, dimid_lat, dimid_lvl, dimid_date/), varid_temperature, deflate_level=9, status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_temperature, "units", values="Kelvin", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_temperature, "long_name", values="air_temperature", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_temperature, "standard_name", values="air temperature at level center", status=status )
!     IF_NOTOK_RETURN(status=1)

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

!     ! write variables
!     call MDF_Put_Var( fid, varid_elapsed, (/(days_elapsed)/), status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_dec, (/(ddate_avg)/), status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_cal_int, (/(i,i=1,size(idate_avg))/), status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_date_int, idate_avg, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_lon, ncout(region)%lonMean, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_lat, ncout(region)%latMean, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_lvl, (/(i,i=1,lm(region))/), status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_bnd, (/(i,i=1,lm(region)+1)/), status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_tracer, reshape(ncout(region)%field, (/imHalo,jmHalo,lm(region),ntracetloc,1/)), status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_tracernm, names, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_press, ncout(region)%pressure, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_gph, ncout(region)%gph, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_orography, ncout(region)%orography, status)
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_temperature, ncout(region)%temperature, status)
!     IF_NOTOK_RETURN(status=1)

!     ! close file
!     call MDF_Close( fid, status )
!     IF_NOTOK_RETURN(status=1)

!     write (gol, '("[user_output_mix_write]  Done writing arrays and closing output file.")'); call goPr

!     ncfile(region)%orography = 0.0
!     ncfile(region)%temperature = 0.0
!     ncfile(region)%gph = 0.0
!     ncfile(region)%ps = 0.0
!     ncfile(region)%field = 0.0
!     ncfile(region)%wda = 0.0
!     ncfile(region)%idate_start = ncfile(region)%idate_end ! end of interval becomes start of next interval  

!     nullify(zoomed)

!     deallocate(ncout(region)%orography)
!     deallocate(ncout(region)%temperature)
!     deallocate(ncout(region)%pressure)
!     deallocate(ncout(region)%gph)
!     deallocate(ncout(region)%field)
!     deallocate(ncout(region)%lonMean)
!     deallocate(ncout(region)%latMean)
!     deallocate(ncout(region)%lonEdge)
!     deallocate(ncout(region)%latEdge)

    ! ok
    status = 0

  end subroutine user_output_mix_write

end module user_output_mix