/[MITgcm]/MITgcm/model/src/solve_for_pressure.F
ViewVC logotype

Diff of /MITgcm/model/src/solve_for_pressure.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by cnh, Mon Jun 8 21:43:02 1998 UTC revision 1.64 by jmc, Sun Aug 24 21:40:18 2008 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )  C     !ROUTINE: SOLVE_FOR_PRESSURE
9  C     /==========================================================\  C     !INTERFACE:
10  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |        SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
 C     | o Controls inversion of two and/or three-dimensional     |  
 C     |   elliptic problems for the pressure field.              |  
 C     \==========================================================/  
11    
12    C     !DESCRIPTION: \bv
13    C     *==========================================================*
14    C     | SUBROUTINE SOLVE_FOR_PRESSURE
15    C     | o Controls inversion of two and/or three-dimensional
16    C     |   elliptic problems for the pressure field.
17    C     *==========================================================*
18    C     \ev
19    
20    C     !USES:
21          IMPLICIT NONE
22  C     == Global variables  C     == Global variables
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26    #include "GRID.h"
27    #include "SURFACE.h"
28    #include "FFIELDS.h"
29  #include "DYNVARS.h"  #include "DYNVARS.h"
30  #include "CG2D.h"  #include "SOLVE_FOR_PRESSURE.h"
31    #ifdef ALLOW_NONHYDROSTATIC
32    #include "SOLVE_FOR_PRESSURE3D.h"
33    #include "NH_VARS.h"
34    #endif
35    #ifdef ALLOW_CD_CODE
36    #include "CD_CODE_VARS.h"
37    #endif
38    #ifdef ALLOW_OBCS
39    #include "OBCS.h"
40    #endif
41    
42    C     === Functions ====
43          LOGICAL  DIFFERENT_MULTIPLE
44          EXTERNAL DIFFERENT_MULTIPLE
45    
46    C     !INPUT/OUTPUT PARAMETERS:
47  C     == Routine arguments ==  C     == Routine arguments ==
48  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myTime :: Current time in simulation
49    C     myIter :: Current iteration number in simulation
50    C     myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
51          _RL myTime
52          INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
54    
55  C     Local variables  C     !LOCAL VARIABLES:
56        INTEGER i,j,bi,bj  C     == Local variables ==
57          INTEGER i,j,k,bi,bj
58          _RL firstResidual,lastResidual
59          _RL tmpFac
60          _RL sumEmP, tileEmP
61          LOGICAL putPmEinXvector
62          INTEGER numIters, ks
63          CHARACTER*10 sufx
64          CHARACTER*(MAX_LEN_MBUF) msgBuf
65    #ifdef ALLOW_NONHYDROSTATIC
66          INTEGER kp1
67          _RL     wFacKm, wFacKp
68          LOGICAL zeroPsNH
69          _RL     uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70          _RL     vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71    #else
72          _RL     cg3d_b(1)
73    #endif
74    CEOP
75    
76    #ifdef ALLOW_NONHYDROSTATIC
77    c       zeroPsNH = .FALSE.
78            zeroPsNH = exactConserv
79    #else
80            cg3d_b(1) = 0.
81    #endif
82    
83    C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
84    C anelastic (always Z-coordinate):
85    C     1) assume that rhoFacF(1)=1 (and ksurf == 1);
86    C        (this reduces the number of lines of code to modify)
87    C     2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
88    C        (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
89    C       => 2 factors cancel in elliptic eq. for Phi_s ,
90    C       but 1rst factor(a) remains in RHS cg2d_b.
91    
92  #ifdef ALLOW_CD  C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
93  C--   Save previous solution.  C     instead of simply etaN ; This can speed-up the solver convergence in
94    C     the case where |Global_mean_PmE| is large.
95          putPmEinXvector = .FALSE.
96    c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
97    
98    C--   Save previous solution & Initialise Vector solution and source term :
99          sumEmP = 0.
100        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
101         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
102          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
103           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
104            cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)  #ifdef ALLOW_CD_CODE
105              etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
106    #endif
107              cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
108              cg2d_b(i,j,bi,bj) = 0.
109           ENDDO           ENDDO
110          ENDDO          ENDDO
111            IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
112             tmpFac = freeSurfFac*mass2rUnit
113             IF (exactConserv)
114         &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
115             DO j=1,sNy
116              DO i=1,sNx
117               cg2d_b(i,j,bi,bj) =
118         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
119              ENDDO
120             ENDDO
121            ENDIF
122            IF ( putPmEinXvector ) THEN
123             tileEmP = 0.
124             DO j=1,sNy
125              DO i=1,sNx
126                tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
127         &                                       *maskH(i,j,bi,bj)
128              ENDDO
129             ENDDO
130             sumEmP = sumEmP + tileEmP
131            ENDIF
132         ENDDO         ENDDO
133        ENDDO        ENDDO
134          IF ( putPmEinXvector ) THEN
135            _GLOBAL_SUM_R8( sumEmP, myThid )
136          ENDIF
137    
138          DO bj=myByLo(myThid),myByHi(myThid)
139           DO bi=myBxLo(myThid),myBxHi(myThid)
140            IF ( putPmEinXvector ) THEN
141              tmpFac = 0.
142              IF (globalArea.GT.0.) tmpFac =
143         &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
144              DO j=1,sNy
145               DO i=1,sNx
146                cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
147         &                        - tmpFac*Bo_surf(i,j,bi,bj)
148               ENDDO
149              ENDDO
150            ENDIF
151    C- RHS: similar to the divergence of the vertically integrated mass transport:
152    C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
153            DO k=Nr,1,-1
154             CALL CALC_DIV_GHAT(
155         I                       bi,bj,k,
156         U                       cg2d_b, cg3d_b,
157         I                       myThid )
158            ENDDO
159           ENDDO
160          ENDDO
161    
162    C--   Add source term arising from w=d/dt (p_s + p_nh)
163          DO bj=myByLo(myThid),myByHi(myThid)
164           DO bi=myBxLo(myThid),myBxHi(myThid)
165    #ifdef ALLOW_NONHYDROSTATIC
166            IF ( use3Dsolver .AND. zeroPsNH ) THEN
167             DO j=1,sNy
168              DO i=1,sNx
169               ks = ksurfC(i,j,bi,bj)
170               IF ( ks.LE.Nr ) THEN
171                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
172         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
173         &         /deltaTMom/deltaTfreesurf
174         &         * etaH(i,j,bi,bj)
175                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
176         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
177         &         /deltaTMom/deltaTfreesurf
178         &         * etaH(i,j,bi,bj)
179               ENDIF
180              ENDDO
181             ENDDO
182            ELSEIF ( use3Dsolver ) THEN
183             DO j=1,sNy
184              DO i=1,sNx
185               ks = ksurfC(i,j,bi,bj)
186               IF ( ks.LE.Nr ) THEN
187                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
188         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
189         &         /deltaTMom/deltaTfreesurf
190         &         *( etaN(i,j,bi,bj)
191         &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
192                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
193         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
194         &         /deltaTMom/deltaTfreesurf
195         &         *( etaN(i,j,bi,bj)
196         &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
197               ENDIF
198              ENDDO
199             ENDDO
200            ELSEIF ( exactConserv ) THEN
201    #else
202            IF ( exactConserv ) THEN
203    #endif /* ALLOW_NONHYDROSTATIC */
204             DO j=1,sNy
205              DO i=1,sNx
206               ks = ksurfC(i,j,bi,bj)
207               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
208         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
209         &         /deltaTMom/deltaTfreesurf
210         &         * etaH(i,j,bi,bj)
211              ENDDO
212             ENDDO
213            ELSE
214             DO j=1,sNy
215              DO i=1,sNx
216               ks = ksurfC(i,j,bi,bj)
217               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
218         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
219         &         /deltaTMom/deltaTfreesurf
220         &         * etaN(i,j,bi,bj)
221              ENDDO
222             ENDDO
223            ENDIF
224    
225    #ifdef ALLOW_OBCS
226            IF (useOBCS) THEN
227             DO i=1,sNx
228    C Northern boundary
229              IF (OB_Jn(i,bi,bj).NE.0) THEN
230               cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
231               cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
232              ENDIF
233    C Southern boundary
234              IF (OB_Js(i,bi,bj).NE.0) THEN
235               cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
236               cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
237              ENDIF
238             ENDDO
239             DO j=1,sNy
240    C Eastern boundary
241              IF (OB_Ie(j,bi,bj).NE.0) THEN
242               cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
243               cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
244              ENDIF
245    C Western boundary
246              IF (OB_Iw(j,bi,bj).NE.0) THEN
247               cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
248               cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
249              ENDIF
250             ENDDO
251            ENDIF
252    #endif /* ALLOW_OBCS */
253    C-    end bi,bj loops
254           ENDDO
255          ENDDO
256    
257    #ifdef ALLOW_DEBUG
258          IF ( debugLevel .GE. debLevB ) THEN
259           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
260         &                        myThid)
261          ENDIF
262  #endif  #endif
263          IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
264           WRITE(sufx,'(I10.10)') myIter
265           CALL WRITE_FLD_XY_RS( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
266          ENDIF
267    
268  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
269  C--   gradient solver.  C--   gradient solver.
270  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
271          firstResidual=0.
272          lastResidual=0.
273          numIters=cg2dMaxIters
274    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
275    #ifdef ALLOW_CG2D_NSA
276    C--   Call the not-self-adjoint version of cg2d
277          CALL CG2D_NSA(
278         U           cg2d_b,
279         U           cg2d_x,
280         O           firstResidual,
281         O           lastResidual,
282         U           numIters,
283         I           myThid )
284    #else /* not ALLOW_CG2D_NSA = default */
285        CALL CG2D(        CALL CG2D(
286         U           cg2d_b,
287         U           cg2d_x,
288         O           firstResidual,
289         O           lastResidual,
290         U           numIters,
291       I           myThid )       I           myThid )
292    #endif /* ALLOW_CG2D_NSA */
293          _EXCH_XY_R8(cg2d_x, myThid )
294    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
295    
296    #ifdef ALLOW_DEBUG
297          IF ( debugLevel .GE. debLevB ) THEN
298           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
299         &                        myThid)
300          ENDIF
301    #endif
302    
303    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
304          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
305         &   ) THEN
306           IF ( debugLevel .GE. debLevA ) THEN
307            _BEGIN_MASTER( myThid )
308            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
309            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
310            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
311            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
312            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
313            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
314            _END_MASTER( myThid )
315           ENDIF
316          ENDIF
317    
318    C--   Transfert the 2D-solution to "etaN" :
319          DO bj=myByLo(myThid),myByHi(myThid)
320           DO bi=myBxLo(myThid),myBxHi(myThid)
321            DO j=1-OLy,sNy+OLy
322             DO i=1-OLx,sNx+OLx
323              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
324             ENDDO
325            ENDDO
326           ENDDO
327          ENDDO
328    
329    #ifdef ALLOW_NONHYDROSTATIC
330          IF ( use3Dsolver ) THEN
331    
332  C--   Solve for a three-dimensional pressure term (NH or IGW or both ).  C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
333  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
334  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
335  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
336             DO j=1,sNy+1
337              DO i=1,sNx+1
338               uf(i,j)=-_recip_dxC(i,j,bi,bj)*
339         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
340               vf(i,j)=-_recip_dyC(i,j,bi,bj)*
341         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
342              ENDDO
343             ENDDO
344    
345    #ifdef ALLOW_OBCS
346             IF (useOBCS) THEN
347              DO i=1,sNx+1
348    C Northern boundary
349              IF (OB_Jn(i,bi,bj).NE.0) THEN
350               vf(i,OB_Jn(i,bi,bj))=0.
351              ENDIF
352    C Southern boundary
353              IF (OB_Js(i,bi,bj).NE.0) THEN
354               vf(i,OB_Js(i,bi,bj)+1)=0.
355              ENDIF
356              ENDDO
357              DO j=1,sNy+1
358    C Eastern boundary
359              IF (OB_Ie(j,bi,bj).NE.0) THEN
360               uf(OB_Ie(j,bi,bj),j)=0.
361              ENDIF
362    C Western boundary
363              IF (OB_Iw(j,bi,bj).NE.0) THEN
364               uf(OB_Iw(j,bi,bj)+1,J)=0.
365              ENDIF
366              ENDDO
367             ENDIF
368    #endif /* ALLOW_OBCS */
369    
370             IF ( usingZCoords ) THEN
371    C-       Z coordinate: assume surface @ level k=1
372               tmpFac = freeSurfFac*deepFac2F(1)
373             ELSE
374    C-       Other than Z coordinate: no assumption on surface level index
375               tmpFac = 0.
376               DO j=1,sNy
377                DO i=1,sNx
378                  ks = ksurfC(i,j,bi,bj)
379                  IF ( ks.LE.Nr ) THEN
380                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
381         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
382         &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom
383                  ENDIF
384                ENDDO
385               ENDDO
386             ENDIF
387             k=1
388             kp1 = MIN(k+1,Nr)
389             wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
390             IF (k.GE.Nr) wFacKp = 0.
391             DO j=1,sNy
392              DO i=1,sNx
393                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
394         &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
395         &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
396         &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
397         &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
398         &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
399         &         -wVel(i,j,kp1,bi,bj)*wFacKp
400         &        )*_rA(i,j,bi,bj)/deltaTmom
401              ENDDO
402             ENDDO
403             DO k=2,Nr
404              kp1 = MIN(k+1,Nr)
405    C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
406    C        both appear in wVel term, but at 2 different levels
407              wFacKm = deepFac2F( k )*rhoFacF( k )
408              wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
409              IF (k.GE.Nr) wFacKp = 0.
410              DO j=1,sNy
411               DO i=1,sNx
412                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
413         &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
414         &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
415         &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
416         &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
417         &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
418         &         -wVel(i,j,kp1,bi,bj)*wFacKp
419         &        )*_rA(i,j,bi,bj)/deltaTmom
420    
421               ENDDO
422              ENDDO
423             ENDDO
424    
425    #ifdef ALLOW_OBCS
426             IF (useOBCS) THEN
427              DO k=1,Nr
428              DO i=1,sNx
429    C Northern boundary
430              IF (OB_Jn(i,bi,bj).NE.0) THEN
431               cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj)=0.
432              ENDIF
433    C Southern boundary
434              IF (OB_Js(i,bi,bj).NE.0) THEN
435               cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj)=0.
436              ENDIF
437              ENDDO
438              DO j=1,sNy
439    C Eastern boundary
440              IF (OB_Ie(j,bi,bj).NE.0) THEN
441               cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj)=0.
442              ENDIF
443    C Western boundary
444              IF (OB_Iw(j,bi,bj).NE.0) THEN
445               cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj)=0.
446              ENDIF
447              ENDDO
448              ENDDO
449             ENDIF
450    #endif /* ALLOW_OBCS */
451    C-    end bi,bj loops
452            ENDDO
453           ENDDO
454    
455          firstResidual=0.
456          lastResidual=0.
457          numIters=cg3dMaxIters
458          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
459          CALL CG3D(
460         U           cg3d_b,
461         U           phi_nh,
462         O           firstResidual,
463         O           lastResidual,
464         U           numIters,
465         I           myThid )
466          _EXCH_XYZ_R8(phi_nh, myThid )
467          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
468    
469          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
470         &   ) THEN
471           IF ( debugLevel .GE. debLevA ) THEN
472            _BEGIN_MASTER( myThid )
473            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
474            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
475            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
476            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
477            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
478            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
479            _END_MASTER( myThid )
480           ENDIF
481          ENDIF
482    
483    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
484          IF ( zeroPsNH ) THEN
485           DO bj=myByLo(myThid),myByHi(myThid)
486            DO bi=myBxLo(myThid),myBxHi(myThid)
487    
488             IF ( usingZCoords ) THEN
489    C-       Z coordinate: assume surface @ level k=1
490              DO k=2,Nr
491               DO j=1-OLy,sNy+OLy
492                DO i=1-OLx,sNx+OLx
493                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
494         &                           - phi_nh(i,j,1,bi,bj)
495                ENDDO
496               ENDDO
497              ENDDO
498              DO j=1-OLy,sNy+OLy
499               DO i=1-OLx,sNx+OLx
500                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
501         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
502                 phi_nh(i,j,1,bi,bj) = 0.
503               ENDDO
504              ENDDO
505             ELSE
506    C-       Other than Z coordinate: no assumption on surface level index
507              DO j=1-OLy,sNy+OLy
508               DO i=1-OLx,sNx+OLx
509                ks = ksurfC(i,j,bi,bj)
510                IF ( ks.LE.Nr ) THEN
511                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
512         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
513                 DO k=Nr,1,-1
514                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
515         &                            - phi_nh(i,j,ks,bi,bj)
516                 ENDDO
517                ENDIF
518               ENDDO
519              ENDDO
520             ENDIF
521    
522            ENDDO
523           ENDDO
524          ENDIF
525    
526          ENDIF
527    #endif /* ALLOW_NONHYDROSTATIC */
528    
529    #ifdef ALLOW_SHOWFLOPS
530          CALL SHOWFLOPS_INSOLVE( myThid)
531    #endif
532    
533        RETURN        RETURN
534        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.64
  ViewVC Help
Powered by ViewVC 1.1.22