!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module grids_util_mod
  !-----------------------------------------------------------------------
  !                   GNU General Public License                        
  !                                                                      
  ! This program is free software; you can redistribute it and/or modify it and  
  ! are expected to follow the terms of the GNU General Public License  
  ! as published by the Free Software Foundation; either version 2 of   
  ! the License, or (at your option) any later version.                 
  !                                                                      
  ! MOM is distributed in the hope that it will be useful, but WITHOUT    
  ! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  
  ! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    
  ! License for more details.                                           
  !                                                                      
  ! For the full text of the GNU General Public License,                
  ! write to: Free Software Foundation, Inc.,                           
  !           675 Mass Ave, Cambridge, MA 02139, USA.                   
  ! or see:   http://www.gnu.org/licenses/gpl.html                      
  !-----------------------------------------------------------------------
  !
  ! <CONTACT EMAIL= "Zhi.Liang@noaa.gov">Z. Liang </CONTACT>
  ! <REVIEWER EMAIL="Stephen.Griffies@noaa.gov">S. M. Griffies</REVIEWER>

  ! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>

  !<OVERVIEW>
  ! <TT>grids_util_mod</TT> contains some public interface used by several modules in
  !  generate_ocean_grid package. 
  !</OVERVIEW>

  use fms_mod,       only : stdout, string
  use mpp_mod,       only : mpp_pe, mpp_npes, mpp_root_pe, mpp_error, mpp_chksum, FATAL, NOTE 
  use mpp_io_mod,    only : mpp_open, MPP_OVERWR, MPP_NETCDF, MPP_SINGLE, axistype, fieldtype
  use mpp_io_mod,    only : mpp_write_meta, mpp_write, mpp_close
  use constants_mod, only : PI

  implicit none
  private

  public :: make_axis, gcell, get_file_unit, write_field_meta, write_field_data, set_grid

  !--- public interface ------------------------------------------------
  ! <INTERFACE NAME="write_field_data">
  !   <OVERVIEW>
  !     Write data to corresponding grid file. 
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     For the purpose of generating higher resolution grid. There is a
  !     2 GB limit of each grid file. The interface will allow the grid information
  !     be stored at multiple files.
  !   </DESCRIPTION>
  !   <IN NAME="filename" TYPE= "character(len=*)">
  !    The name of the grid file to be generated. If the grid file is 
  !    over 2 GB limit, it will break into several files with file name filename,
  !    filename2, filename3 ....
  !   </IN>
  !   <IN NAME="fieldname" TYPE= "character(len=*)">
  !    name of the field to be written into the file filename
  !   </IN>
  !   <IN NAME="fielddata" TYPE= "real, dimension(:,:) or dimension(:,:,:)">
  !    data of fieldname to be written to the file filename.
  !   </IN>
  ! </INTERFACE>
  !
  interface write_field_data
     module procedure  write_field_data_2d
     module procedure  write_field_data_3d
  end interface write_field_data

  real, dimension(:), allocatable :: xt, yt, xc, yc
  integer                         :: ni, nj               ! grid size
  logical                         :: is_root_pe = .false. ! to indicate if current pe is root pe.

  real, parameter :: max_file_size = 4294967295.  !The maximum size in bytes for any one file. 
                                                  ! With NetCDF3, this  should be 2 Gb or less. 
                                                  ! the reason using real number is the maximum integer.
  real, parameter    :: tolr = 1.0e-4

  integer, parameter :: max_fields = 300    ! can be increased if needed
  integer, parameter :: max_files = 50      ! can be increased if needed
  integer            :: num_files = 0       ! if >= 0, means a series file will be created.
  integer            :: num_fields = 0      ! number of fields written to the file.


  character(len=128)   :: opened_files(0:max_files) = ''
  integer              :: files_unit(0:max_files) = 0
  type(axistype),save  :: axis_xt(0:max_files), axis_yt(0:max_files)
  type(axistype),save  :: axis_xc(0:max_files), axis_yc(0:max_files)
  type(axistype),save  :: axis_v(0:max_files)
  type(fieldtype),save :: flds(max_fields)
  character(len=128)   :: flds_name(max_fields)
  integer              :: file_num(max_fields) = -1

contains

   !#######################################################################
  ! <SUBROUTINE NAME="gcell">
  !   <OVERVIEW>
  !     grid cell construction.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     A domain is composed of one or more regions: Build "num" T cells with resolution 
  !     "deltat(n) n=1,num" within the domain composed of regions bounded by "bounds".
  !     Also construct "num" C-cells of resolution "deltau(n) n=1,num" with the relation 
  !     between T and U cells given by: deltat(n) = 0.5*(deltau(n-1) + deltau(n)).
  !     Resolution may be constant or smoothly varying within each region AND there must 
  !     be an integral number of grid cells within each region. The domain is the sum of all regions.
  !   </DESCRIPTION>
  !   <TEMPLATE>
  !     call gcell (maxlen, n_bounds, bounds, d_bounds, nbpts, num, deltat, deltau, stretch)
  !   </TEMPLATE>
  !   <IN NAME="maxlen" TYPE="integer" >
  !     maximum length of "deltat" and "deltau"
  !   </IN>
  !   <IN NAME="n_bounds" TYPE="integer" >
  !     number of bounds needed to define the regions
  !   </IN>
  !   <IN NAME="bounds" TYPE="real, dimension(n_bounds)" >
  !     latitude, longitude, or depth at each bound
  !   </IN>
  !   <IN NAME="d_bounds" TYPE="real, dimension(n_bounds)" >
  !     delta (resolution) at each of the "bounds"
  !   </IN>
  !   <IN NAME="nbpts" TYPE="integer" >
  !     number of extra boundary cells to add to the domain. (usually one at the beginning and end)
  !   </IN>
  !   <IN NAME="stretch" TYPE="real" >
  !     stretching factor for last region (should only be used in the vertical) to provide
  !     increased stretching of grid points. "stretch" = 1.0 gives no increased stretching.
  !     "stretch" = 1.2 gives increased stretching...etc
  !   </IN>
  !   <IN NAME="debug" TYPE="logical, optional" >
  !      flag that controls standard output.
  !   </IN>
  !   <OUT NAME="num" TYPE="integer" >
  !     total number of grid cells within the domain
  !   </OUT>
  !   <OUT NAME="deltat" TYPE="real, dimension(1-nbpts:maxlen)" >
  !     resolution of T grid cells: n=1,num
  !   </OUT>
  !   <OUT NAME="deltau" TYPE="real, dimension(1-nbpts:maxlen)" >
  !     resolution of C grid cells: n=1,num
  !   </OUT>
  ! </SUBROUTINE>
  subroutine gcell (region, bnd1, bnd2, dbnd1, dbnd2, num, ntotal, maxlen, deltau)
    integer, intent(in)             :: region
    real,    intent(in)             :: bnd1, bnd2, dbnd1, dbnd2
    integer, intent(out)            :: num
    integer, intent(in)             :: ntotal, maxlen
    real, dimension(1:maxlen), intent(out) :: deltau

    real    :: avg_res, chg_res, wid, an, err, d_new, del
    integer :: m, i

    avg_res = 0.5*(dbnd1 + dbnd2)
    chg_res = dbnd2 - dbnd1
    wid = abs(bnd2 - bnd1)
    an  = wid/avg_res
    m   = nint(an)
    err = wid - avg_res * m

    if( abs(err) > tolr ) then
       write(stdout(),*) "==>Error: non integral number of cells in region #", region, &
            ", average resolution within region =",avg_res,                            &
            ", this implies ",an," grid cells, Change grid specifications ",           &
            "via table, Here is some help..."
       if(m>1) then
          d_new = (2.0*wid)/(m-1) - dbnd1;
          write(stdout(),*) "Note: to get ",m-1," grid cells within region ",region, &
               ", change resolution from ",dbnd2," to ",d_new
       endif

       d_new = (2.0*wid)/m - dbnd1
       write(stdout(),*)"Note: to get ",m," grid cells within region ",region, &
            ", change resolution from ",dbnd2," to ",d_new      
       d_new = (2.0*wid)/(m+1) - dbnd1
       write(stdout(),*)"Note: to get ",m+1," grid cells within region ",region, &
            ", change resolution from ",dbnd2," to ",d_new

       call mpp_error(FATAL);
    endif

    if(ntotal+m > maxlen) call mpp_error(FATAL, "grids_util_mod: maxlen exceeded in gcell, "// &
                                      "increase size of maxlen in hgrid_mod or vgrid_mod")

    ! Calculate resolution of corner cells: "deltau"
    ! T grid points will be centered in these cells

    do i = 1, m
       del = avg_res - 0.5*chg_res*cos((PI/m)*(i-0.5))
       deltau(i) = del
    enddo

    num = m

    return

  end subroutine gcell

!///////////////////////////////////////////////////////////
!/
!/ Make_axis
!/ define horizontal grid resolution and axis data
!/
!///////////////////////////////////////////////////////////

subroutine make_axis(cart, maxlen, num_regions,bounds,delta,corners,centers, &
                     num, square_grid, extend_square_grid )
    character(len=1),       intent(in) :: cart
    integer,                intent(in) :: maxlen
    integer,             intent(inout) :: num_regions
    real,  dimension(:), intent(inout) :: bounds, delta
    real, dimension(0:), intent(inout) :: corners, centers
    integer,             intent(inout) :: num
    logical,                intent(in) :: square_grid, extend_square_grid

  integer:: i, n, m, ii
  real   :: delta_corners(maxlen)

  if(cart == 'Y') then  
     if(square_grid) then
        call iso_grid (maxlen, num_regions, bounds, delta, num, corners, centers, extend_square_grid)
        return
     endif
   endif

   !--- bounds should increase monotonically
   do n = 2, num_regions
      if( bounds(n-1) .gt. bounds(n)) then
         do m=1,num_regions
            write (stdout(),'(i3,f10.5)') m, bounds(m)
         end do
         call mpp_error(FATAL, &
              ' grids_util_mod: longitude boundaries with cart = '// cart //&
              'do not increase monotonically ')
      endif
   enddo

  num = 0;
  do n = 2, num_regions
    write(stdout(),*) " region # ",n," going from ",bounds(n-1)," (res=",delta(n-1), &
           ") to ",bounds(n)," (res=",delta(n),")"
    
    call gcell (n, bounds(n-1), bounds(n), delta(n-1), delta(n), m, num, maxlen, delta_corners)

    !*********************************************************************
    !   Build the grid points on a "B" grid. The T and C
    !   cells are staggered in the horizontal but at the same level in
    !   the vertical. However, the W cells (for vertical
    !   advection velocities at the bottoms of the C and T cells) are
    !   staggered in the vertical with respect to C and T cells.
    !   (no extra boundary point at the start)
    !*********************************************************************/
    corners(num) = bounds(n-1)

    do i = 1, m
      ii = num+i
      corners(ii) = corners(ii-1) + delta_corners(i)
    enddo
    num = num+m
  enddo

  corners(num) = bounds(num_regions)
  corners(num+1) = 2*corners(num) - corners(num-1)
  centers(0) = bounds(1) - 0.5*delta(1)

  do i = 1, num+1
     centers(i) = 2.0*corners(i-1) - centers(i-1)
  enddo 

  return   

end subroutine make_axis

  !#######################################################################
  subroutine iso_grid (maxlen, nylats, y_lat, dy_lat, nj, corners, centers, extend_square_grid)
    ! compute latitude resolution of grid cells to match the convergence
    ! of meridians.
    ! 
    ! inputs:
    !
    ! maxlen   = maximum length of "dytdeg0" and "dyudeg0"
    ! nylats   = should equal 2 (defines one region)
    ! dy_lat   = latitudinal resolution of grid cell on equator
    ! y_lat(1) = southern boundary of the domain (it will be adjusted
    !            to fit an integral number of cells)
    ! y_lat(2) = northern boundary of the domain (it will be adjusted
    !            to fit an integral number of cells)
    !
    ! outputs:
    !
    ! nj      = number of grid cells
    ! yt0     = latitude of point within T cell (degrees)
    ! yu0     = latitude of point within U cell (degrees)

    integer,                   intent(in)    :: maxlen, nylats
    integer,                   intent(out)   :: nj
    real, dimension(nylats),   intent(inout) :: y_lat
    real, dimension(nylats),   intent(in)    :: dy_lat
    real, dimension(0:maxlen), intent(out)   :: corners, centers
    logical,                   intent(in)    :: extend_square_grid
    !--- local variables -------------------------------------------------
    integer,                    parameter :: jmaxlat = 1000, lenjmax = 2*jmaxlat+1
    real, dimension(-jmaxlat-1:jmaxlat+1) :: dusq, dtsq, usq, tsq 
    real,                    dimension(2) :: y_bound, dy_bound
    integer                               :: n, n1, n2, nps, npn, j1, jmts, jmtn, j
    real                                  :: pole1, pole2, wid, stretch, D2R
    !---------------------------------------------------------------------

    D2R = PI/180.

    if (nylats .gt. 2) then
       call mpp_error(FATAL,'grids_util_mod: when hgrid_nml "square_grid"=true, hgrid_nml "nylats"= ' &
                 //trim(string(nylats))//' should eqaul 2')
    endif

    if (dy_lat(1) .ne. dy_lat(2)) then
       call mpp_error(FATAL,'grids_util_mod: nml dy_lat(1)= '//trim(string(dy_lat(1)))// &
             ' must equal dy_lat(2)= '//trim(string(dy_lat(2)))//' when hgrid_nml "square_grid"=true')
    endif
    if ((abs(y_lat(1)) > 89.9999) .or. (abs(y_lat(2)) > 89.9999)) then
       call mpp_error(FATAL,'grids_util_mod: Cannot specify hgrid_nml "ylat" = 90 deg ' &
                     //'when hgrid_nml "square_grid"=true ')
    endif

    ! build a square grid

    usq(0)  = 0.0
    dusq(0) = dy_lat(1)
    tsq(1)  = 0.5*dusq(0)
    tsq(0)  = -tsq(1)
    dtsq(1) = dusq(0)*cos(tsq(1)*D2R)
    dtsq(0) = dtsq(1)
    do n=1,jmaxlat
       dusq(n)   = 2.0*dtsq(n) - dusq(n-1)
       usq(n)    = tsq(n) + 0.5*dusq(n)
       if (tsq(n) .lt. 90.0) then
          tsq(n+1)  = tsq(n) + dusq(n)
          dtsq(n+1) = dusq(0)*cos(tsq(n+1)*D2R)
       else
          tsq(n+1)  = tsq(n)
          dtsq(n+1) = 0.0
       endif
    enddo
    do n=-1,-jmaxlat,-1
       dusq(n) = dusq(-n)
       usq(n)  = -usq(-n)
       dtsq(n) = dtsq(-(n-1)) 
       tsq(n)  = -tsq(-(n-1))
    enddo

    ! pick out cells between bounding latitudes

    n1 = -jmaxlat
    n2 = jmaxlat
    do n=0,jmaxlat
       if (usq(n) > y_lat(2)) then
          n2 = n
          exit
       endif
    enddo
    do n=0,-jmaxlat,-1
       if (usq(n) < y_lat(1)) then
          n1 = n+1
          exit
       endif
    enddo

    ! re-define bounding latitudes to match square boundaries
    y_lat(1) = usq(n1-1)
    y_lat(2) = usq(n2)  

    do n=n1,n2
       write(stdout(),*) 'n=',n,' tsq=',tsq(n),', dyt=',dtsq(n)
    enddo
    if (n1 .eq. -jmaxlat .or. n2 .eq. jmaxlat) then
       call mpp_error(FATAL,'hgrid_mod: Need to increase jmaxlat to reach the max latitude')
    endif

    if (extend_square_grid) then

       ! extend square grid in southern hemisphere to South pole

       ! set south pole (pole1) to skip over land in antarctica
       pole1 = -80.0
       stretch = 1.0
       wid = (y_lat(1) - pole1)
       nps = nint(wid/dusq(n1)) - 1
       y_bound(1) = pole1
       y_bound(2) = y_lat(1)
       dy_bound(1) = 2.0*wid/nps - dusq(n1)
       dy_bound(2) = dusq(n1)
       j1 = n1-nps
       call gcell(1, y_bound(1), y_bound(2), dy_bound(1), dy_bound(2), jmts, 0, maxlen, dtsq(j1) ) 
       dusq(j1) = 2*dtsq(j1)-dy_bound(1)
       do n = 1, jmts-1
          dusq(j1+n) = 2*dtsq(j1+n)-dusq(j1+n-1)
       enddo

       ! construct latitude of T and U grid points
       do n=n1-1, n1-jmts+1, -1
          tsq(n) = tsq(n+1) - dusq(n)
          usq(n) = tsq(n+1) - 0.5*dusq(n)
          write(stdout(),*)'extended n=',n,' usq=',usq(n),'tsq=',tsq(n), 'dyt=',dusq(n)
       enddo

       n1 = n1 - jmts + 1

       ! extend square grid to North pole

       pole2 = 90.0
       npn = nint((pole2 - y_lat(2))/dusq(n2)) - 1
       wid = (pole2 - y_lat(2))
       y_bound(1) = y_lat(2)
       y_bound(2) = pole2
       dy_bound(1) = dusq(n2)
       dy_bound(2) = 2.0*wid/npn - dusq(n2)

       call gcell(1, y_bound(1), y_bound(2), dy_bound(1), dy_bound(2), jmtn, 0, maxlen, dtsq(n2+1) ) 
       dusq(n2+1) = 2*dtsq(n2+1)-dy_bound(1)

       do n = 1, jmtn-1
          dusq(n2+n) = 2*dtsq(n2+n)-dusq(n2+n-1)
       enddo

       ! construct latitude of T and U grid points
       do n = n2+1, n2+jmtn-1
          tsq(n) = tsq(n-1) + dusq(n)
          usq(n) = tsq(n) + 0.5*dusq(n)
          write(stdout(),*)'extended n=',n,' usq=',usq(n),'tsq=',tsq(n), 'dyt=',dusq(n)
       enddo
       n2 = n2 + jmtn - 1
       ! re-define bounding latitudes to match extended boundaries

       y_lat(1) = pole1
       y_lat(2) = pole2 
    endif
    nj = n2 - n1 + 1

    do j=1,nj
       corners(j)     = usq(j+n1-1)
       centers(j)     = tsq(j+n1-1)
    enddo
    
    corners(0) = 2*centers(1) - corners(1)
    centers(0) = 2*corners(0) - centers(1)

    centers(nj+1) = 2*corners(nj)-centers(nj)
    corners(nj+1) = 2*centers(nj+1)- corners(nj)

  end subroutine iso_grid


  !#######################################################################
  ! This function is only for global meta and vgrid. So get_file_unit should 
  ! always linked to file opened_file(0)
  ! <FUNCTION NAME="get_file_unit">
  !   <OVERVIEW>
  !    returns the io unit corresponding to filename.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !    If the file filename is already open, return the io unit of this 
  !    opened file. Otherwise will open the file and return the io unit.
  !   </DESCRIPTION>
  !   <TEMPLATE>
  !     get_file_unit(filename)
  !   </TEMPLATE>
  !   <IN NAME="filename" TYPE= "character(len=*)">
  !    The name of the grid file to be generated.
  !   </IN>

  function get_file_unit(filename)   
    character(len=*), intent(in) :: filename
    integer                      :: get_file_unit

    if(.not. is_root_pe) return

    if(trim(opened_files(0)) == trim(filename) ) then 
          get_file_unit = files_unit(0)
          return
    endif

    !--- if file is not opened, open the file

    call mpp_open(get_file_unit, trim(filename), MPP_OVERWR,MPP_NETCDF, &
          threading=MPP_SINGLE,  fileset=MPP_SINGLE )

   files_unit(num_files) = get_file_unit
   opened_files(num_files) = trim(filename)

  end function get_file_unit
  ! </FUNCTION>
  !#####################################################################
  ! <SUBROUTINE NAME="write_field_meta">
  !   <OVERVIEW>
  !    Write meta data of a field to a netcdf file.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !    It will check if the grid file will over the 2 GB limit. If do, will open 
  !    a new file with name filename? (? is 1, 2, 3 ....) and write axis metadata
  !    to the new file.
  !   </DESCRIPTION>
  !   <TEMPLATE>
  !    call write_field_meta(filename, fieldname, units, field_longname, fielddim, x_pos, y_pos)
  !   </TEMPLATE>
  !   <IN NAME="filename" TYPE= "character(len=*)">
  !    The name of the grid file to be generated. If the grid file is 
  !    over 2 GB limit, it will break into several files with file name filename,
  !    filename1, filename2, filename3 ....
  !   </IN>
  !   <IN NAME="fieldname" TYPE= "character(len=*)">
  !    name of the field to be written into the file filename
  !   </IN>
  !   <IN NAME="units" TYPE= "character(len=*)">
  !    units of field fieldname.
  !   </IN>
  !   <IN NAME="field_longname" TYPE= "character(len=*)">
  !    longname of fielname.
  !   </IN>
  !   <IN NAME="fielddim" TYPE= "integer">
  !    Indicate the dimension of fieldname. fielddim should be either 2 or 3. 
  !   </IN>
  !   <IN NAME="x_pos, y_pos" TYPE= "character(len=1)">
  !    To indicate the cell position. its value can be "T" or "C".
  !   </IN>
  subroutine write_field_meta(filename, fieldname, units, field_longname, fielddim, x_pos, y_pos)
    character(len=*),           intent(in) :: filename, fieldname
    character(len=*),           intent(in) :: units, field_longname
    integer,                    intent(in) :: fielddim
    character(len=1), intent(in), optional :: x_pos, y_pos

    character(len=128) :: curr_file, basename
    logical            :: is_first_field 
    integer            :: unit, fieldsize, length
    type(axistype)     :: axis_x, axis_y
    logical, save      :: do_write_axis = .true.
    real, save         :: size_in_file = 20000.0

    if(.not. is_root_pe) return

    curr_file = filename
    num_fields = num_fields + 1
    flds_name(num_fields) = trim(fieldname)

    if(num_fields ==1) is_first_field = .true.

    if(size_in_file .gt. max_file_size) then   ! 
       num_files = num_files + 1
       if(num_files .gt. max_files) call mpp_error(FATAL, &
           'grids_util_mod: number of files is over max_files, increase max_files ')

    endif

    if(num_files .gt. 0) then
       length = len_trim(filename)
       if(filename(length-2:length) == '.nc') then
          basename = filename(1:length-3)
       else
          basename = filename(1:length)
       endif
       if(num_files < 10)  then   ! num_files should be less than 100
          write(curr_file, '(a,i1,a)')trim(basename), num_files, '.nc'
       else
          write(curr_file, '(a,i2,a)')trim(basename), num_files, '.nc'
       endif
    endif

    if(trim(opened_files(num_files)) == curr_file) then  ! the file is already opened
       unit = files_unit(num_files)
    else   ! need to open the file and write out axis meta.
       call mpp_open(unit, trim(curr_file), MPP_OVERWR,MPP_NETCDF, &
            threading=MPP_SINGLE,  fileset=MPP_SINGLE )

       opened_files(num_files) = trim(curr_file)
       files_unit(num_files) = unit
       do_write_axis = .true.
    endif
    file_num(num_fields) = num_files
    if(do_write_axis) then
       !--- write axis meta -------------------------------------------
       call mpp_write_meta(unit, axis_xt(num_files),'grid_x_T','degree_east','Nominal Longitude of T-cell center', &
            cartesian ='X', data = xt )
       call mpp_write_meta(unit, axis_yt(num_files),'grid_y_T','degree_north','Nominal Latitude of T-cell center', &
            cartesian ='Y', data = yt )
       call mpp_write_meta(unit, axis_xc(num_files),'grid_x_C','degree_east','Nominal Longitude of C-cell center', &
            cartesian ='X', data = xc )
       call mpp_write_meta(unit, axis_yc(num_files),'grid_y_C','degree_north','Nominal Latitude of C-cell center', &
            cartesian ='Y', data = yc )
       call mpp_write_meta(unit, axis_v(num_files), 'vertex', 'none ',        &
            'Vertex position from southwest couterclockwise', data = (/1.,2.,3.,4./) )   
       do_write_axis = .false.
       !--- open a new file, set size_in_file to 20000
       size_in_file = 20000.

    endif

    if(present(x_pos)) then
      if (x_pos .ne. 'T' .and. x_pos .ne. 'C')  &
        call mpp_error(FATAL,'grids_util_mod: x_pos should be either "T" or "C" ')
    endif

    if(present(y_pos)) then
        if (y_pos .ne. 'T' .and. y_pos .ne. 'C') &
        call mpp_error(FATAL,'grids_util_mod: x_pos should be either "T" or "C" ')
    endif

    axis_x = axis_xt(num_files)
    if(present(x_pos))then
       if( x_pos == 'C') axis_x = axis_xc(num_files)
    endif
    axis_y = axis_yt(num_files)
    if(present(y_pos))then
       if(y_pos == 'C') axis_y = axis_yc(num_files)
    endif
    !--- write out field meta
    if(fielddim==2) then
       fieldsize = 8*ni*nj
       call mpp_write_meta(unit, flds(num_fields), (/axis_x, axis_y/), fieldname, units, field_longname, pack=1) 
    else if(fielddim==3) then
       fieldsize = 32*ni*nj
       call mpp_write_meta(unit, flds(num_fields), (/axis_x, axis_y,axis_v(num_files)/), &
            fieldname, units, field_longname, pack=1) 
    else
       call mpp_error(FATAL, 'grids_util_mod: dimension of field '//trim(fieldname)//' should be either 2 or 3')
    endif

    size_in_file = size_in_file + fieldsize


  end subroutine write_field_meta
  !  </SUBROUTINE>
  !#####################################################################
  !  <SUBROUTINE NAME="write_field_data_2d" INTERFACE="write_field_data">
  !  <IN NAME="filename" TYPE= "character(len=*)"> </IN>
  !   <IN NAME="fieldname" TYPE= "character(len=*)"> </IN>
  !   <IN NAME="fielddata" TYPE= "real, dimension(:,:)"> </IN>
  !<PUBLICROUTINE INTERFACE="write_field_data">
  subroutine write_field_data_2d(filename, fieldname, fielddata)
  !</PUBLICROUTINE>
    character(len=*),     intent(in) :: filename, fieldname
    real, dimension(:,:), intent(in) :: fielddata             
    integer                          :: n, k
    logical                          :: is_first_field, is_last_field

    if(.not. is_root_pe) return

    is_first_field = .false.
    is_last_field = .false.
    do n = 1, num_fields
       if(trim(flds_name(n)) == trim(fieldname)) then
          k = file_num(n)
          if(n ==1) then
             is_first_field = .true.
          else
             if(file_num(n) .ne. file_num(n-1)) then
                is_first_field = .true.
             endif
          endif
          if(n == max_fields) then
             is_last_field = .true.
          else
             if(file_num(n) .ne. file_num(n+1)) is_last_field = .true.
          endif

          ! if current field is the first field in a file, need to writ out axis data
          if(is_first_field ) then 
             call mpp_write(files_unit(k), axis_xt(k) )
             call mpp_write(files_unit(k), axis_yt(k) )
             call mpp_write(files_unit(k), axis_xc(k) )
             call mpp_write(files_unit(k), axis_yc(k) )
             call mpp_write(files_unit(k), axis_v(k) )
          endif

          !--- data should be already on global domain
          if(size(fielddata,1) .ne. ni .or. size(fielddata,2) .ne. nj ) then
             call mpp_error(FATAL,'grids_util_mod: data of field '//trim(fieldname)//' is not on global domain')
          endif

          call mpp_write(files_unit(k), flds(n), fielddata)
          if(is_last_field) call mpp_close(files_unit(k))

          return
       endif
    enddo

    !--- if fieldname is not found in the array flds_name, abort the program
    call mpp_error(FATAL,'grids_util_mod: file '//trim(filename)// &
         ' does not contain meta data of field '//trim(fieldname))

  end subroutine write_field_data_2d
  !  </SUBROUTINE>

  !#####################################################################
  !  <SUBROUTINE NAME="write_field_data_3d" INTERFACE="write_field_data">
  !   <IN NAME="fielddata" TYPE= "real, dimension(:,:,:)"> </IN>
  subroutine write_field_data_3d(filename, fieldname, fielddata)
    character(len=*),       intent(in) :: filename, fieldname
    real, dimension(:,:,:), intent(in) :: fielddata             
    integer                            :: n, k
    logical                            :: is_first_field, is_last_field

    if(.not. is_root_pe) return

    is_first_field = .false.
    is_last_field = .false.
    do n = 1, num_fields
       if(trim(flds_name(n)) == trim(fieldname)) then
             k = file_num(n)
          if(n ==1) then
             is_first_field = .true.
          else
             if(file_num(n) .ne. file_num(n-1)) then
                is_first_field = .true.
             endif
          endif 
        if(n == max_fields) then
           is_last_field = .true.
        else
           if(file_num(n) .ne. file_num(n+1)) is_last_field = .true.
        endif
           
          ! if current field is the first field in a file, need to writ out axis data
          if(is_first_field) then 
             call mpp_write(files_unit(k), axis_xt(k) )
             call mpp_write(files_unit(k), axis_xt(k) )
             call mpp_write(files_unit(k), axis_xc(k) )
             call mpp_write(files_unit(k), axis_yc(k) )
             call mpp_write(files_unit(k), axis_v(k) )
          endif
          call mpp_write(files_unit(k), flds(n), fielddata)
          if(is_last_field) call mpp_close(files_unit(k))

          return
        endif
     enddo

     !--- if fieldname is not found in the array flds_name, abort the program
     call mpp_error(FATAL,'grids_util_mod: file '//trim(filename)// &
          ' does not contain meta data of field '//trim(fieldname))

  end subroutine write_field_data_3d
  !  </SUBROUTINE>

  !#####################################################################
  ! <SUBROUTINE NAME="set_grid">

  !   <OVERVIEW>
  !    set the axis grid information.
  !   </OVERVIEW>     
  !   <TEMPLATE>
  !    call set_grid(grid_xt, grid_yt, grid_xc, grid_yc)
  !   </TEMPLATE>
  !   <IN NAME="grid_xt, grid_yt" TYPE="real, dimension(:)">
  !    longitude and latitude of the T-cell grid.
  !   </IN>
  !   <IN NAME="grid_xc, grid_yc" TYPE="real, dimension(:)">
  !    longitude and latitude of the C-cell grid.
  !   </IN>
  subroutine set_grid(grid_xt, grid_yt, grid_xc, grid_yc)
    real, dimension(:), intent(in) :: grid_xt, grid_yt, grid_xc, grid_yc


    ni = size(grid_xt(:))
    nj = size(grid_yt(:))

    allocate(xt(ni), yt(nj), xc(ni), yc(nj))
    xt = grid_xt
    yt = grid_yt
    xc = grid_xc
    yc = grid_yc 

    if(mpp_pe() == mpp_root_pe()) is_root_pe = .true.

  end subroutine set_grid
  !  </SUBROUTINE>
  !#####################################################################

end module grids_util_mod