module module_write_nemsio

  use esmf
  use nemsio_module
  use module_fv3_io_def, only : write_nemsioflip, write_fsyncflag

  implicit none

  include 'mpif.h'

  private
  logical :: first_nemsio_call
  integer :: im,jm,lm, idate(7),nmeta, nsoil,ncld, idrt, ntrac
  integer :: mype, ntasks, mpi_comm, nbdl
  logical :: hydrostatic
  real(kind=ESMF_KIND_R4) :: varr4
  integer,dimension(200,100)      :: nfldlev
  character(16),dimension(3000,5) :: recname,reclevtyp
  integer,dimension(3000,5)       :: reclev

  integer, dimension(:), allocatable   :: nrec
  integer, dimension(:), allocatable   :: idsl, idvc,idvm
  integer, dimension(:), allocatable   :: fieldcount
  integer, dimension(:), allocatable   :: idisp, irecv
!
  integer,dimension(:), allocatable :: nmetavari,nmetavarc, nmetavarr4,nmetavarr8
  integer,dimension(:), allocatable :: nmetaaryi,nmetaaryc, nmetaaryr4,nmetaaryr8
  character(16),dimension(:,:),allocatable :: variname, varcname, varr4name, varr8name
  character(16),dimension(:,:),allocatable :: aryiname
  integer, dimension(:,:),   allocatable   :: varival, aryilen
  integer, dimension(:,:,:), allocatable   :: aryival
  real(4), dimension(:,:),   allocatable   :: varr4val
  real(8), dimension(:,:),   allocatable   :: varr8val
  character(16), dimension(:,:), allocatable   :: varcval
  logical, dimension(:), allocatable :: extrameta
!
  real(4), dimension(:,:,:), allocatable  :: vcoord
  real(4), dimension(:), allocatable      :: lon1d, lat1d

  logical         :: first_set

  public nemsio_first_call, write_nemsio

  contains

  subroutine nemsio_first_call(fieldbundle, imo, jmo, &
             wrt_mype, wrt_ntasks, wrt_mpi_comm, wrt_nbdl, mybdl,   &
             inidate, lat, lon, rc)
    type(ESMF_FieldBundle), intent(in)     :: fieldbundle
    integer, intent(in)                    :: imo, jmo
    integer, intent(in)                    :: wrt_mype, wrt_ntasks, wrt_mpi_comm
    integer, intent(in)                    :: wrt_nbdl, mybdl
    integer, intent(in)                    :: inidate(7)
    real(8), intent(in)                    :: lat(:), lon(:)
    integer, optional,intent(out)          :: rc

!** local vars
    integer i,j, nfld, nidx, nlen
    integer fieldDimCount,gridDimCount
    character(100) :: fieldname
    type(ESMF_GRID)            :: wrtgrid
    type(ESMF_TypeKind_Flag)   :: typekind

    integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound
    type(ESMF_Field), allocatable  :: fcstField(:)
    real(ESMF_KIND_R8), dimension(:,:), pointer   :: lonPtr, latPtr

!-------------------------------------------------------------------
!
    im    = imo
    jm    = jmo
    nmeta = 8
    idrt  = 4
    idate(1:7) = inidate(1:7)
    mype  = wrt_mype
    ntasks= wrt_ntasks
    nbdl  = wrt_nbdl
    mpi_comm = wrt_mpi_comm
    if(.not.allocated(idisp)) allocate(idisp(ntasks),irecv(ntasks))
    if(.not.allocated(fieldCount)) allocate(fieldCount(nbdl))
    if(.not.allocated(nrec)) allocate(nrec(nbdl))
    if(.not.allocated(idsl)) then
      allocate(idsl(nbdl),idvc(nbdl),idvm(nbdl))
      idsl=-9999; idvc=-9999; idvm=-9999; first_set=.true.
      nsoil = 4
      ntrac = 3
      ncld  = 1
    endif
!
!** get attibute info from fieldbundle
    call get_global_attr(fieldbundle, mybdl, rc=rc)

!** get meta info from fieldbundle
    call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount(mybdl), &
         grid=wrtGrid, rc=rc)
    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
               line=__LINE__, &
               file=__FILE__)) &
               return  ! bail out
!
    call ESMF_GridGet(wrtgrid, dimCount=gridDimCount, rc=rc)
    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
      line=__LINE__, &
      file=__FILE__)) &
      return  ! bail out
!
!** get lat lon:
    if(.not.allocated(lon1d)) allocate(lon1d(imo*jmo))
    if(.not.allocated(lat1d)) allocate(lat1d(imo*jmo))
    do j=1,jmo
      do i=1,imo
        lon1d((j-1)*imo+i) = lon(i)
        lat1d((j-1)*imo+i) = lat(j)
      enddo
    enddo
!
    allocate(fcstField(fieldCount(mybdl)))
    call ESMF_FieldBundleGet(fieldbundle, fieldList=fcstField,     &
         itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc)
    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
      line=__LINE__, &
      file=__FILE__)) &
      return  ! bail out

    nrec(mybdl)=0
    do i=1,fieldcount(mybdl)

       call ESMF_FieldGet(fcstField(i), typekind=typekind, &
            dimCount=fieldDimCount, grid=wrtGrid, name=fieldname, rc=rc)
       if( index(trim(fieldname),"vector") >0) cycle
       if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
         line=__LINE__, &
         file=__FILE__)) &
         return  ! bail out
       if (fieldDimCount > gridDimCount) then
         allocate(ungriddedLBound(fieldDimCount-gridDimCount))
         allocate(ungriddedUBound(fieldDimCount-gridDimCount))
         call ESMF_FieldGet(fcstField(i), ungriddedLBound=ungriddedLBound, &
                    ungriddedUBound=ungriddedUBound, rc=rc)
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
           line=__LINE__, &
           file=__FILE__)) &
           return  ! bail out
         nfldlev(i,mybdl) = ungriddedUBound(fieldDimCount-gridDimCount) - &
                        ungriddedLBound(fieldDimCount-gridDimCount) + 1
         nrec(mybdl) = nrec(mybdl) + nfldlev(i,mybdl)
         deallocate(ungriddedLBound)
         deallocate(ungriddedUBound)
       else if(fieldDimCount == 2) then
         nfldlev(i,mybdl) = 1
         nrec(mybdl) = nrec(mybdl) + 1
       endif

    enddo
!
    nfld = 1
    do i=1,fieldcount(mybdl)
      call ESMF_FieldGet(fcstField(i),name=fieldName,rc=rc)
      if( nfldlev(i,mybdl) == 1) then
        reclevtyp(nfld,mybdl) = "sfc"
        reclev(nfld,mybdl)    = 1
        recname(nfld,mybdl)   = trim(fieldName)

        nlen = len(trim(fieldName))
        if( nlen > 3) then
          if( fieldName(nlen-2:nlen) == 'sfc') recname(nfld,mybdl) = fieldName(1:nlen-3)
          if( fieldName(nlen-2:nlen) == 'nst') recname(nfld,mybdl) = fieldName(1:nlen-3)
          if (fieldName(nlen-2:nlen)== '10m') then
            reclevtyp(nfld,mybdl) = "10 m above gnd"
            if( fieldName(nlen-2:nlen) == "10m") then
              recname(nfld,mybdl) = fieldName(1:nlen-3)
              if (trim(fieldName) == "f10m") recname(nfld,mybdl) = trim(fieldName)
            endif
          else if (fieldName(nlen-1:nlen)== '2m') then
            reclevtyp(nfld,mybdl) = "2 m above gnd"
            if( fieldName(nlen-1:nlen) == '2m') then
              recname(nfld,mybdl) = fieldName(1:nlen-2)
            endif
          else if( index(trim(fieldName),'soil')>0 ) then
            if( index(trim(fieldName),'soilt')>0) then
               recname(nfld,mybdl) = 'tmp'
            else if( index(trim(fieldName),'soilm')>0 .and.  &
              trim(reclevtyp(nfld,mybdl))=='sfc') then
               recname(nfld,mybdl) = 'soilm'
               reclevtyp(nfld,mybdl) = '0-200 cm down'
            else
              recname(nfld,mybdl)   = fieldName(1:nlen-1)
            endif
            if(fieldName(nlen:nlen)=='1') then
              reclevtyp(nfld,mybdl) = '0-10 cm down'
            else if (fieldName(nlen:nlen)=='2') then
              reclevtyp(nfld,mybdl) = '10-40 cm down'
            else if (fieldName(nlen:nlen)=='3') then
              reclevtyp(nfld,mybdl) = '40-100 cm down'
            else if (fieldName(nlen:nlen)=='4') then
              reclevtyp(nfld,mybdl) = '100-200 cm down'
            endif
          else if( fieldName(nlen-2:nlen) =='clm' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'atmos col'
          else if( fieldName(nlen-2:nlen) =='toa' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'nom. top'
          else if( fieldName(nlen-2:nlen) =='lcl' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'low cld lay'
          else if( fieldName(nlen-2:nlen) =='mcl' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'mid cld lay'
          else if( fieldName(nlen-2:nlen) =='hcl' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'high cld lay'
          else if( fieldName(nlen-2:nlen) =='lct' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'low cld top'
          else if( fieldName(nlen-2:nlen) =='mct' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'mid cld top'
          else if( fieldName(nlen-2:nlen) =='hct' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'high cld top'
          else if( fieldName(nlen-2:nlen) =='lcb' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'low cld bot'
          else if( fieldName(nlen-2:nlen) =='mcb' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'mid cld bot'
          else if( fieldName(nlen-2:nlen) =='hcb' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-3)
            reclevtyp(nfld,mybdl) = 'high cld bot'
          endif
        endif
        if( nlen > 5) then
          if( fieldName(nlen-4:nlen) =='cnvcl' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-5)
            reclevtyp(nfld,mybdl) = 'convect-cld laye'
          else if( fieldName(nlen-5:nlen) =='cnvclt' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-6)
            reclevtyp(nfld,mybdl) = 'convect-cld top'
          else if( fieldName(nlen-5:nlen) =='cnvclb' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-6)
            reclevtyp(nfld,mybdl) = 'convect-cld bot'
          else if( fieldName(nlen-4:nlen) =='bndcl' ) then
            recname(nfld,mybdl)   = fieldName(1:nlen-5)
            reclevtyp(nfld,mybdl) = 'bndary-layer cld'

          endif
        endif
        nfld = nfld + 1
      else
        lm   = nfldlev(i,mybdl)
        if(first_set) then
          idsl(mybdl) = 1
          idvc(mybdl) = 2
          idvm(mybdl) = 1
          first_set = .false.
        endif
        do j = 1,nfldlev(i,mybdl)
          recname(nfld,mybdl)   = trim(fieldName)
          if(trim(fieldName) == "preshy" .or. trim(fieldName) == "presnh") then
            recname(nfld,mybdl) = "pres"
          endif
          reclevtyp(nfld,mybdl) = "mid layer"
          reclev(nfld,mybdl)    = j
          nfld = nfld + 1
        enddo
      endif
    enddo
!
  end subroutine nemsio_first_call

!----------------------------------------------------------------------------------------
  subroutine write_nemsio(fieldbundle, filename, nf_hours, &
             nf_minutes, nf_seconds, nfsecond_num, nfsecond_den, mybdl, rc)
!
    type(ESMF_FieldBundle), intent(in)     :: fieldbundle
    character(*), intent(in)               :: filename
    integer, intent(in)                    :: nf_hours, nf_minutes, nf_seconds
    integer, intent(in)                    :: nfsecond_num, nfsecond_den
    integer, intent(in)                    :: mybdl
    integer, optional,intent(out)          :: rc
!
!** local vars
    integer i,ii,j,m,n,k, k1,k2,k3,jrec, nfseconds, nofsync, FFSYNC
    integer istart, iend, jstart, jend, kstart, kend, nlen
    real    fhour
    logical OPENED
    real(4),dimension(:),allocatable    :: tmp
    real(4),dimension(:,:),allocatable  :: arrayr4
    real(4),dimension(:,:),pointer      :: arrayr42d
    real(8),dimension(:,:),pointer      :: arrayr82d
    real(4),dimension(:,:,:),pointer    :: arrayr43d
    real(8),dimension(:,:,:),pointer    :: arrayr83d
    type(ESMF_Field), allocatable  :: fcstField(:)
    type(ESMF_TypeKind_Flag)       :: typekind
    type(nemsio_gfile) ::  nemsiofile
    character(128) :: fieldname
!

!** init nemsio
    call nemsio_init(iret=rc)
!
!**  OPEN NEMSIO FILE
!
    if (mype == 0) then

      nfseconds = nf_seconds*nfsecond_den + nfsecond_num

      print *,'in write_nemsio,bf nemsio_open, filename=',trim(filename),    &
      'idate=',idate,'nfour=',NF_HOURS,NF_MINUTES,NF_SECONDS, 'mybdl=',mybdl,          &
      'dim=',im,jm,lm,'nmeta=',nmeta,'idrt=',idrt,'nsoil=',nsoil,                      &
      'ntrac=',ntrac,'nrec=',nrec(mybdl),'extrameta=',extrameta(mybdl),                &
      'vcoord=',vcoord(1:5,1,1),'nfhours=',nf_hours,nf_minutes,nfseconds,nfsecond_den, &
      'idsl=',idsl(mybdl),'idvc=',idvc(mybdl),idvm(mybdl)
!      'nmetavari=',nmetavari(mybdl),'nmetavarc=',nmetavarc(mybdl)
!    if(nmetavari(mybdl)>0) print *,'in write_nemsio,bf nemsio_open,nmetavari=', &
!      nmetavari(mybdl),'varival=',trim(variname(1,mybdl)),varival(1,mybdl)
!    if(nmetavarc(mybdl)>0) print *,'in write_nemsio,bf nemsio_open,nmetavarc=', &
!      nmetavarc(mybdl),'varcval=',trim(varcname(1,mybdl)),varcval(1,mybdl)
!    if(nmetaaryi(mybdl)>0) print *,'in write_nemsio,bf nemsio_open,nmetaaryi=', &
!      nmetaaryi(mybdl),'aryival=',trim(aryiname(1,mybdl)),aryilen(1,mybdl),aryival(1:3,1,mybdl)
 
      call nemsio_open(nemsiofile,trim(FILENAME),'write',rc,             &
        modelname="FV3GFS", gdatatype="bin4",                            &
        idate=idate,nfhour=nf_hours, nfminute=nf_minutes,                &
        nfsecondn=nfseconds, nfsecondd=nfsecond_den,                     &
        dimx=im,dimy=jm,dimz=lm, nmeta=nmeta,idrt=idrt,                  &
        nsoil=nsoil,ntrac=ntrac,nrec=nrec(mybdl), ncldt=ncld,            &
        idsl=idsl(mybdl),idvc=idvc(mybdl), idvm=idvm(mybdl),             &
        vcoord=vcoord, lon=lon1d,lat=lat1d,                              &
        extrameta=extrameta(mybdl),recname=RECNAME(1:nrec(mybdl),mybdl), &
        reclevtyp=RECLEVTYP(1:nrec(mybdl),mybdl),                        &
        reclev=RECLEV(1:nrec(mybdl),mybdl),                              &
        nmetavari=nmetavari(mybdl), nmetavarr=nmetavarr4(mybdl),         &
        nmetavarc=nmetavarc(mybdl), nmetaaryi=nmetaaryi(mybdl),          &
        variname=variname(1:nmetavari(mybdl),mybdl),                     &
        varival=varival(1:nmetavari(mybdl),mybdl),                       &
        varrname=varr4name(1:nmetavarr4(mybdl),mybdl),                   &
        varrval=varr4val(1:nmetavarr4(mybdl),mybdl),                     &
        varcname=varcname(1:nmetavarc(mybdl),mybdl),                     &
        varcval=varcval(1:nmetavarc(mybdl),mybdl),                       &
        aryiname=aryiname(1:nmetaaryi(mybdl),mybdl),                     &  
        aryilen=aryilen(1:nmetaaryi(mybdl),mybdl),                       & 
        aryival=aryival(1:maxval(aryilen(1:nmetaaryi(mybdl),mybdl)),     &
          1:nmetaaryi(mybdl),mybdl) )
       if(rc/=0) print *,'nemsio_open, file=',trim(filename),' iret=',rc

    endif
!
!** collect data to first pe and write out data
!
    allocate(arrayr4(im,jm),tmp(im*jm))
    allocate(fcstField(fieldCount(mybdl)))
    call ESMF_FieldBundleGet(fieldbundle, fieldList=fcstField,     &
         itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc)
    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
        line=__LINE__, file=__FILE__)) return  ! bail out

    jrec = 1
    do i=1, fieldcount(mybdl)
!
       call ESMF_FieldGet(fcstField(i),typekind=typekind, name=fieldname, rc=rc)

       if( nfldlev(i,mybdl) == 1) then
         if( typekind == ESMF_TYPEKIND_R4) then
           call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr42d, rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__, file=__FILE__)) return  ! bail out
           istart = lbound(arrayr42d,1)
           iend   = ubound(arrayr42d,1)
           jstart = lbound(arrayr42d,2)
           jend   = ubound(arrayr42d,2)
           nlen   = (iend-istart+1) * (jend-jstart+1)
         elseif( typekind == ESMF_TYPEKIND_R8) then
           call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr82d, rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__, file=__FILE__)) return  ! bail out
           istart = lbound(arrayr82d,1)
           iend   = ubound(arrayr82d,1)
           jstart = lbound(arrayr82d,2)
           jend   = ubound(arrayr82d,2)
           nlen   = (iend-istart+1) * (jend-jstart+1)
           allocate( arrayr42d(istart:iend,jstart:jend))
!$omp parallel do default(none) shared(arrayr42d,arrayr82d,jstart,jend,istart,iend) private(m,n)
           do n=jstart,jend
             do m=istart,iend
               arrayr42d(m,n) = arrayr82d(m,n)
             enddo
           enddo
         endif
! send data to task 0
         call mpi_gather(nlen, 1, MPI_INTEGER, irecv(:), 1, MPI_INTEGER, 0, mpi_comm, rc)
         if(mype == 0) then
           idisp(1) = 0
           do n=1,ntasks-1
             idisp(n+1) = idisp(n) + irecv(n)
           enddo
!           if(mype==0) print *,' collect data, idisp=',idisp(:)
!           if(mype==0) print *,' collect data, irecv=',irecv(:)
         endif
!         if( trim(recname(jrec,mybdl))=="HGTsfc" .and. trim(recname(jrec,mybdl))=="sfc") then
!           print *,'in write nemsio,fb=',i,' write jrec=',jrec,' val=',maxval(arrayr42d(istart:iend,jstart:jend)),  &
!             minval(arrayr42d(istart:iend,jstart:jend)),maxloc(arrayr42d(istart:iend,jstart:jend)), &
!             minloc(arrayr42d(istart:iend,jstart:jend))
!         endif
         call mpi_gatherv(arrayr42d,nlen,MPI_REAL, arrayr4,irecv,idisp(:), MPI_REAL, &
              0, mpi_comm, rc)
         if(mype==0) then
!           print *,'in write nemsio, value=',maxval(arrayr4(1:im,1:jm)), &
!              minval(arrayr4(1:im,1:jm)),maxloc(arrayr4(1:im,1:jm)),minloc(arrayr4(1:im,1:jm))
           tmp = reshape(arrayr4, (/im*jm/))
           call nemsio_writerec(nemsiofile, jrec, tmp, iret=rc)
!           print *,'in write nemsio,fb=',i,' write jrec=',jrec,'fld is',  &
!             trim(recname(jrec,mybdl)), 'rc=', &
!             rc, 'value=',maxval(tmp),minval(tmp),maxloc(tmp),minloc(tmp)
         endif
         jrec = jrec + 1

      elseif (nfldlev(i,mybdl) > 1) then

         if( typekind == ESMF_TYPEKIND_R4) then
           call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr43d, rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__, file=__FILE__)) return  ! bail out

           istart = lbound(arrayr43d,1)
           iend   = ubound(arrayr43d,1)
           jstart = lbound(arrayr43d,2)
           jend   = ubound(arrayr43d,2)
           kstart = lbound(arrayr43d,3)
           kend   = ubound(arrayr43d,3)
           nlen   = (iend-istart+1) * (jend-jstart+1)
         elseif( typekind == ESMF_TYPEKIND_R8) then
           call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr83d, rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__, file=__FILE__)) return  ! bail out
           istart = lbound(arrayr83d,1)
           iend   = ubound(arrayr83d,1)
           jstart = lbound(arrayr83d,2)
           jend   = ubound(arrayr83d,2)
           kstart = lbound(arrayr83d,3)
           kend   = ubound(arrayr83d,3)
           nlen   = (iend-istart+1) * (jend-jstart+1)
         endif

         ! send data to task 0
         call mpi_gather(nlen, 1, MPI_INTEGER, irecv, 1, MPI_INTEGER, 0, mpi_comm, rc)
         if(mype == 0) then
           idisp(1) = 0
           do n=1,ntasks-1
             idisp(n+1) = idisp(n) + irecv(n)
           enddo
         endif
! write out all levels
         allocate(arrayr42d(istart:iend,jstart:jend))
!         do k=kstart,kend
         if ( write_nemsioflip ) then
           k1=kend   ; k2=kstart ; k3=-1
         else
           k1=kstart ; k2=kend   ; k3=1
         endif
         do k=k1,k2,k3
           if (typekind == ESMF_TYPEKIND_R4) then
!$omp parallel do default(none) shared(arrayr42d,arrayr43d,jstart,jend,istart,iend,k) private(m,n)
             do n=jstart,jend
               do m=istart,iend
                  arrayr42d(m,n) = arrayr43d(m,n,k)    
               enddo
             enddo
           elseif (typekind == ESMF_TYPEKIND_R8) then
!$omp parallel do default(none) shared(arrayr42d,arrayr83d,jstart,jend,istart,iend,k) private(m,n)
             do n=jstart,jend
               do m=istart,iend
                  arrayr42d(m,n) = arrayr83d(m,n,k)    
               enddo
             enddo
           endif
!
           call mpi_gatherv(arrayr42d, nlen, MPI_REAL, arrayr4, irecv,idisp, MPI_REAL, &
                            0, mpi_comm, rc)
           if(mype == 0) then
             tmp = reshape(arrayr4, (/im*jm/))
             call nemsio_writerec(nemsiofile, jrec, tmp, iret=rc)
             jrec = jrec + 1
           endif
         enddo
         deallocate(arrayr42d)
!
      endif
    enddo
!
    deallocate(tmp)
    deallocate(arrayr4)
    deallocate(fcstField)
!
!** close nemsio file
    call nemsio_close(nemsiofile, iret=rc)
!
    call nemsio_finalize()
!
!** ffsync
    if( write_fsyncflag ) then
      do n=751,900
        inquire(n,opened=opened)
        if(.not.opened)then
          nofsync=n
          exit
        endif
      enddo
!
      open(unit=nofsync, file=trim(FILENAME) )
        rc=FFSYNC(nofsync)
        if (rc.ne.0) then
          print *,"Error returned from ffsync, file=", trim(FILENAME),"rc=",RC
        ENDIF
      close(nofsync)
!
    endif
!
  end subroutine write_nemsio

!----------------------------------------------------------------------------------------

  subroutine write_nemaio_final()

!**
    if(allocated(lon1d))  deallocate(lon1d)
    if(allocated(lat1d))  deallocate(lat1d)
    if(allocated(vcoord)) deallocate(vcoord)
    if(allocated(idsl))   deallocate(idsl, idvc,idvm)
    deallocate(irecv)
    deallocate(idisp)
    deallocate(fieldcount)

  end subroutine write_nemaio_final
!
!----------------------------------------------------------------------------------------

  subroutine get_global_attr(fldbundle, mybdl, rc)
    type(ESMF_FieldBundle), intent(in) :: fldbundle
    integer, intent(in)                :: mybdl
    integer, intent(out)               :: rc

! local variable
   integer i,j,k,n,kz, attcount
   integer ni,naryi,nr4,nr8, nc
   integer aklen
   character(80) attName, hydrostatics, fldname
   type(ESMF_TypeKind_Flag)                :: typekind
   real(4), dimension(:), allocatable :: ak4,bk4
   real(8), dimension(:), allocatable :: ak8,bk8
!
! look at the field bundle attributes
    call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
                           attnestflag=ESMF_ATTNEST_OFF, Count=attcount, rc=rc)
    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
      line=__LINE__, &
      file=__FILE__)) &
      return  ! bail out

! first loop over all the attributes to find the count for integer attr, real
! attr, etc
    j=1
    k=1
    if (.not. allocated(nmetavari))  then
      allocate(nmetavari(nbdl),nmetavarr4(nbdl),nmetavarr8(nbdl),nmetavarc(nbdl))
      allocate(nmetaaryi(nbdl),nmetaaryr4(nbdl),nmetaaryr8(nbdl),nmetaaryc(nbdl))
      allocate(extrameta(nbdl))
    endif
    nmetavari(mybdl)=0; nmetavarr4(mybdl)=0; nmetavarr8(mybdl)=0; nmetavarc(mybdl)=0
    nmetaaryi(mybdl)=0; nmetaaryr4(mybdl)=0; nmetaaryr8(mybdl)=0; nmetaaryc(mybdl)=0
    aklen=0.
    do i=1, attCount

      call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
        attnestflag=ESMF_ATTNEST_OFF, attributeIndex=i,name=attName, typekind=typekind, &
        itemCount=n, rc=rc)
      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
        line=__LINE__, &
        file=__FILE__)) &
        return  ! bail out

! add this attribute to the list of transfers
      if (typekind==ESMF_TYPEKIND_CHARACTER) then
        if( n == 1) then
          nmetavarc(mybdl) = nmetavarc(mybdl) + 1
        else if (n > 1) then
          nmetaaryc(mybdl) = nmetaaryc(mybdl) + 1
        endif
      else if (typekind==ESMF_TYPEKIND_I4) then
        if( n == 1) then
          nmetavari(mybdl) = nmetavari(mybdl) + 1
        else if (n > 1) then
          nmetaaryi(mybdl) = nmetaaryi(mybdl) + 1
        endif
      else if (typekind==ESMF_TYPEKIND_R4) then
        if( n == 1) then
          nmetavarr4(mybdl) = nmetavarr4(mybdl) + 1
        else if (n > 1) then
          if( trim(attName) == "ak" .or. trim(attName) == "bk") then
            aklen = n
          else
            nmetaaryr4(mybdl) = nmetaaryr4(mybdl) + 1
          endif
        endif
      else if (typekind==ESMF_TYPEKIND_R8) then
        if( n == 1) then
          nmetavarr8(mybdl) = nmetavarr8(mybdl) + 1
        else if (n > 1) then
          if( trim(attName) == "ak" .or. trim(attName) == "bk") then
            aklen = n
          else
            nmetaaryr8(mybdl) = nmetaaryr8(mybdl) + 1
          endif
        endif
      endif
    enddo
    if(mybdl==2) nmetaaryi(mybdl) = nmetaaryi(mybdl) + 1
!    print *,'in get _global_attr, nmetavarc=',nmetavarc(mybdl),'nmetaaryc=',nmetaaryc(mybdl), &
!      'nmetavari=',nmetavari(mybdl),'nmetaaryi=',nmetaaryi(mybdl),'nmetavarr4=',nmetavarr4(mybdl), &
!      'nmetavarr8=',nmetavarr8(mybdl) , 'aklen=',aklen
!
! get value:
    if (nmetavari(mybdl) > 0) then
      if(.not.allocated(variname))   allocate(variname(100,nbdl),varival(100,nbdl))
    endif
    if (nmetavarr4(mybdl) > 0) then
      if(.not.allocated(varr4name))  allocate(varr4name(100,nbdl),varr4val(100,nbdl))
    endif
    if (nmetavarr8(mybdl) > 0) then
      if(.not.allocated(varr8name))  allocate(varr8name(100,nbdl),varr8val(100,nbdl))
    endif
    nmetavarc(mybdl) = nmetavarc(mybdl) + 2
    if (nmetavarc(mybdl) > 0)  then
      if(.not.allocated(varcname))  allocate(varcname(100,nbdl),varcval(100,nbdl))
    endif
    if (nmetaaryi(mybdl) > 0)  then
      if(.not.allocated(aryiname))  then
        allocate(aryiname(20,nbdl),aryilen(20,nbdl))
        allocate(aryival(2000,20,nbdl))
      endif
    endif
    if (aklen > 0)  then
      if(.not.allocated(vcoord))  allocate(vcoord(aklen,3,2))
      vcoord = 0.
    endif
!
    ni=0; nr4=0; nr8=0; nc=0; naryi=0
    do i=1, attCount

      call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
        attnestflag=ESMF_ATTNEST_OFF, attributeIndex=i, name=attName, &
        typekind=typekind, itemCount=n,  rc=rc)
       print *,'in write nemsio fist get att, att=',trim(attName),'n=',n

      if (typekind==ESMF_TYPEKIND_I4 ) then
        if(n==1) then
         ni = ni + 1
         variname(ni,mybdl) = trim(attName)
         call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
           name=trim(variname(ni,mybdl)), value=varival(ni,mybdl), rc=rc)
         if (trim(variname(ni,mybdl)) == 'ncnsto') ntrac=varival(ni,mybdl)
         if (trim(variname(ni,mybdl)) == 'ncld')   ncld=varival(ni,mybdl)
         if (trim(variname(ni,mybdl)) == 'nsoil')  nsoil=varival(ni,mybdl)
        else
         naryi = naryi + 1
         aryiname(naryi,mybdl) = trim(attName)
         aryilen(naryi,mybdl)  = n
         call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
              name=trim(attName), valueList=aryival(1:n,naryi,mybdl),rc=rc)
        endif
      else if (typekind==ESMF_TYPEKIND_R4) then
        if(n==1) then
         nr4 = nr4 + 1
         varr4name(nr4,mybdl) = trim(attName)
         call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
           name=trim(varr4name(nr4,mybdl)), value=varr4val(nr4,mybdl), rc=rc)
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
           line=__LINE__, &
           file=__FILE__)) &
           return  ! bail out
        else
          if(trim(attName) =="ak") then
            allocate(ak4(n))
            call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
              name=trim(attName), valueList=ak4, rc=rc)
            if( write_nemsioflip ) then
              kz=size(ak4)
              do k=1,kz
                vcoord(k,1,1) = ak4(kz-k+1)
              enddo
            else
              vcoord(:,1,1) = ak4(:)
            endif
            deallocate(ak4)
          else if(trim(attName) =="bk") then
            allocate(bk4(n))
            call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
              name=trim(attName), valueList=bk4, rc=rc)
            if( write_nemsioflip ) then
              kz=size(bk4)
              do k=1,kz
                vcoord(k,2,1) = bk4(kz-k+1)
              enddo
            else
              vcoord(:,2,1) = bk4(:)
            endif
            deallocate(bk4)
          endif
        endif
      else if (typekind==ESMF_TYPEKIND_R8) then
        if(n==1) then
         nr8 = nr8 + 1
         varr8name(nr8,mybdl) = trim(attName)
         call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
           name=trim(varr8name(nr8,mybdl)), value=varr8val(nr8,mybdl), rc=rc)
        else
          if(trim(attName) =="ak") then
            allocate(ak8(n))
            call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
              name=trim(attName), valueList=ak8, rc=rc)
            if( write_nemsioflip ) then
              kz=size(ak8)
              do k=1,kz
                vcoord(k,1,1) = ak8(kz-k+1)
              enddo
            else
              vcoord(:,1,1) = ak8(:)
            endif
            deallocate(ak8)
          else if(trim(attName) =="bk") then
            allocate(bk8(n))
            call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
              name=trim(attName), valueList=bk8, rc=rc)
            if( write_nemsioflip ) then
              kz=size(bk8)
              do k=1,kz
                vcoord(k,2,1) = bk8(kz-k+1)
              enddo
            else
              vcoord(:,2,1) = bk8(:)
            endif
            deallocate(bk8)
          endif
        endif
      else if (typekind==ESMF_TYPEKIND_CHARACTER) then
        if(n==1) then
         nc = nc + 1
         varcname(nc,mybdl) = trim(attName)
         call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", &
           name=trim(varcname(nc,mybdl)), value=varcval(nc,mybdl), rc=rc)
        endif
      endif

!      if(nmetavari(mybdl)>0) print *,'variname=',variname(1,mybdl),'varival=',varival(1,mybdl)
!      if(nmetavarc(mybdl)>0) print *,'varcname=',varcname(1,mybdl),'varcval=',varcval(1,mybdl)
!
    enddo
!
    if(write_nemsioflip) then
      nc = nc + 1
      varcname(nc,mybdl) = 'y-direction'
      varcval(nc,mybdl)  = 'north2south'
      nc = nc + 1
      varcname(nc,mybdl) = 'z-direction'
      varcval(nc,mybdl)  = 'bottom2top'
    else
      nc = nc + 1
      varcname(nc,mybdl) = 'y-direction'
      varcval(nc,mybdl)  = 'south2north'
      nc = nc + 1
      varcname(nc,mybdl) = 'z-direction'
      varcval(nc,mybdl)  = 'top2bottom'
    endif
!
!output lpl
    if( mybdl == 2 ) then
      naryi = naryi + 1
      aryiname(naryi,mybdl) = 'lpl'
      aryilen(naryi,mybdl)  = jm
      aryival(1:jm,naryi,mybdl) = im
    endif

    if( nmetavari(mybdl)>0 .or. nmetavarc(mybdl)>0 .or. nmetavarr4(mybdl) >0 .or. nmetavarr8(mybdl)>0 &
      .or. nmetaaryi(mybdl)>0) then
      extrameta(mybdl) = .true.
    else
      extrameta(mybdl) = .false.
    endif
!
  end subroutine get_global_attr
!
!----------------------------------------------------------------------------------------

end module module_write_nemsio