C $Header: /u/gcmpack/MITgcm/pkg/mnc/mnc_cw_cvars.F,v 1.16 2011/05/23 01:08:22 jmc Exp $
C $Name:  $

#include "MNC_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 1
C     !ROUTINE: MNC_CW_WRITE_CVAR

C     !INTERFACE:
      SUBROUTINE MNC_CW_WRITE_CVAR(
     I     fname,
     I     cvname,
     I     fid,
     I     did,
     I     bi, bj,
     I     myThid )

C     !DESCRIPTION:
C     Write a CF-convention coordinate variable (a vector).

C     !USES:
      implicit none
#include "MNC_COMMON.h"
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_SIZE.h"
#include "W2_EXCH2_TOPOLOGY.h"
#endif
#include "netcdf.inc"

C     Functions
      integer IFNBLNK, ILNBLNK

C     !INPUT PARAMETERS:
      character*(*) fname
      character*(*) cvname
      integer fid, did, bi,bj
      integer myThid
CEOP

C     !LOCAL VARIABLES:
      integer i, vid, nnf, nnl, doit, err
      integer nids, cv_did(1), xtmin,ytmin
      character*(MAX_LEN_MBUF) msgbuf
      integer cv_start(1), cv_count(1)
      _RS rtmp(sNx + 2*OLx + sNy + 2*OLy + Nr)
C     variables for text attributes
      integer MAX_LEN_NAME, ia
      PARAMETER ( MAX_LEN_NAME = 128 )
      character*(MAX_LEN_NAME) units, long_name, positive

      DO i=1,MAX_LEN_NAME
         units(i:i)     = ' '
         long_name(i:i) = ' '
         positive(i:i)  = ' '
      ENDDO

      nnf = IFNBLNK(cvname)
      nnl = ILNBLNK(cvname)

      xtmin = 0
      ytmin = 0
#ifdef ALLOW_EXCH2
      xtmin = exch2_tbasex(W2_myTileList(bi,bj))
      ytmin = exch2_tbasey(W2_myTileList(bi,bj))
#else
      IF ( .NOT. useCubedSphereExchange ) THEN
C     make sure for a non-cubed-sphere curvi-linear grid,
C     that the X/Y coordinate variables are monotonous
C     bi+(myXGlobalLo-1)/sNx is the i-tile number
C     bj+(myYGlobalLo-1)/sNy is the j-tile number
       xtmin = sNx * ( bi+(myXGlobalLo-1)/sNx - 1 )
       ytmin = sNy * ( bj+(myYGlobalLo-1)/sNy - 1 )
      ENDIF
#endif
      doit = 1
      nids = 1
      cv_did(1)= did

C     Check all the coordinate variables that we know about
      IF (cvname(nnf:nnl) .EQ. 'X') THEN

        cv_start(1) = 1
        cv_count(1) = sNx
#ifdef ALLOW_EXCH2
        DO i = cv_start(1),cv_count(1)
         rtmp(i) = xtmin + i
        ENDDO
#else
        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = xtmin + i
         ENDDO
        ELSE
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = xC(i,1,bi,bj)
         ENDDO
        ENDIF
#endif
        IF ( usingCartesianGrid ) THEN
         long_name = 'X-coordinate of cell center'
         units     = 'meters'
        ELSEIF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         long_name = 'i-index of cell center'
         units     = 'none'
        ELSEIF ( usingSphericalPolarGrid ) THEN
         long_name = 'longitude of cell center'
         units     = 'degrees_east'
        ELSEIF ( usingCylindricalGrid ) THEN
         long_name = 'polar angle coordinate of cell center'
         units     = 'degrees'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown horizontal grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Xp1') THEN

        cv_start(1) = 1
        cv_count(1) = sNx + 1
#ifdef ALLOW_EXCH2
        DO i = cv_start(1),cv_count(1)
         rtmp(i) = xtmin + i
        ENDDO
#else
        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = xtmin + i
         ENDDO
        ELSE
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = xG(i,1,bi,bj)
         ENDDO
        ENDIF
#endif
        IF ( usingCartesianGrid ) THEN
         long_name = 'X-Coordinate of cell corner'
         units     = 'meters'
        ELSEIF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         long_name = 'i-index of cell corner'
         units     = 'none'
        ELSEIF ( usingSphericalPolarGrid ) THEN
         long_name = 'longitude of cell corner'
         units     = 'degrees_east'
        ELSEIF ( usingCylindricalGrid ) THEN
         long_name = 'polar angle  of cell corner'
         units     = 'degrees'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown horizontal grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Xwh') THEN

        cv_start(1) = 1
        cv_count(1) = sNx + 2*OLx
#ifdef ALLOW_EXCH2
        DO i = cv_start(1),cv_count(1)
         rtmp(i) = xtmin + i
        ENDDO
#else
        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = xtmin - OLx + i
         ENDDO
        ELSE
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = xC(i,1,bi,bj)
CML????          rtmp(i) = xC(i-Olx,1,bi,bj)
         ENDDO
        ENDIF
#endif
        IF ( usingCartesianGrid ) THEN
         long_name = 'X-Coordinate of cell center including overlaps'
         units     = 'meters'
        ELSEIF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         long_name = 'i-index of cell center including overlaps'
         units     = 'none'
        ELSEIF ( usingSphericalPolarGrid ) THEN
         long_name = 'longitude of cell center including overlaps'
         units     = 'degrees_east'
        ELSEIF ( usingCylindricalGrid ) THEN
         long_name =
     &        'polar angle coordinate of cell center including overlaps'
         units     = 'degrees'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown horizontal grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Y') THEN

        cv_start(1) = 1
        cv_count(1) = sNy
#ifdef ALLOW_EXCH2
        DO i = cv_start(1),cv_count(1)
         rtmp(i) = ytmin + i
        ENDDO
#else
        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = ytmin + i
         ENDDO
        ELSE
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = yC(1,i,bi,bj)
         ENDDO
        ENDIF
#endif
        IF ( usingCartesianGrid ) THEN
         long_name = 'Y-Coordinate of cell center'
         units     = 'meters'
        ELSEIF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         long_name = 'j-index of cell center'
         units     = 'none'
        ELSEIF ( usingSphericalPolarGrid ) THEN
         long_name = 'latitude of cell center'
         units     = 'degrees_north'
        ELSEIF ( usingCylindricalGrid ) THEN
         long_name = 'radial coordinate of cell center'
         units     = 'meters'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown horizontal grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Yp1') THEN

        cv_start(1) = 1
        cv_count(1) = sNy + 1
#ifdef ALLOW_EXCH2
        DO i = cv_start(1),cv_count(1)
         rtmp(i) = ytmin + i
        ENDDO
#else
        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = ytmin + i
         ENDDO
        ELSE
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = yG(1,i,bi,bj)
         ENDDO
        ENDIF
#endif
        IF ( usingCartesianGrid ) THEN
         long_name = 'Y-Coordinate of cell corner'
         units     = 'meters'
        ELSEIF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         long_name = 'j-index of cell corner'
         units     = 'none'
        ELSEIF ( usingSphericalPolarGrid ) THEN
         long_name = 'latitude of cell corner'
         units     = 'degrees_north'
        ELSEIF ( usingCylindricalGrid ) THEN
         long_name = 'radial coordinate of cell corner'
         units     = 'meters'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown horizontal grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Ywh') THEN

        cv_start(1) = 1
        cv_count(1) = sNy + 2*OLy
#ifdef ALLOW_EXCH2
        DO i = cv_start(1),cv_count(1)
         rtmp(i) = ytmin + i
        ENDDO
#else
        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = ytmin - OLy + i
         ENDDO
        ELSE
         DO i = cv_start(1),cv_count(1)
          rtmp(i) = yC(1,i-OLy,bi,bj)
         ENDDO
        ENDIF
#endif
        IF ( usingCartesianGrid ) THEN
         long_name = 'Y-Coordinate of cell center including overlaps'
         units     = 'meters'
        ELSEIF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
         long_name = 'j-index of cell center including overlaps'
         units     = 'none'
        ELSEIF ( usingSphericalPolarGrid ) THEN
         long_name = 'latitude of cell center including overlaps'
         units     = 'degrees_north'
        ELSEIF ( usingCylindricalGrid ) THEN
         long_name =
     &        'radial coordinate of cell center including overlaps'
         units     = 'meters'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown horizontal grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Z') THEN

        cv_start(1) = 1
        cv_count(1) = Nr
        DO i = cv_start(1),cv_count(1)
          rtmp(i) = rC(i)
        ENDDO
C
        long_name = 'vertical coordinate of cell center'
        IF ( usingZCoords ) THEN
         units     = 'meters'
         positive  = 'up'
        ELSEIF ( usingPCoords ) THEN
         units     = 'pascal'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown vertical grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Zp1') THEN

        cv_start(1) = 1
        cv_count(1) = Nr + 1
        DO i = cv_start(1),cv_count(1)
          rtmp(i) = rF(i)
        ENDDO
C
        long_name = 'vertical coordinate of cell interface'
        IF ( usingZCoords ) THEN
         units     = 'meters'
         positive  = 'up'
        ELSEIF ( usingPCoords ) THEN
         units     = 'pascal'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown vertical grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Zu') THEN

        cv_start(1) = 1
        cv_count(1) = Nr
        DO i = cv_start(1),cv_count(1)
          rtmp(i) = rF(i + 1)
        ENDDO
C
        IF ( usingZCoords ) THEN
         long_name = 'vertical coordinate of lower cell interface'
         units     = 'meters'
         positive  = 'up'
        ELSEIF ( usingPCoords ) THEN
         long_name = 'vertical coordinate of upper cell interface'
         units     = 'pascal'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown vertical grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Zl') THEN

        cv_start(1) = 1
        cv_count(1) = Nr
        DO i = cv_start(1),cv_count(1)
          rtmp(i) = rF(i)
        ENDDO
C
        IF ( usingZCoords ) THEN
         long_name = 'vertical coordinate of upper cell interface'
         units     = 'meters'
         positive  = 'up'
        ELSEIF ( usingPCoords ) THEN
         long_name = 'vertical coordinate of lower cell interface'
         units     = 'pascal'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown vertical grid!'
        ENDIF

      ELSEIF (cvname(nnf:nnl) .EQ. 'Zm1') THEN

        cv_start(1) = 1
        cv_count(1) = Nr - 1
        DO i = cv_start(1),cv_count(1)
          rtmp(i) = rF(i + 1)
        ENDDO
C
        IF ( usingZCoords ) THEN
         long_name = 'vertical coordinate of lower cell interface'
         units     = 'meters'
         positive  = 'up'
        ELSEIF ( usingPCoords ) THEN
         long_name = 'vertical coordinate of upper cell interface'
         units     = 'pascal'
        ELSE
C       unknown grid type
         print *, 'S/R MNC_CW_CVARS: Ooops, unknown vertical grid!'
        ENDIF

      ELSE

        doit = 0

      ENDIF

      IF ( doit .EQ. 1 ) THEN

        CALL MNC_FILE_REDEF(fname, myThid)
#ifdef REAL4_IS_SLOW
        err = NF_DEF_VAR(fid, cvname, NF_DOUBLE,
     &       nids, cv_did, vid)
#else
        err = NF_DEF_VAR(fid, cvname, NF_FLOAT,
     &       nids, cv_did, vid)
#endif /* REAL4_IS_SLOW */
        i = ILNBLNK( fname )
        write(msgbuf,'(5a)') 'defining coordinate variable ''',
     &       cvname(nnf:nnl), ''' in file ''', fname(1:i), ''''
        CALL MNC_HANDLE_ERR(err, msgbuf, myThid)
C     add attributes if set
        ia = ILNBLNK(long_name)
        IF ( ia .GT. 0 ) THEN
         err = NF_PUT_ATT_TEXT(fid, vid, 'long_name', ia, long_name)
         write(msgbuf,'(5a)')
     &     'adding attribute ''long_name'' to coordinate variable ''',
     &     cvname(nnf:nnl), ''' in file ''', fname(1:i), ''''
         CALL MNC_HANDLE_ERR(err, msgbuf, myThid)
        ENDIF
        ia = ILNBLNK(units)
        IF ( ia .GT. 0 ) THEN
         err = NF_PUT_ATT_TEXT(fid, vid, 'units', ia, units)
         write(msgbuf,'(5a)')
     &     'adding attribute ''units'' to coordinate variable ''',
     &     cvname(nnf:nnl), ''' in file ''', fname(1:i), ''''
         CALL MNC_HANDLE_ERR(err, msgbuf, myThid)
        ENDIF
        ia = ILNBLNK(positive)
        IF ( ia .GT. 0 ) THEN
         err = NF_PUT_ATT_TEXT(fid, vid, 'positive', ia, positive)
         write(msgbuf,'(5a)')
     &     'adding attribute ''positive'' to coordinate variable ''',
     &     cvname(nnf:nnl), ''' in file ''', fname(1:i), ''''
         CALL MNC_HANDLE_ERR(err, msgbuf, myThid)
        ENDIF
C
        CALL MNC_FILE_ENDDEF(fname, myThid)
#ifdef REAL4_IS_SLOW
        err = NF_PUT_VARA_DOUBLE(fid, vid,
     &       cv_start, cv_count, rtmp)
#else
        err = NF_PUT_VARA_REAL(fid, vid,
     &       cv_start, cv_count, rtmp)
#endif /* REAL4_IS_SLOW */
        write(msgbuf,'(5a)') 'writing coordinate variable ''',
     &       cvname(nnf:nnl), ''' in file ''', fname(1:i), ''''
        CALL MNC_HANDLE_ERR(err, msgbuf, myThid)

      ENDIF

      RETURN
      END

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