/[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.17 by jmc, Tue Mar 6 16:57:10 2001 UTC revision 1.58 by jmc, Tue Dec 5 05:25:08 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #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)
11  C     | o Controls inversion of two and/or three-dimensional     |  
12  C     |   elliptic problems for the pressure field.              |  C     !DESCRIPTION: \bv
13  C     \==========================================================/  C     *==========================================================*
14        IMPLICIT NONE  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"
 #include "DYNVARS.h"  
26  #include "GRID.h"  #include "GRID.h"
27  #include "SURFACE.h"  #include "SURFACE.h"
28    #include "FFIELDS.h"
29    #include "DYNVARS.h"
30    #include "SOLVE_FOR_PRESSURE.h"
31  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
32  #include "CG3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
33  #include "GW.h"  #include "NH_VARS.h"
34    #endif
35    #ifdef ALLOW_CD_CODE
36    #include "CD_CODE_VARS.h"
37  #endif  #endif
38  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
39  #include "OBCS.h"  #include "OBCS.h"
40  #endif  #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  C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector  C     == Local variables ==
 C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  
57        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
58        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL firstResidual,lastResidual
61        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL tmpFac
62          _RL sumEmP, tileEmP
63          LOGICAL putPmEinXvector
64          INTEGER numIters, ks
65          CHARACTER*(MAX_LEN_MBUF) msgBuf
66    #ifdef ALLOW_NONHYDROSTATIC
67          INTEGER kp1
68          _RL     wFacKm, wFacKp
69          LOGICAL zeroPsNH
70    #endif
71    CEOP
72    
73    #ifdef TIME_PER_TIMESTEP_SFP
74    CCE107 common block for per timestep timing
75    C     !TIMING VARIABLES
76    C     == Timing variables ==
77          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
78          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
79    #endif
80    #ifdef USE_PAPI_FLOPS_SFP
81    CCE107 common block for PAPI summary performance
82    #include <fpapi.h>
83          INTEGER*8 flpops, instr
84          INTEGER check
85          REAL*4 real_time, proc_time, mflops, ipc
86          COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
87    #else
88    #ifdef USE_PCL_FLOPS_SFP
89    CCE107 common block for PCL summary performance
90    #include <pclh.f>
91          INTEGER pcl_counter_list(5), flags, nevents, res, ipcl
92          INTEGER*8 i_result(5), descr
93          REAL*8 fp_result(5)
94          COMMON /pclvars/ i_result, descr, fp_result, pcl_counter_list,
95         $     flags, nevents
96          INTEGER nmaxevents
97          PARAMETER (nmaxevents = 61)
98          CHARACTER*22 pcl_counter_name(0:nmaxevents-1)
99          COMMON /pclnames/ pcl_counter_name
100    #endif
101    #endif
102    
103    #ifdef ALLOW_NONHYDROSTATIC
104    c       zeroPsNH = .FALSE.
105            zeroPsNH = exactConserv
106    #endif
107    
108    C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
109    C anelastic (always Z-coordinate):
110    C     1) assume that rhoFacF(1)=1 (and ksurf == 1);
111    C        (this reduces the number of lines of code to modify)
112    C     2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
113    C        (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
114    C       => 2 factors cancel in elliptic eq. for Phi_s ,
115    C       but 1rst factor(a) remains in RHS cg2d_b.
116    
117    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
118    C     instead of simply etaN ; This can speed-up the solver convergence in
119    C     the case where |Global_mean_PmE| is large.
120          putPmEinXvector = .FALSE.
121    c     putPmEinXvector = useRealFreshWaterFlux
122    
123  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
124          sumEmP = 0.
125        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
126         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
127          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
128           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
129  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
130            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
131  #endif  #endif
132            cg2d_x(i,j,bi,bj) = etaN(i,j,bi,bj)            cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
133            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
 #ifdef USE_NATURAL_BCS  
      &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*  
      &       EmPmR(I,J,bi,bj)/deltaTMom  
 #endif  
134           ENDDO           ENDDO
135          ENDDO          ENDDO
136            IF (useRealFreshWaterFlux) THEN
137             tmpFac = freeSurfFac*convertEmP2rUnit
138             IF (exactConserv)
139         &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
140             DO j=1,sNy
141              DO i=1,sNx
142               cg2d_b(i,j,bi,bj) =
143         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
144              ENDDO
145             ENDDO
146            ENDIF
147            IF ( putPmEinXvector ) THEN
148             tileEmP = 0.
149             DO j=1,sNy
150              DO i=1,sNx
151                tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
152         &                                       *maskH(i,j,bi,bj)
153              ENDDO
154             ENDDO
155             sumEmP = sumEmP + tileEmP
156            ENDIF
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159          IF ( putPmEinXvector ) THEN
160            _GLOBAL_SUM_R8( sumEmP, myThid )
161          ENDIF
162    
163        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
164         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
165            IF ( putPmEinXvector ) THEN
166              tmpFac = 0.
167              IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
168         &                          *convertEmP2rUnit*sumEmP/globalArea
169              DO j=1,sNy
170               DO i=1,sNx
171                cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
172         &                        - tmpFac*Bo_surf(i,j,bi,bj)
173               ENDDO
174              ENDDO
175            ENDIF
176    C- RHS: similar to the divergence of the vertically integrated mass transport:
177    C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
178          DO K=Nr,1,-1          DO K=Nr,1,-1
179           DO j=1,sNy+1           DO j=1,sNy+1
180            DO i=1,sNx+1            DO i=1,sNx+1
181             uf(i,j) = _dyG(i,j,bi,bj)             uf(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
182       &      *drF(k)*_hFacW(i,j,k,bi,bj)       &               *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
183             vf(i,j) = _dxG(i,j,bi,bj)             vf(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
184       &      *drF(k)*_hFacS(i,j,k,bi,bj)       &               *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
185            ENDDO            ENDDO
186           ENDDO           ENDDO
187           CALL CALC_DIV_GHAT(           CALL CALC_DIV_GHAT(
# Line 84  C--   Add source term arising from w=d/d Line 197  C--   Add source term arising from w=d/d
197        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
198         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
199  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
200          DO j=1,sNy          IF ( use3Dsolver .AND. zeroPsNH ) THEN
201           DO i=1,sNx           DO j=1,sNy
202            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            DO i=1,sNx
203       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(             ks = ksurfC(i,j,bi,bj)
204       &          -cg2d_x(I,J,bi,bj)             IF ( ks.LE.Nr ) THEN
205       &          -cg3d_x(I,J,1,bi,bj)              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
206       &        )/deltaTMom/deltaTMom       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
207            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)       &         /deltaTMom/deltaTfreesurf
208       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &         * etaH(i,j,bi,bj)
209       &         -cg2d_x(I,J,bi,bj)              cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
210       &         -cg3d_x(I,J,1,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
211       &       )/deltaTMom/deltaTMom       &         /deltaTMom/deltaTfreesurf
212         &         * etaH(i,j,bi,bj)
213               ENDIF
214              ENDDO
215           ENDDO           ENDDO
216          ENDDO          ELSEIF ( use3Dsolver ) THEN
217             DO j=1,sNy
218              DO i=1,sNx
219               ks = ksurfC(i,j,bi,bj)
220               IF ( ks.LE.Nr ) THEN
221                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
222         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
223         &         /deltaTMom/deltaTfreesurf
224         &         *( etaN(i,j,bi,bj)
225         &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
226                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
227         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
228         &         /deltaTMom/deltaTfreesurf
229         &         *( etaN(i,j,bi,bj)
230         &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
231               ENDIF
232              ENDDO
233             ENDDO
234            ELSEIF ( exactConserv ) THEN
235  #else  #else
236          DO j=1,sNy          IF ( exactConserv ) THEN
237           DO i=1,sNx  #endif /* ALLOW_NONHYDROSTATIC */
238            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)           DO j=1,sNy
239       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(            DO i=1,sNx
240       &          -cg2d_x(I,J,bi,bj)             ks = ksurfC(i,j,bi,bj)
241       &        )/deltaTMom/deltaTMom             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
242         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
243         &         /deltaTMom/deltaTfreesurf
244         &         * etaH(i,j,bi,bj)
245              ENDDO
246           ENDDO           ENDDO
247          ENDDO          ELSE
248  #endif           DO j=1,sNy
249              DO i=1,sNx
250               ks = ksurfC(i,j,bi,bj)
251               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
252         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
253         &         /deltaTMom/deltaTfreesurf
254         &         * etaN(i,j,bi,bj)
255              ENDDO
256             ENDDO
257            ENDIF
258    
259  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
260          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 115  C--   Add source term arising from w=d/d Line 262  C--   Add source term arising from w=d/d
262  C Northern boundary  C Northern boundary
263            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(I,bi,bj).NE.0) THEN
264             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
265               cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
266            ENDIF            ENDIF
267  C Southern boundary  C Southern boundary
268            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(I,bi,bj).NE.0) THEN
269             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
270               cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
271            ENDIF            ENDIF
272           ENDDO           ENDDO
273           DO j=1,sNy           DO j=1,sNy
274  C Eastern boundary  C Eastern boundary
275            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
276             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
277               cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
278            ENDIF            ENDIF
279  C Western boundary  C Western boundary
280            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
281             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
282               cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
283            ENDIF            ENDIF
284           ENDDO           ENDDO
285          ENDIF          ENDIF
286  #endif  #endif /* ALLOW_OBCS */
287    C-    end bi,bj loops
288         ENDDO         ENDDO
289        ENDDO        ENDDO
290    
291    #ifdef ALLOW_DEBUG
292          IF ( debugLevel .GE. debLevB ) THEN
293           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
294         &                        myThid)
295          ENDIF
296    #endif
297    
298  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
299  C--   gradient solver.  C--   gradient solver.
300  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
301          firstResidual=0.
302          lastResidual=0.
303          numIters=cg2dMaxIters
304    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
305    #ifdef ALLOW_CG2D_NSA
306    C--   Call the not-self-adjoint version of cg2d
307          CALL CG2D_NSA(
308         U           cg2d_b,
309         U           cg2d_x,
310         O           firstResidual,
311         O           lastResidual,
312         U           numIters,
313         I           myThid )
314    #else /* not ALLOW_CG2D_NSA = default */
315        CALL CG2D(        CALL CG2D(
316       I           cg2d_b,       U           cg2d_b,
317       U           cg2d_x,       U           cg2d_x,
318         O           firstResidual,
319         O           lastResidual,
320         U           numIters,
321       I           myThid )       I           myThid )
322    #endif /* ALLOW_CG2D_NSA */
323        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
324    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
325    
326    #ifdef ALLOW_DEBUG
327          IF ( debugLevel .GE. debLevB ) THEN
328           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
329         &                        myThid)
330          ENDIF
331    #endif
332    
333    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
334          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
335         &   ) THEN
336           IF ( debugLevel .GE. debLevA ) THEN
337            _BEGIN_MASTER( myThid )
338            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
339            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
340            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
341            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
342            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
343            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
344            _END_MASTER( myThid )
345           ENDIF
346          ENDIF
347    
348  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
349        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
350         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
351          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
352           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
353            etaN(i,j,bi,bj) = cg2d_x(i,j,bi,bj)            etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
354           ENDDO           ENDDO
355          ENDDO          ENDDO
356         ENDDO         ENDDO
357        ENDDO        ENDDO
358    
359  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
360        IF ( nonHydrostatic ) THEN        IF ( use3Dsolver ) THEN
361    
362  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 ).
363  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
# Line 167  C     see CG3D.h for the interface to th Line 365  C     see CG3D.h for the interface to th
365          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
366           DO j=1,sNy+1           DO j=1,sNy+1
367            DO i=1,sNx+1            DO i=1,sNx+1
368             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
369       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
370             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
371       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
372            ENDDO            ENDDO
373           ENDDO           ENDDO
# Line 197  C Western boundary Line 395  C Western boundary
395            ENDIF            ENDIF
396            ENDDO            ENDDO
397           ENDIF           ENDIF
398  #endif  #endif /* ALLOW_OBCS */
399    
400             IF ( usingZCoords ) THEN
401    C-       Z coordinate: assume surface @ level k=1
402               tmpFac = freeSurfFac*deepFac2F(1)
403             ELSE
404    C-       Other than Z coordinate: no assumption on surface level index
405               tmpFac = 0.
406               DO j=1,sNy
407                DO i=1,sNx
408                  ks = ksurfC(i,j,bi,bj)
409                  IF ( ks.LE.Nr ) THEN
410                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
411         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
412         &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom
413                  ENDIF
414                ENDDO
415               ENDDO
416             ENDIF
417           K=1           K=1
418             kp1 = MIN(k+1,Nr)
419             wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
420             IF (k.GE.Nr) wFacKp = 0.
421           DO j=1,sNy           DO j=1,sNy
422            DO i=1,sNx            DO i=1,sNx
423             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
424       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(K)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
425       &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)       &       -drF(K)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
426       &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)       &       +drF(K)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
427       &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )       &       -drF(K)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
428       &       +(       &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
429       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*wFacKp
430       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
      &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(  
      &          +cg2d_x(I,J,bi,bj)  
      &        )/deltaTMom/deltaTMom  
431            ENDDO            ENDDO
432           ENDDO           ENDDO
433           DO K=2,Nr-1           DO K=2,Nr
434              kp1 = MIN(k+1,Nr)
435    C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
436    C        both appear in wVel term, but at 2 different levels
437              wFacKm = deepFac2F( k )*rhoFacF( k )
438              wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
439              IF (k.GE.Nr) wFacKp = 0.
440            DO j=1,sNy            DO j=1,sNy
441             DO i=1,sNx             DO i=1,sNx
442              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
443       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(K)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
444       &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)       &       -drF(K)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
445       &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)       &       +drF(K)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
446       &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )       &       -drF(K)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
447       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
448       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*wFacKp
449       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
450    
451             ENDDO             ENDDO
452            ENDDO            ENDDO
453           ENDDO           ENDDO
          K=Nr  
          DO j=1,sNy  
           DO i=1,sNx  
             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)  
      &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)  
      &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)  
      &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)  
      &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )  
      &       +( wVel(i,j,k  ,bi,bj)  
      &        )*_rA(i,j,bi,bj)/deltaTmom  
   
           ENDDO  
          ENDDO  
454    
455  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
456           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 269  C Western boundary Line 477  C Western boundary
477            ENDDO            ENDDO
478            ENDDO            ENDDO
479           ENDIF           ENDIF
480  #endif  #endif /* ALLOW_OBCS */
481    C-    end bi,bj loops
482            ENDDO
483           ENDDO
484    
485          ENDDO ! bi        firstResidual=0.
486         ENDDO ! bj        lastResidual=0.
487          numIters=cg3dMaxIters
488          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
489          CALL CG3D(
490         U           cg3d_b,
491         U           phi_nh,
492         O           firstResidual,
493         O           lastResidual,
494         U           numIters,
495         I           myThid )
496          _EXCH_XYZ_R8(phi_nh, myThid )
497          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
498    
499         CALL CG3D( myThid )        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
500         _EXCH_XYZ_R8(cg3d_x, myThid )       &   ) THEN
501           IF ( debugLevel .GE. debLevA ) THEN
502            _BEGIN_MASTER( myThid )
503            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
504            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
505            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
506            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
507            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
508            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
509            _END_MASTER( myThid )
510           ENDIF
511          ENDIF
512    
513    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
514          IF ( zeroPsNH ) THEN
515           DO bj=myByLo(myThid),myByHi(myThid)
516            DO bi=myBxLo(myThid),myBxHi(myThid)
517    
518             IF ( usingZCoords ) THEN
519    C-       Z coordinate: assume surface @ level k=1
520              DO k=2,Nr
521               DO j=1-OLy,sNy+OLy
522                DO i=1-OLx,sNx+OLx
523                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
524         &                           - phi_nh(i,j,1,bi,bj)
525                ENDDO
526               ENDDO
527              ENDDO
528              DO j=1-OLy,sNy+OLy
529               DO i=1-OLx,sNx+OLx
530                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
531         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
532                 phi_nh(i,j,1,bi,bj) = 0.
533               ENDDO
534              ENDDO
535             ELSE
536    C-       Other than Z coordinate: no assumption on surface level index
537              DO j=1-OLy,sNy+OLy
538               DO i=1-OLx,sNx+OLx
539                ks = ksurfC(i,j,bi,bj)
540                IF ( ks.LE.Nr ) THEN
541                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
542         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
543                 DO k=Nr,1,-1
544                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
545         &                            - phi_nh(i,j,ks,bi,bj)
546                 ENDDO
547                ENDIF
548               ENDDO
549              ENDDO
550             ENDIF
551    
552            ENDDO
553           ENDDO
554        ENDIF        ENDIF
 #endif  
555    
556          ENDIF
557    #endif /* ALLOW_NONHYDROSTATIC */
558    
559    #ifdef TIME_PER_TIMESTEP_SFP
560    CCE107 Time per timestep information
561          _BEGIN_MASTER( myThid )
562          CALL TIMER_GET_TIME( utnew, stnew, wtnew )
563    C Only output timing information after the 1st timestep
564          IF ( wtold .NE. 0.0D0 ) THEN
565            WRITE(msgBuf,'(A34,3F10.6)')
566         $        'User, system and wallclock time:', utnew - utold,
567         $        stnew - stold, wtnew - wtold
568             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
569          ENDIF
570          utold = utnew
571          stold = stnew
572          wtold = wtnew
573          _END_MASTER( myThid )
574    #endif
575    #ifdef USE_PAPI_FLOPS_SFP
576    CCE107 PAPI summary performance
577          _BEGIN_MASTER( myThid )
578    #ifdef USE_FLIPS
579          call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
580    #else
581          call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
582    #endif
583          WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
584         $     'Mflop/s during this timestep:', mflops, ' ', mflops
585         $     *proc_time/(real_time + 1E-36)
586          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
587    #ifdef PAPI_VERSION
588          call PAPIF_ipc(real_time, proc_time, instr, ipc, check)
589          WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
590         $     'IPC during this timestep:', ipc, ' ', ipc*proc_time
591         $     /(real_time + 1E-36)
592          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
593    #endif
594          _END_MASTER( myThid )
595    #else
596    #ifdef USE_PCL_FLOPS_SFP
597    CCE107 PCL summary performance
598          _BEGIN_MASTER( myThid )
599          PCLstop(descr, i_result, fp_result, nevents)
600          do ipcl = 1, nevents
601             WRITE(msgBuf,'(A22,A26,F10.6)'),
602         $        pcl_counter_name(pcl_counter_list(ipcl)),
603         $        'during this timestep:', fp_results(ipcl)
604             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
605          enddo
606          PCLstart(descr, pcl_counter_list, nevents, flags)
607          _END_MASTER( myThid )
608    #endif
609    #endif
610        RETURN        RETURN
611        END        END
612    
613    #ifdef TIME_PER_TIMESTEP_SFP
614    CCE107 Initialization of common block for per timestep timing
615          BLOCK DATA settimers
616    C     !TIMING VARIABLES
617    C     == Timing variables ==
618          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
619          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
620          DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
621          END
622    #endif
623    #ifdef USE_PAPI_FLOPS_SFP
624    CCE107 Initialization of common block for PAPI summary performance
625          BLOCK DATA setpapis
626          INTEGER*8 flpops, instr
627          REAL real_time, proc_time, mflops, ipc
628          COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
629          DATA flpops, instr, real_time, proc_time, mflops, ipc /2*0,4*0.E0/
630          END
631    #endif
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.58
  ViewVC Help
Powered by ViewVC 1.1.22