--- MITgcm/pkg/diagnostics/diagnostics_utils.F 2004/02/26 22:20:36 1.5 +++ MITgcm/pkg/diagnostics/diagnostics_utils.F 2008/02/05 15:31:19 1.25 @@ -1,319 +1,402 @@ - subroutine getdiag (myThid,lev,ipoint,undef,qtmp) -C*********************************************************************** -C PURPOSE +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagnostics_utils.F,v 1.25 2008/02/05 15:31:19 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 DIAGS_MK_UNITS (Function) +C-- o DIAGS_MK_TITLE (Function) +C-- o DIAGNOSTICS_GET_POINTERS + +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 INPUT: -C lev ..... Diagnostic LEVEL -C ipoint ..... DIAGNOSTIC NUMBER FROM MENU -C undef ..... UNDEFINED VALUE -C bi ..... X-direction process(or) number -C bj ..... Y-direction process(or) number -C -C OUTPUT: -C qtmp ..... AVERAGED DIAGNOSTIC QUANTITY -C -C*********************************************************************** - implicit none +C !USES: + IMPLICIT NONE #include "EEPARAMS.h" -#include "CPP_OPTIONS.h" #include "SIZE.h" +#include "DIAGNOSTICS_SIZE.h" +#include "DIAGNOSTICS.h" -#ifdef ALLOW_FIZHI -#include "fizhi_SIZE.h" -#else - integer Nrphys - parameter (Nrphys=1) -#endif - -#include "diagnostics_SIZE.h" -#include "diagnostics.h" - - integer myThid,lev,ipoint +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 - _RL qtmp(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nr+Nrphys,Nsx,Nsy) + INTEGER ndId, mate, ip, im + INTEGER bi,bj, myThid - _RL factor - integer i,j,ipnt,klev - integer bi,bj +C !OUTPUT PARAMETERS: +C qtmp ..... AVERAGED DIAGNOSTIC QUANTITY + _RL qtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) +CEOP - if (ipoint.lt.1) go to 999 +C !LOCAL VARIABLES: + _RL factor + INTEGER i, j, ipnt,ipCt + INTEGER lev, levCt, klev - klev = kdiag(ipoint) - if(klev.ge.lev) then - ipnt = idiag(ipoint) + lev - 1 - factor = 1.0 - if(ndiag(ipoint).ne.0) factor = 1.0/ndiag(ipoint) - - do bj=myByLo(myThid), myByHi(myThid) - do bi=myBxLo(myThid), myBxHi(myThid) - - do j = 1,sNy - do i = 1,sNx - if( qdiag(i,j,ipnt,bi,bj).ne.undef ) then - qtmp(i,j,lev,bi,bj) = qdiag(i,j,ipnt,bi,bj)*factor - else - qtmp(i,j,lev,bi,bj) = undef - endif - enddo - enddo - - enddo - enddo - - endif - - 999 return - end - - subroutine getdiag2 (myThid,lev,ipoint,undef,qtmp) -C*********************************************************************** -C PURPOSE -C Retrieve averaged model diagnostic -C INPUT: -C lev ..... Diagnostic LEVEL -C ipoint ..... DIAGNOSTIC NUMBER FROM MENU -C undef ..... UNDEFINED VALUE -C -C OUTPUT: -C qtmp ..... AVERAGED DIAGNOSTIC QUANTITY -C -C*********************************************************************** - implicit none + 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 + + 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 +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 -#include "EEPARAMS.h" -#include "CPP_OPTIONS.h" -#include "SIZE.h" + ENDIF + ENDIF + ENDIF -#ifdef ALLOW_FIZHI -#include "fizhi_SIZE.h" -#else - integer Nrphys - parameter (Nrphys=1) -#endif + RETURN + END -#include "diagnostics_SIZE.h" -#include "diagnostics.h" +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - integer myThid,lev,ipoint - _RL undef - _RL qtmp(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nr+Nrphys,Nsx,Nsy) +CBOP 0 +C !ROUTINE: DIAGNOSTICS_COUNT +C !INTERFACE: + SUBROUTINE DIAGNOSTICS_COUNT (chardiag, + I biArg, bjArg, myThid) - integer i,j,ipnt,klev - integer bi,bj +C !DESCRIPTION: +C*********************************************************************** +C routine to increment the diagnostic counter only +C*********************************************************************** +C !USES: + IMPLICIT NONE - if (ipoint.lt.1) go to 999 +C == Global variables === +#include "EEPARAMS.h" +#include "SIZE.h" +#include "DIAGNOSTICS_SIZE.h" +#include "DIAGNOSTICS.h" - klev = kdiag(ipoint) - if(klev.ge.lev) then - ipnt = idiag(ipoint) + lev - 1 - - do bj=myByLo(myThid), myByHi(myThid) - do bi=myBxLo(myThid), myBxHi(myThid) - - do j = 1,sNy - do i = 1,sNx - if( qdiag(i,j,ipnt,bi,bj).ne.undef ) then - qtmp(i,j,lev,bi,bj) = qdiag(i,j,ipnt,bi,bj) - else - qtmp(i,j,lev,bi,bj) = undef - endif - enddo - enddo - - enddo - enddo - - endif - - 999 return - end - subroutine clrindx (myThid,listnum) +C !INPUT PARAMETERS: C*********************************************************************** -C -C PURPOSE -C DRIVER TO CLEAR DIAGNOSTICS SPECIFIED IN DIAGNOSTIC INDEX LIST -C -C ARGUMENT DESCRIPTION -C listnum .... diagnostics list number -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 +c CHARACTER*(MAX_LEN_MBUF) msgBuf + +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,1,1).GE.0) THEN +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 + bi = MIN(biArg,nSx) + bj = MIN(bjArg,nSy) + 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-|--+----| - implicit none +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" -#include "CPP_OPTIONS.h" -#include "SIZE.h" -#ifdef ALLOW_FIZHI -#include "fizhi_SIZE.h" -#else - integer Nrphys - parameter (Nrphys=1) -#endif - -#include "diagnostics_SIZE.h" -#include "diagnostics.h" - - integer myThid, listnum - - integer m, n - character*8 parms1 - character*1 parse1(8) - character*3 mate_index - integer mate - - equivalence ( parms1 , parse1(1) ) - equivalence ( mate_index , parse1(6) ) - - do n=1,nfields(listnum) - do m=1,ndiagt - if( flds(n,listnum).eq.cdiag(m) .and. idiag(m).ne.0 ) then - call clrdiag (myThid,m) - -c Check for Counter Diagnostic -c ---------------------------- - parms1 = gdiag(m) - if( parse1(5).eq.'C' ) then - read (mate_index,100) mate - call clrdiag (myThid,mate) - endif - endif - enddo - enddo - - 100 format(i3) - RETURN - END - - - subroutine clrdiag (myThid,index) -C*********************************************************************** -C PURPOSE -C ZERO OUT MODEL DIAGNOSTIC ARRAY ELEMENTS -C*********************************************************************** - - implicit none +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" -#include "CPP_OPTIONS.h" -#include "SIZE.h" -#ifdef ALLOW_FIZHI -#include "fizhi_SIZE.h" -#else - integer Nrphys - parameter (Nrphys=1) -#endif - -#include "diagnostics_SIZE.h" -#include "diagnostics.h" - - integer myThid, index - - integer bi,bj - integer i,j,k - -C ********************************************************************** -C **** SET DIAGNOSTIC AND COUNTER TO ZERO **** -C ********************************************************************** - - do bj=myByLo(myThid), myByHi(myThid) - do bi=myBxLo(myThid), myBxHi(myThid) - do k = 1,kdiag(index) - do j = 1,sNy - do i = 1,sNx - qdiag(i,j,idiag(index)+k-1,bi,bj) = 0.0 - enddo - enddo - enddo - enddo - enddo +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 - ndiag(index) = 0 + RETURN + END - return - end +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - subroutine setdiag (myThid,num,ndiagmx) -C*********************************************************************** -C -C PURPOSE -C SET POINTER LOCATIONS, NAMES, LEVELS and TITLES FOR DIAGNOSTIC NUM -C -C*********************************************************************** +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 *================================================================* - implicit none -#include "CPP_OPTIONS.h" +C !USES: + IMPLICIT NONE +#include "EEPARAMS.h" #include "SIZE.h" +#include "DIAGNOSTICS_SIZE.h" +#include "DIAGNOSTICS.h" -#ifdef ALLOW_FIZHI -#include "fizhi_SIZE.h" -#else - integer Nrphys - parameter (Nrphys=1) -#endif - -#include "diagnostics_SIZE.h" -#include "diagnostics.h" - - integer num,myThid,ndiagmx - integer ipointer - - DATA IPOINTER / 1 / - - character*8 parms1 - character*1 parse1(8) - character*3 mate_index - integer mate - - equivalence ( parms1 , parse1(1) ) - equivalence ( mate_index , parse1(6) ) - -C ********************************************************************** -C **** SET POINTERS FOR DIAGNOSTIC NUM **** -C ********************************************************************** - - parms1 = gdiag(num) - - IF( IDIAG(NUM).EQ.0 ) THEN - if(ndiagmx+kdiag(num).gt.numdiags) then - write(6,4000)num,cdiag(num) - else - IDIAG(NUM) = IPOINTER - IPOINTER = IPOINTER + KDIAG(NUM) - ndiagmx = ndiagmx + KDIAG(NUM) - if(myThid.eq.1) WRITE(6,2000)KDIAG(NUM),NUM,CDIAG(NUM),ndiagmx - endif - ELSE - if(myThid.eq.1) WRITE(6,3000) NUM, CDIAG(NUM) - ENDIF +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 -c Check for Counter Diagnostic -c ---------------------------- - if( parse1(5).eq.'C') then - read (mate_index,100) mate - - IF( IDIAG(mate).EQ.0 ) THEN - if(ndiagmx+kdiag(num).gt.numdiags) then - write(6,5000)num,cdiag(num) - else - IDIAG(mate) = IPOINTER - IPOINTER = IPOINTER + KDIAG(mate) - ndiagmx = ndiagmx + KDIAG(mate) - if(myThid.eq.1)WRITE(6,2000)KDIAG(mate),mate,CDIAG(mate),ndiagmx - endif ELSE - if(myThid.eq.1) WRITE(6,3000) mate, CDIAG(mate) + STOP 'DIAGNOSTICS_GET_POINTERS: invalid listId number' ENDIF - endif RETURN - - 100 format(i3) - 2000 FORMAT(1X,'Allocating ',I2,' Level(s) for Diagnostic # ',I3, - . ' (',A8,'), Total Number of Diagnostics: ',I5) - 3000 FORMAT(1X,'Diagnostic # ',I3,' (',A8,') has already been set') - 4000 FORMAT(1X,'Unable to allocate space for Diagnostic # ',I3, - . ' (',A8,')') - 5000 FORMAT(1X,'Unable to allocate space for Counter Diagnostic # ', - . I3,' (',A8,')',' WARNING - Diag will not accumulate properly') END