model/src/find_hyd_press_1d.F

759f945c40 Jean*0001 #include "CPP_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: FIND_HYD_PRESS_1D
                0005 C     !INTERFACE:
                0006       SUBROUTINE FIND_HYD_PRESS_1D(
                0007      O                    pCen, pInt,
                0008      U                    rhoCen,
                0009      I                    tCen, sCen, maxResid,
                0010      I                    belowCritNb, maxIterNb, myThid )
                0011
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | S/R FIND_HYD_PRESS_1D
                0015 C     | o Over one column, find pressure in hydrostatic balance
                0016 C     |   with updated density (which is function of pressure)
                0017 C     *==========================================================*
                0018 C     \ev
                0019
                0020 C     !USES:
                0021       IMPLICIT NONE
                0022 C     == Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027
                0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     pCen      :: Pressure vertical profile at grid-cell center
                0030 C     pInt      :: Pressure vertical profile at grid-cell interface
                0031 C     rhoCen    :: Density vertical profile
                0032 C     tCen      :: Potential temperature vertical profile
                0033 C     sCen      :: Salinity/water-vapor vertical profile
                0034 C     maxResid  :: Maximum density difference criteria
                0035 C   belowCritNb :: Required number of consecutive iter below maxResid
                0036 C     maxIterNb :: Maximum number of iterations to perform
                0037 C     myThid    :: my Thread Id number
                0038       _RL     pCen(Nr), pInt(Nr+1)
                0039       _RL     rhoCen(Nr)
                0040       _RL     tCen(Nr), sCen(Nr)
                0041       _RL     maxResid
                0042       INTEGER belowCritNb, maxIterNb
                0043       INTEGER myThid
                0044
                0045 C     !LOCAL VARIABLES:
                0046 C     msgBuf     :: Informational/error message buffer
                0047 C     k, n       :: loop counter
                0048 C     nIter      :: iteration counter
                0049 C     nUnderCrit :: count number of iter below density Diff Criteria
                0050 C     searchForP :: Continue to iterate/search for P(rho(P))
                0051 C     rhoLoc     :: new density profile (computed using new Pres)
                0052 C     diffMax    :: Maximum density difference (new minus previous)
                0053 C     diffRMS    :: Root Mean Squared density difference
                0054       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0055       INTEGER k, n
                0056       INTEGER nIter, nUnderCrit
                0057       LOGICAL searchForP
                0058       _RL     rhoLoc(Nr)
                0059       _RL     dRlocM, dRlocP, dRho
                0060       _RL     diffMax, diffRMS
                0061 CEOP
                0062
                0063 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0064
                0065       WRITE(msgBuf,'(A)') ' '
                0066       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0067      &                    SQUEEZE_RIGHT, myThid )
                0068       WRITE(msgBuf,'(2A,I5,A)') 'FIND_HYD_PRESS_1D: ',
                0069      &  'Start to iterate (MaxIter=', maxIterNb,
                0070      &   ' ) until P(rho(P))'
                0071       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0072      &                    SQUEEZE_RIGHT, myThid )
                0073       WRITE(msgBuf,'(2A,I3,A,1PE13.6)') 'FIND_HYD_PRESS_1D: ',
                0074      &   ' converges ; critera (x', belowCritNb,
                0075      &   ') on Rho diff=', maxResid
                0076       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0077      &                    SQUEEZE_RIGHT, myThid )
                0078
2ad79bdf32 Jean*0079       pInt(1)  = top_Pres
759f945c40 Jean*0080       searchForP = .TRUE.
                0081       nUnderCrit = 0
                0082       DO n=1,maxIterNb
                0083        IF ( searchForP ) THEN
                0084         nIter = n
                0085
                0086 C--   Integrate Hydrostatic pressure
                0087         IF ( integr_GeoPot.EQ.1 ) THEN
                0088          DO k=1,Nr
                0089           pCen(k)   = pInt(k)
                0090      &              + rhoCen(k)*gravity*gravFacC(k)*drF(k)*halfRL
                0091           pInt(k+1) = pInt(k)
                0092      &              + rhoCen(k)*gravity*gravFacC(k)*drF(k)
                0093          ENDDO
                0094         ELSE
                0095          DO k=1,Nr
                0096           dRlocM = halfRL*drC(k)
                0097           dRlocP = halfRL*drC(k+1)
                0098           IF (k.EQ.1)  dRlocM = drC(k)
                0099           IF (k.EQ.Nr) dRlocP = drC(k+1)
                0100           pCen(k)   = pInt(k)
                0101      &              + rhoCen(k)*gravity*gravFacF(k)*dRlocM
                0102           pInt(k+1) = pCen(k)
                0103      &              + rhoCen(k)*gravity*gravFacF(k+1)*dRlocP
                0104          ENDDO
                0105         ENDIF
                0106
                0107 C--   Compute new density
                0108         DO k=1,Nr
                0109           CALL FIND_RHO_SCALAR(
                0110      I                          tCen(k), sCen(k), pCen(k),
                0111      O                          rhoLoc(k), myThid )
                0112         ENDDO
                0113
                0114 C--   Test for convergence:
                0115         diffRMS = 0.
                0116         diffMax = 0.
                0117         DO k=1,Nr
                0118           dRho = rhoLoc(k)-rhoCen(k)
                0119           IF ( ABS(dRho) .GT. ABS(diffMax) ) diffMax = dRho
                0120           diffRMS = diffRMS + dRho*dRho
                0121         ENDDO
                0122         diffRMS = diffRMS/DFLOAT(Nr)
                0123         IF ( diffRMS.GT.0. ) diffRMS = SQRT(diffRMS)
                0124         IF ( ABS(diffMax).LE.maxResid ) THEN
                0125           nUnderCrit = nUnderCrit + 1
                0126         ELSE
                0127           nUnderCrit = 0
                0128         ENDIF
                0129 C-    Double criteria: stop if perfect convergence or if below
                0130 C      criteria for at least "belowCritNb" iterations
                0131         searchForP = ( nUnderCrit.LT.belowCritNb )
                0132      &         .AND. ( diffMax.NE.zeroRL )
                0133
                0134         IF ( debugLevel.GE.debLevB ) THEN
                0135          WRITE(msgBuf,'(A,I5,2(A,1PE20.12))')
                0136      &    ' iter', nIter,', RMS-diff=', diffRMS,', Max-diff=', diffMax
                0137          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0138      &                       SQUEEZE_RIGHT, myThid )
                0139         ENDIF
                0140
                0141         IF ( searchForP .AND. nIter.LT.maxIterNb ) THEN
                0142 C--   Update rhoCen for new iteration
                0143          DO k=1,Nr
                0144           rhoCen(k) = rhoLoc(k)
                0145          ENDDO
                0146         ENDIF
                0147
                0148        ENDIF
                0149       ENDDO
                0150
                0151       IF ( searchForP ) THEN
                0152         WRITE(msgBuf,'(2A,I5,A,I3,A)') 'FIND_HYD_PRESS_1D: ',
                0153      &    'No convergence after', nIter,
                0154      &    ' iters (nUnderCrit=', nUnderCrit, ' )'
                0155         CALL PRINT_ERROR( msgBuf, myThid )
                0156         STOP 'ABNORMAL END: S/R FIND_HYD_PRESS_1D'
                0157       ELSE
                0158         WRITE(msgBuf,'(2A,I5,A,I3,A)') 'FIND_HYD_PRESS_1D: ',
                0159      &    'converged after', nIter,
                0160      &    ' iters (nUnderCrit=', nUnderCrit, ' )'
                0161         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0162      &                      SQUEEZE_RIGHT, myThid )
                0163 C--   Update rhoCen with new solution
                0164         DO k=1,Nr
                0165           rhoCen(k) = rhoLoc(k)
                0166         ENDDO
                0167       ENDIF
                0168
                0169       RETURN
                0170       END