!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#ifdef WRITER
      integer function nf_fwrite (A, ncid, varid, record, type)
#else
      integer function nf_fread  (A, ncid, varid, record, type)
# include "cppdefs.h"
#endif
!
! Read/Write a floating point array from/into an input/output
! NetCDF file.
!
! Arguments: A       real array of standard horizontal dimensions
!                             which is to be read or written.
!            ncid    NetCDF ID of in the file.
!            varid   variable ID of that variable in NetCDF file.
!            record  record number.
!            type    type of the grid (RHO-, U, V, W, PSI etc.)
!
! Because significant portion of the code calculates starting and
! stopping indices for the subarray (which are exactly the same for
! both read and write operations) the code has been unified to
! guarantee that both functions are exactly adjoint.
!
#if defined MUSTANG
      USE comMUSTANG, only: nk_nivsed_out
#endif
#if defined BIOLOGY && defined PISCES && defined key_sediment
      USE par_sed, only: jpksed
#endif
      USE buffer, only: buff, max_buffer_N_size
      implicit none
#include "param.h"
      real A(GLOBAL_2D_ARRAY,max_buffer_N_size)

      integer ncid,  type,   vert_type,  imin, imax,  start(4)
     &      , varid, record, horiz_type, jmin, jmax,  count(4)
     &      , i,j,k, shift,  ierr

#if defined MPI
# include "scalars.h"
#endif
#include "netcdf.inc"
                                      ! Decode grid type into
      vert_type=type/4                ! vertical and horizontal
      horiz_type=type-4*vert_type     ! grid types, then calculate
      jmin=horiz_type/2               ! starting indices indices
      imin=horiz_type-2*jmin          ! in horizontal directions.

      start(1)=1
      start(2)=1
#if defined AGRIF && defined AGRIF_ADAPTIVE && !defined WRITER
      if (.Not.Agrif_Root()) then
        start(1)=1+(Agrif_ix()-1)*Agrif_irhox()
        start(2)=1+(Agrif_iy()-1)*Agrif_irhoy()
      endif
#endif

#ifdef MPI

# ifdef PARALLEL_FILES
!#  ifdef EW_PERIODIC
!      imin=1
!      imax=Lm
!#  else
      if (ii.gt.0) imin=1
      if (ii.eq.NP_XI-1) then
        imax=Lmmpi+1
      else
        imax=Lmmpi
      endif
!#  endif
!#  ifdef NS_PERIODIC
!      jmin=1
!      jmax=Mm
!#  else
      if (jj.gt.0) jmin=1
      if (jj.eq.NP_ETA-1) then
        jmax=Mmmpi+1
      else
        jmax=Mmmpi
      endif
!#  endif
# else
      if (ii.gt.0) then
        start(1)=1-imin+iminmpi
        imin=1
      endif
      if (ii.eq.NP_XI-1) then
        imax=Lmmpi+1
      else
        imax=Lmmpi
      endif
      if (jj.gt.0) then
        start(2)=1-jmin+jminmpi
        jmin=1
      endif
      if (jj.eq.NP_ETA-1) then
        jmax=Mmmpi+1
      else
        jmax=Mmmpi
      endif
# endif
#else
      imax=Lm+1
      jmax=Mm+1
#endif


!! If MPI domain voisin is set to null we have to expand
!! the zone read in netcdf file.

      count(1)=imax-imin+1
      count(2)=jmax-jmin+1

c**   write(stdout,'(1x,A,i4,1x,A,i2,2(3x,A,I2,2x,A,I3,2x,A,I3))')
c**  &      'NF_READ/WRITE: mynode=',mynode,'horiz_grid',horiz_type,
c**  &                         'ii=',ii, 'imin=',imin, 'imax=',imax,
C**  &                         'jj=',jj, 'jmin=',jmin, 'jmax=',jmax


      if (vert_type.eq.0) then    ! Sort out vertical grids.
        count(3)=1                !--> 2D variables
        start(3)=record
      elseif (vert_type.eq.1) then
        count(3)=N
        count(4)=1                !--> 3D RHO-grid
        start(3)=1
        start(4)=record
      elseif (vert_type.eq.2) then
        count(3)=N+1
        count(4)=1                !--> 3D W-grid
        start(3)=1
        start(4)=record
#ifdef SEDIMENT
      elseif (vert_type.eq.3) then
        count(3)=NLAY
        count(4)=1                !--> 3D BED-grid
        start(3)=1
        start(4)=record
#endif
#ifdef MUSTANG
      elseif (vert_type.eq.3) then
        count(3)=nk_nivsed_out
        count(4)=1                !--> 3D BED-grid
        start(3)=1
        start(4)=record
#endif
#ifdef ABL1D
      elseif (vert_type.eq.4) then
        count(3)=N_abl
        count(4)=1                !--> 3D ABL-grid
        start(3)=1
        start(4)=record
#endif
#if defined BIOLOGY && defined PISCES && defined key_sediment
      elseif (vert_type.eq.5) then
        count(3)=jpksed
        count(4)=1                !--> 3D DIAGENITIC-grid
        start(3)=1
        start(4)=record
#endif
      else
#ifndef WRITER
        write(*,'(/1x,2A,I3/)') 'NF_FREAD ERROR: ',
     &                'illegal grid type', type
        nf_fread=nf_noerr+1
#else
        write(*,'(/1x,2A,I4/)') 'NF_FWRITE ERROR: ',
     &                'illegal grid type', type
        nf_fwrite=nf_noerr+1
#endif
        return
      endif

!
! Read/Write array from the disk.
!===== ===== ===== ==== === =====
!
#ifndef WRITER
# ifdef ABL1D
      if (vert_type.eq.4) then
        ierr=nf_get_vara_FTYPE (ncid, varid, start, count, buff(1))
      else
        ierr=nf_get_vara_FTYPE (ncid, varid, start, count, buff(1))
      endif
# else
      ierr=nf_get_vara_FTYPE (ncid, varid, start, count, buff(1))
# endif
      nf_fread=ierr
      if (ierr .ne. nf_noerr) then
        write(*,'(/1x,2A,I5,1x,A,I4/)') 'NF_FREAD ERROR: ',
     &             'nf_get_vara netCDF error code =', ierr
     &              MYID
        return
      endif
#endif
      do k=1,count(3)
        do j=jmin,jmax
          shift=1-imin+count(1)*(j-jmin+(k-1)*count(2))
          do i=imin,imax
#ifdef WRITER
# ifdef ABL1D
            if (vert_type.eq.4) then
              buff(i+shift)=A(i,j,k)
            else
              buff(i+shift)=A(i,j,k)
            endif
# else
            buff(i+shift)=A(i,j,k)
# endif
#else
# ifdef ABL1D
            if (vert_type.eq.4) then
              A(i,j,k)=buff(i+shift)
            else
              A(i,j,k)=buff(i+shift)
            endif
# else
            A(i,j,k)=buff(i+shift)
# endif
#endif
          enddo
        enddo
      enddo
#ifndef WRITER
!$acc update if(compute_on_device)
!$acc&  device(A(imin:imax,jmin:jmax,1:count(3)))
#endif
#ifdef WRITER
# ifdef ABL1D
      if (vert_type.eq.4) then
        ierr=nf_put_vara_FTYPE (ncid, varid, start, count, buff(1))
      else
        ierr=nf_put_vara_FTYPE (ncid, varid, start, count, buff(1))
      endif
# else
      ierr=nf_put_vara_FTYPE (ncid, varid, start, count, buff(1))
# endif
      nf_fwrite=ierr
      if (ierr.ne.nf_noerr) then
        write(*,'(/1x,2A,I5,1x,A,I4/)') 'NF_FWRITE ERROR: ',
     &              'nf_put_vara netCDF error code =', ierr
     &               MYID
      endif
#endif

!
! Exchange periodic and computational margins (reader only).
!
#ifdef MPI
# define LOCALLM Lmmpi
# define LOCALMM Mmmpi
#else
# define LOCALLM Lm
# define LOCALMM Mm
#endif
#ifndef WRITER
# if defined EW_PERIODIC || defined NS_PERIODIC  || defined MPI
      if (horiz_type.eq.0 .and. vert_type.eq.0) then
        call exchange_r2d_tile (1,LOCALLM,1,LOCALMM, A)
      elseif (horiz_type.eq.1 .and. vert_type.eq.0) then
        call exchange_u2d_tile (1,LOCALLM,1,LOCALMM, A)
      elseif (horiz_type.eq.2 .and. vert_type.eq.0) then
        call exchange_v2d_tile (1,LOCALLM,1,LOCALMM, A)
      elseif (horiz_type.eq.3 .and. vert_type.eq.0) then
        call exchange_p2d_tile (1,LOCALLM,1,LOCALMM, A)
#  ifdef SOLVE3D
      elseif (horiz_type.eq.0 .and. vert_type.eq.1) then
#   ifdef THREE_GHOST_POINTS_TS
        call exchange_r3d_3pts_tile (1,LOCALLM,1,LOCALMM, A)
#   else
        call exchange_r3d_tile (1,LOCALLM,1,LOCALMM, A)
#   endif
      elseif (horiz_type.eq.1 .and. vert_type.eq.1) then
#   ifdef THREE_GHOST_POINTS_UV
        call exchange_u3d_3pts_tile (1,LOCALLM,1,LOCALMM, A)
#   else
        call exchange_u3d_tile (1,LOCALLM,1,LOCALMM, A)
#   endif
      elseif (horiz_type.eq.2 .and. vert_type.eq.1) then
#   ifdef THREE_GHOST_POINTS_UV
        call exchange_v3d_3pts_tile (1,LOCALLM,1,LOCALMM, A)
#   else
        call exchange_v3d_tile (1,LOCALLM,1,LOCALMM, A)
#   endif
      elseif (horiz_type.eq.3 .and. vert_type.eq.1) then
        call exchange_p3d_tile (1,LOCALLM,1,LOCALMM, A)

      elseif (horiz_type.eq.0 .and. vert_type.eq.2) then
        call exchange_w3d_tile (1,LOCALLM,1,LOCALMM, A)
#   ifdef ABL1D
      elseif (horiz_type.eq.0 .and. vert_type.eq.4) then
      do k=1, N_abl
        call exchange_r2d_tile (1,LOCALLM,1,LOCALMM, A(:,:,k))
      enddo
#   endif
CC      elseif (horiz_type.eq.0 .and. vert_type.eq.3) then
CC        call exchange_b3d_tile (1,LOCALLM,1,LOCALMM, A)
#  endif
      endif
!$acc update if(compute_on_device) host(A(:,:,1:count(3)))
# endif
#endif
      return
      end

#ifndef WRITER
# define WRITER
# include "nf_fread.F"
#endif
