Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:39:06 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
3e5de6a370 Jean*0001 #include "DIAG_OPTIONS.h"
                0002 
858abfd78a Jean*0003 C--   File diagstats_fill.F:
                0004 C--    Contents:
                0005 C--    o DIAGSTATS_FILL
                0006 C--    o DIAGSTATS_TO_RL
                0007 
                0008 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0009 
3e5de6a370 Jean*0010 CBOP
                0011 C     !ROUTINE: DIAGSTATS_FILL
                0012 C     !INTERFACE:
8c4f953ef4 Jean*0013       SUBROUTINE DIAGSTATS_FILL(
858abfd78a Jean*0014      I               inpFldRL, fracFldRL,
                0015 #ifndef REAL4_IS_SLOW
                0016      I               inpFldRS, fracFldRS,
                0017 #endif
                0018      I               scaleFact, power, arrType, nLevFract,
662c087d92 Jean*0019      I               ndId, kInQSd, region2fill, kLev, nLevs,
c9169fbe09 Jean*0020      I               bibjFlg, biArg, bjArg, myThid )
3e5de6a370 Jean*0021 
                0022 C     !DESCRIPTION:
                0023 C***********************************************************************
8c4f953ef4 Jean*0024 C   compute statistics over 1 tile
                0025 C   and increment the diagnostics array
662c087d92 Jean*0026 C     using a scaling factor & square option (power=2),
8c4f953ef4 Jean*0027 C     and with the option to use a fraction-weight (assumed
                0028 C         to be the counter-mate of the current diagnostics)
3e5de6a370 Jean*0029 C***********************************************************************
                0030 C     !USES:
                0031       IMPLICIT NONE
                0032 
                0033 C     == Global variables ===
                0034 #include "EEPARAMS.h"
                0035 #include "SIZE.h"
                0036 #include "DIAGNOSTICS_SIZE.h"
                0037 #include "DIAGNOSTICS.h"
                0038 
                0039 C     !INPUT PARAMETERS:
                0040 C***********************************************************************
                0041 C  Arguments Description
                0042 C  ----------------------
858abfd78a Jean*0043 C     inpFldRL  :: Field to increment diagnostics array (arrType=0,1)
                0044 C     fracFldRL :: fraction used for weighted average diagnostics (arrType=0,2)
                0045 C     inpFldRS  :: Field to increment diagnostics array (arrType=2,3)
                0046 C     fracFldRS :: fraction used for weighted average diagnostics (arrType=1,3)
662c087d92 Jean*0047 C     scaleFact :: scaling factor
                0048 C     power     :: option to fill-in with the field square (power=2)
858abfd78a Jean*0049 C     arrType   :: select which array & fraction (RL/RS) to process:
                0050 C                  0: both RL ; 1: inpRL & fracRS ; 2: inpRS,fracRL ; 3: both RS
                0051 C     nLevFract :: number of levels of the fraction field, =0: do not use fraction
662c087d92 Jean*0052 C     ndId      :: Diagnostics Id Number (in available diag list) of diag to process
                0053 C     kInQSd    :: Pointer to the slot in qSdiag to fill
                0054 C   region2fill :: array, indicates whether to compute statistics over region
3e5de6a370 Jean*0055 C                   "j" (if region2fill(j)=1) or not (if region2fill(j)=0)
662c087d92 Jean*0056 C     kLev      :: Integer flag for vertical levels:
3e5de6a370 Jean*0057 C                  > 0 (any integer): WHICH single level to increment in qSdiag.
                0058 C                  0,-1 to increment "nLevs" levels in qSdiag,
8c4f953ef4 Jean*0059 C                  0 : fill-in in the same order as the input array
3e5de6a370 Jean*0060 C                  -1: fill-in in reverse order.
662c087d92 Jean*0061 C     nLevs     :: indicates Number of levels of the input field array
3e5de6a370 Jean*0062 C                  (whether to fill-in all the levels (kLev<1) or just one (kLev>0))
c9169fbe09 Jean*0063 C     bibjFlg   :: Integer flag to indicate instructions for bi bj loop
3e5de6a370 Jean*0064 C                  0 indicates that the bi-bj loop must be done here
                0065 C                  1 indicates that the bi-bj loop is done OUTSIDE
                0066 C                  2 indicates that the bi-bj loop is done OUTSIDE
                0067 C                     AND that we have been sent a local array (with overlap regions)
                0068 C                  3 indicates that the bi-bj loop is done OUTSIDE
                0069 C                     AND that we have been sent a local array
                0070 C                     AND that the array has no overlap region (interior only)
c9169fbe09 Jean*0071 C                  NOTE - bibjFlg can be NEGATIVE to indicate not to increment counter
                0072 C     biArg     :: X-direction tile number - used for bibjFlg=1-3
                0073 C     bjArg     :: Y-direction tile number - used for bibjFlg=1-3
662c087d92 Jean*0074 C     myThid    :: my thread Id number
3e5de6a370 Jean*0075 C***********************************************************************
                0076 C                  NOTE: User beware! If a local (1 tile only) array
c9169fbe09 Jean*0077 C                        is sent here, bibjFlg MUST NOT be set to 0
3e5de6a370 Jean*0078 C                        or there will be out of bounds problems!
                0079 C***********************************************************************
858abfd78a Jean*0080       _RL inpFldRL(*)
                0081       _RL fracFldRL(*)
                0082 #ifndef REAL4_IS_SLOW
                0083       _RS inpFldRS(*)
                0084       _RS fracFldRS(*)
                0085 #endif
8c4f953ef4 Jean*0086       _RL scaleFact
662c087d92 Jean*0087       INTEGER power
858abfd78a Jean*0088       INTEGER arrType
8c4f953ef4 Jean*0089       INTEGER nLevFract
3e5de6a370 Jean*0090       INTEGER ndId, kInQSd
                0091       INTEGER region2fill(0:nRegions)
c9169fbe09 Jean*0092       INTEGER kLev, nLevs, bibjFlg, biArg, bjArg
3e5de6a370 Jean*0093       INTEGER myThid
                0094 CEOP
                0095 
                0096 C     !LOCAL VARIABLES:
                0097 C ===============
662c087d92 Jean*0098 C     useFract  :: flag to increment (or not) with fraction-weigted inpFld
8c4f953ef4 Jean*0099       LOGICAL useFract
                0100       INTEGER sizF
3e5de6a370 Jean*0101       INTEGER sizI1,sizI2,sizJ1,sizJ2
                0102       INTEGER sizTx,sizTy
                0103       INTEGER iRun, jRun, k, bi, bj
                0104       INTEGER kFirst, kLast
                0105       INTEGER kd, kd0, ksgn, kStore
                0106       CHARACTER*8 parms1
                0107       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0108       INTEGER km, km0
858abfd78a Jean*0109 #ifndef REAL4_IS_SLOW
                0110       _RL tmpFldRL( sNx+1,sNy+1)
                0111       _RL tmpFracRL(sNx+1,sNy+1)
                0112 #endif
3e5de6a370 Jean*0113 
                0114 C If-sequence to see if we are a valid and an active diagnostic
                0115 c     IF ( ndId.NE.0 .AND. kInQSd.NE.0 ) THEN
                0116 
                0117 C-      select range for 1rst & 2nd indices to accumulate
8c4f953ef4 Jean*0118 C         depending on variable location on C-grid,
3e5de6a370 Jean*0119         parms1 = gdiag(ndId)(1:8)
                0120         IF ( parms1(2:2).EQ.'Z' ) THEN
                0121          iRun = sNx+1
                0122          jRun = sNy+1
                0123 c       ELSEIF ( parms1(2:2).EQ.'U' ) THEN
                0124 c        iRun = sNx+1
                0125 c        jRun = sNy
                0126 c       ELSEIF ( parms1(2:2).EQ.'V' ) THEN
                0127 c        iRun = sNx
                0128 c        jRun = sNy+1
                0129         ELSE
                0130          iRun = sNx
                0131          jRun = sNy
                0132         ENDIF
                0133 
                0134 C-      Dimension of the input array:
c9169fbe09 Jean*0135         IF (ABS(bibjFlg).EQ.3) THEN
3e5de6a370 Jean*0136           sizI1 = 1
                0137           sizI2 = sNx
                0138           sizJ1 = 1
                0139           sizJ2 = sNy
                0140           iRun = sNx
                0141           jRun = sNy
                0142         ELSE
                0143           sizI1 = 1-OLx
                0144           sizI2 = sNx+OLx
                0145           sizJ1 = 1-OLy
                0146           sizJ2 = sNy+OLy
                0147         ENDIF
c9169fbe09 Jean*0148         IF (ABS(bibjFlg).GE.2) THEN
3e5de6a370 Jean*0149          sizTx = 1
                0150          sizTy = 1
                0151         ELSE
                0152          sizTx = nSx
                0153          sizTy = nSy
                0154         ENDIF
                0155 C-      Which part of inpFld to add : k = 3rd index,
                0156 C         and do the loop >> do k=kFirst,kLast <<
                0157         IF (kLev.LE.0) THEN
                0158           kFirst = 1
                0159           kLast  = nLevs
                0160         ELSEIF ( nLevs.EQ.1 ) THEN
                0161           kFirst = 1
                0162           kLast  = 1
                0163         ELSEIF ( kLev.LE.nLevs ) THEN
                0164           kFirst = kLev
                0165           kLast  = kLev
                0166         ELSE
                0167           STOP 'ABNORMAL END in DIAGSTATS_FILL: kLev > nLevs > 0'
                0168         ENDIF
8c4f953ef4 Jean*0169 C-      Which part of qSdiag to update: kd = 3rd index,
3e5de6a370 Jean*0170 C         and do the loop >> do k=kFirst,kLast ; kd = kd0 + k*ksgn <<
8c4f953ef4 Jean*0171 C  1rst try this: for the mask: km = km0 + k*ksgn so that kd= km + kInQSd - 1
3e5de6a370 Jean*0172         IF ( kLev.EQ.-1 ) THEN
                0173           ksgn = -1
                0174           kd0 = kInQSd + nLevs
                0175           km0 = 1 + nLevs
                0176         ELSEIF ( kLev.EQ.0 ) THEN
                0177           ksgn = 1
                0178           kd0 = kInQSd - 1
                0179           km0 = 0
                0180         ELSE
                0181           ksgn = 0
                0182           kd0 = kInQSd + kLev - 1
                0183           km0 = kLev
                0184         ENDIF
8c4f953ef4 Jean*0185 C-      Set fraction-weight option :
                0186         useFract = nLevFract.GT.0
                0187         IF ( useFract ) THEN
                0188           sizF = nLevFract
                0189         ELSE
                0190           sizF = 1
                0191         ENDIF
3e5de6a370 Jean*0192 
                0193 C-      Check for consistency with Nb of levels reserved in storage array
                0194         kStore = kd0 + MAX(ksgn*kFirst,ksgn*kLast) - kInQSd + 1
                0195         IF ( kStore.GT.kdiag(ndId) ) THEN
                0196          _BEGIN_MASTER(myThid)
e129400813 Jean*0197           WRITE(msgBuf,'(2A,I4,A)') 'DIAGSTATS_FILL: ',
3e5de6a370 Jean*0198      &     'exceed Nb of levels(=',kdiag(ndId),' ) reserved '
                0199           CALL PRINT_ERROR( msgBuf , myThid )
e129400813 Jean*0200           WRITE(msgBuf,'(2A,I6,2A)') 'DIAGSTATS_FILL: ',
3e5de6a370 Jean*0201      &     'for Diagnostics #', ndId, ' : ', cdiag(ndId)
                0202           CALL PRINT_ERROR( msgBuf , myThid )
                0203           WRITE(msgBuf,'(2A,2I4,I3)') 'calling DIAGSTATS_FILL ',
                0204      I     'with kLev,nLevs,bibjFlg=', kLev,nLevs,bibjFlg
                0205           CALL PRINT_ERROR( msgBuf , myThid )
                0206           WRITE(msgBuf,'(2A,I6,A)') 'DIAGSTATS_FILL: ',
                0207      I     '==> trying to store up to ', kStore, ' levels'
                0208           CALL PRINT_ERROR( msgBuf , myThid )
                0209           STOP 'ABNORMAL END: S/R DIAGSTATS_FILL'
                0210          _END_MASTER(myThid)
                0211         ENDIF
                0212 
858abfd78a Jean*0213 #ifndef REAL4_IS_SLOW
                0214       IF ( arrType.EQ.0 .OR. ( arrType.EQ.1 .AND. .NOT.useFract ) ) THEN
                0215 #endif
c9169fbe09 Jean*0216         IF ( bibjFlg.EQ.0 ) THEN
3e5de6a370 Jean*0217          DO bj=myByLo(myThid), myByHi(myThid)
                0218           DO bi=myBxLo(myThid), myBxHi(myThid)
                0219            DO k = kFirst,kLast
                0220             kd = kd0 + ksgn*k
                0221             km = km0 + ksgn*k
                0222             CALL DIAGSTATS_LOCAL(
                0223      U                  qSdiag(0,0,kd,bi,bj),
858abfd78a Jean*0224      I                  inpFldRL, fracFldRL,
662c087d92 Jean*0225      I                  scaleFact, power, useFract, sizF,
3e5de6a370 Jean*0226      I                  sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
                0227      I                  iRun,jRun,k,bi,bj,
c9169fbe09 Jean*0228      I                  km, bi, bj, bibjFlg, region2fill,
3e5de6a370 Jean*0229      I                  ndId, gdiag(ndId), myThid )
                0230            ENDDO
                0231           ENDDO
                0232          ENDDO
                0233         ELSE
                0234           bi = MIN(biArg,sizTx)
                0235           bj = MIN(bjArg,sizTy)
                0236           DO k = kFirst,kLast
                0237             kd = kd0 + ksgn*k
                0238             km = km0 + ksgn*k
                0239             CALL DIAGSTATS_LOCAL(
                0240      U                  qSdiag(0,0,kd,biArg,bjArg),
858abfd78a Jean*0241      I                  inpFldRL, fracFldRL,
662c087d92 Jean*0242      I                  scaleFact, power, useFract, sizF,
3e5de6a370 Jean*0243      I                  sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
                0244      I                  iRun,jRun,k,bi,bj,
c9169fbe09 Jean*0245      I                  km, biArg, bjArg, bibjFlg, region2fill,
3e5de6a370 Jean*0246      I                  ndId, gdiag(ndId), myThid )
                0247           ENDDO
                0248         ENDIF
                0249 
858abfd78a Jean*0250 #ifndef REAL4_IS_SLOW
                0251       ELSE
c9169fbe09 Jean*0252         IF ( bibjFlg.EQ.0 ) THEN
858abfd78a Jean*0253          DO bj=myByLo(myThid), myByHi(myThid)
                0254           DO bi=myBxLo(myThid), myBxHi(myThid)
                0255            DO k = kFirst,kLast
                0256             kd = kd0 + ksgn*k
                0257             km = km0 + ksgn*k
                0258             CALL DIAGSTATS_TO_RL(
                0259      I                  inpFldRL, fracFldRL, inpFldRS, fracFldRS,
                0260      O                  tmpFldRL, tmpFracRL,
                0261      I                  arrType, useFract, sizF,
                0262      I                  sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
                0263      I                  iRun,jRun,k,bi,bj, myThid )
                0264             CALL DIAGSTATS_LOCAL(
                0265      U                  qSdiag(0,0,kd,bi,bj),
                0266      I                  tmpFldRL, tmpFracRL,
                0267      I                  scaleFact, power, useFract, 1,
                0268      I                  1, iRun, 1, jRun, 1, 1, 1,
                0269      I                  iRun, jRun, 1, 1, 1,
c9169fbe09 Jean*0270      I                  km, bi, bj, bibjFlg, region2fill,
858abfd78a Jean*0271      I                  ndId, gdiag(ndId), myThid )
                0272            ENDDO
                0273           ENDDO
                0274          ENDDO
                0275         ELSE
                0276           bi = MIN(biArg,sizTx)
                0277           bj = MIN(bjArg,sizTy)
                0278           DO k = kFirst,kLast
                0279             kd = kd0 + ksgn*k
                0280             km = km0 + ksgn*k
                0281             CALL DIAGSTATS_TO_RL(
                0282      I                  inpFldRL, fracFldRL, inpFldRS, fracFldRS,
                0283      O                  tmpFldRL, tmpFracRL,
                0284      I                  arrType, useFract, sizF,
                0285      I                  sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
                0286      I                  iRun,jRun,k,bi,bj, myThid )
                0287             CALL DIAGSTATS_LOCAL(
                0288      U                  qSdiag(0,0,kd,biArg,bjArg),
                0289      I                  tmpFldRL, tmpFracRL,
                0290      I                  scaleFact, power, useFract, 1,
                0291      I                  1, iRun, 1, jRun, 1, 1, 1,
                0292      I                  iRun, jRun, 1, 1, 1,
c9169fbe09 Jean*0293      I                  km, biArg, bjArg, bibjFlg, region2fill,
858abfd78a Jean*0294      I                  ndId, gdiag(ndId), myThid )
                0295           ENDDO
                0296         ENDIF
                0297       ENDIF
                0298 #endif /* ndef REAL4_IS_SLOW */
                0299 
3e5de6a370 Jean*0300 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0301 c     ELSE
                0302 
                0303 c     ENDIF
                0304 
8c4f953ef4 Jean*0305       RETURN
3e5de6a370 Jean*0306       END
858abfd78a Jean*0307 
                0308 #ifndef REAL4_IS_SLOW
                0309 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0310 
                0311 CBOP
                0312 C     !ROUTINE: DIAGSTATS_TO_RL
                0313 C     !INTERFACE:
                0314       SUBROUTINE DIAGSTATS_TO_RL(
                0315      I                  inpFldRL, inpFrcRL, inpFldRS, inpFrcRS,
                0316      O                  outFldRL, outFrcRL,
                0317      I                  arrType, useFract, sizF,
                0318      I                  sizI1,sizI2,sizJ1,sizJ2,sizK,sizTx,sizTy,
                0319      I                  iRun,jRun,kIn,biIn,bjIn,
                0320      I                  myThid )
                0321 
                0322 C     !DESCRIPTION:
                0323 C     Do a local copy with conversion to RL type array
                0324 
                0325 C     !USES:
                0326       IMPLICIT NONE
                0327 
                0328 #include "EEPARAMS.h"
                0329 #include "SIZE.h"
                0330 
                0331 C     !INPUT/OUTPUT PARAMETERS:
                0332 C     == Routine Arguments ==
                0333 C     inpFldRL    :: input field array    to convert (arrType=0,1)
                0334 C     inpFrcRL    :: input fraction array to convert (arrType=0,2)
                0335 C     inpFldRS    :: input field array    to convert (arrType=2,3)
                0336 C     inpFrcRS    :: input fraction array to convert (arrType=1,3)
                0337 C     outFldRL    :: output field array
                0338 C     outFrcRL    :: output fraction array
                0339 C     arrType     :: select which array & fraction (RL/RS) to process:
                0340 C                    0: both RL ; 1: fldRL & frcRS ; 2: fldRS,frcRL ; 3: both RS
                0341 C     useFract    :: if True, process fraction-weight
                0342 C     sizF        :: size of inpFrc array: 3rd  dimension
                0343 C     sizI1,sizI2 :: size of inpFld array: 1rst index range (min,max)
                0344 C     sizJ1,sizJ2 :: size of inpFld array: 2nd  index range (min,max)
                0345 C     sizK        :: size of inpFld array: 3rd  dimension
                0346 C     sizTx,sizTy :: size of inpFld array: tile dimensions
                0347 C     iRun,jRun   :: range of 1rst & 2nd index
                0348 C     kIn         :: level index of inpFld array to process
                0349 C     biIn,bjIn   :: tile indices of inpFld array to process
                0350 C     myThid      :: my Thread Id number
                0351       INTEGER sizI1,sizI2,sizJ1,sizJ2
                0352       INTEGER sizF,sizK,sizTx,sizTy
                0353       INTEGER iRun, jRun, kIn, biIn, bjIn
                0354       _RL     inpFldRL(sizI1:sizI2,sizJ1:sizJ2,sizK,sizTx,sizTy)
                0355       _RL     inpFrcRL(sizI1:sizI2,sizJ1:sizJ2,sizF,sizTx,sizTy)
                0356       _RS     inpFldRS(sizI1:sizI2,sizJ1:sizJ2,sizK,sizTx,sizTy)
                0357       _RS     inpFrcRS(sizI1:sizI2,sizJ1:sizJ2,sizF,sizTx,sizTy)
                0358       _RL     outFldRL(1:iRun,1:jRun)
                0359       _RL     outFrcRL(1:iRun,1:jRun)
                0360       INTEGER arrType
                0361       LOGICAL useFract
                0362       INTEGER myThid
                0363 CEOP
                0364 
                0365 C     !LOCAL VARIABLES:
                0366 C     i,j    :: loop indices
                0367       INTEGER i, j, kFr
                0368 
                0369       IF ( arrType.LE.1 ) THEN
                0370         DO j=1,jRun
                0371          DO i=1,iRun
                0372            outFldRL(i,j) = inpFldRL(i,j,kIn,biIn,bjIn)
                0373          ENDDO
                0374         ENDDO
                0375       ELSE
                0376         DO j=1,jRun
                0377          DO i=1,iRun
                0378            outFldRL(i,j) = inpFldRS(i,j,kIn,biIn,bjIn)
                0379          ENDDO
                0380         ENDDO
                0381       ENDIF
                0382 
                0383       IF ( useFract ) THEN
                0384        kFr = MIN(kIn,sizF)
                0385        IF ( arrType.EQ.0 .OR. arrType.EQ.2 ) THEN
                0386         DO j=1,jRun
                0387          DO i=1,iRun
                0388            outFrcRL(i,j) = inpFrcRL(i,j,kFr,biIn,bjIn)
                0389          ENDDO
                0390         ENDDO
                0391        ELSE
                0392         DO j=1,jRun
                0393          DO i=1,iRun
                0394            outFrcRL(i,j) = inpFrcRS(i,j,kFr,biIn,bjIn)
                0395          ENDDO
                0396         ENDDO
                0397        ENDIF
                0398       ENDIF
                0399 
                0400       RETURN
                0401       END
                0402 #endif /* ndef REAL4_IS_SLOW */