C $Header: /u/gcmpack/MITgcm/model/src/find_hyd_press_1d.F,v 1.2 2016/04/04 21:29:00 jmc Exp $ C $Name: $ #include "CPP_OPTIONS.h" CBOP C !ROUTINE: FIND_HYD_PRESS_1D C !INTERFACE: SUBROUTINE FIND_HYD_PRESS_1D( O pCen, pInt, U rhoCen, I tCen, sCen, maxResid, I belowCritNb, maxIterNb, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R FIND_HYD_PRESS_1D C | o Over one column, find pressure in hydrostatic balance C | with updated density (which is function of pressure) C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" C !INPUT/OUTPUT PARAMETERS: C pCen :: Pressure vertical profile at grid-cell center C pInt :: Pressure vertical profile at grid-cell interface C rhoCen :: Density vertical profile C tCen :: Potential temperature vertical profile C sCen :: Salinity/water-vapor vertical profile C maxResid :: Maximum density difference criteria C belowCritNb :: Required number of consecutive iter below maxResid C maxIterNb :: Maximum number of iterations to perform C myThid :: my Thread Id number _RL pCen(Nr), pInt(Nr+1) _RL rhoCen(Nr) _RL tCen(Nr), sCen(Nr) _RL maxResid INTEGER belowCritNb, maxIterNb INTEGER myThid C !LOCAL VARIABLES: C msgBuf :: Informational/error message buffer C k, n :: loop counter C nIter :: iteration counter C nUnderCrit :: count number of iter below density Diff Criteria C searchForP :: Continue to iterate/search for P(rho(P)) C rhoLoc :: new density profile (computed using new Pres) C diffMax :: Maximum density difference (new minus previous) C diffRMS :: Root Mean Squared density difference CHARACTER*(MAX_LEN_MBUF) msgBuf INTEGER k, n INTEGER nIter, nUnderCrit LOGICAL searchForP _RL rhoLoc(Nr) _RL dRlocM, dRlocP, dRho _RL diffMax, diffRMS CEOP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| WRITE(msgBuf,'(A)') ' ' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,'(2A,I5,A)') 'FIND_HYD_PRESS_1D: ', & 'Start to iterate (MaxIter=', maxIterNb, & ' ) until P(rho(P))' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,'(2A,I3,A,1PE13.6)') 'FIND_HYD_PRESS_1D: ', & ' converges ; critera (x', belowCritNb, & ') on Rho diff=', maxResid CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) pInt(1) = top_Pres searchForP = .TRUE. nUnderCrit = 0 DO n=1,maxIterNb IF ( searchForP ) THEN nIter = n C-- Integrate Hydrostatic pressure IF ( integr_GeoPot.EQ.1 ) THEN DO k=1,Nr pCen(k) = pInt(k) & + rhoCen(k)*gravity*gravFacC(k)*drF(k)*halfRL pInt(k+1) = pInt(k) & + rhoCen(k)*gravity*gravFacC(k)*drF(k) ENDDO ELSE DO k=1,Nr dRlocM = halfRL*drC(k) dRlocP = halfRL*drC(k+1) IF (k.EQ.1) dRlocM = drC(k) IF (k.EQ.Nr) dRlocP = drC(k+1) pCen(k) = pInt(k) & + rhoCen(k)*gravity*gravFacF(k)*dRlocM pInt(k+1) = pCen(k) & + rhoCen(k)*gravity*gravFacF(k+1)*dRlocP ENDDO ENDIF C-- Compute new density DO k=1,Nr CALL FIND_RHO_SCALAR( I tCen(k), sCen(k), pCen(k), O rhoLoc(k), myThid ) ENDDO C-- Test for convergence: diffRMS = 0. diffMax = 0. DO k=1,Nr dRho = rhoLoc(k)-rhoCen(k) IF ( ABS(dRho) .GT. ABS(diffMax) ) diffMax = dRho diffRMS = diffRMS + dRho*dRho ENDDO diffRMS = diffRMS/DFLOAT(Nr) IF ( diffRMS.GT.0. ) diffRMS = SQRT(diffRMS) IF ( ABS(diffMax).LE.maxResid ) THEN nUnderCrit = nUnderCrit + 1 ELSE nUnderCrit = 0 ENDIF C- Double criteria: stop if perfect convergence or if below C criteria for at least "belowCritNb" iterations searchForP = ( nUnderCrit.LT.belowCritNb ) & .AND. ( diffMax.NE.zeroRL ) IF ( debugLevel.GE.debLevB ) THEN WRITE(msgBuf,'(A,I5,2(A,1PE20.12))') & ' iter', nIter,', RMS-diff=', diffRMS,', Max-diff=', diffMax CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) ENDIF IF ( searchForP .AND. nIter.LT.maxIterNb ) THEN C-- Update rhoCen for new iteration DO k=1,Nr rhoCen(k) = rhoLoc(k) ENDDO ENDIF ENDIF ENDDO IF ( searchForP ) THEN WRITE(msgBuf,'(2A,I5,A,I3,A)') 'FIND_HYD_PRESS_1D: ', & 'No convergence after', nIter, & ' iters (nUnderCrit=', nUnderCrit, ' )' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R FIND_HYD_PRESS_1D' ELSE WRITE(msgBuf,'(2A,I5,A,I3,A)') 'FIND_HYD_PRESS_1D: ', & 'converged after', nIter, & ' iters (nUnderCrit=', nUnderCrit, ' )' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) C-- Update rhoCen with new solution DO k=1,Nr rhoCen(k) = rhoLoc(k) ENDDO ENDIF RETURN END