| 1 |
C $Header: /u/gcmpack/MITgcm/model/src/find_hyd_press_1d.F,v 1.1 2016/03/10 20:54:57 jmc Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "CPP_OPTIONS.h" |
| 5 |
|
| 6 |
CBOP |
| 7 |
C !ROUTINE: FIND_HYD_PRESS_1D |
| 8 |
C !INTERFACE: |
| 9 |
SUBROUTINE FIND_HYD_PRESS_1D( |
| 10 |
O pCen, pInt, |
| 11 |
U rhoCen, |
| 12 |
I tCen, sCen, maxResid, |
| 13 |
I belowCritNb, maxIterNb, myThid ) |
| 14 |
|
| 15 |
C !DESCRIPTION: \bv |
| 16 |
C *==========================================================* |
| 17 |
C | S/R FIND_HYD_PRESS_1D |
| 18 |
C | o Over one column, find pressure in hydrostatic balance |
| 19 |
C | with updated density (which is function of pressure) |
| 20 |
C *==========================================================* |
| 21 |
C \ev |
| 22 |
|
| 23 |
C !USES: |
| 24 |
IMPLICIT NONE |
| 25 |
C == Global variables === |
| 26 |
#include "SIZE.h" |
| 27 |
#include "EEPARAMS.h" |
| 28 |
#include "PARAMS.h" |
| 29 |
#include "GRID.h" |
| 30 |
|
| 31 |
C !INPUT/OUTPUT PARAMETERS: |
| 32 |
C pCen :: Pressure vertical profile at grid-cell center |
| 33 |
C pInt :: Pressure vertical profile at grid-cell interface |
| 34 |
C rhoCen :: Density vertical profile |
| 35 |
C tCen :: Potential temperature vertical profile |
| 36 |
C sCen :: Salinity/water-vapor vertical profile |
| 37 |
C maxResid :: Maximum density difference criteria |
| 38 |
C belowCritNb :: Required number of consecutive iter below maxResid |
| 39 |
C maxIterNb :: Maximum number of iterations to perform |
| 40 |
C myThid :: my Thread Id number |
| 41 |
_RL pCen(Nr), pInt(Nr+1) |
| 42 |
_RL rhoCen(Nr) |
| 43 |
_RL tCen(Nr), sCen(Nr) |
| 44 |
_RL maxResid |
| 45 |
INTEGER belowCritNb, maxIterNb |
| 46 |
INTEGER myThid |
| 47 |
|
| 48 |
C !LOCAL VARIABLES: |
| 49 |
C msgBuf :: Informational/error message buffer |
| 50 |
C k, n :: loop counter |
| 51 |
C nIter :: iteration counter |
| 52 |
C nUnderCrit :: count number of iter below density Diff Criteria |
| 53 |
C searchForP :: Continue to iterate/search for P(rho(P)) |
| 54 |
C rhoLoc :: new density profile (computed using new Pres) |
| 55 |
C diffMax :: Maximum density difference (new minus previous) |
| 56 |
C diffRMS :: Root Mean Squared density difference |
| 57 |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
| 58 |
INTEGER k, n |
| 59 |
INTEGER nIter, nUnderCrit |
| 60 |
LOGICAL searchForP |
| 61 |
_RL rhoLoc(Nr) |
| 62 |
_RL dRlocM, dRlocP, dRho |
| 63 |
_RL diffMax, diffRMS |
| 64 |
CEOP |
| 65 |
|
| 66 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 67 |
|
| 68 |
WRITE(msgBuf,'(A)') ' ' |
| 69 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
| 70 |
& SQUEEZE_RIGHT, myThid ) |
| 71 |
WRITE(msgBuf,'(2A,I5,A)') 'FIND_HYD_PRESS_1D: ', |
| 72 |
& 'Start to iterate (MaxIter=', maxIterNb, |
| 73 |
& ' ) until P(rho(P))' |
| 74 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
| 75 |
& SQUEEZE_RIGHT, myThid ) |
| 76 |
WRITE(msgBuf,'(2A,I3,A,1PE13.6)') 'FIND_HYD_PRESS_1D: ', |
| 77 |
& ' converges ; critera (x', belowCritNb, |
| 78 |
& ') on Rho diff=', maxResid |
| 79 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
| 80 |
& SQUEEZE_RIGHT, myThid ) |
| 81 |
|
| 82 |
pInt(1) = top_Pres |
| 83 |
searchForP = .TRUE. |
| 84 |
nUnderCrit = 0 |
| 85 |
DO n=1,maxIterNb |
| 86 |
IF ( searchForP ) THEN |
| 87 |
nIter = n |
| 88 |
|
| 89 |
C-- Integrate Hydrostatic pressure |
| 90 |
IF ( integr_GeoPot.EQ.1 ) THEN |
| 91 |
DO k=1,Nr |
| 92 |
pCen(k) = pInt(k) |
| 93 |
& + rhoCen(k)*gravity*gravFacC(k)*drF(k)*halfRL |
| 94 |
pInt(k+1) = pInt(k) |
| 95 |
& + rhoCen(k)*gravity*gravFacC(k)*drF(k) |
| 96 |
ENDDO |
| 97 |
ELSE |
| 98 |
DO k=1,Nr |
| 99 |
dRlocM = halfRL*drC(k) |
| 100 |
dRlocP = halfRL*drC(k+1) |
| 101 |
IF (k.EQ.1) dRlocM = drC(k) |
| 102 |
IF (k.EQ.Nr) dRlocP = drC(k+1) |
| 103 |
pCen(k) = pInt(k) |
| 104 |
& + rhoCen(k)*gravity*gravFacF(k)*dRlocM |
| 105 |
pInt(k+1) = pCen(k) |
| 106 |
& + rhoCen(k)*gravity*gravFacF(k+1)*dRlocP |
| 107 |
ENDDO |
| 108 |
ENDIF |
| 109 |
|
| 110 |
C-- Compute new density |
| 111 |
DO k=1,Nr |
| 112 |
CALL FIND_RHO_SCALAR( |
| 113 |
I tCen(k), sCen(k), pCen(k), |
| 114 |
O rhoLoc(k), myThid ) |
| 115 |
ENDDO |
| 116 |
|
| 117 |
C-- Test for convergence: |
| 118 |
diffRMS = 0. |
| 119 |
diffMax = 0. |
| 120 |
DO k=1,Nr |
| 121 |
dRho = rhoLoc(k)-rhoCen(k) |
| 122 |
IF ( ABS(dRho) .GT. ABS(diffMax) ) diffMax = dRho |
| 123 |
diffRMS = diffRMS + dRho*dRho |
| 124 |
ENDDO |
| 125 |
diffRMS = diffRMS/DFLOAT(Nr) |
| 126 |
IF ( diffRMS.GT.0. ) diffRMS = SQRT(diffRMS) |
| 127 |
IF ( ABS(diffMax).LE.maxResid ) THEN |
| 128 |
nUnderCrit = nUnderCrit + 1 |
| 129 |
ELSE |
| 130 |
nUnderCrit = 0 |
| 131 |
ENDIF |
| 132 |
C- Double criteria: stop if perfect convergence or if below |
| 133 |
C criteria for at least "belowCritNb" iterations |
| 134 |
searchForP = ( nUnderCrit.LT.belowCritNb ) |
| 135 |
& .AND. ( diffMax.NE.zeroRL ) |
| 136 |
|
| 137 |
IF ( debugLevel.GE.debLevB ) THEN |
| 138 |
WRITE(msgBuf,'(A,I5,2(A,1PE20.12))') |
| 139 |
& ' iter', nIter,', RMS-diff=', diffRMS,', Max-diff=', diffMax |
| 140 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
| 141 |
& SQUEEZE_RIGHT, myThid ) |
| 142 |
ENDIF |
| 143 |
|
| 144 |
IF ( searchForP .AND. nIter.LT.maxIterNb ) THEN |
| 145 |
C-- Update rhoCen for new iteration |
| 146 |
DO k=1,Nr |
| 147 |
rhoCen(k) = rhoLoc(k) |
| 148 |
ENDDO |
| 149 |
ENDIF |
| 150 |
|
| 151 |
ENDIF |
| 152 |
ENDDO |
| 153 |
|
| 154 |
IF ( searchForP ) THEN |
| 155 |
WRITE(msgBuf,'(2A,I5,A,I3,A)') 'FIND_HYD_PRESS_1D: ', |
| 156 |
& 'No convergence after', nIter, |
| 157 |
& ' iters (nUnderCrit=', nUnderCrit, ' )' |
| 158 |
CALL PRINT_ERROR( msgBuf, myThid ) |
| 159 |
STOP 'ABNORMAL END: S/R FIND_HYD_PRESS_1D' |
| 160 |
ELSE |
| 161 |
WRITE(msgBuf,'(2A,I5,A,I3,A)') 'FIND_HYD_PRESS_1D: ', |
| 162 |
& 'converged after', nIter, |
| 163 |
& ' iters (nUnderCrit=', nUnderCrit, ' )' |
| 164 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
| 165 |
& SQUEEZE_RIGHT, myThid ) |
| 166 |
C-- Update rhoCen with new solution |
| 167 |
DO k=1,Nr |
| 168 |
rhoCen(k) = rhoLoc(k) |
| 169 |
ENDDO |
| 170 |
ENDIF |
| 171 |
|
| 172 |
RETURN |
| 173 |
END |
Parent Directory
| 
Revision Log
| 
Revision Graph