C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.64 2017/07/23 00:24:18 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, myTime, myIter, 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 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,jLen INTEGER ioUnit CHARACTER*(MAX_LEN_FNAM) fn CHARACTER*(10) 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) IF ( rwSuffixType.EQ.0 ) THEN WRITE(suff,'(I10.10)') myIter ELSE CALL RW_GET_SUFFIX( suff, myTime, myIter, myThid ) ENDIF iLen = ILNBLNK(fnames(listId)) WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff IF ( diag_mdsio.AND.(diagMdsDir.NE.' ') ) THEN jLen = ILNBLNK(diagMdsDir) WRITE( fn, '(5A)' ) diagMdsDir(1:jLen),'/', & fnames(listId)(1:iLen),'.',suff 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 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 = ABS(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-|--+----|