C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mdsio/mdsio_write_field.F,v 1.1 2006/12/29 05:41:27 jmc Exp $
C $Name:  $

#include "MDSIO_OPTIONS.h"

CBOP
C !ROUTINE: MDS_WRITE_FIELD
C !INTERFACE:
      SUBROUTINE MDS_WRITE_FIELD(
     I   fName,
     I   filePrec,
     I   globalFile,
     I   useCurrentDir,
     I   arrType,
     I   zSize,nNz,
     I   arr,
     I   jrecord,
     I   myIter,
     I   myThid )

C !DESCRIPTION:
C Arguments:
C
C fName     (string)  :: base name for file to write
C filePrec  (integer) :: number of bits per word in file (32 or 64)
C globalFile (logical):: selects between writing a global or tiled file
C useCurrentDir(logic):: always write to the current directory (even if
C                        "mdsioLocalDir" is set)
C arrType   (char(2)) :: declaration of "arr": either "RS" or "RL"
C zSize     (integer) :: size of third dimension: normally either 1 or Nr
C nNz       (integer) :: number of vertical levels to write
C arr       ( RS/RL ) :: array to write, arr(:,:,zSize,:,:)
C irecord   (integer) :: record number to write
C myIter    (integer) :: time step number
C myThid    (integer) :: thread identifier
C
C MDS_WRITE_FIELD creates either a file of the form "fName.data" and
C "fName.meta" if the logical flag "globalFile" is set true. Otherwise
C it creates MDS tiled files of the form "fName.xxx.yyy.data" and
C "fName.xxx.yyy.meta". If jrecord > 0, a meta-file is created.
C Currently, the meta-files are not read because it is difficult
C to parse files in fortran. We should read meta information before
C adding records to an existing multi-record file.
C The precision of the file is decsribed by filePrec, set either
C to floatPrec32 or floatPrec64. The precision or declaration of
C the array argument must be consistently described by the char*(2)
C string arrType, either "RS" or "RL". nNz allows for both 2-D and
C 3-D arrays to be handled. nNz=1 implies a 2-D model field and
C nNz=Nr implies a 3-D model field. irecord=|jrecord| is the record number
C to be written and must be >= 1. NOTE: It is currently assumed that
C the highest record number in the file was the last record written.
C Nor is there a consistency check between the routine arguments and file.
C ie. If your write record 2 after record 4 the meta information
C will record the number of records to be 2. This, again, is because
C we have read the meta information. To be fixed.
C
C Created: 03/16/99 adcroft@mit.edu
C Changed: 01/06/02 menemenlis@jpl.nasa.gov
C          added useSingleCpuIO hack
C changed:  1/23/04 afe@ocean.mit.edu
C          added exch2 handling -- yes, the globalfile logic is nuts
CEOP

C !USES:
      IMPLICIT NONE
C Global variables / common blocks
#include "SIZE.h"
#include "EEPARAMS.h"
#include "EESUPPORT.h"
#include "PARAMS.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_TOPOLOGY.h"
#include "W2_EXCH2_PARAMS.h"
#endif /* ALLOW_EXCH2 */
#include "MDSIO_SCPU.h"

C !INPUT PARAMETERS:
      CHARACTER*(*) fName
      INTEGER filePrec
      LOGICAL globalFile
      LOGICAL useCurrentDir
      CHARACTER*(2) arrType
      INTEGER zSize, nNz
cph(
cph      Real arr(*)
      _RL arr(1-oLx:sNx+oLx,1-oLy:sNy+oLy,zSize,nSx,nSy)
cph)
      INTEGER jrecord
      INTEGER myIter
      INTEGER myThid
C !OUTPUT PARAMETERS:

C !FUNCTIONS
      INTEGER  ILNBLNK
      INTEGER  MDS_RECLEN
      LOGICAL  MASTER_CPU_IO
      EXTERNAL ILNBLNK
      EXTERNAL MDS_RECLEN
      EXTERNAL MASTER_CPU_IO

C !LOCAL VARIABLES:
      CHARACTER*(MAX_LEN_FNAM) dataFName,metaFName,pfName
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      LOGICAL fileIsOpen
      LOGICAL iAmDoingIO
      LOGICAL writeMetaF
      INTEGER irecord
      INTEGER iG,jG,irec,bi,bj,i,j,k,dUnit,IL,pIL
      INTEGER dimList(3,3),nDims
      INTEGER x_size,y_size,length_of_rec
#if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
      INTEGER iG_IO,jG_IO,npe
      PARAMETER ( x_size = exch2_domain_nxt * sNx )
      PARAMETER ( y_size = exch2_domain_nyt * sNy )
#else
      PARAMETER ( x_size = Nx )
      PARAMETER ( y_size = Ny )
#endif
      Real*4 r4seg(sNx)
      Real*8 r8seg(sNx)
      Real*4 xy_buffer_r4(x_size,y_size)
      Real*8 xy_buffer_r8(x_size,y_size)
      Real*8 globalBuf(Nx,Ny)
#ifdef ALLOW_EXCH2
c     INTEGER tGy,tGx,tNy,tNx,tn
      INTEGER tGy,tGx,    tNx,tn
#endif /* ALLOW_EXCH2 */
      INTEGER tNy

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C Assume nothing
      fileIsOpen = .FALSE.
      IL  = ILNBLNK( fName )
      pIL = ILNBLNK( mdsioLocalDir )
      irecord = ABS(jrecord)
      writeMetaF = jrecord.GT.0

C Only do I/O if I am the master thread (and mpi process 0 IF useSingleCpuIO):
      iAmDoingIO = MASTER_CPU_IO(myThid)

C Only do I/O if I am the master thread
      IF ( iAmDoingIO ) THEN

C Record number must be >= 1
        IF (irecord .LT. 1) THEN
         WRITE(msgBuf,'(A,I9.8)')
     &     ' MDS_WRITE_FIELD: argument irecord = ',irecord
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT , myThid)
         WRITE(msgBuf,'(A)')
     &     ' MDS_WRITE_FIELD: invalid value for irecord'
         CALL PRINT_ERROR( msgBuf, myThid )
         STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
        ENDIF

C Assign special directory
        IF ( useCurrentDir .OR. pIL.EQ.0 ) THEN
         pfName = fName
        ELSE
         WRITE(pfName,'(2A)') mdsioLocalDir(1:pIL), fName(1:IL)
        ENDIF
        pIL=ILNBLNK( pfName )

C Assign a free unit number as the I/O channel for this routine
        CALL MDSFINDUNIT( dUnit, myThid )

C- endif iAmDoingIO
      ENDIF

C If option globalFile is desired but does not work or if
C globalFile is too slow, then try using single-CPU I/O.
      IF (useSingleCpuIO) THEN

C Master thread of process 0, only, opens a global file
       IF ( iAmDoingIO ) THEN
         WRITE(dataFName,'(2a)') fName(1:IL),'.data'
         length_of_rec=MDS_RECLEN(filePrec,x_size*y_size,myThid)
         IF (irecord .EQ. 1) THEN
          OPEN( dUnit, file=dataFName, status=_NEW_STATUS,
     &        access='direct', recl=length_of_rec )
         ELSE
          OPEN( dUnit, file=dataFName, status=_OLD_STATUS,
     &        access='direct', recl=length_of_rec )
         ENDIF
       ENDIF

C Gather array and WRITE it to file, one vertical level at a time
       DO k=1,nNz
C-      copy from arr(level=k) to 2-D "local":
        IF ( arrType.EQ.'RS' ) THEN
          CALL MDS_PASStoRS(sharedLocalBuf,arr,k,zSize,.FALSE.,myThid)
        ELSEIF ( arrType.EQ.'RL' ) THEN
          CALL MDS_PASStoRL(sharedLocalBuf,arr,k,zSize,.FALSE.,myThid)
        ELSE
          WRITE(msgBuf,'(A)')
     &         ' MDS_WRITE_FIELD: illegal value for arrType'
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
        ENDIF
        CALL GATHER_2D( globalBuf, sharedLocalBuf, myThid )

        IF ( iAmDoingIO ) THEN
          irec=k+nNz*(irecord-1)
          IF (filePrec .EQ. precFloat32) THEN
#if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
           DO J=1,y_size
            DO I=1,x_size
             xy_buffer_r4(I,J) = 0.0
            ENDDO
           ENDDO
           bj=1
           DO npe=1,nPx*nPy
            DO bi=1,nSx
             DO J=1,sNy
              DO I=1,sNx
#ifdef ALLOW_USE_MPI
               iG=mpi_myXGlobalLo(npe)-1+(bi-1)*sNx+i
               jG=mpi_myYGlobalLo(npe)-1+(bj-1)*sNy+j
#else
               iG= myXGlobalLo-1+(bi-1)*sNx+i
               jG= myYGlobalLo-1+(bj-1)*sNy+j
#endif /* ALLOW_USE_MPI */
               iG_IO=exch2_txGlobalo(W2_mpi_myTileList(npe,bi))+i-1
               jG_IO=exch2_tyGlobalo(W2_mpi_myTileList(npe,bi))+j-1
               xy_buffer_r4(iG_IO,jG_IO) = globalBuf(iG,jG)
              ENDDO
             ENDDO
            ENDDO
           ENDDO
#else /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
           DO J=1,Ny
            DO I=1,Nx
             xy_buffer_r4(I,J) = globalBuf(I,J)
            ENDDO
           ENDDO
#endif /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
#ifdef _BYTESWAPIO
           CALL MDS_BYTESWAPR4( x_size*y_size, xy_buffer_r4 )
#endif
           WRITE(dUnit,rec=irec) xy_buffer_r4
          ELSEIF (filePrec .EQ. precFloat64) THEN
#if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
           DO J=1,y_size
            DO I=1,x_size
             xy_buffer_r8(I,J) = 0.0
            ENDDO
           ENDDO
           bj=1
           DO npe=1,nPx*nPy
            DO bi=1,nSx
             DO J=1,sNy
              DO I=1,sNx
#ifdef ALLOW_USE_MPI
               iG=mpi_myXGlobalLo(npe)-1+(bi-1)*sNx+i
               jG=mpi_myYGlobalLo(npe)-1+(bj-1)*sNy+j
#else
               iG= myXGlobalLo-1+(bi-1)*sNx+i
               jG= myYGlobalLo-1+(bj-1)*sNy+j
#endif /* ALLOW_USE_MPI */
               iG_IO=exch2_txGlobalo(W2_mpi_myTileList(npe,bi))+i-1
               jG_IO=exch2_tyGlobalo(W2_mpi_myTileList(npe,bi))+j-1
               xy_buffer_r8(iG_IO,jG_IO) = globalBuf(iG,jG)
              ENDDO
             ENDDO
            ENDDO
           ENDDO
#else /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
           DO J=1,Ny
            DO I=1,Nx
             xy_buffer_r8(I,J) = globalBuf(I,J)
            ENDDO
           ENDDO
#endif /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
#ifdef _BYTESWAPIO
           CALL MDS_BYTESWAPR8( x_size*y_size, xy_buffer_r8 )
#endif
           WRITE(dUnit,rec=irec) xy_buffer_r8
          ELSE
           WRITE(msgBuf,'(A)')
     &       ' MDS_WRITE_FIELD: illegal value for filePrec'
           CALL PRINT_ERROR( msgBuf, myThid )
           STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
          ENDIF
        ENDIF
       ENDDO

C Close data-file
       IF ( iAmDoingIO ) THEN
         CLOSE( dUnit )
       ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C---  else .NOT.useSingleCpuIO
      ELSE

C Only do I/O if I am the master thread
       IF ( iAmDoingIO ) THEN

C If we are writing to a global file then we open it here
        IF (globalFile) THEN
         WRITE(dataFName,'(2a)') fName(1:IL),'.data'
         IF (irecord .EQ. 1) THEN
          length_of_rec=MDS_RECLEN( filePrec, sNx, myThid )
          OPEN( dUnit, file=dataFName, status=_NEW_STATUS,
     &            access='direct', recl=length_of_rec )
          fileIsOpen=.TRUE.
         ELSE
          length_of_rec=MDS_RECLEN( filePrec, sNx, myThid )
          OPEN( dUnit, file=dataFName, status=_OLD_STATUS,
     &            access='direct', recl=length_of_rec )
          fileIsOpen=.TRUE.
         ENDIF
        ENDIF

C Loop over all tiles
        DO bj=1,nSy
         DO bi=1,nSx
C If we are writing to a tiled MDS file then we open each one here
          IF (.NOT. globalFile) THEN
           iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
           jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
           WRITE(dataFName,'(2A,I3.3,A,I3.3,A)')
     &              pfName(1:pIL),'.',iG,'.',jG,'.data'
           IF (irecord .EQ. 1) THEN
            length_of_rec=MDS_RECLEN( filePrec, sNx, myThid )
            OPEN( dUnit, file=dataFName, status=_NEW_STATUS,
     &            access='direct', recl=length_of_rec )
            fileIsOpen=.TRUE.
           ELSE
            length_of_rec=MDS_RECLEN( filePrec, sNx, myThid )
            OPEN( dUnit, file=dataFName, status=_OLD_STATUS,
     &            access='direct', recl=length_of_rec )
            fileIsOpen=.TRUE.
           ENDIF
          ENDIF
          IF (fileIsOpen) THEN
           tNy = sNy
#ifdef ALLOW_EXCH2
           tn = W2_myTileList(bi)
           tGy = exch2_tyGlobalo(tn)
           tGx = exch2_txGlobalo(tn)
           tNy = exch2_tNy(tn)
           tNx = exch2_tNx(tn)
#endif /* ALLOW_EXCH2 */
           DO k=1,nNz
            DO j=1,tNy
             IF (globalFile) THEN
#ifdef ALLOW_EXCH2
              irec = 1 + (tGx-1)/tNx
     &                 + ( j-1 + tGy-1 )*exch2_domain_nxt
     &                 + ( k-1 + (irecord-1)*nNz
     &                   )*tNy*exch2_domain_nyt*exch2_domain_nxt
#else /* ALLOW_EXCH2 */
              iG = myXGlobalLo-1 + (bi-1)*sNx
              jG = myYGlobalLo-1 + (bj-1)*sNy
              irec= 1 + INT(iG/sNx) + nSx*nPx*(jG+j-1)
     &                + nSx*nPx*Ny*(k-1)
     &                + nSx*nPx*Ny*nNz*(irecord-1)
#endif /* ALLOW_EXCH2 */
             ELSE
              iG = 0
              jG = 0
              irec=j + sNy*(k-1) + sNy*nNz*(irecord-1)
             ENDIF
             IF (filePrec .EQ. precFloat32) THEN
              IF (arrType .EQ. 'RS') THEN
               CALL MDS_SEG4toRS( j,bi,bj,k,zSize, r4seg,.FALSE.,arr )
              ELSEIF (arrType .EQ. 'RL') THEN
               CALL MDS_SEG4toRL( j,bi,bj,k,zSize, r4seg,.FALSE.,arr )
              ELSE
               WRITE(msgBuf,'(A)')
     &           ' MDS_WRITE_FIELD: illegal value for arrType'
               CALL PRINT_ERROR( msgBuf, myThid )
               STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
              ENDIF
#ifdef _BYTESWAPIO
              CALL MDS_BYTESWAPR4( sNx, r4seg )
#endif
              WRITE(dUnit,rec=irec) r4seg
             ELSEIF (filePrec .EQ. precFloat64) THEN
              IF (arrType .EQ. 'RS') THEN
               CALL MDS_SEG8toRS( j,bi,bj,k,zSize, r8seg,.FALSE.,arr )
              ELSEIF (arrType .EQ. 'RL') THEN
               CALL MDS_SEG8toRL( j,bi,bj,k,zSize, r8seg,.FALSE.,arr )
              ELSE
               WRITE(msgBuf,'(A)')
     &           ' MDS_WRITE_FIELD: illegal value for arrType'
               CALL PRINT_ERROR( msgBuf, myThid )
               STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
              ENDIF
#ifdef _BYTESWAPIO
              CALL MDS_BYTESWAPR8( sNx, r8seg )
#endif
              WRITE(dUnit,rec=irec) r8seg
             ELSE
              WRITE(msgBuf,'(A)')
     &          ' MDS_WRITE_FIELD: illegal value for filePrec'
              CALL PRINT_ERROR( msgBuf, myThid )
              STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
             ENDIF
C End of j loop
            ENDDO
C End of k loop
           ENDDO
          ELSE
C fileIsOpen=F
           WRITE(msgBuf,'(A)')
     &       ' MDS_WRITE_FIELD: I should never get to this point'
           CALL PRINT_ERROR( msgBuf, myThid )
           STOP 'ABNORMAL END: S/R MDS_WRITE_FIELD'
          ENDIF
C If we were writing to a tiled MDS file then we close it here
          IF (fileIsOpen .AND. (.NOT. globalFile)) THEN
           CLOSE( dUnit )
           fileIsOpen = .FALSE.
          ENDIF
C Create meta-file for each tile if we are tiling
          IF ( .NOT.globalFile .AND. writeMetaF ) THEN
           iG=bi+(myXGlobalLo-1)/sNx
           jG=bj+(myYGlobalLo-1)/sNy
           WRITE(metaFname,'(2A,I3.3,A,I3.3,A)')
     &              pfName(1:pIL),'.',iG,'.',jG,'.meta'
#if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
           tn = W2_myTileList(bi)
           dimList(1,1)=x_size
           dimList(2,1)=exch2_txGlobalo(tn)
           dimList(3,1)=exch2_txGlobalo(tn)+sNx-1
           dimList(1,2)=y_size
           dimList(2,2)=exch2_tyGlobalo(tn)
           dimList(3,2)=exch2_tyGlobalo(tn)+sNy-1
#else /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
C- jmc: if MISSING_TILE_IO, keep meta files unchanged
C       to stay consistent with global file structure
           dimList(1,1)=Nx
           dimList(2,1)=myXGlobalLo+(bi-1)*sNx
           dimList(3,1)=myXGlobalLo+bi*sNx-1
           dimList(1,2)=Ny
           dimList(2,2)=myYGlobalLo+(bj-1)*sNy
           dimList(3,2)=myYGlobalLo+bj*sNy-1
#endif /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
           dimList(1,3)=nNz
           dimList(2,3)=1
           dimList(3,3)=nNz
           nDims=3
           IF ( nNz.EQ.1 ) nDims=2
           CALL MDS_WRITE_META(
     I              metaFName, dataFName, the_run_name, ' ',
     I              filePrec, nDims, dimList, 0,  ' ',
     I              0, UNSET_RL, irecord, myIter, myThid )
          ENDIF
C End of bi,bj loops
         ENDDO
        ENDDO

C If global file was opened then close it
        IF (fileIsOpen .AND. globalFile) THEN
         CLOSE( dUnit )
         fileIsOpen = .FALSE.
        ENDIF

C- endif iAmDoingIO
       ENDIF

C     if useSingleCpuIO / else / end
      ENDIF

C Create meta-file for the global-file (also if useSingleCpuIO)
      IF ( writeMetaF .AND. iAmDoingIO .AND.
     &    (globalFile .OR. useSingleCpuIO) ) THEN
         WRITE(metaFName,'(2A)') fName(1:IL),'.meta'
         dimList(1,1)=x_size
         dimList(2,1)=1
         dimList(3,1)=x_size
         dimList(1,2)=y_size
         dimList(2,2)=1
         dimList(3,2)=y_size
         dimList(1,3)=nNz
         dimList(2,3)=1
         dimList(3,3)=nNz
         ndims=3
         IF ( nNz.EQ.1 ) ndims=2
         CALL MDS_WRITE_META(
     I              metaFName, dataFName, the_run_name, ' ',
     I              filePrec, nDims, dimList, 0,  ' ',
     I              0, UNSET_RL, irecord, myIter, myThid )
c    I              metaFName, dataFName, the_run_name, titleLine,
c    I              filePrec, nDims, dimList, nFlds,  fldList,
c    I              nTimRec, timList, irecord, myIter, myThid )
      ENDIF

C To be safe, make other processes wait for I/O completion
      _BARRIER

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      RETURN
      END