--- MITgcm/pkg/diagnostics/diagnostics_out.F 2005/06/26 16:51:49 1.15 +++ MITgcm/pkg/diagnostics/diagnostics_out.F 2008/05/22 09:53:21 1.38 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.15 2005/06/26 16:51:49 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.38 2008/05/22 09:53:21 mlosch Exp $ C $Name: $ #include "DIAG_OPTIONS.h" @@ -26,13 +26,8 @@ #include "DIAGNOSTICS_SIZE.h" #include "DIAGNOSTICS.h" -#ifdef ALLOW_FIZHI -#include "fizhi_SIZE.h" -#else - INTEGER Nrphys - PARAMETER (Nrphys=0) -#endif - + INTEGER NrMax + PARAMETER( NrMax = numLevels ) C !INPUT PARAMETERS: C listId :: Diagnostics list number being written @@ -45,19 +40,27 @@ C !LOCAL VARIABLES: C i,j,k :: loop indices +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 - INTEGER i, j, k +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 + _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy) + + INTEGER i, j, k, lm, klev INTEGER bi, bj INTEGER md, ndId, ip, im INTEGER mate, mVec - CHARACTER*8 parms1 - CHARACTER*3 mate_index - _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy) + CHARACTER*10 gcode _RL undef, getcon + _RL tmpLev EXTERNAL getcon INTEGER ILNBLNK EXTERNAL ILNBLNK @@ -67,11 +70,13 @@ CHARACTER*(MAX_LEN_FNAM) fn CHARACTER*(MAX_LEN_MBUF) suff CHARACTER*(MAX_LEN_MBUF) msgBuf + INTEGER prec, nRec +#ifdef ALLOW_MDSIO LOGICAL glf +#endif #ifdef ALLOW_MNC INTEGER ii CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn - CHARACTER*(5) ctmp INTEGER CW_DIMS, NLEN PARAMETER ( CW_DIMS = 10 ) PARAMETER ( NLEN = 80 ) @@ -79,20 +84,39 @@ CHARACTER*(NLEN) dn(CW_DIMS) CHARACTER*(NLEN) d_cw_name CHARACTER*(NLEN) dn_blnk - _RS ztmp(Nr+Nrphys) +#ifdef DIAG_MNC_COORD_NEEDSWORK + CHARACTER*(5) ctmp + _RS ztmp(NrMax) +#endif + REAL*8 misvalLoc + REAL*8 misval_r8(2) + REAL*4 misval_r4(2) + INTEGER misvalIntLoc, misval_int(2) #endif /* ALLOW_MNC */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| ioUnit= standardMessageUnit undef = getcon('UNDEF') - glf = globalFiles WRITE(suff,'(I10.10)') myIter ilen = ILNBLNK(fnames(listId)) WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10) #ifdef ALLOW_MNC IF (useMNC .AND. diag_mnc) THEN +#ifdef DIAGNOSTICS_MISSING_VALUE +C Handle missing value attribute (land points) + misvalLoc = undef +C Defaults to UNSET_I + misvalIntLoc = UNSET_I + DO ii=1,2 +C misval_r4(ii) = UNSET_FLOAT4 +C misval_r8(ii) = UNSET_FLOAT8 + misval_r4(ii) = misvalLoc + misval_r8(ii) = misvalLoc + misval_int(ii) = UNSET_I + ENDDO +#endif /* DIAGNOSTICS_MISSING_VALUE */ DO i = 1,MAX_LEN_FNAM diag_mnc_bn(i:i) = ' ' ENDDO @@ -105,6 +129,11 @@ CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid) CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid) CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid) + CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',myIter,myThid) + +C NOTE: at some point it would be a good idea to add a time_bounds +C variable that has dimension (2,T) and clearly denotes the +C beginning and ending times for each diagnostics period dn(1)(1:NLEN) = dn_blnk(1:NLEN) WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId) @@ -112,20 +141,38 @@ ib(1) = 1 ie(1) = nlevels(listId) - CALL MNC_CW_ADD_GNAME('diag_levels', 1, + CALL MNC_CW_ADD_GNAME('diag_levels', 1, & dim, dn, ib, ie, myThid) - CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels', + CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels', & 0,0, myThid) CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description', & 'Idicies of vertical levels within the source arrays', & myThid) - +#ifdef DIAGNOSTICS_MISSING_VALUE + CALL MNC_CW_VATTR_MISSING('diag_levels', 0, + I misval_r8, misval_r4, misval_int, + I myThid ) +#endif /* DIAGNOSTICS_MISSING_VALUE */ + CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0, & 'diag_levels', levs(1,listId), myThid) CALL MNC_CW_DEL_VNAME('diag_levels', myThid) CALL MNC_CW_DEL_GNAME('diag_levels', myThid) +#ifdef DIAG_MNC_COORD_NEEDSWORK +C This part has been placed in an #ifdef because, as its currently +C written, it will only work with variables defined on a dynamics +C grid. As we start using diagnostics for physics grids, ice +C levels, land levels, etc. the different vertical coordinate +C dimensions will have to be taken into account. + +C 20051021 JMC & EH3 : We need to extend this so that a few +C variables each defined on different grids do not have the same +C vertical dimension names so we should be using a pattern such +C as: Z[uml]td000000 where the 't' is the type as specified by +C gdiag(10) + C Now define: Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx ctmp(1:5) = 'mul ' DO i = 1,3 @@ -137,8 +184,8 @@ C The following three ztmp() loops should eventually be modified C to reflect the fractional nature of levs(j,l) -- they should C do something like: -C ztmp(j) = rC(INT(FLOOR(levs(j,l)))) -C + ( rC(INT(FLOOR(levs(j,l)))) +C ztmp(j) = rC(INT(FLOOR(levs(j,l)))) +C + ( rC(INT(FLOOR(levs(j,l)))) C + rC(INT(CEIL(levs(j,l)))) ) C / ( levs(j,l) - FLOOR(levs(j,l)) ) C for averaged levels. @@ -164,22 +211,43 @@ & 'Dimensional coordinate value at the lower point', & myThid) ENDIF +#ifdef DIAGNOSTICS_MISSING_VALUE + CALL MNC_CW_VATTR_MISSING(dn(1), 0, + I misval_r8, misval_r4, misval_int, + I myThid ) +#endif /* DIAGNOSTICS_MISSING_VALUE */ CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid) CALL MNC_CW_DEL_VNAME(dn(1), myThid) CALL MNC_CW_DEL_GNAME(dn(1), myThid) ENDDO +#endif /* DIAG_MNC_COORD_NEEDSWORK */ ENDIF #endif /* ALLOW_MNC */ +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + DO md = 1,nfields(listId) ndId = jdiag(md,listId) - parms1 = gdiag(ndId)(1:8) - IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN + gcode = gdiag(ndId)(1:10) + mate = 0 + mVec = 0 + IF ( gcode(5:5).EQ.'C' ) THEN +C- Check for Mate of a Counter Diagnostic + mate = hdiag(ndId) + 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 : + DO lm=1,averageCycle(listId) - ip = ABS(idiag(md,listId)) + ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1) im = mdiag(md,listId) + IF (mate.GT.0) im = im + kdiag(mate)*(lm-1) + IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1) + IF ( ndiag(ip,1,1).EQ.0 ) THEN C- Empty diagnostics case : @@ -188,14 +256,20 @@ & '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter CALL PRINT_MESSAGE( msgBuf, errorMessageUnit, & SQUEEZE_RIGHT, myThid) - WRITE(msgBuf,'(A,I4,3A,I3,2A)') + 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) - WRITE(msgBuf,'(A,I2,A)') - & '- WARNING - has not been filled (ndiag=', - & ndiag(ip,1,1), ' )' + 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)') @@ -218,34 +292,22 @@ ELSE C- diagnostics is not empty : - IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)') + IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN + WRITE(ioUnit,'(A,I6,3A,I8,2A)') & ' Computing Diagnostic # ', ndId, ' ', cdiag(ndId), & ' Counter:',ndiag(ip,1,1),' Parms: ',gdiag(ndId) - - IF ( parms1(5:5).EQ.'C' ) THEN -C Check for Mate of a Counter Diagnostic -C -------------------------------------- - mate_index = parms1(6:8) - READ (mate_index,'(I3)') mate - IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,2A)') + IF ( mate.GT.0 ) THEN + WRITE(ioUnit,'(3A,I6,2A)') & ' use Counter Mate for ', cdiag(ndId), & ' Diagnostic # ',mate, ' ', cdiag(mate) - - ELSE - mate = 0 - -C Check for Mate of a Vector Diagnostic -C ------------------------------------- - IF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN - mate_index = parms1(6:8) - READ (mate_index,'(I3)') mVec + ELSEIF ( mVec.GT.0 ) THEN IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN - IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)') + WRITE(ioUnit,'(3A,I6,3A)') & ' Vector Mate for ', cdiag(ndId), & ' Diagnostic # ',mVec, ' ', cdiag(mVec), & ' exists ' ELSE - IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)') + WRITE(ioUnit,'(3A,I6,3A)') & ' Vector Mate for ', cdiag(ndId), & ' Diagnostic # ',mVec, ' ', cdiag(mVec), & ' not enabled' @@ -253,38 +315,71 @@ ENDIF ENDIF - 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) + IF ( fflags(listId)(2:2).EQ.'P' ) THEN +C- get all the levels (for vertical interpolation) + DO bj = myByLo(myThid), myByHi(myThid) + DO bi = myBxLo(myThid), myBxHi(myThid) + DO k = 1,kdiag(ndId) + tmpLev = k + CALL GETDIAG( + I tmpLev,undef, + O qtmp1(1-OLx,1-OLy,k,bi,bj), + I ndId,mate,ip,im,bi,bj,myThid) + ENDDO + ENDDO ENDDO - ENDDO - ENDDO + ELSE +C- get only selected levels: + 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 + ENDDO + ENDDO + ENDIF C- end of empty diag / not empty block ENDIF -#ifdef ALLOW_MDSIO -C Prepare for mdsio optionality - IF (diag_mdsio) THEN - IF (fflags(listId)(1:1) .EQ. ' ') THEN -C This is the old default behavior - CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL', - & Nr+Nrphys,nlevels(listId),qtmp1,md,myIter,myThid) - ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN -C Force it to be 32-bit precision - CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL', - & Nr+Nrphys,nlevels(listId),qtmp1,md,myIter,myThid) - ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN -C Force it to be 64-bit precision - CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL', - & Nr+Nrphys,nlevels(listId),qtmp1,md,myIter,myThid) - ENDIF +C----------------------------------------------------------------------- +C Check to see if we need to 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, + I undef, 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 + +C-- Ready to write field "md", element "lm" in averageCycle(listId) + +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 a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R + CALL WRITE_REC_LEV_RL( + I fn, prec, + I NrMax, 1, nlevels(listId), + I qtmp1, -nRec, myIter, myThid ) ENDIF -#endif /* ALLOW_MDSIO */ #ifdef ALLOW_MNC IF (useMNC .AND. diag_mnc) THEN @@ -310,7 +405,7 @@ dim(2) = sNy + 2*OLy ib(1) = OLx + 1 ib(2) = OLy + 1 - IF (gdiag(ndId)(2:2) .EQ. 'M') THEN + IF (gdiag(ndId)(2:2) .EQ. 'M') THEN dn(1)(1:2) = 'X' ie(1) = OLx + sNx dn(2)(1:2) = 'Y' @@ -331,7 +426,7 @@ dn(2)(1:3) = 'Yp1' ie(2) = OLy + sNy + 1 ENDIF - + C Z is special since it varies WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId) IF ( (gdiag(ndId)(10:10) .EQ. 'R') @@ -346,7 +441,7 @@ & .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId) ENDIF - dim(3) = Nr+Nrphys + dim(3) = NrMax ib(3) = 1 ie(3) = nlevels(listId) @@ -356,24 +451,76 @@ ib(4) = 1 ie(4) = 1 - CALL MNC_CW_ADD_GNAME(d_cw_name, 4, + CALL MNC_CW_ADD_GNAME(d_cw_name, 4, & dim, dn, ib, ie, myThid) - CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name, + CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name, & 4,5, myThid) CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description', & tdiag(ndId),myThid) CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units', & udiag(ndId),myThid) - IF ((fflags(listId)(1:1) .EQ. ' ') +#ifdef DIAGNOSTICS_MISSING_VALUE +C Handle missing value attribute (land points) + IF ( misvalFlt(listId) .NE. UNSET_RL ) THEN + misvalLoc = misvalFlt(listId) + ELSE + misvalLoc = undef + ENDIF +C Defaults to UNSET_I + misvalIntLoc = misvalInt(listId) + DO ii=1,2 +C misval_r4(ii) = UNSET_FLOAT4 +C misval_r8(ii) = UNSET_FLOAT8 + misval_r4(ii) = misvalLoc + misval_r8(ii) = misvalLoc + misval_int(ii) = UNSET_I + ENDDO +C Missing values only for scalar diagnostics at mass points (so far) + IF ( gdiag(ndId)(1:2) .EQ. 'SM' ) THEN +C assign missing values and set flag for adding the netCDF atttibute + CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 2, + I misval_r8, misval_r4, misval_int, + I myThid ) +C and now use the missing values for masking out the land points + DO bj = myByLo(myThid), myByHi(myThid) + DO bi = myBxLo(myThid), myBxHi(myThid) + DO k = 1,nlevels(listId) + klev = NINT(levs(k,listId)) + DO j = 1-OLy,sNy+OLy + DO i = 1-OLx,sNx+OLx + IF ( _hFacC(I,J,klev,bi,bj) .EQ. 0. ) + & qtmp1(i,j,k,bi,bj) = misvalLoc + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ELSE +C suppress the missing value attribute (iflag = 0) +C Note: We have to call the following subroutine for each mnc that has +C been created "on the fly" by mnc_cw_add_vname and will be deleted +C by mnc_cw_del_vname, because all of these variables use the same +C identifier so that mnc_cw_vfmv(indv) needs to be overwritten for +C each of these variables + CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 0, + I misval_r8, misval_r4, misval_int, + I myThid ) + ENDIF +#endif /* DIAGNOSTICS_MISSING_VALUE */ + + IF ( ( (writeBinaryPrec .EQ. precFloat32) + & .AND. (fflags(listId)(1:1) .NE. 'D') + & .AND. (fflags(listId)(1:1) .NE. 'R') ) & .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0, & cdiag(ndId), qtmp1, myThid) - ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN + ELSEIF ( (writeBinaryPrec .EQ. precFloat64) + & .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0, & cdiag(ndId), qtmp1, myThid) ENDIF - + CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid) CALL MNC_CW_DEL_GNAME(d_cw_name, myThid) @@ -382,10 +529,25 @@ ENDIF #endif /* ALLOW_MNC */ + ENDDO C-- end of Processing Fld # md ENDIF ENDDO +#ifdef ALLOW_MDSIO + IF (diag_mdsio) THEN +C- Note: temporary: since it's 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 = nfields(listId)*averageCycle(listId) + CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE., + & 0, 0, nlevels(listId), ' ', + & nfields(listId), flds(1,listId), 1, myTime, + & nRec, myIter, myThid) + ENDIF +#endif /* ALLOW_MDSIO */ + RETURN END