C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_utils.F,v 1.36 2017/07/23 00:24:18 jmc Exp $
C $Name: $
#include "DIAG_OPTIONS.h"
C-- File diagnostics_utils.F: General purpose support routines
C-- Contents:
C-- o DIAGNOSTICS_COUNT
C-- o DIAGNOSTICS_GET_DIAG
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-- o DIAGS_RENAMED (Function)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C !ROUTINE: DIAGNOSTICS_COUNT
C !INTERFACE:
SUBROUTINE DIAGNOSTICS_COUNT( diagName,
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 diagName :: name of diagnostic 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 diagName
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 ( diagName.EQ.flds(m,n) .AND. idiag(m,n).GT.0 ) THEN
ipt = idiag(m,n)
IF (ndiag(ipt,bi,bj).GE.0) THEN
ndId = ABS(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_DIAG
C !INTERFACE:
SUBROUTINE DIAGNOSTICS_GET_DIAG(
I kl, undefRL,
O qtmp,
I ndId, mate, ip, im, bi, bj, myThid )
C !DESCRIPTION:
C Retrieve time-averaged (or snap-shot) diagnostic field
C !USES:
IMPLICIT NONE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"
C !INPUT PARAMETERS:
C kl :: level selection: >0 : single selected lev ; =0 : all kdiag levels
C undefRL :: undefined "_RL" type value
C ndId :: diagnostic Id number (in available diagnostics list)
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
INTEGER kl
_RL undefRL
INTEGER ndId, mate, ip, im
INTEGER bi, bj, myThid
C !OUTPUT PARAMETERS:
C qtmp :: time-averaged (or snap-shot) diagnostic field
_RL qtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,*)
CEOP
C !LOCAL VARIABLES:
_RL factor
INTEGER i, j, ipnt, ipCt
INTEGER k, kd, km, kLev
IF (ndId.GE.1) THEN
kLev = kdiag(ndId)
IF ( kl.GE.1 .AND. kl.LE.kLev ) THEN
kLev = 1
ELSEIF ( kl.NE.0 ) THEN
kLev = 0
ENDIF
DO k = 1,kLev
kd = k
IF ( kl.GE.1 ) kd = kl
IF ( mate.EQ.0 ) THEN
C- No counter diagnostics => average = Sum / ndiag :
ipnt = ip + kd - 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. undefRL ) THEN
qtmp(i,j,k) = qdiag(i,j,ipnt,bi,bj)*factor
ELSE
qtmp(i,j,k) = undefRL
ENDIF
ENDDO
ENDDO
#else /* ALLOW_FIZHI */
DO j = 1,sNy+1
DO i = 1,sNx+1
qtmp(i,j,k) = qdiag(i,j,ipnt,bi,bj)*factor
ENDDO
ENDDO
#endif /* ALLOW_FIZHI */
ELSE
C- With counter diagnostics => average = Sum / counter:
ipnt = ip + kd - 1
km = MIN(kd,kdiag(mate))
ipCt = im + km - 1
DO j = 1,sNy+1
DO i = 1,sNx+1
IF ( qdiag(i,j,ipCt,bi,bj) .NE. 0. ) THEN
qtmp(i,j,k) = qdiag(i,j,ipnt,bi,bj)
& / qdiag(i,j,ipCt,bi,bj)
ELSE
qtmp(i,j,k) = undefRL
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDIF
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 = ABS(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 = ABS(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 = ABS(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-|--+----|
_BEGIN_MASTER( myThid)
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 ( diag_pkgStatus.NE.ready2setDiags ) THEN
CALL DIAGNOSTICS_STATUS_ERROR( 'DIAGNOSTICS_SETKLEV',
& ' ', diagName, ready2setDiags, myThid )
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
_END_MASTER( myThid)
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,nbc
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
nbc = j
IF ( nbc.GT.16 ) THEN
C- try to reduce length by changing m^2 & m^3 to m2 & m3:
DIAGS_MK_UNITS = ' '
j = 0
DO i=1,n
IF ( diagUnitsInPieces(i:i) .NE. ' ' ) THEN
IF ( j.GE.1 .AND. nbc.GT.16 .AND.
& diagUnitsInPieces(i:i).EQ.'^' ) THEN
IF ( diagUnitsInPieces(i-1:i-1).EQ.'m' ) THEN
nbc = nbc - 1
ELSE
j = j+1
IF ( j.LE.16 ) DIAGS_MK_UNITS(j:j) = diagUnitsInPieces(i:i)
ENDIF
ELSE
j = j+1
IF ( j.LE.16 ) DIAGS_MK_UNITS(j:j) = diagUnitsInPieces(i:i)
ENDIF
ENDIF
ENDDO
ENDIF
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
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C !ROUTINE: DIAGS_RENAMED
C !INTERFACE:
CHARACTER*8 FUNCTION DIAGS_RENAMED(
I diagName, myThid )
C !DESCRIPTION:
C *==========================================================*
C | FUNCTION DIAGS_RENAMED
C | o In case of an old diagnostics name,
C | provides the corresponding new name
C *==========================================================*
C !USES:
IMPLICIT NONE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"
C !INPUT PARAMETERS:
C diagName :: name of diagnostic to rename (or not)
C myThid :: my Thread Id number
CHARACTER*8 diagName
INTEGER myThid
CEOP
C !LOCAL VARIABLES:
CHARACTER*8 newName
CHARACTER*(MAX_LEN_MBUF) msgBuf
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
newName = blkName
IF ( useSEAICE ) THEN
IF ( diagName .EQ. 'SIfu ' ) newName = 'oceTAUX '
IF ( diagName .EQ. 'SIfv ' ) newName = 'oceTAUY '
IF ( diagName .EQ. 'SIuwind ' ) newName = 'EXFuwind'
IF ( diagName .EQ. 'SIvwind ' ) newName = 'EXFvwind'
IF ( diagName .EQ. 'SIsigI ' ) newName = 'SIsig1 '
IF ( diagName .EQ. 'SIsigII ' ) newName = 'SIsig2 '
ENDIF
IF ( newName.EQ.blkName ) THEN
DIAGS_RENAMED = diagName
ELSE
DIAGS_RENAMED = newName
WRITE(msgBuf,'(6A)') '** WARNING ** (DIAGS_RENAMED):',
& ' diagnostics "', diagName, '" replaced by "', newName, '"'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT , myThid )
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , myThid )
ENDIF
RETURN
END