pkg/diagnostics/diagnostics_utils.F

C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_utils.F,v 1.28 2009/01/25 17:00:20 jmc Exp $
C $Name:  $

#include "DIAG_OPTIONS.h"

C--   File diagnostics_utils.F: General purpose support routines
C--    Contents:
C--    o GETDIAG
C--    o DIAGNOSTICS_COUNT
C--    o DIAGNOSTICS_GET_POINTERS
C--    o DIAGNOSTICS_SETKLEV
C--    o DIAGS_GET_PARMS_I (Function)
C--    o DIAGS_MK_UNITS (Function)
C--    o DIAGS_MK_TITLE (Function)

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

C     !INTERFACE:
      SUBROUTINE GETDIAG(
     I                    levreal, undef,
     O                    qtmp,
     I                    ndId, mate, ip, im, bi, bj, myThid )

C     !DESCRIPTION:
C     Retrieve averaged model diagnostic

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

C     !INPUT PARAMETERS:
C     levreal :: Diagnostic LEVEL
C     undef   :: UNDEFINED VALUE
C     ndId    :: DIAGNOSTIC NUMBER FROM MENU
C     mate    :: counter DIAGNOSTIC NUMBER if any ; 0 otherwise
C     ip      :: pointer to storage array location for diag.
C     im      :: pointer to storage array location for mate
C     bi      :: X-direction tile number
C     bj      :: Y-direction tile number
C     myThid  :: my thread Id number
      _RL levreal
      _RL undef
      INTEGER ndId, mate, ip, im
      INTEGER bi,bj, myThid

C     !OUTPUT PARAMETERS:
C     qtmp    ..... AVERAGED DIAGNOSTIC QUANTITY
      _RL qtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

C     !LOCAL VARIABLES:
      _RL factor
      INTEGER i, j, ipnt,ipCt
      INTEGER lev, levCt, klev

      IF (ndId.GE.1) THEN
       lev = NINT(levreal)
       klev = kdiag(ndId)
       IF (lev.LE.klev) THEN

        IF ( mate.EQ.0 ) THEN
C-      No counter diagnostics => average = Sum / ndiag :

          ipnt = ip + lev - 1
          factor = FLOAT(ndiag(ip,bi,bj))
          IF (ndiag(ip,bi,bj).NE.0) factor = 1. _d 0 / factor

#ifdef ALLOW_FIZHI
          DO j = 1,sNy+1
            DO i = 1,sNx+1
              IF ( qdiag(i,j,ipnt,bi,bj) .LE. undef ) THEN
                qtmp(i,j) = qdiag(i,j,ipnt,bi,bj)*factor
              ELSE
                qtmp(i,j) = undef
              ENDIF
            ENDDO
          ENDDO
#else /* ALLOW_FIZHI */
          DO j = 1,sNy+1
            DO i = 1,sNx+1
              qtmp(i,j) = qdiag(i,j,ipnt,bi,bj)*factor
            ENDDO
          ENDDO
#endif /* ALLOW_FIZHI */

        ELSE
C-      With counter diagnostics => average = Sum / counter:

          ipnt = ip + lev - 1
          levCt= MIN(lev,kdiag(mate))
          ipCt = im + levCt - 1
          DO j = 1,sNy+1
            DO i = 1,sNx+1
              IF ( qdiag(i,j,ipCt,bi,bj) .NE. 0. ) THEN
                qtmp(i,j) = qdiag(i,j,ipnt,bi,bj)
     &                    / qdiag(i,j,ipCt,bi,bj)
              ELSE
                qtmp(i,j) = undef
              ENDIF
            ENDDO
          ENDDO

        ENDIF
       ENDIF
      ENDIF

      RETURN
      END

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

CBOP 0
C     !ROUTINE: DIAGNOSTICS_COUNT
C     !INTERFACE:
      SUBROUTINE DIAGNOSTICS_COUNT (chardiag,
     I                              biArg, bjArg, myThid)

C     !DESCRIPTION:
C***********************************************************************
C   routine to increment the diagnostic counter only
C***********************************************************************
C     !USES:
      IMPLICIT NONE

C     == Global variables ===
#include "EEPARAMS.h"
#include "SIZE.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"

C     !INPUT PARAMETERS:
C***********************************************************************
C  Arguments Description
C  ----------------------
C     chardiag :: Character expression for diag to increment the counter
C     biArg    :: X-direction tile number, or 0 if called outside bi,bj loops
C     bjArg    :: Y-direction tile number, or 0 if called outside bi,bj loops
C     myThid   :: my thread Id number
C***********************************************************************
      CHARACTER*8 chardiag
      INTEGER biArg, bjArg
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
C ===============
      INTEGER m, n
      INTEGER bi, bj
      INTEGER ipt, ndId
c     CHARACTER*(MAX_LEN_MBUF) msgBuf

      IF ( biArg.EQ.0 .AND. bjArg.EQ.0 ) THEN
        bi = myBxLo(myThid)
        bj = myByLo(myThid)
      ELSE
        bi = MIN(biArg,nSx)
        bj = MIN(bjArg,nSy)
      ENDIF

C--   Run through list of active diagnostics to find which counter
C     to increment (needs to be a valid & active diagnostic-counter)
      DO n=1,nlists
       DO m=1,nActive(n)
        IF ( chardiag.EQ.flds(m,n) .AND. idiag(m,n).GT.0 ) THEN
         ipt = idiag(m,n)
         IF (ndiag(ipt,bi,bj).GE.0) THEN
          ndId = jdiag(m,n)
          ipt = ipt + pdiag(n,bi,bj)*kdiag(ndId)
C-    Increment the counter for the diagnostic
          IF ( biArg.EQ.0 .AND. bjArg.EQ.0 ) THEN
           DO bj=myByLo(myThid), myByHi(myThid)
            DO bi=myBxLo(myThid), myBxHi(myThid)
             ndiag(ipt,bi,bj) = ndiag(ipt,bi,bj) + 1
            ENDDO
           ENDDO
          ELSE
             ndiag(ipt,bi,bj) = ndiag(ipt,bi,bj) + 1
          ENDIF
C-    Increment is done
         ENDIF
        ENDIF
       ENDDO
      ENDDO

      RETURN
      END

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

CBOP 0
C     !ROUTINE: DIAGNOSTICS_GET_POINTERS
C     !INTERFACE:
      SUBROUTINE DIAGNOSTICS_GET_POINTERS(
     I                       diagName, listId,
     O                       ndId, ip,
     I                       myThid )

C     !DESCRIPTION:
C     *================================================================*
C     | o Returns the diagnostic Id number and diagnostic
C     |   pointer to storage array for a specified diagnostic.
C     *================================================================*
C     | Note: A diagnostics field can be stored multiple times
C     |       (for different output frequency,phase, ...).
C     | operates in 2 ways:
C     | o listId =0 => find 1 diagnostics Id & pointer which name matches.
C     | o listId >0 => find the unique diagnostic Id & pointer with
C     |      the right name and same output time as "listId" output-list
C     | o return ip=0 if did not find the right diagnostic;
C     |   (ndId <>0 if diagnostic exist but output time does not match)
C     *================================================================*

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

C     !INPUT PARAMETERS:
C     diagName :: diagnostic identificator name (8 characters long)
C     listId   :: list number that specify the output frequency
C     myThid   :: my Thread Id number
C     !OUTPUT PARAMETERS:
C     ndId     :: diagnostics  Id number (in available diagnostics list)
C     ip       :: diagnostics  pointer to storage array


      CHARACTER*8 diagName
      INTEGER listId
      INTEGER ndId, ip
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
      INTEGER n,m

      ip   = 0
      ndId = 0

      IF ( listId.LE.0 ) THEN
C--   select the 1rst one which name matches:

C-    search for this diag. in the active 2D/3D diagnostics list
        DO n=1,nlists
         DO m=1,nActive(n)
           IF ( ip.EQ.0 .AND. diagName.EQ.flds(m,n)
     &                  .AND. idiag(m,n).NE.0 ) THEN
            ip   = ABS(idiag(m,n))
            ndId = jdiag(m,n)
           ENDIF
         ENDDO
        ENDDO

      ELSEIF ( listId.LE.nlists ) THEN
C--   select the unique diagnostic with output-time identical to listId

C-    search for this diag. in the active 2D/3D diagnostics list
        DO n=1,nlists
         IF ( ip.EQ.0
     &        .AND. freq(n) .EQ. freq(listId)
     &        .AND. phase(n).EQ.phase(listId)
     &        .AND. averageFreq(n) .EQ.averageFreq(listId)
     &        .AND. averagePhase(n).EQ.averagePhase(listId)
     &        .AND. averageCycle(n).EQ.averageCycle(listId)
     &      ) THEN
          DO m=1,nActive(n)
           IF ( ip.EQ.0 .AND. diagName.EQ.flds(m,n)
     &                  .AND. idiag(m,n).NE.0 ) THEN
            ip   = ABS(idiag(m,n))
            ndId = jdiag(m,n)
           ENDIF
          ENDDO
         ELSEIF ( ip.EQ.0 ) THEN
          DO m=1,nActive(n)
           IF ( ip.EQ.0 .AND. diagName.EQ.flds(m,n)
     &                  .AND. idiag(m,n).NE.0 ) THEN
            ndId = jdiag(m,n)
           ENDIF
          ENDDO
         ENDIF
        ENDDO

      ELSE
        STOP 'DIAGNOSTICS_GET_POINTERS: invalid listId number'
      ENDIF

      RETURN
      END

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

CBOP 0
C     !ROUTINE: DIAGNOSTICS_SETKLEV

C     !INTERFACE:
      SUBROUTINE DIAGNOSTICS_SETKLEV(
     I                                diagName, nLevDiag, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | S/R DIAGNOSTICS_SETKLEV
C     | o Define explicitly the number of level (stored in kdiag)
C     |   of a diagnostic field. For most diagnostics, the number
C     |   of levels is derived (in S/R SET_LEVELS) from gdiag(10)
C     |   but occasionally one may want to set it explicitly.
C     *==========================================================*

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

C     !INPUT PARAMETERS:
C     diagName  :: diagnostic identificator name (8 characters long)
C     nLevDiag  :: number of level to set for this diagnostics field
C     myThid    :: my Thread Id number
      CHARACTER*8  diagName
      INTEGER nLevDiag
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER n, ndId

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

C--   Check if this S/R is called from the right place ;
C     needs to be after DIAGNOSTICS_INIT_EARLY and before DIAGNOSTICS_INIT_FIXED
      IF ( .NOT.settingDiags ) THEN
        WRITE(msgBuf,'(4A,I5)') 'DIAGNOSTICS_SETKLEV: ',
     &     'diagName="', diagName, '" , nLevDiag=', nLevDiag
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SETKLEV: ',
     &     '<== called from the WRONG place, i.e.'
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SETKLEV: ',
     &     'outside diagnostics setting section = from'
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SETKLEV: ',
     &     '   Diag_INIT_EARLY down to Diag_INIT_FIXED'
        CALL PRINT_ERROR( msgBuf, myThid )
        STOP 'ABNORMAL END: S/R DIAGNOSTICS_SETKLEV'
      ENDIF

C--   Find this diagnostics in the list of available diag.
      ndId = 0
      DO n = 1,ndiagt
        IF ( diagName.EQ.cdiag(n) ) THEN
          ndId = n
        ENDIF
      ENDDO
      IF ( ndId.EQ.0 ) THEN
        WRITE(msgBuf,'(4A)') 'DIAGNOSTICS_SETKLEV: ',
     &     'diagName="', diagName, '" not known.'
        CALL PRINT_ERROR( msgBuf, myThid )
        STOP 'ABNORMAL END: S/R DIAGNOSTICS_SETKLEV'
      ENDIF

C-    Optional level number diagnostics (X): set number of levels
      IF ( kdiag(ndId).EQ.0
     &   .AND. gdiag(ndId)(10:10).EQ.'X' ) THEN
        kdiag(ndId) = nLevDiag
      ELSEIF ( kdiag(ndId).EQ.nLevDiag
     &   .AND. gdiag(ndId)(10:10).EQ.'X' ) THEN
C-    level number already set to same value: send warning
        WRITE(msgBuf,'(4A,I5)') '** WARNING ** DIAGNOSTICS_SETKLEV: ',
     &     'diagName="', diagName, '" , nLevDiag=', nLevDiag
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT , myThid )
        WRITE(msgBuf,'(2A,I5,A)')'** WARNING ** DIAGNOSTICS_SETKLEV:',
     &     ' level Nb (=', kdiag(ndId), ') already set.'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT , myThid )
      ELSEIF ( gdiag(ndId)(10:10).EQ.'X' ) THEN
C-    level number already set to a different value: do not reset but stop
        WRITE(msgBuf,'(4A,I5)') 'DIAGNOSTICS_SETKLEV: ',
     &     'diagName="', diagName, '" , nLevDiag=', nLevDiag
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(2A,I5,3A)') 'DIAGNOSTICS_SETKLEV: ',
     &     'level Nb already set to', kdiag(ndId), ' => STOP'
        CALL PRINT_ERROR( msgBuf, myThid )
      ELSE
C-    for now, do nothing but just send a warning
        WRITE(msgBuf,'(4A,I5)') '** WARNING ** DIAGNOSTICS_SETKLEV: ',
     &     'diagName="', diagName, '" , nLevDiag=', nLevDiag
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT , myThid )
        WRITE(msgBuf,'(2A,I5,3A)') '** WARNING ** will set level Nb',
     &     ' from diagCode(ndId=', ndId, ')="', gdiag(ndId)(1:10), '"'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT , myThid )
        WRITE(msgBuf,'(4A)') '** WARNING ** DIAGNOSTICS_SETKLEV',
     &     '("', diagName, '") <== Ignore this call.'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT , myThid )
      ENDIF

      RETURN
      END

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

CBOP 0
C     !ROUTINE: DIAGS_GET_PARMS_I

C     !INTERFACE:
      INTEGER FUNCTION DIAGS_GET_PARMS_I(
     I                            parName, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | FUNCTION DIAGS_GET_PARMS_I
C     | o Return the value of integer parameter
C     |   from one of the DIAGNOSTICS.h common blocs
C     *==========================================================*

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

C     !INPUT PARAMETERS:
C     parName   :: string used to identify which parameter to get
C     myThid    :: my Thread Id number
      CHARACTER*(*) parName
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER n

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

      n = LEN(parName)
c     write(0,'(3A,I4)')
c    &  'DIAGS_GET_PARMS_I: parName="',parName,'" , length=',n

      IF ( parName.EQ.'LAST_DIAG_ID' ) THEN
         DIAGS_GET_PARMS_I = ndiagt
      ELSE
         WRITE(msgBuf,'(4A)') 'DIAGS_GET_PARMS_I: ',
     &    ' parName="', parName, '" not known.'
         CALL PRINT_ERROR( msgBuf, myThid )
         STOP 'ABNORMAL END: S/R DIAGS_GET_PARMS_I'
      ENDIF

      RETURN
      END

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

CBOP 0
C     !ROUTINE: DIAGS_MK_UNITS

C     !INTERFACE:
      CHARACTER*16 FUNCTION DIAGS_MK_UNITS(
     I                            diagUnitsInPieces, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | FUNCTION DIAGS_MK_UNITS
C     | o Return the diagnostic units string (16c) removing
C     |   blanks from the input string
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
#include "EEPARAMS.h"

C     !INPUT PARAMETERS:
C     diagUnitsInPieces :: string for diagnostic units: in several
C                          pieces, with blanks in between
C     myThid            ::  my thread Id number
      CHARACTER*(*) diagUnitsInPieces
      INTEGER      myThid
CEOP

C     !LOCAL VARIABLES:
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER i,j,n

      DIAGS_MK_UNITS = '                '
      n = LEN(diagUnitsInPieces)

      j = 0
      DO i=1,n
       IF (diagUnitsInPieces(i:i) .NE. ' ' ) THEN
         j = j+1
         IF ( j.LE.16 ) DIAGS_MK_UNITS(j:j) = diagUnitsInPieces(i:i)
       ENDIF
      ENDDO

      IF ( j.GT.16 ) THEN
         WRITE(msgBuf,'(2A,I4,A)') '** WARNING ** ',
     &   'DIAGS_MK_UNITS: too long (',j,' >16) input string'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &       SQUEEZE_RIGHT , myThid)
         WRITE(msgBuf,'(3A)') '** WARNING ** ',
     &   'DIAGS_MK_UNITS: input=', diagUnitsInPieces
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &       SQUEEZE_RIGHT , myThid)
      ENDIF

      RETURN
      END

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

CBOP 0
C     !ROUTINE: DIAGS_MK_TITLE

C     !INTERFACE:
      CHARACTER*80 FUNCTION DIAGS_MK_TITLE(
     I                            diagTitleInPieces, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | FUNCTION DIAGS_MK_TITLE
C     | o Return the diagnostic title string (80c) removing
C     |   consecutive blanks from the input string
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
#include "EEPARAMS.h"

C     !INPUT PARAMETERS:
C     diagTitleInPieces :: string for diagnostic units: in several
C                          pieces, with blanks in between
C     myThid            ::  my Thread Id number
      CHARACTER*(*) diagTitleInPieces
      INTEGER      myThid
CEOP

C     !LOCAL VARIABLES:
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      LOGICAL flag
      INTEGER i,j,n

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

      DIAGS_MK_TITLE = '                                        '
     &               //'                                        '
      n = LEN(diagTitleInPieces)

      j = 0
      flag = .FALSE.
      DO i=1,n
       IF (diagTitleInPieces(i:i) .NE. ' ' ) THEN
         IF ( flag ) THEN
           j = j+1
           IF (j.LE.80) DIAGS_MK_TITLE(j:j) = ' '
         ENDIF
         j = j+1
         IF ( j.LE.80 ) DIAGS_MK_TITLE(j:j) = diagTitleInPieces(i:i)
         flag = .FALSE.
       ELSE
         flag = j.GE.1
       ENDIF
      ENDDO

      IF ( j.GT.80 ) THEN
         WRITE(msgBuf,'(2A,I4,A)') '** WARNING ** ',
     &   'DIAGS_MK_TITLE: too long (',j,' >80) input string'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &       SQUEEZE_RIGHT , myThid)
         WRITE(msgBuf,'(3A)') '** WARNING ** ',
     &   'DIAGS_MK_TITLE: input=', diagTitleInPieces
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &       SQUEEZE_RIGHT , myThid)
      ENDIF

      RETURN
      END
1	jmc	1.29	C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_utils.F,v 1.28 2009/01/25 17:00:20 jmc Exp $
2	edhill	1.8	C $Name: $
3
4	edhill	1.9	#include "DIAG_OPTIONS.h"
5
6	jmc	1.25	C-- File diagnostics_utils.F: General purpose support routines
7			C-- Contents:
8			C-- o GETDIAG
9			C-- o DIAGNOSTICS_COUNT
10	jmc	1.26	C-- o DIAGNOSTICS_GET_POINTERS
11	jmc	1.29	C-- o DIAGNOSTICS_SETKLEV
12	jmc	1.26	C-- o DIAGS_GET_PARMS_I (Function)
13	jmc	1.25	C-- o DIAGS_MK_UNITS (Function)
14			C-- o DIAGS_MK_TITLE (Function)
15
16	edhill	1.12	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
17			CBOP 0
18			C !ROUTINE: GETDIAG
19
20			C !INTERFACE:
21	jmc	1.18	SUBROUTINE GETDIAG(
22			I levreal, undef,
23			O qtmp,
24	jmc	1.21	I ndId, mate, ip, im, bi, bj, myThid )
25	edhill	1.12
26			C !DESCRIPTION:
27	edhill	1.13	C Retrieve averaged model diagnostic
28	jmc	1.15
29	edhill	1.12	C !USES:
30	jmc	1.18	IMPLICIT NONE
31	molod	1.3	#include "EEPARAMS.h"
32	molod	1.1	#include "SIZE.h"
33	jmc	1.15	#include "DIAGNOSTICS_SIZE.h"
34			#include "DIAGNOSTICS.h"
35	molod	1.3
36	jmc	1.18	C !INPUT PARAMETERS:
37	jmc	1.21	C levreal :: Diagnostic LEVEL
38			C undef :: UNDEFINED VALUE
39			C ndId :: DIAGNOSTIC NUMBER FROM MENU
40			C mate :: counter DIAGNOSTIC NUMBER if any ; 0 otherwise
41			C ip :: pointer to storage array location for diag.
42			C im :: pointer to storage array location for mate
43			C bi :: X-direction tile number
44			C bj :: Y-direction tile number
45			C myThid :: my thread Id number
46	jmc	1.15	_RL levreal
47	molod	1.3	_RL undef
48	jmc	1.21	INTEGER ndId, mate, ip, im
49	jmc	1.18	INTEGER bi,bj, myThid
50	jmc	1.15
51	jmc	1.18	C !OUTPUT PARAMETERS:
52	edhill	1.12	C qtmp ..... AVERAGED DIAGNOSTIC QUANTITY
53	jmc	1.18	_RL qtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54			CEOP
55	molod	1.11
56	jmc	1.18	C !LOCAL VARIABLES:
57	molod	1.3	_RL factor
58	jmc	1.18	INTEGER i, j, ipnt,ipCt
59			INTEGER lev, levCt, klev
60	molod	1.11
61	jmc	1.21	IF (ndId.GE.1) THEN
62	jmc	1.15	lev = NINT(levreal)
63	jmc	1.21	klev = kdiag(ndId)
64	jmc	1.18	IF (lev.LE.klev) THEN
65	jmc	1.15
66	jmc	1.18	IF ( mate.EQ.0 ) THEN
67			C- No counter diagnostics => average = Sum / ndiag :
68	jmc	1.15
69	jmc	1.21	ipnt = ip + lev - 1
70			factor = FLOAT(ndiag(ip,bi,bj))
71			IF (ndiag(ip,bi,bj).NE.0) factor = 1. _d 0 / factor
72	jmc	1.18
73	jmc	1.27	#ifdef ALLOW_FIZHI
74	jmc	1.18	DO j = 1,sNy+1
75			DO i = 1,sNx+1
76			IF ( qdiag(i,j,ipnt,bi,bj) .LE. undef ) THEN
77			qtmp(i,j) = qdiag(i,j,ipnt,bi,bj)*factor
78			ELSE
79			qtmp(i,j) = undef
80			ENDIF
81			ENDDO
82			ENDDO
83	jmc	1.27	#else /* ALLOW_FIZHI */
84			DO j = 1,sNy+1
85			DO i = 1,sNx+1
86			qtmp(i,j) = qdiag(i,j,ipnt,bi,bj)*factor
87			ENDDO
88			ENDDO
89			#endif /* ALLOW_FIZHI */
90	jmc	1.18
91			ELSE
92			C- With counter diagnostics => average = Sum / counter:
93
94	jmc	1.21	ipnt = ip + lev - 1
95	jmc	1.18	levCt= MIN(lev,kdiag(mate))
96	jmc	1.21	ipCt = im + levCt - 1
97	jmc	1.18	DO j = 1,sNy+1
98			DO i = 1,sNx+1
99			IF ( qdiag(i,j,ipCt,bi,bj) .NE. 0. ) THEN
100			qtmp(i,j) = qdiag(i,j,ipnt,bi,bj)
101			& / qdiag(i,j,ipCt,bi,bj)
102			ELSE
103			qtmp(i,j) = undef
104			ENDIF
105			ENDDO
106			ENDDO
107	molod	1.1
108	jmc	1.18	ENDIF
109			ENDIF
110			ENDIF
111	molod	1.3
112	jmc	1.15	RETURN
113			END
114	edhill	1.12
115			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
116
117	jmc	1.15	CBOP 0
118	jmc	1.19	C !ROUTINE: DIAGNOSTICS_COUNT
119			C !INTERFACE:
120	jmc	1.21	SUBROUTINE DIAGNOSTICS_COUNT (chardiag,
121			I biArg, bjArg, myThid)
122	jmc	1.19
123			C !DESCRIPTION:
124			C***********************************************************************
125			C routine to increment the diagnostic counter only
126			C***********************************************************************
127			C !USES:
128			IMPLICIT NONE
129
130			C == Global variables ===
131			#include "EEPARAMS.h"
132			#include "SIZE.h"
133			#include "DIAGNOSTICS_SIZE.h"
134			#include "DIAGNOSTICS.h"
135
136			C !INPUT PARAMETERS:
137			C***********************************************************************
138			C Arguments Description
139			C ----------------------
140			C chardiag :: Character expression for diag to increment the counter
141			C biArg :: X-direction tile number, or 0 if called outside bi,bj loops
142			C bjArg :: Y-direction tile number, or 0 if called outside bi,bj loops
143			C myThid :: my thread Id number
144			C***********************************************************************
145			CHARACTER*8 chardiag
146			INTEGER biArg, bjArg
147			INTEGER myThid
148			CEOP
149
150			C !LOCAL VARIABLES:
151			C ===============
152	jmc	1.21	INTEGER m, n
153			INTEGER bi, bj
154	jmc	1.28	INTEGER ipt, ndId
155	jmc	1.19	c CHARACTER*(MAX_LEN_MBUF) msgBuf
156
157	jmc	1.28	IF ( biArg.EQ.0 .AND. bjArg.EQ.0 ) THEN
158			bi = myBxLo(myThid)
159			bj = myByLo(myThid)
160			ELSE
161			bi = MIN(biArg,nSx)
162			bj = MIN(bjArg,nSy)
163			ENDIF
164
165	jmc	1.21	C-- Run through list of active diagnostics to find which counter
166			C to increment (needs to be a valid & active diagnostic-counter)
167	jmc	1.19	DO n=1,nlists
168			DO m=1,nActive(n)
169	jmc	1.21	IF ( chardiag.EQ.flds(m,n) .AND. idiag(m,n).GT.0 ) THEN
170			ipt = idiag(m,n)
171	jmc	1.28	IF (ndiag(ipt,bi,bj).GE.0) THEN
172			ndId = jdiag(m,n)
173			ipt = ipt + pdiag(n,bi,bj)*kdiag(ndId)
174	jmc	1.21	C- Increment the counter for the diagnostic
175			IF ( biArg.EQ.0 .AND. bjArg.EQ.0 ) THEN
176			DO bj=myByLo(myThid), myByHi(myThid)
177			DO bi=myBxLo(myThid), myBxHi(myThid)
178			ndiag(ipt,bi,bj) = ndiag(ipt,bi,bj) + 1
179			ENDDO
180			ENDDO
181			ELSE
182			ndiag(ipt,bi,bj) = ndiag(ipt,bi,bj) + 1
183			ENDIF
184			C- Increment is done
185			ENDIF
186	jmc	1.19	ENDIF
187			ENDDO
188			ENDDO
189
190	jmc	1.21	RETURN
191	jmc	1.19	END
192
193			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
194
195			CBOP 0
196	jmc	1.26	C !ROUTINE: DIAGNOSTICS_GET_POINTERS
197			C !INTERFACE:
198			SUBROUTINE DIAGNOSTICS_GET_POINTERS(
199			I diagName, listId,
200			O ndId, ip,
201			I myThid )
202
203			C !DESCRIPTION:
204			C ================================================================
205			C \| o Returns the diagnostic Id number and diagnostic
206			C \| pointer to storage array for a specified diagnostic.
207			C ================================================================
208			C \| Note: A diagnostics field can be stored multiple times
209			C \| (for different output frequency,phase, ...).
210			C \| operates in 2 ways:
211			C \| o listId =0 => find 1 diagnostics Id & pointer which name matches.
212			C \| o listId >0 => find the unique diagnostic Id & pointer with
213			C \| the right name and same output time as "listId" output-list
214			C \| o return ip=0 if did not find the right diagnostic;
215			C \| (ndId <>0 if diagnostic exist but output time does not match)
216			C ================================================================
217
218			C !USES:
219			IMPLICIT NONE
220			#include "EEPARAMS.h"
221			#include "SIZE.h"
222			#include "DIAGNOSTICS_SIZE.h"
223			#include "DIAGNOSTICS.h"
224
225			C !INPUT PARAMETERS:
226			C diagName :: diagnostic identificator name (8 characters long)
227			C listId :: list number that specify the output frequency
228			C myThid :: my Thread Id number
229			C !OUTPUT PARAMETERS:
230			C ndId :: diagnostics Id number (in available diagnostics list)
231			C ip :: diagnostics pointer to storage array
232
233
234			CHARACTER*8 diagName
235			INTEGER listId
236			INTEGER ndId, ip
237			INTEGER myThid
238			CEOP
239
240			C !LOCAL VARIABLES:
241			INTEGER n,m
242
243			ip = 0
244			ndId = 0
245
246			IF ( listId.LE.0 ) THEN
247			C-- select the 1rst one which name matches:
248
249			C- search for this diag. in the active 2D/3D diagnostics list
250			DO n=1,nlists
251			DO m=1,nActive(n)
252			IF ( ip.EQ.0 .AND. diagName.EQ.flds(m,n)
253			& .AND. idiag(m,n).NE.0 ) THEN
254			ip = ABS(idiag(m,n))
255			ndId = jdiag(m,n)
256			ENDIF
257			ENDDO
258			ENDDO
259
260			ELSEIF ( listId.LE.nlists ) THEN
261			C-- select the unique diagnostic with output-time identical to listId
262
263			C- search for this diag. in the active 2D/3D diagnostics list
264			DO n=1,nlists
265			IF ( ip.EQ.0
266			& .AND. freq(n) .EQ. freq(listId)
267			& .AND. phase(n).EQ.phase(listId)
268			& .AND. averageFreq(n) .EQ.averageFreq(listId)
269			& .AND. averagePhase(n).EQ.averagePhase(listId)
270			& .AND. averageCycle(n).EQ.averageCycle(listId)
271			& ) THEN
272			DO m=1,nActive(n)
273			IF ( ip.EQ.0 .AND. diagName.EQ.flds(m,n)
274			& .AND. idiag(m,n).NE.0 ) THEN
275			ip = ABS(idiag(m,n))
276			ndId = jdiag(m,n)
277			ENDIF
278			ENDDO
279			ELSEIF ( ip.EQ.0 ) THEN
280			DO m=1,nActive(n)
281			IF ( ip.EQ.0 .AND. diagName.EQ.flds(m,n)
282			& .AND. idiag(m,n).NE.0 ) THEN
283			ndId = jdiag(m,n)
284			ENDIF
285			ENDDO
286			ENDIF
287			ENDDO
288
289			ELSE
290			STOP 'DIAGNOSTICS_GET_POINTERS: invalid listId number'
291			ENDIF
292
293			RETURN
294			END
295
296			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
297
298			CBOP 0
299	jmc	1.29	C !ROUTINE: DIAGNOSTICS_SETKLEV
300
301			C !INTERFACE:
302			SUBROUTINE DIAGNOSTICS_SETKLEV(
303			I diagName, nLevDiag, myThid )
304
305			C !DESCRIPTION:
306			C ==========================================================
307			C \| S/R DIAGNOSTICS_SETKLEV
308			C \| o Define explicitly the number of level (stored in kdiag)
309			C \| of a diagnostic field. For most diagnostics, the number
310			C \| of levels is derived (in S/R SET_LEVELS) from gdiag(10)
311			C \| but occasionally one may want to set it explicitly.
312			C ==========================================================
313
314			C !USES:
315			IMPLICIT NONE
316			#include "EEPARAMS.h"
317			#include "SIZE.h"
318			#include "DIAGNOSTICS_SIZE.h"
319			#include "DIAGNOSTICS.h"
320
321			C !INPUT PARAMETERS:
322			C diagName :: diagnostic identificator name (8 characters long)
323			C nLevDiag :: number of level to set for this diagnostics field
324			C myThid :: my Thread Id number
325			CHARACTER*8 diagName
326			INTEGER nLevDiag
327			INTEGER myThid
328			CEOP
329
330			C !LOCAL VARIABLES:
331			CHARACTER*(MAX_LEN_MBUF) msgBuf
332			INTEGER n, ndId
333
334			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
335
336			C-- Check if this S/R is called from the right place ;
337			C needs to be after DIAGNOSTICS_INIT_EARLY and before DIAGNOSTICS_INIT_FIXED
338			IF ( .NOT.settingDiags ) THEN
339			WRITE(msgBuf,'(4A,I5)') 'DIAGNOSTICS_SETKLEV: ',
340			& 'diagName="', diagName, '" , nLevDiag=', nLevDiag
341			CALL PRINT_ERROR( msgBuf, myThid )
342			WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SETKLEV: ',
343			& '<== called from the WRONG place, i.e.'
344			CALL PRINT_ERROR( msgBuf, myThid )
345			WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SETKLEV: ',
346			& 'outside diagnostics setting section = from'
347			CALL PRINT_ERROR( msgBuf, myThid )
348			WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SETKLEV: ',
349			& ' Diag_INIT_EARLY down to Diag_INIT_FIXED'
350			CALL PRINT_ERROR( msgBuf, myThid )
351			STOP 'ABNORMAL END: S/R DIAGNOSTICS_SETKLEV'
352			ENDIF
353
354			C-- Find this diagnostics in the list of available diag.
355			ndId = 0
356			DO n = 1,ndiagt
357			IF ( diagName.EQ.cdiag(n) ) THEN
358			ndId = n
359			ENDIF
360			ENDDO
361			IF ( ndId.EQ.0 ) THEN
362			WRITE(msgBuf,'(4A)') 'DIAGNOSTICS_SETKLEV: ',
363			& 'diagName="', diagName, '" not known.'
364			CALL PRINT_ERROR( msgBuf, myThid )
365			STOP 'ABNORMAL END: S/R DIAGNOSTICS_SETKLEV'
366			ENDIF
367
368			C- Optional level number diagnostics (X): set number of levels
369			IF ( kdiag(ndId).EQ.0
370			& .AND. gdiag(ndId)(10:10).EQ.'X' ) THEN
371			kdiag(ndId) = nLevDiag
372			ELSEIF ( kdiag(ndId).EQ.nLevDiag
373			& .AND. gdiag(ndId)(10:10).EQ.'X' ) THEN
374			C- level number already set to same value: send warning
375			WRITE(msgBuf,'(4A,I5)') ' WARNING DIAGNOSTICS_SETKLEV: ',
376			& 'diagName="', diagName, '" , nLevDiag=', nLevDiag
377			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
378			& SQUEEZE_RIGHT , myThid )
379			WRITE(msgBuf,'(2A,I5,A)')' WARNING DIAGNOSTICS_SETKLEV:',
380			& ' level Nb (=', kdiag(ndId), ') already set.'
381			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
382			& SQUEEZE_RIGHT , myThid )
383			ELSEIF ( gdiag(ndId)(10:10).EQ.'X' ) THEN
384			C- level number already set to a different value: do not reset but stop
385			WRITE(msgBuf,'(4A,I5)') 'DIAGNOSTICS_SETKLEV: ',
386			& 'diagName="', diagName, '" , nLevDiag=', nLevDiag
387			CALL PRINT_ERROR( msgBuf, myThid )
388			WRITE(msgBuf,'(2A,I5,3A)') 'DIAGNOSTICS_SETKLEV: ',
389			& 'level Nb already set to', kdiag(ndId), ' => STOP'
390			CALL PRINT_ERROR( msgBuf, myThid )
391			ELSE
392			C- for now, do nothing but just send a warning
393			WRITE(msgBuf,'(4A,I5)') ' WARNING DIAGNOSTICS_SETKLEV: ',
394			& 'diagName="', diagName, '" , nLevDiag=', nLevDiag
395			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
396			& SQUEEZE_RIGHT , myThid )
397			WRITE(msgBuf,'(2A,I5,3A)') ' WARNING will set level Nb',
398			& ' from diagCode(ndId=', ndId, ')="', gdiag(ndId)(1:10), '"'
399			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
400			& SQUEEZE_RIGHT , myThid )
401			WRITE(msgBuf,'(4A)') ' WARNING DIAGNOSTICS_SETKLEV',
402			& '("', diagName, '") <== Ignore this call.'
403			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
404			& SQUEEZE_RIGHT , myThid )
405			ENDIF
406
407			RETURN
408			END
409
410			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
411
412			CBOP 0
413	jmc	1.26	C !ROUTINE: DIAGS_GET_PARMS_I
414
415			C !INTERFACE:
416			INTEGER FUNCTION DIAGS_GET_PARMS_I(
417			I parName, myThid )
418
419			C !DESCRIPTION:
420			C ==========================================================
421			C \| FUNCTION DIAGS_GET_PARMS_I
422			C \| o Return the value of integer parameter
423			C \| from one of the DIAGNOSTICS.h common blocs
424			C ==========================================================
425
426			C !USES:
427			IMPLICIT NONE
428			#include "EEPARAMS.h"
429			#include "SIZE.h"
430			#include "DIAGNOSTICS_SIZE.h"
431			#include "DIAGNOSTICS.h"
432
433			C !INPUT PARAMETERS:
434			C parName :: string used to identify which parameter to get
435			C myThid :: my Thread Id number
436			CHARACTER() parName
437			INTEGER myThid
438			CEOP
439
440			C !LOCAL VARIABLES:
441			CHARACTER*(MAX_LEN_MBUF) msgBuf
442			INTEGER n
443
444			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
445
446			n = LEN(parName)
447			c write(0,'(3A,I4)')
448			c & 'DIAGS_GET_PARMS_I: parName="',parName,'" , length=',n
449
450			IF ( parName.EQ.'LAST_DIAG_ID' ) THEN
451			DIAGS_GET_PARMS_I = ndiagt
452			ELSE
453			WRITE(msgBuf,'(4A)') 'DIAGS_GET_PARMS_I: ',
454			& ' parName="', parName, '" not known.'
455			CALL PRINT_ERROR( msgBuf, myThid )
456			STOP 'ABNORMAL END: S/R DIAGS_GET_PARMS_I'
457			ENDIF
458
459			RETURN
460			END
461
462			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
463
464			CBOP 0
465	jmc	1.17	C !ROUTINE: DIAGS_MK_UNITS
466
467			C !INTERFACE:
468	jmc	1.21	CHARACTER*16 FUNCTION DIAGS_MK_UNITS(
469	jmc	1.17	I diagUnitsInPieces, myThid )
470
471			C !DESCRIPTION:
472			C ==========================================================
473			C \| FUNCTION DIAGS_MK_UNITS
474	jmc	1.21	C \| o Return the diagnostic units string (16c) removing
475	jmc	1.17	C \| blanks from the input string
476			C ==========================================================
477
478			C !USES:
479			IMPLICIT NONE
480			#include "EEPARAMS.h"
481
482			C !INPUT PARAMETERS:
483	jmc	1.21	C diagUnitsInPieces :: string for diagnostic units: in several
484	jmc	1.17	C pieces, with blanks in between
485			C myThid :: my thread Id number
486			CHARACTER() diagUnitsInPieces
487			INTEGER myThid
488			CEOP
489
490			C !LOCAL VARIABLES:
491			CHARACTER*(MAX_LEN_MBUF) msgBuf
492			INTEGER i,j,n
493
494	jmc	1.29	DIAGS_MK_UNITS = ' '
495	jmc	1.17	n = LEN(diagUnitsInPieces)
496	jmc	1.21
497	jmc	1.17	j = 0
498			DO i=1,n
499			IF (diagUnitsInPieces(i:i) .NE. ' ' ) THEN
500			j = j+1
501			IF ( j.LE.16 ) DIAGS_MK_UNITS(j:j) = diagUnitsInPieces(i:i)
502			ENDIF
503			ENDDO
504
505			IF ( j.GT.16 ) THEN
506	jmc	1.29	WRITE(msgBuf,'(2A,I4,A)') ' WARNING ',
507	jmc	1.17	& 'DIAGS_MK_UNITS: too long (',j,' >16) input string'
508			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
509			& SQUEEZE_RIGHT , myThid)
510	jmc	1.29	WRITE(msgBuf,'(3A)') ' WARNING ',
511	jmc	1.17	& 'DIAGS_MK_UNITS: input=', diagUnitsInPieces
512			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
513			& SQUEEZE_RIGHT , myThid)
514			ENDIF
515
516			RETURN
517			END
518	jmc	1.23
519			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
520
521			CBOP 0
522			C !ROUTINE: DIAGS_MK_TITLE
523
524			C !INTERFACE:
525			CHARACTER*80 FUNCTION DIAGS_MK_TITLE(
526			I diagTitleInPieces, myThid )
527
528			C !DESCRIPTION:
529			C ==========================================================
530			C \| FUNCTION DIAGS_MK_TITLE
531			C \| o Return the diagnostic title string (80c) removing
532			C \| consecutive blanks from the input string
533			C ==========================================================
534
535			C !USES:
536			IMPLICIT NONE
537			#include "EEPARAMS.h"
538
539			C !INPUT PARAMETERS:
540			C diagTitleInPieces :: string for diagnostic units: in several
541			C pieces, with blanks in between
542			C myThid :: my Thread Id number
543			CHARACTER() diagTitleInPieces
544			INTEGER myThid
545			CEOP
546
547			C !LOCAL VARIABLES:
548			CHARACTER*(MAX_LEN_MBUF) msgBuf
549			LOGICAL flag
550			INTEGER i,j,n
551
552	molod	1.22	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
553	jmc	1.23
554			DIAGS_MK_TITLE = ' '
555			& //' '
556			n = LEN(diagTitleInPieces)
557
558			j = 0
559			flag = .FALSE.
560			DO i=1,n
561			IF (diagTitleInPieces(i:i) .NE. ' ' ) THEN
562			IF ( flag ) THEN
563			j = j+1
564			IF (j.LE.80) DIAGS_MK_TITLE(j:j) = ' '
565			ENDIF
566			j = j+1
567			IF ( j.LE.80 ) DIAGS_MK_TITLE(j:j) = diagTitleInPieces(i:i)
568			flag = .FALSE.
569			ELSE
570			flag = j.GE.1
571			ENDIF
572			ENDDO
573
574			IF ( j.GT.80 ) THEN
575	jmc	1.29	WRITE(msgBuf,'(2A,I4,A)') ' WARNING ',
576	jmc	1.23	& 'DIAGS_MK_TITLE: too long (',j,' >80) input string'
577			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
578			& SQUEEZE_RIGHT , myThid)
579	jmc	1.29	WRITE(msgBuf,'(3A)') ' WARNING ',
580	jmc	1.23	& 'DIAGS_MK_TITLE: input=', diagTitleInPieces
581			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
582			& SQUEEZE_RIGHT , myThid)
583			ENDIF
584
585			RETURN
586			END