!======================================================================
! 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
!======================================================================
!
#include "cppdefs.h"
#ifdef PARALLEL_FILES
      program partit 
!
! Generic netCDF partitioning tool: reads netCDF files corresonding
! to the whole physical grid and prepares multiple files which hold
! data corresponding to different subdomains. These files can be
! then read in parallel by different MPI processes. 
!
! Usage:  partit NP_XI NP_ETA ncname1 ... ncnameN
! ------  where NP_XI  number of subdomains along XI-direction
!               NP_ETA number of subdomains along ETA-direction
!               ncname1 ... ncnameN  names of netCDF files
!
! Non-partitionable objects of netCDF files, such as scalar variables
! and attributes (both global and attributes to variables) are copied
! redundantly into the partitioned files, while partitionable array
! data is subdivided into subdomains and distributed among the
! partitioned files in such a manner that all files contain
! individual data without any overlap or redundantly stored data.
!
! The partitioning algorithm works as follows: The partitionable
! dimensions ('xi_rho', 'xi_u', 'eta_rho' and 'eta_v') are identified
! by name, then then their values are read and compared in pairs to
! detect if any of the directions have periodicity. It is assumed
! that ghost points corresponding to physical boundaries are stored
! in the file, but computational margins (including periodic margins) 
! are not. Consequently, if xi_rho and xi_u are equal to each other, 
! XI-direction is periodic, and if they differ by one, it is not.
! ETA-direction is treated similarly. Once periodicity type is
! determined, the internal number of internal points in each
! direction (i.e. excluding ghost points corresponding to physical
! boundaries) id divided by the number of subdomains in that
! direction and then physical boundary points are attached to
! subdomains which are adjacent to the boundaries. This results in
! slightly different dimension sizes of netCDF files corresponding
! to diffeent subdomains.
!
! Once all dimensions are sorted out, data corresponding to
! subdomains is extracted from the source file and copied into
! partial files.
!
      implicit none
      integer stdout, max_buff_size, maxdims, maxvars, maxnodes
      parameter (stdout=6,           max_buff_size=3000*2000*100,
     &           maxdims=40,         maxvars=128,     maxnodes=4000)
      character*120 ncname0, ncname(0:maxnodes-1), string 
      character*32 dimname(maxdims), varname(maxvars)
      logical XiPeriodic, EtaPeriodic, part_switch(maxvars),
     &                                           series(maxvars)
      real*8 buff(max_buff_size)
      integer narg,  NP_XI, NNODES,  xi_rho,  id_xi_rho,  id_xi_psi,
     &        arg,  NP_ETA,   node,  xi_u,    id_xi_u,    id_xi_v,
     &        iargc, subXI,     ii,  eta_rho, id_eta_rho, id_eta_psi,
     &        ierr, subETA,     jj,  eta_v,   id_eta_v,   id_eta_u,
     &        ndims, nvars, ngatts,  tsize,   unlimdimid, varatts,
     &        i,j,k,  lstr,   lvar,  lenstr,  rec,  size, ncid0, 
     &  ncid(0:maxnodes-1), dimid(maxdims),   dimsize(maxdims),
     &      varid(maxvars), vartype(maxvars), vardims(maxvars),
     &      start(maxdims), count(maxdims),   start1(maxdims),
     &      dimids(maxdims,maxvars),          ibuff(maxdims)
      common /partit_main/  buff,  ncid,  dimid,  dimsize, varid,
     &    vartype, vardims, start, count, start1, dimids,  ibuff
      integer LLm, MMm
      integer chunk_size_X, margin_X
      integer chunk_size_E, margin_E      
      integer istrmpi, jstrmpi
      integer iendmpi, jendmpi
      integer Lmmpi, Mmmpi
      
#include "netcdf.inc"
!
! Check how many arguments are given, complain about the error,
! if too few, otherwise extract NP_X and NP_E from the first two
! arguments.
!
      narg=iargc()
      if (narg .lt. 3) then
        write(stdout,'(/1x,A,1x,A/32x,A/)') 'Usage of partit',
     &    'should be:', 'partit NP_X NP_E ncname1 ... ncnameN'
        stop
      endif

      call getarg(1,string)
      lvar=lenstr(string)
      NP_XI=0
      do i=1,lvar
        j=ichar(string(i:i))-48
        if (j.ge.0 .and. j.le.9) then
          NP_XI=10*NP_XI+j
        else
          write(stdout,'(/8x,2(A,1x),2A/)') 'ERROR: illegal first',
     &     'argument', string(1:lvar), ', must be an integer number.'
          stop
        endif
      enddo

      call getarg(2,string)
      lvar=lenstr(string)
      NP_ETA=0
      do i=1,lvar
        j=ichar(string(i:i))-48
        if (j.ge.0 .and. j.le.9) then
          NP_ETA=10*NP_ETA+j
        else
          write(stdout,'(/8x,2(A,1x),2A/)') 'ERROR: illegal second',
     &      'argument', string(1:lvar), ', must be an integer number.'
          stop
        endif
      enddo
      NNODES=NP_XI*NP_ETA
      if (NNODES.gt.maxnodes) then
        write(stdout,'(/8x,A,1x,A,I3,1x,A/15x,A/)') 'ERROR:',
     &      'requested number of nodes',NNODES,'exceeds limit.',
     &      'Increase parameter maxnodes in file "partit.F".'
        stop
      endif
      write(stdout,'(/4x,2(4x,A,I3)/)') 'NP_XI =',  NP_XI,
     &                                  'NP_ETA =', NP_ETA
!
! Process netCDF files: open, determine if it is already
! a partitioned file, then make general inquiry. Complain
! about error if it already partitioned, or if number of
! variables and/or dimensions exceeds specified limits.
!
      do arg=3,narg
        call getarg(arg,ncname0)
        lstr=lenstr(ncname0)
        ierr=nf_open (ncname0(1:lstr), nf_nowrite, ncid0)
        if (ierr .ne. nf_noerr) then
          write(stdout,'(/8x,A,1x,A,1x,A,A/)') 'ERROR: Cannot',
     &                  'open netCDF file', ncname0(1:lstr),'.'
          goto 97     !--> next file
        endif
        ierr=nf_inq_att (ncid0, nf_global, 'partition', i,j)
        if (ierr .eq. nf_noerr) then
          write(stdout,'(/8x,3A/17x,A/)') 'WARNING: netCDF file ''',
     &                 ncname0(1:lstr), ''' is already partitioned',
     &                 'file. It cannot be partitioned any further.'
          goto 97     !--> next file
        endif

        write(stdout,'(8x,3A)') 'Processing netCDF file ''',
     &                               ncname0(1:lstr), '''...'

        ierr=nf_inq (ncid0, ndims, nvars, ngatts, unlimdimid)
        if (ierr .ne. nf_noerr) then
           write(stdout,'(/8x,A,1x,A/15x,A,1x,A,A/)') 'ERROR:',
     &        'Cannot determine number of dimensions, variables',
     &        'and attributes in netCDF file',ncname0(1:lstr),'.'
          goto 97     !--> next file
        elseif (ndims .gt. maxdims) then
          write(stdout,'(/8x,A,I4,1x,4A/15x,A,1x,A/)')
     &        'ERROR: number of dimensions', ndims,  'in netCDF',
     &        'file ''', ncname0(1:lstr), '''', 'exceeds limit.',
     &        'Increase parameter maxdims in file "partit.F".'
          goto 97     !--> next file
         elseif (nvars .gt. maxvars) then
          write(stdout,'(/8x,A,I4,1x,4A/15x,A,1x,A/)')
     &        'ERROR: number of variables',  nvars,  'in netCDF',
     &        'file ''', ncname0(1:lstr), '''', 'exceeds limit.',
     &        'Increase parameter maxvars in file "partit.F".'
          goto 97     !--> next file
        endif
!
! Sort out dimensions: For each dimension find and save its name and
! size. Then check whether all partitionable dimensions (identified
! by names 'xi_rho', 'xi_u', 'eta_rho' and 'eta_v')  are present and
! save their IDs and sizes. 
!
        tsize=1      ! <-- default value.
        do i=1,ndims
          ierr=nf_inq_dimname (ncid0, i, dimname(i))
          if (ierr .ne. nf_noerr) then
             write(stdout,'(/8x,A,I3/15x,3A/)')
     &           'ERROR: Cannot determine name for dimension ID =',
     &            i,  'in netCDF file ''',  ncname0(1:lstr),  '''.' 
             goto 97     !--> next file
          endif
          ierr=nf_inq_dimlen  (ncid0, i, dimsize(i))
          if (ierr .ne. nf_noerr) then
             lvar=lenstr(dimname(i))
             write(stdout,'(/8x,A,A,A/15x,3A/)')
     &            'ERROR: Cannot determine length of dimension ''',
     &             dimname(i)(1:lvar),  '''',  'in netCDF file ''',
     &                                       ncname0(1:lstr), '''.'
             goto 97     !--> next file
          endif
          if (i.eq. unlimdimid) then
            tsize=dimsize(i)
            dimsize(i)=nf_unlimited
          endif
        enddo
!
! Determine IDs and sizes of partitionable dimensions, 'xi_rho',
! 'xi_u', 'eta_rho' and 'eta_v'. Also save IDs of obsolete dimensions
! 'xi_psi', 'xi_v', 'eta_psi' and and 'eta_u'. These are used to
! readress obsolete dimensions according to the rules:
!
        id_xi_rho=0                     ! xi_psi  --> xi_u
        id_xi_u=0                       ! xi_v    --> xi_rho
        id_eta_rho=0                    ! eta_psi --> eta_v
        id_eta_v=0                      ! eta_u   --> eta_rho 
        id_xi_psi=0
        id_xi_v=0
        id_eta_psi=0
        id_eta_u=0
        do i=1,ndims
          lvar=lenstr(dimname(i))
          if (lvar.eq.6 .and. dimname(i)(1:lvar).eq.'xi_rho') then
            id_xi_rho=i
            xi_rho=dimsize(i)
          elseif (lvar.eq.4 .and. dimname(i)(1:lvar).eq.'xi_u') then
            id_xi_u=i
            xi_u=dimsize(i)
          elseif (lvar.eq.7.and.dimname(i)(1:lvar).eq.'eta_rho') then
            id_eta_rho=i
            eta_rho=dimsize(i)
          elseif (lvar.eq.5 .and. dimname(i)(1:lvar).eq.'eta_v') then
            id_eta_v=i
            eta_v=dimsize(i)
          elseif (lvar.eq.6 .and.dimname(i)(1:lvar).eq.'xi_psi') then
            id_xi_psi=i
          elseif (lvar.eq.4 .and. dimname(i)(1:lvar).eq.'xi_v') then
            id_xi_v=i
          elseif (lvar.eq.7.and.dimname(i)(1:lvar).eq.'eta_psi') then
            id_eta_psi=i
          elseif (lvar.eq.5 .and. dimname(i)(1:lvar).eq.'eta_u') then
            id_eta_u=i
          endif
c**       write(*,'(I3,1x,A,T16,I3)') i,dimname(i)(1:lvar),dimsize(i)
        enddo
        if (id_xi_rho.eq.0  .or. id_xi_u.eq.0 .or.
     &      id_eta_rho.eq.0 .or. id_eta_v.eq.0) then
          write(stdout,'(/8x,2A/15x,3A/)') 'ERROR: not all ',
     &            'partitionable dimensions are found',
     &            'in netCDF file ''', ncname0(1:lstr), '''.'
          goto 97     !--> next file 
        endif
!
! Determine subdomain dimensions. Here "subXI" and "subETA" are
! the nimbers of internal grid points within each subdomains (that
! is, excluding physical boundary points and computational margins). 
! The number of internal points in either direction for the the
! whole computational domain must be divisible by NP_XI and NP_ETA  
! respectively. If it cannot be divided, complain about the error
! and exit.
!
        LLm = xi_rho - 2
        MMm = eta_rho - 2
        if (xi_rho .eq. xi_u) then
          XiPeriodic=.true.
          subXI=xi_rho/NP_XI
          if (subXI*NP_XI .ne. xi_rho) then
            write(stdout,'(/8x,A,1x,A,I4,1x,A/15x,A,I3,1x,A/)')
     &         'ERROR: Cannot partition XI-direction:', 'xi_rho =',
     &          xi_rho,  'is',  'not divisible by NP_XI =',  NP_XI,
     &         'in XI-periodic case.'

            goto 97
          endif
        elseif (xi_rho .eq. xi_u+1) then 
          XiPeriodic=.false.
          subXI=(xi_rho-2)/NP_XI
          chunk_size_X=(LLm+NP_XI-1)/NP_XI
          margin_X=(NP_XI*chunk_size_X-LLm)/2
          
C          if (subXI*NP_XI .ne. xi_rho-2) then
C            write(stdout,'(/8x,A,1x,A,I4,1x,A/15x,A,I3,1x,A/)')
C     &         'ERROR: Cannot partition XI-direction:', 'xi_rho-2 =',
C     &          xi_rho-2,  'is',  'not divisible by NP_XI =',  NP_XI,
C     &         'in nonperiodic XI case.'
C            goto 97
C          endif
        else
          write(stdout,'(/8x,2A/)') 'ERROR: inconsistent ',
     &                'dimensions ''xi_rho'' and ''xi_u''.' 
          goto 97     !--> next file
        endif

        if (eta_rho .eq. eta_v) then
          EtaPeriodic=.true.
          subETA=eta_rho/NP_ETA
      
          if (subETA*NP_ETA .ne. eta_rho) then
            write(stdout,'(/8x,A,1x,A,I4,1x,A/15x,A,I3,1x,A/)')
     &         'ERROR: Cannot partition ETA-direction:', 'eta_rho =',
     &          eta_rho,  'is',  'not divisible by NP_ETA =', NP_ETA,
     &         'in ETA-periodic case.'

            goto 97
          endif
        elseif (eta_rho .eq. eta_v+1) then
          EtaPeriodic=.false.
          subETA=(eta_rho-2)/NP_ETA
          
          chunk_size_E=(MMm+NP_ETA-1)/NP_ETA
          margin_E=(NP_ETA*chunk_size_E-MMm)/2          

C          if (subETA*NP_ETA .ne. eta_rho-2) then
C            write(stdout,'(/8x,A,1x,A,I4,1x,A/15x,A,I3,1x,A/)')
C     &        'ERROR: Cannot partition ETA-direction:','eta_rho-2 =',
C     &         eta_rho-2, 'is',  'not divisible by NP_ETA =', NP_ETA,
C     &        'in nonperiodic ETA case.'
C            goto 97
C          endif
        else
          write(stdout,'(/8x,A,1x,A/)') 'ERROR: inconsistent',
     &                 'dimensions ''eta_rho'' and ''eta_v''.'
          goto 97     !--> next file
        endif
!
! Create partitioned files.
! ====== =========== ======
!
        do node=0,NNODES-1
          lstr=lenstr(ncname0)
          ncname(node)=ncname0
          ierr=0
          call insert_node (ncname(node), lstr, node, NNODES, ierr)
          if (ierr. ne. 0) goto 97     !--> next file
!          ierr=nf_create (ncname(node)(1:lstr),nf_clobber,ncid(node))
          ierr=nf_create(ncname(node)(1:lstr),
     &                   nf_64bit_offset,ncid(node))
          if (ierr .eq. nf_noerr) then
            write(stdout,'(12x,3A)') 'Created partitioned file ''',
     &                                  ncname(node)(1:lstr), '''.' 
          else
            write(stdout,'(/8x,A,1x,3A/)') 'ERROR: cannot create',
     &              'netCDF file ''', ncname(node)(1:lstr), '''.'
            goto 97     !--> next file
          endif
!
! Define dimensions of partitioned files.
!
          jj=node/NP_XI
          ii=node-jj*NP_XI
          
      istrmpi=1+ii*chunk_size_X-margin_X
      iendmpi=istrmpi+chunk_size_X-1
      istrmpi=max(istrmpi,1)
      iendmpi=min(iendmpi,LLm)

      jstrmpi=1+jj*chunk_size_E-margin_E
      jendmpi=jstrmpi+chunk_size_E-1
      jstrmpi=max(jstrmpi,1)
      jendmpi=min(jendmpi,MMm)     
 
      Lmmpi=iendmpi-istrmpi+1
      Mmmpi=jendmpi-jstrmpi+1
      
          do i=1,ndims
            size=dimsize(i)
            if (i .eq. id_xi_rho) then
              size=Lmmpi
              if (.not.XiPeriodic) then
                if (ii.eq.0      ) size=size+1
                if (ii.eq.NP_XI-1) size=size+1
              endif
            elseif (i .eq. id_xi_u) then
              size=Lmmpi
              if (.not.XiPeriodic) then
                if (ii.eq.NP_XI-1) size=size+1 
              endif
            elseif (i .eq. id_eta_rho) then
              size=Mmmpi
              if (.not.EtaPeriodic) then
                if (jj.eq.0       ) size=size+1 
                if (jj.eq.NP_ETA-1) size=size+1 
              endif
            elseif (i .eq. id_eta_v) then
              size=Mmmpi
              if (.not.EtaPeriodic) then
                if (jj.eq.NP_ETA-1) size=size+1 
              endif
            endif
            if (i.ne.id_xi_psi   .and.  i.ne.id_xi_v  .and. 
     &          i.ne.id_eta_psi  .and.  i.ne.id_eta_u)  then
              lvar=lenstr(dimname(i))
              ierr=nf_def_dim (ncid(node), dimname(i)(1:lvar),
     &                                         size, dimid(i))
              if (ierr .ne. nf_noerr) then
                write(stdout,'(/8x,4A/15x,A,I4,A)') 'ERROR: ',
     &            'Cannot define dimension ''', dimname(i)(1:lvar),
     &            '''.',    'netCDF ettor status =',       ierr, '.'
              endif
c**           write(*,'(2I3,4x,2I3,I4,1x,A)') i, dimid(i), ii,
c**  &                            jj, size, dimname(i)(1:lvar)
            else
              dimid(i)=0
            endif
          enddo
!
! WARNING!!! ...After this moment array dimid(1:ndims) contains
! the set of NEW dimension IDs. Since the four dimensions, 'xi_psi',
! 'eta_psi', 'xi_v' and 'eta_u' have been eliminated, dimid(i) does
! not correspond to the set of dimension IDs of the original file
! [which would be just dimid(i)=i], but it is rather different.
! Array dimid(1:ndims) will be used later to remap old dimension
! IDs into new ones, see the remapping procedure approximately 80
! lines below.
!
! Put global attribute 'partition' which identifies subdomain
! within the processor grid individually for each file.
!
          ibuff(1)=ii
          ibuff(2)=jj
          ibuff(3)=NP_XI
          ibuff(4)=NP_ETA
          ierr=nf_put_att_int (ncid(node), nf_global, 'partition',
     &                                           nf_int, 4, ibuff)
        enddo
!
! Copy global attributes
!
        do i=1,ngatts 
          ierr=nf_inq_attname (ncid0, nf_global, i, string) 
          if (ierr. eq. nf_noerr) then
            lvar=lenstr(string)
            do node=0,NNODES-1
              ierr=nf_copy_att (ncid0, nf_global, string(1:lvar),
     &                                     ncid(node), nf_global)
              if (ierr. ne. nf_noerr) then
                lstr=lenstr(ncname(node))
                write(stdout,'(/8x,4A/15x,3A/)')  'ERROR: Cannot ',
     &            'copy global attribute ''', string(1:lvar), '''',
     &            'into netCDF file ''', ncname(node)(1:lstr),'''.'
                goto 97
              endif
            enddo
          else
            lstr=lenstr(ncname(0))
            write(stdout,'(/8x,2A,I3/15x,3A/)') 'ERROR: Cannot ',
     &         'determine mame of global attribute with ID =', i,
     &         'from netCDF file ''',    ncname0(1:lstr),   '''.'
            goto 97
          endif
        enddo
!
! Define variables and their attributes.
!
        do i=1,nvars
          ierr=nf_inq_var (ncid0,   i, varname(i),  vartype(i),
     &                      vardims(i), dimids(1,i),   varatts)
!
! Readress obsolete dimensions, if any:
!
          do j=1,vardims(i)
            if (dimids(j,i).eq.id_xi_psi) then
              dimids(j,i)=id_xi_u
            elseif (dimids(j,i).eq.id_xi_v) then
              dimids(j,i)=id_xi_rho
            elseif (dimids(j,i).eq.id_eta_psi) then
              dimids(j,i)=id_eta_v
            elseif (dimids(j,i).eq.id_eta_u) then
              dimids(j,i)=id_eta_rho
            endif
          enddo
!
! Determine whether partitionable dimensions or unlimited dimension
! are present for this variable. 
!
          series(i)=.false.
          part_switch(i)=.false.
          do j=1,vardims(i)
            if (dimids(j,i).eq.id_xi_rho .or.
     &          dimids(j,i).eq.id_xi_u    .or.
     &          dimids(j,i).eq.id_eta_rho .or.
     &          dimids(j,i).eq.id_eta_v) then
              part_switch(i)=.true. 
            elseif (dimids(j,i).eq.unlimdimid) then
              series(i)=.true. 
            endif
          enddo
!
! WARNING: Since dimids(1:vardims(i),i) contains dimension IDs
! corresponding to the set of IDs of the ORIGINAL file, and since
! some of the original dimensions were eliminated (merged), the
! set of dimension IDs in the NEW definitions is obtained by 
! inverse mapping of dimids(j,i) onto ibuff(j) using dimid(k) as
! a mapping array.  
!
          do j=1,vardims(i)
            do k=1,ndims
              if (dimids(j,i).eq.k) ibuff(j)=dimid(k)
            enddo
          enddo
c**       write(*,*) 'old_dimids:', (dimids(j,i),j=1,vardims(i))
c**       write(*,*) 'new_dimids:',    (ibuff(j),j=1,vardims(i))

          lvar=lenstr(varname(i))
          do node=0,NNODES-1
            ierr=nf_def_var (ncid(node), varname(i)(1:lvar),
     &              vartype(i), vardims(i), ibuff, varid(i))
          enddo
c**       write(stdout,'(I3,1x,A,T20,I3,1x,L1,1x,L1)') i,
c**  &                        varname(i)(1:lvar), vardims(i),
c**  &                            part_switch(i),  series(i)
          do j=1,varatts
            ierr=nf_inq_attname (ncid0, varid(i), j, string)
            lvar=lenstr(string)
            do node=0,NNODES-1
            ierr=nf_copy_att (ncid0, i, string(1:lvar),
     &                            ncid(node), varid(i))
            enddo
          enddo
        enddo
!
! Leave definition mode
!
        do node=0,NNODES-1
c          ierr=nf_set_fill (ncid(node), nf_nofill, i)
          ierr=nf_enddef(ncid(node))
        enddo
!
! Transfer variables into newly created files. 
!
        do rec=1,tsize
          if (tsize.gt.1) write(stdout,'(16x,A,I4,A)')
     &                 'Processing record', rec, '...'
          do i=1,nvars
            if (series(i) .or. rec.eq.1) then
              if (.not.part_switch(i) .and. .not.series(i)) then
!
! Scalar (zero-dimensional) variables:
!
                if (vartype(i) .eq. nf_char) then
                  ierr=nf_get_var_text (ncid0, i, buff)
                elseif (vartype(i) .eq. nf_int) then
                  ierr=nf_get_var_int    (ncid0, i, buff)
                elseif (vartype(i) .eq. nf_real) then
                  ierr=nf_get_var_real   (ncid0, i, buff)
                elseif (vartype(i) .eq. nf_double) then
                  ierr=nf_get_var_double (ncid0, i, buff)
                else
                  lvar=lenstr(varname(i))
                  write(stdout,'(/8x,4A/)') 'ERROR: scalar variable',
     &              ' ''', varname(i)(1:lvar), ''' has unknown type.'
                  goto 97
                endif
                if (ierr .eq. nf_noerr) then
                  do node=0,NNODES-1
                    if (vartype(i) .eq. nf_char) then
                     ierr=nf_put_var_text (ncid(node),varid(i),buff)
                    elseif (vartype(i) .eq. nf_int) then
                     ierr=nf_put_var_int   (ncid(node),varid(i),buff)
                    elseif (vartype(i) .eq. nf_real) then
                     ierr=nf_put_var_real  (ncid(node),varid(i),buff)
                    elseif (vartype(i) .eq. nf_double) then
                     ierr=nf_put_var_double(ncid(node),varid(i),buff)
                    endif
                    if (ierr .ne. nf_noerr) then
                      lvar=lenstr(varname(i))
                      lstr=lenstr(ncname(node))
                      write(stdout,'(/8x,3A/15x,3A,I4,A/)')
     &                    'ERROR: Cannot write scalar variable ''',
     &                     varname(i)(1:lvar), ''' into netCDF file',
     &                    '''',  ncname(node)(1:lstr), 
     &                    '''.     netCDF error code =',  ierr,   '.'
                      goto 97
                    endif
                  enddo
                else
                  lvar=lenstr(varname(i))
                  write(stdout,'(/8x,4A/)') 'ERROR: Cannot read ',
     &             'scalar variable ''', varname(i)(1:lvar), '''.'
                  goto 97
                endif
              elseif (.not.part_switch(i)) then
!
! Non-partitionable array.
!
                size=1
                do j=1,vardims(i) 
                  if (dimids(j,i).eq.unlimdimid) then
                    start(j)=rec
                    count(j)=1
                  else
                    start(j)=1
                    count(j)=dimsize(dimids(j,i))
                  endif 
                  size=size*count(j)
                enddo
                if (vartype(i) .eq. nf_char) then
                  size=size*1
                elseif (vartype(i) .eq. nf_int) then
                  size=size*4
                elseif (vartype(i) .eq. nf_real) then
                  size=size*4
                elseif (vartype(i) .eq. nf_double) then
                  size=size*8
                else
                  lvar=lenstr(varname(i))
                  write(stdout,'(/8x,3A/)') 'ERROR: variable ''',
     &                 varname(i)(1:lvar), ''' has unknown type.'
                  goto 97
                endif
                if (size .gt. 8*max_buff_size) then
                  write(stdout,'(/8x,A,3(/15x,A,I10,1x,A)/)')
     &              'ERROR: unsufficient buffer size in "partit.F":',
     &              'requested:',         size,      'Bytes,',
     &              'available:',   8*max_buff_size, 'Bytes.',
     &              'Increase parameter max_buff_size and recompile.'
                  goto 97
                endif

                if (vartype(i) .eq. nf_char) then
                  ierr=nf_get_vara_text   (ncid0, i, start,
     &                                          count, buff)
                elseif (vartype(i) .eq. nf_int) then
                  ierr=nf_get_vara_int    (ncid0, i, start,
     &                                         count, buff)
                elseif (vartype(i) .eq. nf_real) then
                  ierr=nf_get_vara_real   (ncid0, i, start,
     &                                         count, buff)
                elseif (vartype(i) .eq. nf_double) then
                  ierr=nf_get_vara_double (ncid0, i, start,
     &                                         count, buff)
                endif
                if (ierr .eq. nf_noerr) then
                  do node=0,NNODES-1
                    if (vartype(i) .eq. nf_char) then
                      ierr=nf_put_vara_text   (ncid(node), varid(i),
     &                                           start, count, buff)
                    elseif (vartype(i) .eq. nf_int) then
                      ierr=nf_put_vara_int    (ncid(node), varid(i),
     &                                           start, count, buff)
                    elseif (vartype(i) .eq. nf_real) then
                      ierr=nf_put_vara_real   (ncid(node), varid(i),
     &                                           start, count, buff)
                    elseif (vartype(i) .eq. nf_double) then
                      ierr=nf_put_vara_double (ncid(node), varid(i),
     &                                           start, count, buff)
                    endif
                    if (ierr .ne. nf_noerr) then
                      lvar=lenstr(varname(i))
                      lstr=lenstr(ncname(node))
                      write(stdout,'(/8x,3A,I3/15x,3A,I4,A/)')
     &                  'ERROR: Cannot write variable ''',
     &                   varname(i)(1:lvar),''' for time record',rec,
     &                  'into netCDF file ''',  ncname(node)(1:lstr),
     &                  '''. netCDF error code =', ierr, '.'
                      goto 97
                    endif
                  enddo
                else
                  lstr=lenstr(ncname0)
                  lvar=lenstr(varname(i))
                  write(stdout,'(/8x,4A,I3/15x,3A,I4/)') 'ERROR: ',
     &              'Cannot read variable ''',  varname(i)(1:lvar),
     &              ''' for time record',rec,'from netCDF file ''',
     &              ncname0(1:lstr),'''. netCDF error code =',ierr
                  goto 97
                endif
              elseif (part_switch(i)) then
!
! Partitioned array: 
!
                do node=0,NNODES-1
                  jj=node/NP_XI
                  ii=node-jj*NP_XI
                  
      istrmpi=1+ii*chunk_size_X-margin_X
      iendmpi=istrmpi+chunk_size_X-1
      istrmpi=max(istrmpi,1)
      iendmpi=min(iendmpi,LLm)

      jstrmpi=1+jj*chunk_size_E-margin_E
      jendmpi=jstrmpi+chunk_size_E-1
      jstrmpi=max(jstrmpi,1)
      jendmpi=min(jendmpi,MMm)     
 
      Lmmpi=iendmpi-istrmpi+1
      Mmmpi=jendmpi-jstrmpi+1
      
                  size=1
                  do j=1,vardims(i)
                    if (dimids(j,i).eq.id_xi_rho) then
                      start(j)=istrmpi
                      count(j)=Lmmpi
                      if (.not.XiPeriodic) then
                        if (ii.gt.0      ) start(j)=start(j)+1
                        if (ii.eq.0      ) count(j)=count(j)+1
                        if (ii.eq.NP_XI-1) count(j)=count(j)+1
                      endif
                      start1(j)=1

                    elseif (dimids(j,i).eq.id_xi_u) then
                      start(j)=istrmpi
                      count(j)=Lmmpi
                      if (.not.XiPeriodic) then
                        if (ii.eq.NP_XI-1) count(j)=count(j)+1
                      endif
                      start1(j)=1

                    elseif (dimids(j,i).eq.id_eta_rho) then
                      start(j)=jstrmpi
                      count(j)=Mmmpi
                      if (.not.EtaPeriodic) then
                        if (jj.gt.0       ) start(j)=start(j)+1
                        if (jj.eq.0       ) count(j)=count(j)+1
                        if (jj.eq.NP_ETA-1) count(j)=count(j)+1
                      endif
                      start1(j)=1

                    elseif (dimids(j,i).eq.id_eta_v) then
                      start(j)=jstrmpi
                      count(j)=Mmmpi
                      if (.not.EtaPeriodic) then
                        if (jj.eq.NP_ETA-1) count(j)=count(j)+1
                      endif
                      start1(j)=1

                    elseif (dimids(j,i).eq.unlimdimid) then
                      start(j)=rec
                      count(j)=1
                      start1(j)=rec
                    else
                      start(j)=1
                      count(j)=dimsize(dimids(j,i))
                      start1(j)=1
                    endif
                    size=size*count(j)
                  enddo
c**               write(*,*) 'dimids:', (dimids(j,i),j=1,vardims(i))
c**               write(*,*) ' start:',    (start(j),j=1,vardims(i))
c**               write(*,*) ' count:',    (count(j),j=1,vardims(i))


                  if (vartype(i) .eq. nf_char) then
                    size=size*1
                  elseif (vartype(i) .eq. nf_int) then
                    size=size*4
                  elseif (vartype(i) .eq. nf_real) then
                    size=size*4
                  elseif (vartype(i) .eq. nf_double) then
                    size=size*8
                  else
                    lvar=lenstr(varname(i))
                    write(stdout,'(/8x,4A/)') 'ERROR: variable ''',
     &                   varname(i)(1:lvar), ''' has unknown type.'
                    goto 97
                  endif
                  if (size .gt. 8*max_buff_size) then
                    write(stdout,'(/8x,A,3(/15x,A,I10,1x,A)/)')
     &              'ERROR: unsufficient buffer size in "partit.F":',
     &              'requested:',         size,      'Bytes,',
     &              'available:',   8*max_buff_size, 'Bytes.',
     &              'Increase parameter max_buff_size and recompile.'
                    goto 97
                  endif

                  if (vartype(i) .eq. nf_char) then
                    ierr=nf_get_vara_text   (ncid0, i, start,
     &                                           count, buff)
                  elseif (vartype(i) .eq. nf_int) then
                    ierr=nf_get_vara_int    (ncid0, i, start,
     &                                           count, buff)
                  elseif (vartype(i) .eq. nf_real) then
                    ierr=nf_get_vara_real   (ncid0, i, start,
     &                                           count, buff)
                  elseif (vartype(i) .eq. nf_double) then
                    ierr=nf_get_vara_double (ncid0, i, start,
     &                                           count, buff)
                  endif

                  if (ierr .eq. nf_noerr) then
                    if (vartype(i) .eq. nf_char) then
                      ierr=nf_put_vara_text   (ncid(node), varid(i),
     &                                          start1, count, buff)
                    elseif (vartype(i) .eq. nf_int) then
                      ierr=nf_put_vara_int    (ncid(node), varid(i),
     &                                          start1, count, buff)
                    elseif (vartype(i) .eq. nf_real) then
                      ierr=nf_put_vara_real   (ncid(node), varid(i),
     &                                          start1, count, buff)
                    elseif (vartype(i) .eq. nf_double) then
                      ierr=nf_put_vara_double (ncid(node), varid(i),
     &                                          start1, count, buff)
                    endif
                    if (ierr .ne. nf_noerr) then
                      lvar=lenstr(varname(i))
                      lstr=lenstr(ncname(node))
                      write(stdout,'(/8x,3A,I3/15x,3A,I4,A/)')
     &                  'ERROR: Cannot write partitioned array ''',
     &                   varname(i)(1:lvar),''' for time record',rec,
     &                  'into netCDF file ''',  ncname(node)(1:lstr),
     &                  '''. netCDF error code =', ierr, '.'
                      goto 97
                    endif
                  else
                    lstr=lenstr(ncname0)
                    lvar=lenstr(varname(i))
                    write(stdout,'(/8x,3A,I3/15x,3A,I4,A/)')
     &                  'ERROR: Cannot read partitioned array ''',
     &                   varname(i)(1:lvar),   ''' for time record',
     &                   rec, 'from netCDF file ''',ncname0(1:lstr),
     &                  '''. netCDF error code =',   ierr,    '.'

                    goto 97
                  endif
                enddo       ! <-- node=0,NNODES-1
              endif
            endif       ! <--series(i) .or. rec.eq.1
          enddo       ! <-- i=1,nvars
        enddo       ! <-- rec=1,tsize
!
! Close all netCDF files
!
  97    ierr=nf_close (ncid0)
        do node=0,NNODES-1
         ierr=nf_close (ncid(node))
        enddo
      enddo
      stop
      end

#else
      program partit_empty
      stop
      end
#endif /* PARALLEL_FILES */

