--- MITgcm/pkg/diagnostics/diagnostics_out.F 2011/06/12 13:58:33 1.51 +++ MITgcm/pkg/diagnostics/diagnostics_out.F 2017/01/11 00:22:48 1.62 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.51 2011/06/12 13:58:33 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.62 2017/01/11 00:22:48 gforget Exp $ C $Name: $ #include "DIAG_OPTIONS.h" @@ -41,10 +41,6 @@ C !FUNCTIONS: INTEGER ILNBLNK EXTERNAL ILNBLNK -#ifdef ALLOW_FIZHI - _RL getcon - EXTERNAL getcon -#endif C !LOCAL VARIABLES: C i,j,k :: loop indices @@ -52,55 +48,61 @@ C lm :: loop index (averageCycle) C md :: field number in the list "listId". C ndId :: diagnostics Id number (in available diagnostics list) -C mate :: counter mate 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 :: thread-shared temporary array (needs to be in common block): -C to write a diagnostic field to file, copy it first from (big) -C diagnostic storage qdiag into it. - COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1 +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, ip, im - INTEGER mate, mVec + INTEGER md, ndId, nn, ip, im + INTEGER mate, mDbl, mVec + INTEGER ppFld, isComputed CHARACTER*10 gcode - _RL undef - _RL tmpLev - INTEGER iLen - INTEGER nLevOutp + _RL undefRL + INTEGER nLevOutp, kLev + INTEGER iLen,jLen 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 - INTEGER ll, llMx, jj, jjMx CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn - LOGICAL useMissingValue - REAL*8 misValLoc #endif /* ALLOW_MNC */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C--- set file properties ioUnit= standardMessageUnit - undef = UNSET_RL -#ifdef ALLOW_FIZHI - IF ( useFIZHI ) undef = getcon('UNDEF') -#endif + undefRL = misValFlt(listId) + WRITE(suff,'(I10.10)') myIter iLen = ILNBLNK(fnames(listId)) WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10) + IF ( diag_mdsio.AND.(diagMdsDir.NE.' ') ) THEN + jLen = ILNBLNK(diagMdsDir) + WRITE( fn, '(5A)' ) diagMdsDir(1:jLen),'/', + & fnames(listId)(1:iLen),'.',suff(1:10) + ENDIF C- for now, if integrate vertically, output field has just 1 level: nLevOutp = nlevels(listId) IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1 @@ -126,8 +128,8 @@ timeRec(1) = (timeRec(1)-phase(listId))/freq(listId) i = INT( timeRec(1) ) IF ( timeRec(1).LT.0. ) THEN - tmpLev = FLOAT(i) - IF ( timeRec(1).NE.tmpLev ) i = i - 1 + 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 @@ -137,8 +139,8 @@ timeRec(1) = (timeRec(1)-baseTime)/deltaTClock i = NINT( timeRec(1) ) C if just half way, NINT will return the next time-step: correct this - tmpLev = FLOAT(i) - 0.5 _d 0 - IF ( timeRec(1).EQ.tmpLev ) i = i - 1 + 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 @@ -147,60 +149,67 @@ c timeRec(i) = timeRec(i)/deltaTClock c ENDDO -#ifdef ALLOW_MNC -C-- this is a trick to reverse the order of the loops on md (= field) -C and lm (= averagePeriod): binary output: lm loop inside md loop ; -C mnc ouput: md loop inside lm loop. - IF (useMNC .AND. diag_mnc) THEN - jjMx = averageCycle(listId) - llMx = 1 - ELSE - jjMx = 1 - llMx = averageCycle(listId) - ENDIF - DO jj=1,jjMx +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 - misValLoc = undef - IF ( misvalFlt(listId).NE.UNSET_RL ) - & misValLoc = misvalFlt(listId) CALL DIAGNOSTICS_MNC_SET( - I nLevOutp, listId, jj, - O diag_mnc_bn, useMissingValue, - I misValLoc, myTime, myIter, myThid ) + 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 : -#ifdef ALLOW_MNC - DO ll=1,llMx - lm = jj+ll-1 -#else - DO lm=1,averageCycle(listId) -#endif 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 ( ndiag(ip,1,1).EQ.0 ) THEN + 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)') @@ -242,11 +251,27 @@ ELSE C- diagnostics is not empty : + isComputed = 0 IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN - WRITE(ioUnit,'(A,I6,3A,I8,2A)') + 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), @@ -266,31 +291,52 @@ ENDIF ENDIF - IF ( fflags(listId)(2:2).NE.' ' ) THEN -C- get all the levels (for vertical post-processing) + 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,kdiag(ndId) - tmpLev = k - CALL GETDIAG( - I tmpLev,undef, + 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) + 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 only selected levels: +C- get all the levels (for vertical post-processing) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) - DO k = 1,nlevels(listId) - CALL GETDIAG( - I levs(k,listId),undef, - O qtmp1(1-OLx,1-OLy,k,bi,bj), - I ndId,mate,ip,im,bi,bj,myThid) - ENDDO + 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----------------------------------------------------------------------- @@ -302,8 +348,8 @@ 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, - I undef, myTime, myIter, myThid ) + U qtmp1, qtmp2, + I undefRL, myTime, myIter, myThid ) ELSE WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ', & 'INTERP_VERT not allowed in this config' @@ -316,7 +362,25 @@ CALL DIAGNOSTICS_SUM_LEVELS( I listId, md, ndId, ip, im, lm, U qtmp1, - I undef, myTime, myIter, myThid ) + 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 @@ -326,40 +390,49 @@ C- write to binary file, using MDSIO pkg: IF ( diag_mdsio ) THEN - nRec = lm + (md-1)*averageCycle(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 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, - I diag_mnc_bn, - I useMissingValue, misValLoc, - I qtmp1, - I myTime, myIter, myThid ) + I NrMax, nLevOutp, listId, ndId, mate, + I diag_mnc_bn, qtmp2, + I undefRL, myTime, myIter, myThid ) + ENDIF ENDIF #endif /* ALLOW_MNC */ -C-- end loop on lm (or ll if ALLOW_MNC) counter - ENDDO C-- end of Processing Fld # md ENDIF ENDDO -#ifdef ALLOW_MNC -C-- end loop on jj counter +C-- end loop on lm counter (= averagePeriod) ENDDO -#endif #ifdef ALLOW_MDSIO IF (diag_mdsio) THEN @@ -367,10 +440,11 @@ 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 = nfields(listId)*averageCycle(listId) + nRec = averageCycle(listId)*nfields(listId) CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE., & 0, 0, nLevOutp, ' ', - & nfields(listId), flds(1,listId), nTimRec, timeRec, + & nfields(listId), flds(1,listId), + & nTimRec, timeRec, undefRL, & nRec, myIter, myThid) ENDIF #endif /* ALLOW_MDSIO */