Back to home page

MITgcm

 
 

    


File indexing completed on 2022-03-25 05:09:35 UTC

view on githubraw file Latest commit 7e00d7e8 on 2022-03-24 16:33:13 UTC
051aff16bf Mart*0001 #include "CPP_OPTIONS.h"
                0002 
375581a429 Jean*0003 C--  File remove_mean.F:
                0004 C--   Contents
                0005 C--   o REMOVE_MEAN_RL
                0006 C--   o REMOVE_MEAN_RS
                0007 
                0008 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
051aff16bf Mart*0009 CBOP
                0010 C     !ROUTINE: REMOVE_MEAN_RL
                0011 C     !INTERFACE:
                0012       SUBROUTINE REMOVE_MEAN_RL(
375581a429 Jean*0013      I                myNr,
                0014      U                arrFld,
                0015      I                arrhFac, arrMask, arrArea, arrDr,
                0016      I                arrName, myTime, myIter, myThid )
051aff16bf Mart*0017 C     !DESCRIPTION: \bv
7e00d7e8f9 Jean*0018 C     Correct for global-mean imbalance of global "_RL" array "arrFld"
                0019 C      (i.e., output arrFld has zero mean):
                0020 C     Apply either (if arrDr(1) > 0) uniform correction
                0021 C      or (if arrDr(1) < 0) a correction scaled by "arrhFac" weight.
                0022 C     The global mean imbalance is computed as area (arrArea), mask
                0023 C      (arrMask) and thickness (arrDr) weighted, as well as thickness
                0024 C      factor (arrhFac) weighted (1rst case) or not (2nd case).
                0025 C     Comment: currently this S/R is only applied to 2-D field.
051aff16bf Mart*0026 C     \ev
                0027 
375581a429 Jean*0028 C     !USES:
051aff16bf Mart*0029       IMPLICIT NONE
                0030 C     === Global data ===
                0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
7e00d7e8f9 Jean*0034 #include "GRID.h"
051aff16bf Mart*0035 
375581a429 Jean*0036 C     !INPUT/OUTPUT PARAMETERS:
7e00d7e8f9 Jean*0037 C     myNr      :: third dimension of field to correct
                0038 C     arrFld    :: field array to correct for imbalance
                0039 C     arrhFac   :: same size as arrFld ; thickness factor (1rst case)
                0040 C               :: or correction scaling factor (2nd case)
                0041 C     arrMask   :: 2-D mask for global mean calculation
                0042 C     arrArea   :: 2-D grid-cell area for global mean calculation
                0043 C     arrDr     :: myNr long vector (level thickness)
                0044 C     arrName   :: name of field to correct
                0045 C     myTime    :: Current time in simulation
                0046 C     myIter    :: Current iteration number in simulation
                0047 C     myThid    :: my Thread ID number
051aff16bf Mart*0048       INTEGER myNr
375581a429 Jean*0049       _RL arrFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
051aff16bf Mart*0050       _RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
375581a429 Jean*0051       _RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
051aff16bf Mart*0052       _RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053       _RS arrDr(myNr)
                0054       CHARACTER*(*) arrName
                0055       _RL myTime
375581a429 Jean*0056       INTEGER myIter
051aff16bf Mart*0057       INTEGER myThid
                0058 
375581a429 Jean*0059 C     !LOCAL VARIABLES:
                0060       INTEGER bi, bj, i, j, k
                0061       _RL volTile(nSx,nSy), sumTile(nSx,nSy)
7e00d7e8f9 Jean*0062       _RL tmpVol, volGlob, sumGlob, theMean
                0063       _RL meanCorr
375581a429 Jean*0064       CHARACTER*(MAX_LEN_MBUF) msgBuf
051aff16bf Mart*0065 CEOP
                0066 
                0067       DO bj=myByLo(myThid),myByHi(myThid)
                0068        DO bi=myBxLo(myThid),myBxHi(myThid)
375581a429 Jean*0069         volTile(bi,bj) = 0. _d 0
                0070         sumTile(bi,bj) = 0. _d 0
7e00d7e8f9 Jean*0071         IF ( arrDr(1).GE.zeroRS ) THEN
                0072 C--   do (arrhFac) weighted average to apply uniform correction
                0073          DO k=1,myNr
                0074           DO j=1,sNy
                0075            DO i=1,sNx
                0076             tmpVol = arrArea(i,j,bi,bj)*arrMask(i,j,bi,bj)
                0077      &             * arrhFac(i,j,k,bi,bj)*arrDr(k)
                0078             volTile(bi,bj) = volTile(bi,bj) + tmpVol
                0079             sumTile(bi,bj) = sumTile(bi,bj) + tmpVol*arrFld(i,j,k,bi,bj)
                0080            ENDDO
051aff16bf Mart*0081           ENDDO
                0082          ENDDO
7e00d7e8f9 Jean*0083         ELSE
                0084 C-    do simple average to apply (arrhFac) weighted correction
                0085          DO k=1,myNr
                0086           DO j=1,sNy
                0087            DO i=1,sNx
                0088             tmpVol = arrArea(i,j,bi,bj)*arrMask(i,j,bi,bj)
                0089      &             * ABS(arrDr(k))
                0090             volTile(bi,bj) = volTile(bi,bj)
                0091      &                     + tmpVol*arrhFac(i,j,k,bi,bj)
                0092             sumTile(bi,bj) = sumTile(bi,bj)
                0093      &                     + tmpVol*arrFld(i,j,k,bi,bj)
                0094            ENDDO
                0095           ENDDO
                0096          ENDDO
                0097         ENDIF
051aff16bf Mart*0098        ENDDO
                0099       ENDDO
                0100 
375581a429 Jean*0101       CALL GLOBAL_SUM_TILE_RL( volTile, volGlob, myThid )
                0102       CALL GLOBAL_SUM_TILE_RL( sumTile, sumGlob, myThid )
051aff16bf Mart*0103 
375581a429 Jean*0104       IF ( volGlob.GT.zeroRL ) THEN
7e00d7e8f9 Jean*0105        meanCorr = sumGlob/volGlob
051aff16bf Mart*0106        DO bj=myByLo(myThid),myByHi(myThid)
                0107         DO bi=myBxLo(myThid),myBxHi(myThid)
7e00d7e8f9 Jean*0108          IF ( arrDr(1).GE.zeroRS ) THEN
                0109 C--   apply uniform correction
                0110           DO k=1,myNr
                0111            DO j=1-OLy,sNy+OLy
                0112             DO i=1-OLx,sNx+OLx
                0113              IF ( arrhFac(i,j,k,bi,bj).NE.zeroRS ) THEN
                0114               arrFld(i,j,k,bi,bj) = arrFld(i,j,k,bi,bj) - meanCorr
                0115              ENDIF
                0116             ENDDO
051aff16bf Mart*0117            ENDDO
                0118           ENDDO
7e00d7e8f9 Jean*0119          ELSE
                0120 C-    apply (arrhFac) weighted correction
                0121           DO k=1,myNr
                0122            DO j=1-OLy,sNy+OLy
                0123             DO i=1-OLx,sNx+OLx
                0124              IF ( arrMask(i,j,bi,bj).NE.zeroRS ) THEN
                0125               arrFld(i,j,k,bi,bj) = arrFld(i,j,k,bi,bj)
                0126      &                            - meanCorr*arrhFac(i,j,k,bi,bj)
                0127              ENDIF
                0128             ENDDO
                0129            ENDDO
                0130           ENDDO
                0131 C-    end type correction selection
                0132          ENDIF
051aff16bf Mart*0133         ENDDO
                0134        ENDDO
375581a429 Jean*0135       ELSE
7e00d7e8f9 Jean*0136        meanCorr = 0. _d 0
051aff16bf Mart*0137       ENDIF
                0138 
375581a429 Jean*0139 C     Print the global mean to standard output, this is a measure
                0140 C     of the magnitude of the correction to array arrFld
051aff16bf Mart*0141       IF ( balancePrintMean ) THEN
                0142        _BEGIN_MASTER( myThid )
7e00d7e8f9 Jean*0143        IF ( arrDr(1).GE.zeroRS ) THEN
                0144         theMean = meanCorr
                0145         WRITE(msgBuf,'(3A,1PE21.14,A,I10)')
                0146      &        'REMOVE_MEAN_RL: Global mean of ', arrName, ' = ',
                0147      &        theMean, '  @ it=', myIter
                0148         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0149      &                      SQUEEZE_RIGHT, myThid )
                0150        ELSE
                0151         theMean = 0. _d 0
                0152         IF ( globalArea.GT.zeroRL ) theMean = sumGlob/globalArea
                0153         WRITE(msgBuf,'(3A,2(1PE21.14,A),I10)')
                0154      &        'REMOVE_MEAN_RL: ', arrName, ' Global mean= ',
                0155      &        theMean, ', remove: ', meanCorr, ' @ it=', myIter
                0156         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0157      &                      SQUEEZE_RIGHT, myThid )
                0158        ENDIF
051aff16bf Mart*0159        _END_MASTER( myThid )
                0160       ENDIF
                0161 
                0162       RETURN
                0163       END
                0164 
375581a429 Jean*0165 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
051aff16bf Mart*0166 CBOP
                0167 C     !ROUTINE: REMOVE_MEAN_RS
                0168 C     !INTERFACE:
                0169       SUBROUTINE REMOVE_MEAN_RS(
375581a429 Jean*0170      I                myNr,
                0171      U                arrFld,
                0172      I                arrhFac, arrMask, arrArea, arrDr,
                0173      I                arrName, myTime, myIter, myThid )
051aff16bf Mart*0174 C     !DESCRIPTION: \bv
7e00d7e8f9 Jean*0175 C     Correct for global-mean imbalance of global "_RS" array "arrFld"
                0176 C      (i.e., output arrFld has zero mean):
                0177 C     Apply either (if arrDr(1) > 0) uniform correction
                0178 C      or (if arrDr(1) < 0) a correction scaled by "arrhFac" weight.
                0179 C     The global mean imbalance is computed as area (arrArea), mask
                0180 C      (arrMask) and thickness (arrDr) weighted, as well as thickness
                0181 C      factor (arrhFac) weighted (1rst case) or not (2nd case).
                0182 C     Comment: currently this S/R is only applied to 2-D field.
051aff16bf Mart*0183 C     \ev
                0184 
375581a429 Jean*0185 C     !USES:
051aff16bf Mart*0186       IMPLICIT NONE
                0187 C     === Global data ===
                0188 #include "SIZE.h"
                0189 #include "EEPARAMS.h"
                0190 #include "PARAMS.h"
7e00d7e8f9 Jean*0191 #include "GRID.h"
051aff16bf Mart*0192 
375581a429 Jean*0193 C     !INPUT/OUTPUT PARAMETERS:
7e00d7e8f9 Jean*0194 C     myNr      :: third dimension of field to correct
                0195 C     arrFld    :: field array to correct for imbalance
                0196 C     arrhFac   :: same size as arrFld ; thickness factor (1rst case)
                0197 C               :: or correction scaling factor (2nd case)
                0198 C     arrMask   :: 2-D mask for global mean calculation
                0199 C     arrArea   :: 2-D grid-cell area for global mean calculation
                0200 C     arrDr     :: myNr long vector (level thickness)
                0201 C     arrName   :: name of field to correct
                0202 C     myTime    :: Current time in simulation
                0203 C     myIter    :: Current iteration number in simulation
                0204 C     myThid    :: my Thread ID number
051aff16bf Mart*0205       INTEGER myNr
7e00d7e8f9 Jean*0206       _RS arrFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
051aff16bf Mart*0207       _RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
375581a429 Jean*0208       _RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
051aff16bf Mart*0209       _RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0210       _RS arrDr(myNr)
                0211       CHARACTER*(*) arrName
                0212       _RL myTime
375581a429 Jean*0213       INTEGER myIter
051aff16bf Mart*0214       INTEGER myThid
                0215 
375581a429 Jean*0216 C     !LOCAL VARIABLES:
                0217       INTEGER bi, bj, i, j, k
                0218       _RL volTile(nSx,nSy), sumTile(nSx,nSy)
7e00d7e8f9 Jean*0219       _RL tmpVol, volGlob, sumGlob, theMean
                0220       _RS meanCorr
375581a429 Jean*0221       CHARACTER*(MAX_LEN_MBUF) msgBuf
051aff16bf Mart*0222 CEOP
                0223 
                0224       DO bj=myByLo(myThid),myByHi(myThid)
                0225        DO bi=myBxLo(myThid),myBxHi(myThid)
375581a429 Jean*0226         volTile(bi,bj) = 0. _d 0
                0227         sumTile(bi,bj) = 0. _d 0
7e00d7e8f9 Jean*0228         IF ( arrDr(1).GE.zeroRS ) THEN
                0229 C--   do (arrhFac) weighted average to apply uniform correction
                0230          DO k=1,myNr
                0231           DO j=1,sNy
                0232            DO i=1,sNx
                0233             tmpVol = arrArea(i,j,bi,bj)*arrMask(i,j,bi,bj)
                0234      &             * arrhFac(i,j,k,bi,bj)*arrDr(k)
                0235             volTile(bi,bj) = volTile(bi,bj) + tmpVol
                0236             sumTile(bi,bj) = sumTile(bi,bj) + tmpVol*arrFld(i,j,k,bi,bj)
                0237            ENDDO
051aff16bf Mart*0238           ENDDO
                0239          ENDDO
7e00d7e8f9 Jean*0240         ELSE
                0241 C-    do simple average to apply (arrhFac) weighted correction
                0242          DO k=1,myNr
                0243           DO j=1,sNy
                0244            DO i=1,sNx
                0245             tmpVol = arrArea(i,j,bi,bj)*arrMask(i,j,bi,bj)
                0246      &             * ABS(arrDr(k))
                0247             volTile(bi,bj) = volTile(bi,bj)
                0248      &                     + tmpVol*arrhFac(i,j,k,bi,bj)
                0249             sumTile(bi,bj) = sumTile(bi,bj)
                0250      &                     + tmpVol*arrFld(i,j,k,bi,bj)
                0251            ENDDO
                0252           ENDDO
                0253          ENDDO
                0254         ENDIF
051aff16bf Mart*0255        ENDDO
                0256       ENDDO
                0257 
375581a429 Jean*0258       CALL GLOBAL_SUM_TILE_RL( volTile, volGlob, myThid )
                0259       CALL GLOBAL_SUM_TILE_RL( sumTile, sumGlob, myThid )
051aff16bf Mart*0260 
375581a429 Jean*0261       IF ( volGlob.GT.zeroRL ) THEN
7e00d7e8f9 Jean*0262        meanCorr = sumGlob/volGlob
051aff16bf Mart*0263        DO bj=myByLo(myThid),myByHi(myThid)
                0264         DO bi=myBxLo(myThid),myBxHi(myThid)
7e00d7e8f9 Jean*0265          IF ( arrDr(1).GE.zeroRS ) THEN
                0266 C--   apply uniform correction
                0267           DO k=1,myNr
                0268            DO j=1-OLy,sNy+OLy
                0269             DO i=1-OLx,sNx+OLx
                0270              IF ( arrhFac(i,j,k,bi,bj).NE.zeroRS ) THEN
                0271               arrFld(i,j,k,bi,bj) = arrFld(i,j,k,bi,bj) - meanCorr
                0272              ENDIF
                0273             ENDDO
051aff16bf Mart*0274            ENDDO
                0275           ENDDO
7e00d7e8f9 Jean*0276          ELSE
                0277 C-    apply (arrhFac) weighted correction
                0278           DO k=1,myNr
                0279            DO j=1-OLy,sNy+OLy
                0280             DO i=1-OLx,sNx+OLx
                0281              IF ( arrMask(i,j,bi,bj).NE.zeroRS ) THEN
                0282               arrFld(i,j,k,bi,bj) = arrFld(i,j,k,bi,bj)
                0283      &                            - meanCorr*arrhFac(i,j,k,bi,bj)
                0284              ENDIF
                0285             ENDDO
                0286            ENDDO
                0287           ENDDO
                0288 C-    end type correction selection
                0289          ENDIF
051aff16bf Mart*0290         ENDDO
                0291        ENDDO
375581a429 Jean*0292       ELSE
7e00d7e8f9 Jean*0293        meanCorr = 0. _d 0
051aff16bf Mart*0294       ENDIF
                0295 
375581a429 Jean*0296 C     Print the global mean to standard output, this is a measure
                0297 C     of the magnitude of the correction to array arrFld
051aff16bf Mart*0298       IF ( balancePrintMean ) THEN
                0299        _BEGIN_MASTER( myThid )
7e00d7e8f9 Jean*0300        IF ( arrDr(1).GE.zeroRS ) THEN
                0301         theMean = meanCorr
                0302         WRITE(msgBuf,'(3A,1PE21.14,A,I10)')
                0303      &        'REMOVE_MEAN_RS: Global mean of ', arrName, ' = ',
                0304      &        theMean, '  @ it=', myIter
                0305         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0306      &                      SQUEEZE_RIGHT, myThid )
                0307        ELSE
                0308         theMean = 0. _d 0
                0309         IF ( globalArea.GT.zeroRL ) theMean = sumGlob/globalArea
                0310         WRITE(msgBuf,'(3A,2(1PE21.14,A),I10)')
                0311      &        'REMOVE_MEAN_RS: ', arrName, ' Global mean= ',
                0312      &        theMean, ', remove: ', meanCorr, ' @ it=', myIter
                0313         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0314      &                      SQUEEZE_RIGHT, myThid )
                0315        ENDIF
051aff16bf Mart*0316        _END_MASTER( myThid )
                0317       ENDIF
                0318 
                0319       RETURN
                0320       END