pkg/diagnostics/diagnostics_mnc_out.F

380c427652 Jean*0001 #include "DIAG_OPTIONS.h"
                0002 #undef DIAG_MNC_COORD_NEEDSWORK
                0003
                0004 C--   File diagnostics_mnc_out.F: Routines to write MNC diagnostics output
                0005 C--    Contents:
                0006 C--    o DIAGNOSTICS_MNC_SET
                0007 C--    o DIAGNOSTICS_MNC_OUT
                0008
                0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0010 CBOP 0
                0011 C     !ROUTINE: DIAGNOSTICS_MNC_SET
                0012
                0013 C     !INTERFACE:
                0014       SUBROUTINE DIAGNOSTICS_MNC_SET(
                0015      I     nLevOutp, listId, lm,
9473248f34 Jean*0016      O     diag_mnc_bn,
bd9081423e Jean*0017      I     misValLoc, myTime, myIter, myThid )
380c427652 Jean*0018
                0019 C     !DESCRIPTION:
                0020 C     Set MNC file for writing diagnostics fields.
                0021
                0022 C     !USES:
                0023       IMPLICIT NONE
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
                0028 #include "DIAGNOSTICS_SIZE.h"
                0029 #include "DIAGNOSTICS.h"
                0030
                0031
a22b7a769d Jean*0032 C     !INPUT/OUTPUT PARAMETERS:
                0033 C     nLevOutp        :: number of levels to write in output file
                0034 C     listId          :: Diagnostics list number being written
                0035 C     lm              :: loop index (averageCycle)
                0036 C     diag_mnc_bn     :: NetCDF output file name
                0037 C     misValLoc       :: local Missing Value
                0038 C     myTime          :: current time of simulation (s)
                0039 C     myIter          :: current iteration number
                0040 C     myThid          :: my Thread Id number
380c427652 Jean*0041       INTEGER nLevOutp
                0042       INTEGER listId, lm
                0043       CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
a22b7a769d Jean*0044       _RL     misValLoc
380c427652 Jean*0045       _RL     myTime
                0046       INTEGER myIter, myThid
                0047 CEOP
                0048
                0049 #ifdef ALLOW_MNC
                0050 C     !FUNCTIONS:
                0051       INTEGER ILNBLNK
                0052       EXTERNAL ILNBLNK
                0053
                0054 C     !LOCAL VARIABLES:
                0055       _RL tmpLev
                0056       INTEGER iLen
                0057
                0058 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
                0059       INTEGER ii, klev
                0060       INTEGER CW_DIMS, NLEN
                0061       PARAMETER ( CW_DIMS = 10 )
                0062       PARAMETER ( NLEN    = 80 )
                0063       INTEGER dim(CW_DIMS), ib(CW_DIMS), ie(CW_DIMS)
                0064       CHARACTER*(NLEN) dn(CW_DIMS)
                0065 c     CHARACTER*(NLEN) d_cw_name
                0066 c     CHARACTER*(NLEN) dn_blnk
                0067 #ifdef DIAG_MNC_COORD_NEEDSWORK
a22b7a769d Jean*0068       INTEGER NrMax
                0069       PARAMETER( NrMax = numLevels )
380c427652 Jean*0070       INTEGER i, j
                0071       CHARACTER*(5) ctmp
                0072       _RS ztmp(NrMax)
                0073 #endif
                0074       REAL*8  misval_r8(2)
                0075       REAL*4  misval_r4(2)
                0076       INTEGER misval_int(2)
                0077
                0078 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0079 c     IF (useMNC .AND. diag_mnc) THEN
                0080
                0081 C     Handle missing value attribute (land points)
                0082 C     Defaults to UNSET_I
                0083         DO ii=1,2
                0084           misval_r4(ii)  = misValLoc
                0085           misval_r8(ii)  = misValLoc
                0086           misval_int(ii) = UNSET_I
                0087         ENDDO
                0088 c       DO i = 1,MAX_LEN_FNAM
                0089 c         diag_mnc_bn(i:i) = ' '
                0090 c       ENDDO
                0091 c       DO i = 1,NLEN
                0092 c         dn_blnk(i:i) = ' '
                0093 c       ENDDO
                0094         iLen = ILNBLNK(fnames(listId))
                0095         WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:iLen)
                0096
                0097 C       Update the record dimension by writing the iteration number
                0098         klev = myIter + lm - averageCycle(listId)
                0099         tmpLev = myTime + deltaTClock*( lm - averageCycle(listId) )
                0100         CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
                0101         CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',tmpLev,myThid)
                0102         CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
                0103         CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',klev,myThid)
                0104
                0105 C       NOTE: at some point it would be a good idea to add a time_bounds
                0106 C       variable that has dimension (2,T) and clearly denotes the
                0107 C       beginning and ending times for each diagnostics period
                0108
                0109 c       dn(1)(1:NLEN) = dn_blnk(1:NLEN)
                0110         WRITE(dn(1),'(a,i6.6)') 'Zmd', nLevOutp
                0111         dim(1) = nLevOutp
                0112         ib(1)  = 1
                0113         ie(1)  = nLevOutp
                0114
                0115         CALL MNC_CW_ADD_GNAME('diag_levels', 1,
                0116      &       dim, dn, ib, ie, myThid)
                0117         CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
                0118      &       0,0, myThid)
                0119         CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
                0120      &       'Idicies of vertical levels within the source arrays',
                0121      &       myThid)
                0122 C     suppress the missing value attribute (iflag = 0)
a22b7a769d Jean*0123         CALL MNC_CW_VATTR_MISSING('diag_levels', 0,
380c427652 Jean*0124      I       misval_r8, misval_r4, misval_int, myThid )
                0125
                0126         CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
                0127      &       'diag_levels', levs(1,listId), myThid)
                0128
                0129         CALL MNC_CW_DEL_VNAME('diag_levels', myThid)
                0130         CALL MNC_CW_DEL_GNAME('diag_levels', myThid)
                0131
                0132 #ifdef DIAG_MNC_COORD_NEEDSWORK
                0133 C       This part has been placed in an #ifdef because, as its currently
                0134 C       written, it will only work with variables defined on a dynamics
                0135 C       grid.  As we start using diagnostics for physics grids, ice
                0136 C       levels, land levels, etc. the different vertical coordinate
                0137 C       dimensions will have to be taken into account.
                0138
                0139 C       20051021 JMC & EH3 : We need to extend this so that a few
                0140 C       variables each defined on different grids do not have the same
                0141 C       vertical dimension names so we should be using a pattern such
                0142 C       as: Z[uml]td000000 where the 't' is the type as specified by
                0143 C       gdiag(10)
                0144
                0145 C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
                0146         ctmp(1:5) = 'mul  '
                0147         DO i = 1,3
                0148 c         dn(1)(1:NLEN) = dn_blnk(1:NLEN)
                0149           WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)
                0150           CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)
                0151           CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)
                0152
                0153 C         The following three ztmp() loops should eventually be modified
                0154 C         to reflect the fractional nature of levs(j,l) -- they should
                0155 C         do something like:
                0156 C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))
                0157 C                      + ( rC(INT(FLOOR(levs(j,l))))
                0158 C                          + rC(INT(CEIL(levs(j,l)))) )
                0159 C                        / ( levs(j,l) - FLOOR(levs(j,l)) )
                0160 C         for averaged levels.
                0161           IF (i .EQ. 1) THEN
                0162             DO j = 1,nlevels(listId)
                0163               ztmp(j) = rC(NINT(levs(j,listId)))
                0164             ENDDO
                0165             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
                0166      &           'Dimensional coordinate value at the mid point',
                0167      &           myThid)
                0168           ELSEIF (i .EQ. 2) THEN
                0169             DO j = 1,nlevels(listId)
                0170               ztmp(j) = rF(NINT(levs(j,listId)) + 1)
                0171             ENDDO
                0172             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
                0173      &           'Dimensional coordinate value at the upper point',
                0174      &           myThid)
                0175           ELSEIF (i .EQ. 3) THEN
                0176             DO j = 1,nlevels(listId)
                0177               ztmp(j) = rF(NINT(levs(j,listId)))
                0178             ENDDO
                0179             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
                0180      &           'Dimensional coordinate value at the lower point',
                0181      &           myThid)
                0182           ENDIF
                0183 C     suppress the missing value attribute (iflag = 0)
9473248f34 Jean*0184           IF (useMissingValue)
380c427652 Jean*0185      &          CALL MNC_CW_VATTR_MISSING(dn(1), 0,
                0186      I          misval_r8, misval_r4, misval_int, myThid )
                0187           CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)
                0188           CALL MNC_CW_DEL_VNAME(dn(1), myThid)
                0189           CALL MNC_CW_DEL_GNAME(dn(1), myThid)
                0190         ENDDO
                0191 #endif /*  DIAG_MNC_COORD_NEEDSWORK  */
                0192
                0193 c     ENDIF
                0194 #endif /*  ALLOW_MNC  */
                0195
                0196       RETURN
                0197       END
                0198
                0199 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0200 CBOP 0
                0201 C     !ROUTINE: DIAGNOSTICS_MNC_OUT
                0202
                0203 C     !INTERFACE:
                0204       SUBROUTINE DIAGNOSTICS_MNC_OUT(
a22b7a769d Jean*0205      I     NrMax, nLevOutp, listId, ndId, mate,
9473248f34 Jean*0206      I     diag_mnc_bn, qtmp,
                0207      I     misValLoc, myTime, myIter, myThid )
380c427652 Jean*0208
                0209 C     !DESCRIPTION:
                0210 C     write diagnostics fields to MNC file.
                0211
                0212 C     !USES:
                0213       IMPLICIT NONE
                0214 #include "SIZE.h"
                0215 #include "EEPARAMS.h"
                0216 #include "PARAMS.h"
                0217 #include "GRID.h"
                0218 #include "DIAGNOSTICS_SIZE.h"
                0219 #include "DIAGNOSTICS.h"
                0220
                0221 C     !INPUT PARAMETERS:
a22b7a769d Jean*0222 C     NrMax           :: 3rd dimension of output-field array to write
                0223 C     nLevOutp        :: number of levels to write in output file
                0224 C     listId          :: Diagnostics list number being written
                0225 C     ndId            :: diagnostics Id number (in available diagnostics list)
                0226 C     mate            :: counter diagnostic number if any ; 0 otherwise
                0227 C     diag_mnc_bn     :: NetCDF output file name
                0228 C     qtmp            :: output-field array to write
9473248f34 Jean*0229 C     misValLoc       :: local Missing Value
a22b7a769d Jean*0230 C     myTime          :: current time of simulation (s)
                0231 C     myIter          :: current iteration number
                0232 C     myThid          :: my Thread Id number
380c427652 Jean*0233       INTEGER NrMax
                0234       INTEGER nLevOutp
a22b7a769d Jean*0235       INTEGER listId, ndId, mate
380c427652 Jean*0236       CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
a22b7a769d Jean*0237       _RL     qtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
9473248f34 Jean*0238       _RL     misValLoc
380c427652 Jean*0239       _RL     myTime
                0240       INTEGER myIter, myThid
                0241 CEOP
                0242
a22b7a769d Jean*0243 #ifdef ALLOW_MNC
380c427652 Jean*0244 C     !FUNCTIONS:
                0245 c     INTEGER ILNBLNK
                0246 c     EXTERNAL ILNBLNK
                0247
                0248 C     !LOCAL VARIABLES:
                0249 C     i,j,k :: loop indices
                0250 C     bi,bj :: tile indices
                0251       INTEGER i, j, k
                0252       INTEGER bi, bj
                0253
                0254 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
                0255 c     INTEGER ll, llMx, jj, jjMx
                0256       INTEGER ii, klev
                0257       INTEGER CW_DIMS, NLEN
                0258       PARAMETER ( CW_DIMS = 10 )
                0259       PARAMETER ( NLEN    = 80 )
                0260       INTEGER dim(CW_DIMS), ib(CW_DIMS), ie(CW_DIMS)
                0261       CHARACTER*(NLEN) dn(CW_DIMS)
                0262       CHARACTER*(NLEN) d_cw_name
                0263 c     CHARACTER*(NLEN) dn_blnk
                0264       LOGICAL useMisValForThisDiag
                0265       REAL*8  misval_r8(2)
                0266       REAL*4  misval_r4(2)
                0267       INTEGER misval_int(2)
                0268
                0269 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0270 c     IF (useMNC .AND. diag_mnc) THEN
                0271
                0272         _BEGIN_MASTER( myThid )
                0273
                0274          DO ii = 1,CW_DIMS
                0275 c          d_cw_name(1:NLEN) = dn_blnk(1:NLEN)
                0276 c          dn(ii)(1:NLEN) = dn_blnk(1:NLEN)
                0277            dn(ii) =  ' '
                0278          ENDDO
                0279          DO ii=1,2
                0280            misval_r4(ii)  = misValLoc
                0281            misval_r8(ii)  = misValLoc
                0282            misval_int(ii) = UNSET_I
                0283          ENDDO
                0284
                0285 C        Note that the "d_cw_name" variable is a hack that hides a
                0286 C        subtlety within MNC.  Basically, each MNC-wrapped file is
                0287 C        caching its own concept of what each "grid name" (that is, a
                0288 C        dimension group name) means.  So one cannot re-use the same
                0289 C        "grid" name for different collections of dimensions within a
                0290 C        given file.  By appending the "ndId" values to each name, we
                0291 C        guarantee uniqueness within each MNC-produced file.
                0292          WRITE(d_cw_name,'(a,i6.6)') 'd_cw_',ndId
                0293
                0294 C        XY dimensions
                0295          dim(1)       = sNx + 2*OLx
                0296          dim(2)       = sNy + 2*OLy
                0297          ib(1)        = OLx + 1
                0298          ib(2)        = OLy + 1
                0299          IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
                0300            dn(1)(1:2) = 'X'
                0301            ie(1)      = OLx + sNx
                0302            dn(2)(1:2) = 'Y'
                0303            ie(2)      = OLy + sNy
                0304          ELSEIF (gdiag(ndId)(2:2) .EQ. 'U') THEN
                0305            dn(1)(1:3) = 'Xp1'
                0306            ie(1)      = OLx + sNx + 1
                0307            dn(2)(1:2) = 'Y'
                0308            ie(2)      = OLy + sNy
                0309          ELSEIF (gdiag(ndId)(2:2) .EQ. 'V') THEN
                0310            dn(1)(1:2) = 'X'
                0311            ie(1)      = OLx + sNx
                0312            dn(2)(1:3) = 'Yp1'
                0313            ie(2)      = OLy + sNy + 1
                0314          ELSEIF (gdiag(ndId)(2:2) .EQ. 'Z') THEN
                0315            dn(1)(1:3) = 'Xp1'
                0316            ie(1)      = OLx + sNx + 1
                0317            dn(2)(1:3) = 'Yp1'
                0318            ie(2)      = OLy + sNy + 1
                0319          ENDIF
                0320
                0321 C        Z is special since it varies
                0322          WRITE(dn(3),'(a,i6.6)') 'Zd', nLevOutp
                0323          IF ( (gdiag(ndId)(10:10) .EQ. 'R')
                0324      &        .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
                0325            WRITE(dn(3),'(a,i6.6)') 'Zmd', nLevOutp
                0326          ENDIF
                0327          IF ( (gdiag(ndId)(10:10) .EQ. 'R')
                0328      &        .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
                0329            WRITE(dn(3),'(a,i6.6)') 'Zld', nLevOutp
                0330          ENDIF
                0331          IF ( (gdiag(ndId)(10:10) .EQ. 'R')
                0332      &        .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
                0333            WRITE(dn(3),'(a,i6.6)') 'Zud', nLevOutp
                0334          ENDIF
                0335          dim(3) = NrMax
                0336          ib(3)  = 1
                0337          ie(3)  = nLevOutp
                0338
                0339 C        Time dimension
                0340          dn(4)(1:1) = 'T'
                0341          dim(4) = -1
                0342          ib(4)  = 1
                0343          ie(4)  = 1
                0344
                0345          CALL MNC_CW_ADD_GNAME( d_cw_name, 4,
                0346      &                          dim, dn, ib, ie, myThid )
                0347          CALL MNC_CW_ADD_VNAME( cdiag(ndId), d_cw_name,
                0348      &                          4, 5, myThid )
                0349          CALL MNC_CW_ADD_VATTR_TEXT( cdiag(ndId),'description',
                0350      &                               tdiag(ndId), myThid )
                0351          CALL MNC_CW_ADD_VATTR_TEXT( cdiag(ndId),'units',
                0352      &                               udiag(ndId), myThid )
                0353
a22b7a769d Jean*0354          useMisValForThisDiag = mate.GT.0
                0355 C     Use the missing values for masking out the land points:
                0356 C     only for scalar diagnostics at mass points (so far)
9473248f34 Jean*0357          IF ( useMissingValue.AND.gdiag(ndId)(1:2).EQ.'SM' ) THEN
a22b7a769d Jean*0358            useMisValForThisDiag = .TRUE.
380c427652 Jean*0359 C     note: better to use 2-D mask if kdiag <> Nr or vert.integral
                0360            DO bj = myByLo(myThid), myByHi(myThid)
                0361             DO bi = myBxLo(myThid), myBxHi(myThid)
                0362              DO k = 1,nLevOutp
                0363               klev = NINT(levs(k,listId))
                0364               IF ( fflags(listId)(2:2).EQ.'I' ) kLev = 1
                0365               DO j = 1-OLy,sNy+OLy
                0366                DO i = 1-OLx,sNx+OLx
                0367                 IF ( maskC(i,j,klev,bi,bj) .EQ. 0. )
                0368      &                qtmp(i,j,k,bi,bj) = misValLoc
                0369                ENDDO
                0370               ENDDO
                0371              ENDDO
                0372             ENDDO
                0373            ENDDO
a22b7a769d Jean*0374          ENDIF
                0375          IF ( useMisValForThisDiag ) THEN
                0376 C     assign missing values and set flag for adding the netCDF atttibute
                0377            CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 2,
                0378      I                misval_r8, misval_r4, misval_int, myThid )
380c427652 Jean*0379          ELSE
                0380 C     suppress the missing value attribute (iflag = 0)
                0381 C     Note: We have to call the following subroutine for each mnc that has
                0382 C     been created "on the fly" by mnc_cw_add_vname and will be deleted
                0383 C     by mnc_cw_del_vname, because all of these variables use the same
                0384 C     identifier so that mnc_cw_vfmv(indv) needs to be overwritten for
                0385 C     each of these variables
                0386            CALL MNC_CW_VATTR_MISSING( cdiag(ndId), 0,
                0387      I                 misval_r8, misval_r4, misval_int, myThid )
                0388          ENDIF
                0389
                0390          IF (  ((writeBinaryPrec .EQ. precFloat32).AND.
                0391      &          (fflags(listId)(1:1) .NE. 'D'))
                0392      &         .OR. (fflags(listId)(1:1) .EQ. 'R') ) THEN
                0393            CALL MNC_CW_RL_W( 'R',diag_mnc_bn,0,0,
                0394      &                       cdiag(ndId), qtmp, myThid)
                0395          ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
                0396      &         .OR.  (fflags(listId)(1:1) .EQ. 'D') ) THEN
                0397            CALL MNC_CW_RL_W( 'D',diag_mnc_bn,0,0,
                0398      &                       cdiag(ndId), qtmp, myThid)
                0399          ENDIF
                0400
                0401          CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
                0402          CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
                0403
                0404         _END_MASTER( myThid )
                0405
                0406 c     ENDIF
                0407 #endif /*  ALLOW_MNC  */
                0408
                0409       RETURN
                0410       END