C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.61 2013/02/06 21:25:26 jmc Exp $
C $Name:  $

#include "DIAG_OPTIONS.h"

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

C     !INTERFACE:
      SUBROUTINE DIAGNOSTICS_OUT(
     I     listId,
     I     myTime,
     I     myIter,
     I     myThid )

C     !DESCRIPTION:
C     Write output for diagnostics fields.

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"

      INTEGER NrMax
      PARAMETER( NrMax = numLevels )

C     !INPUT PARAMETERS:
C     listId  :: Diagnostics list number being written
C     myIter  :: current iteration number
C     myTime  :: current time of simulation (s)
C     myThid  :: my Thread Id number
      _RL     myTime
      INTEGER listId, myIter, myThid
CEOP

C     !FUNCTIONS:
      INTEGER ILNBLNK
      EXTERNAL ILNBLNK

C     !LOCAL VARIABLES:
C     i,j,k :: loop indices
C     bi,bj :: tile indices
C     lm    :: loop index (averageCycle)
C     md    :: field number in the list "listId".
C     ndId  :: diagnostics  Id number (in available diagnostics list)
C     ip    :: diagnostics  pointer to storage array
C     im    :: counter-mate pointer to storage array
C     mate  :: counter mate Id number (in available diagnostics list)
C     mDbl  :: processing mate Id number (in case processing requires 2 diags)
C     mVec  :: vector mate Id number
C     ppFld :: post-processed diag or not (=0): =1 stored in qtmp1 ; =2 in qtmp2
C   isComputed :: previous post-processed diag (still available in qtmp)
C     nLevOutp :: number of levels to write in output file
C
C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
C     qtmp1 :: temporary array; used to store a copy of diag. output field.
C     qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
C-  Note: local common block no longer needed.
c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
      _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
      _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)

      INTEGER i, j, k, lm
      INTEGER bi, bj
      INTEGER md, ndId, nn, ip, im
      INTEGER mate, mDbl, mVec
      INTEGER ppFld, isComputed
      CHARACTER*10 gcode
      _RL undefRL
      INTEGER nLevOutp, kLev

      INTEGER iLen
      INTEGER ioUnit
      CHARACTER*(MAX_LEN_FNAM) fn
      CHARACTER*(MAX_LEN_MBUF) suff
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER prec, nRec, nTimRec
      _RL     timeRec(2)
      _RL     tmpLoc
#ifdef ALLOW_MDSIO
      LOGICAL glf
#endif
#ifdef ALLOW_MNC
      CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
#endif /*  ALLOW_MNC  */

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

C---  set file properties
      ioUnit= standardMessageUnit
      undefRL = misValFlt(listId)

      WRITE(suff,'(I10.10)') myIter
      iLen = ILNBLNK(fnames(listId))
      WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
C-    for now, if integrate vertically, output field has just 1 level:
      nLevOutp = nlevels(listId)
      IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1

C--   Set time information:
      IF ( freq(listId).LT.0. ) THEN
C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
        nTimRec = 1
        timeRec(1) = myTime
      ELSE
C-    Time-average: store the 2 edges of the time-averaging interval.
C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
        nTimRec = 2

C-    end of time-averaging interval:
        timeRec(2) = myTime

C-    begining of time-averaging interval:
c       timeRec(1) = myTime - freq(listId)
C     a) find the time of the previous multiple of output freq:
        timeRec(1) = myTime-deltaTClock*0.5 _d 0
        timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
        i = INT( timeRec(1) )
        IF ( timeRec(1).LT.0. ) THEN
          tmpLoc = FLOAT(i)
          IF ( timeRec(1).NE.tmpLoc ) i = i - 1
        ENDIF
        timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
        timeRec(1) = MAX( timeRec(1), startTime )

C     b) round off to nearest multiple of time-step:
        timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
        i = NINT( timeRec(1) )
C     if just half way, NINT will return the next time-step: correct this
        tmpLoc = FLOAT(i) - 0.5 _d 0
        IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
        timeRec(1) = baseTime + deltaTClock*FLOAT(i)
c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
      ENDIF
C--   Convert time to iteration number (debug)
c     DO i=1,nTimRec
c       timeRec(i) = timeRec(i)/deltaTClock
c     ENDDO

C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
      DO lm=1,averageCycle(listId)

#ifdef ALLOW_MNC
       IF (useMNC .AND. diag_mnc) THEN
         CALL DIAGNOSTICS_MNC_SET(
     I                    nLevOutp, listId, lm,
     O                    diag_mnc_bn,
     I                    undefRL, myTime, myIter, myThid )
       ENDIF
#endif /*  ALLOW_MNC  */

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

       isComputed = 0
       DO md = 1,nfields(listId)
        ndId = jdiag(md,listId)
        gcode = gdiag(ndId)(1:10)
        mate = 0
        mVec = 0
        mDbl = 0
        ppFld = 0
        IF ( gcode(5:5).EQ.'C' ) THEN
C-      Check for Mate of a Counter Diagnostic
           mate = hdiag(ndId)
        ELSEIF ( gcode(5:5).EQ.'P' ) THEN
           ppFld = 1
           IF ( gdiag(hdiag(ndId))(5:5).EQ.'P' ) ppFld = 2
C-      Also load the mate (if stored) for Post-Processing
           nn = ndId
           DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
             nn = hdiag(nn)
           ENDDO
           IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
c          write(0,*) ppFld,' ndId=', ndId, nn, mDbl, isComputed
        ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
C-      Check for Mate of a Vector Diagnostic
           mVec = hdiag(ndId)
        ENDIF
        IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
C--     Start processing 1 Fld :

          ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
          im = mdiag(md,listId)
          IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
          IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
          IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)

          IF ( ppFld.EQ.2 .AND. isComputed.EQ.hdiag(ndId) ) THEN
C-        Post-Processed diag from an other Post-Processed diag -and-
C         both of them have just been calculated and are still stored in qtmp:
C         => skip computation and just write qtmp2
            IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
               WRITE(ioUnit,'(A,I6,3A,I6)')
     &         '  get Post-Proc. Diag # ', ndId, '  ', cdiag(ndId),
     &         ' from previous computation of Diag # ', isComputed
            ENDIF
            isComputed = 0
          ELSEIF ( ndiag(ip,1,1).EQ.0 ) THEN
C-        Empty diagnostics case :
            isComputed = 0

            _BEGIN_MASTER( myThid )
            WRITE(msgBuf,'(A,I10)')
     &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                          SQUEEZE_RIGHT, myThid)
            WRITE(msgBuf,'(A,I6,3A,I4,2A)')
     &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
     &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                          SQUEEZE_RIGHT, myThid)
            IF ( averageCycle(listId).GT.1 ) THEN
             WRITE(msgBuf,'(A,2(I3,A))')
     &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
     &                                            ndiag(ip,1,1), ' )'
            ELSE
             WRITE(msgBuf,'(A,2(I3,A))')
     &        '- WARNING -   has not been filled (ndiag=',
     &                                            ndiag(ip,1,1), ' )'
            ENDIF
            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                          SQUEEZE_RIGHT, myThid)
            WRITE(msgBuf,'(A)')
     &       'WARNING DIAGNOSTICS_OUT  => write ZEROS instead'
            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                          SQUEEZE_RIGHT, myThid)
            _END_MASTER( myThid )
            DO bj = myByLo(myThid), myByHi(myThid)
              DO bi = myBxLo(myThid), myBxHi(myThid)
                DO k = 1,nLevOutp
                  DO j = 1-OLy,sNy+OLy
                    DO i = 1-OLx,sNx+OLx
                      qtmp1(i,j,k,bi,bj) = 0. _d 0
                    ENDDO
                  ENDDO
                ENDDO
              ENDDO
            ENDDO

          ELSE
C-        diagnostics is not empty :
            isComputed = 0

            IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
              IF ( ppFld.GE.1 ) THEN
               WRITE(ioUnit,'(A,I6,7A,I8,2A)')
     &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
     &         '   Parms: ',gdiag(ndId)
               IF ( mDbl.EQ.0 ) THEN
                WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
     &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
               ELSE
                WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
     &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
     &          ' and diag: ',
     &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
               ENDIF
              ELSE
               WRITE(ioUnit,'(A,I6,3A,I8,2A)')
     &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
     &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
              ENDIF
              IF ( mate.GT.0 ) THEN
               WRITE(ioUnit,'(3A,I6,2A)')
     &         '       use Counter Mate for  ', cdiag(ndId),
     &         '     Diagnostic # ',mate, '  ', cdiag(mate)
              ELSEIF ( mVec.GT.0 ) THEN
                IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
                 WRITE(ioUnit,'(3A,I6,3A)')
     &             '           Vector  Mate for  ', cdiag(ndId),
     &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
     &             ' exists '
                ELSE
                 WRITE(ioUnit,'(3A,I6,3A)')
     &             '           Vector  Mate for  ', cdiag(ndId),
     &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
     &             ' not enabled'
                ENDIF
              ENDIF
            ENDIF

            IF ( fflags(listId)(2:2).EQ.' ' ) THEN
C-       get only selected levels:
              DO bj = myByLo(myThid), myByHi(myThid)
               DO bi = myBxLo(myThid), myBxHi(myThid)
                DO k = 1,nlevels(listId)
                  kLev = NINT(levs(k,listId))
                  CALL DIAGNOSTICS_GET_DIAG(
     I                         kLev, undefRL,
     O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
     I                         ndId, mate, ip, im, bi, bj, myThid )
                ENDDO
               ENDDO
              ENDDO
              IF ( mDbl.GT.0 ) THEN
               DO bj = myByLo(myThid), myByHi(myThid)
                DO bi = myBxLo(myThid), myBxHi(myThid)
                 DO k = 1,nlevels(listId)
                  kLev = NINT(levs(k,listId))
                  CALL DIAGNOSTICS_GET_DIAG(
     I                         kLev, undefRL,
     O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
     I                         mDbl, 0, im, 0, bi, bj, myThid )
                 ENDDO
                ENDDO
               ENDDO
              ENDIF
            ELSE
C-       get all the levels (for vertical post-processing)
              DO bj = myByLo(myThid), myByHi(myThid)
               DO bi = myBxLo(myThid), myBxHi(myThid)
                  CALL DIAGNOSTICS_GET_DIAG(
     I                         0, undefRL,
     O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
     I                         ndId, mate, ip, im, bi, bj, myThid )
               ENDDO
              ENDDO
              IF ( mDbl.GT.0 ) THEN
               DO bj = myByLo(myThid), myByHi(myThid)
                DO bi = myBxLo(myThid), myBxHi(myThid)
                  CALL DIAGNOSTICS_GET_DIAG(
     I                         0, undefRL,
     O                         qtmp2(1-OLx,1-OLy,1,bi,bj),
     I                         mDbl, 0, im, 0, bi, bj, myThid )
                ENDDO
               ENDDO
              ENDIF
            ENDIF

C-----------------------------------------------------------------------
C--     Apply specific post-processing (e.g., interpolate) before output
C-----------------------------------------------------------------------
            IF ( fflags(listId)(2:2).EQ.'P' ) THEN
C-          Do vertical interpolation:
             IF ( fluidIsAir ) THEN
C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
              CALL DIAGNOSTICS_INTERP_VERT(
     I                         listId, md, ndId, ip, im, lm,
     U                         qtmp1, qtmp2,
     I                         undefRL, myTime, myIter, myThid )
             ELSE
               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
     &           'INTERP_VERT not allowed in this config'
               CALL PRINT_ERROR( msgBuf , myThid )
               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
             ENDIF
            ENDIF
            IF ( fflags(listId)(2:2).EQ.'I' ) THEN
C-          Integrate vertically: for now, output field has just 1 level:
              CALL DIAGNOSTICS_SUM_LEVELS(
     I                         listId, md, ndId, ip, im, lm,
     U                         qtmp1,
     I                         undefRL, myTime, myIter, myThid )
            ENDIF
            IF ( ppFld.GE.1 ) THEN
C-          Do Post-Processing:
             IF ( flds(md,listId).EQ.'PhiVEL  '
     &       .OR. flds(md,listId).EQ.'PsiVEL  '
     &          ) THEN
              CALL DIAGNOSTICS_CALC_PHIVEL(
     I                         listId, md, ndId, ip, im, lm,
     I                         NrMax,
     U                         qtmp1, qtmp2,
     I                         myTime, myIter, myThid )
              isComputed = ndId
             ELSE
               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
     &           'unknown Processing method for diag="',cdiag(ndId),'"'
               CALL PRINT_ERROR( msgBuf , myThid )
               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
             ENDIF
            ENDIF

C--     End of empty diag / not-empty diag block
          ENDIF

C--     Ready to write field "md", element "lm" in averageCycle(listId)

C-        write to binary file, using MDSIO pkg:
          IF ( diag_mdsio ) THEN
c          nRec = lm + (md-1)*averageCycle(listId)
           nRec = md + (lm-1)*nfields(listId)
C         default precision for output files
           prec = writeBinaryPrec
C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
           IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
           IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
           IF ( ppFld.LE.1 ) THEN
            CALL WRITE_REC_LEV_RL(
     I                            fn, prec,
     I                            NrMax, 1, nLevOutp,
     I                            qtmp1, -nRec, myIter, myThid )
           ELSE
            CALL WRITE_REC_LEV_RL(
     I                            fn, prec,
     I                            NrMax, 1, nLevOutp,
     I                            qtmp2, -nRec, myIter, myThid )
           ENDIF
          ENDIF

#ifdef ALLOW_MNC
          IF (useMNC .AND. diag_mnc) THEN
           IF ( ppFld.LE.1 ) THEN
            CALL DIAGNOSTICS_MNC_OUT(
     I                       NrMax, nLevOutp, listId, ndId, mate,
     I                       diag_mnc_bn, qtmp1,
     I                       undefRL, myTime, myIter, myThid )
           ELSE
            CALL DIAGNOSTICS_MNC_OUT(
     I                       NrMax, nLevOutp, listId, ndId, mate,
     I                       diag_mnc_bn, qtmp2,
     I                       undefRL, myTime, myIter, myThid )
           ENDIF
          ENDIF
#endif /*  ALLOW_MNC  */

C--     end of Processing Fld # md
        ENDIF
       ENDDO

C--   end loop on lm counter (= averagePeriod)
      ENDDO

#ifdef ALLOW_MDSIO
      IF (diag_mdsio) THEN
C-    Note: temporary: since it is a pain to add more arguments to
C     all MDSIO S/R, uses instead this specific S/R to write only
C     meta files but with more informations in it.
            glf = globalFiles
            nRec = averageCycle(listId)*nfields(listId)
            CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
     &              0, 0, nLevOutp, ' ',
     &              nfields(listId), flds(1,listId),
     &              nTimRec, timeRec, undefRL,
     &              nRec, myIter, myThid)
      ENDIF
#endif /*  ALLOW_MDSIO  */

      RETURN
      END

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