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 |