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

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

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

Revision 1.51 - (hide annotations) (download)
Tue Dec 20 20:31:28 2005 UTC (18 years, 5 months ago) by jmc
Branch: MAIN
Changes since 1.50: +52 -35 lines 
make 3.D solver compatible with Free-surface at k > 1 (p-coordinate):
 compute & store in commom block solver main diagonal element.
1 jmc 1.51 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.50 2005/12/19 21:10:34 jmc Exp $
2 heimbach 1.21 C $Name: $
3 cnh 1.1
4 edhill 1.39 #include "PACKAGES_CONFIG.h"
5 adcroft 1.5 #include "CPP_OPTIONS.h"
6 cnh 1.1
7 cnh 1.27 CBOP
8     C !ROUTINE: SOLVE_FOR_PRESSURE
9     C !INTERFACE:
10 jmc 1.29 SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
11 cnh 1.27
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 adcroft 1.8 IMPLICIT NONE
22 cnh 1.4 C == Global variables
23     #include "SIZE.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26 adcroft 1.12 #include "GRID.h"
27 jmc 1.17 #include "SURFACE.h"
28 jmc 1.28 #include "FFIELDS.h"
29 jmc 1.48 #include "DYNVARS.h"
30     #include "SOLVE_FOR_PRESSURE.h"
31 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
32 adcroft 1.25 #include "SOLVE_FOR_PRESSURE3D.h"
33 jmc 1.48 #include "NH_VARS.h"
34     #endif
35     #ifdef ALLOW_CD_CODE
36     #include "CD_CODE_VARS.h"
37 adcroft 1.12 #endif
38 adcroft 1.11 #ifdef ALLOW_OBCS
39 adcroft 1.9 #include "OBCS.h"
40 adcroft 1.11 #endif
41 cnh 1.4
42 jmc 1.32 C === Functions ====
43 jmc 1.46 LOGICAL DIFFERENT_MULTIPLE
44     EXTERNAL DIFFERENT_MULTIPLE
45 jmc 1.32
46 cnh 1.27 C !INPUT/OUTPUT PARAMETERS:
47 cnh 1.1 C == Routine arguments ==
48 jmc 1.28 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 jmc 1.29 INTEGER myThid
54 cnh 1.4
55 cnh 1.27 C !LOCAL VARIABLES:
56 adcroft 1.22 C == Local variables ==
57 cnh 1.6 INTEGER i,j,k,bi,bj
58 adcroft 1.9 _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59     _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 adcroft 1.22 _RL firstResidual,lastResidual
61 jmc 1.36 _RL tmpFac
62 jmc 1.47 _RL sumEmP, tileEmP
63     LOGICAL putPmEinXvector
64 adcroft 1.19 INTEGER numIters
65 adcroft 1.25 CHARACTER*(MAX_LEN_MBUF) msgBuf
66 jmc 1.49 #ifdef ALLOW_NONHYDROSTATIC
67 jmc 1.51 INTEGER ks, kp1
68     _RL maskKp1
69 jmc 1.49 LOGICAL zeroPsNH
70     #endif
71 cnh 1.27 CEOP
72 jmc 1.17
73 ce107 1.44 #ifdef TIME_PER_TIMESTEP
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    
81 jmc 1.49 #ifdef ALLOW_NONHYDROSTATIC
82     c zeroPsNH = .FALSE.
83     zeroPsNH = exactConserv
84     #endif
85    
86 jmc 1.47 C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
87     C instead of simply etaN ; This can speed-up the solver convergence in
88     C the case where |Global_mean_PmE| is large.
89     putPmEinXvector = .FALSE.
90     c putPmEinXvector = useRealFreshWaterFlux
91    
92 jmc 1.17 C-- Save previous solution & Initialise Vector solution and source term :
93 jmc 1.47 sumEmP = 0.
94 jmc 1.17 DO bj=myByLo(myThid),myByHi(myThid)
95     DO bi=myBxLo(myThid),myBxHi(myThid)
96     DO j=1-OLy,sNy+OLy
97     DO i=1-OLx,sNx+OLx
98 edhill 1.40 #ifdef ALLOW_CD_CODE
99 jmc 1.17 etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
100 jmc 1.26 #endif
101 jmc 1.18 cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
102 jmc 1.17 cg2d_b(i,j,bi,bj) = 0.
103     ENDDO
104     ENDDO
105 jmc 1.29 IF (useRealFreshWaterFlux) THEN
106 jmc 1.36 tmpFac = freeSurfFac*convertEmP2rUnit
107 mlosch 1.35 IF (exactConserv)
108 jmc 1.36 & tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
109 jmc 1.29 DO j=1,sNy
110     DO i=1,sNx
111     cg2d_b(i,j,bi,bj) =
112     & tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
113     ENDDO
114     ENDDO
115     ENDIF
116 jmc 1.47 IF ( putPmEinXvector ) THEN
117     tileEmP = 0.
118     DO j=1,sNy
119     DO i=1,sNx
120     tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
121     & *maskH(i,j,bi,bj)
122     ENDDO
123     ENDDO
124     sumEmP = sumEmP + tileEmP
125     ENDIF
126 jmc 1.17 ENDDO
127     ENDDO
128 jmc 1.47 IF ( putPmEinXvector ) THEN
129     _GLOBAL_SUM_R8( sumEmP, myThid )
130     ENDIF
131 adcroft 1.12
132     DO bj=myByLo(myThid),myByHi(myThid)
133     DO bi=myBxLo(myThid),myBxHi(myThid)
134 jmc 1.47 IF ( putPmEinXvector ) THEN
135     tmpFac = 0.
136     IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
137     & *convertEmP2rUnit*sumEmP/globalArea
138     DO j=1,sNy
139     DO i=1,sNx
140     cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
141     & - tmpFac*Bo_surf(i,j,bi,bj)
142     ENDDO
143     ENDDO
144     ENDIF
145 adcroft 1.12 DO K=Nr,1,-1
146     DO j=1,sNy+1
147     DO i=1,sNx+1
148     uf(i,j) = _dyG(i,j,bi,bj)
149     & *drF(k)*_hFacW(i,j,k,bi,bj)
150     vf(i,j) = _dxG(i,j,bi,bj)
151     & *drF(k)*_hFacS(i,j,k,bi,bj)
152     ENDDO
153     ENDDO
154     CALL CALC_DIV_GHAT(
155     I bi,bj,1,sNx,1,sNy,K,
156     I uf,vf,
157 jmc 1.17 U cg2d_b,
158 adcroft 1.12 I myThid)
159     ENDDO
160     ENDDO
161     ENDDO
162 cnh 1.4
163 adcroft 1.12 C-- Add source term arising from w=d/dt (p_s + p_nh)
164     DO bj=myByLo(myThid),myByHi(myThid)
165     DO bi=myBxLo(myThid),myBxHi(myThid)
166 adcroft 1.13 #ifdef ALLOW_NONHYDROSTATIC
167 jmc 1.49 IF ( nonHydrostatic .AND. zeroPsNH ) THEN
168     DO j=1,sNy
169     DO i=1,sNx
170 jmc 1.51 ks = ksurfC(i,j,bi,bj)
171     IF ( ks.LE.Nr ) THEN
172     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
173 jmc 1.49 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
174     & * etaH(i,j,bi,bj)
175 jmc 1.51 cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
176 jmc 1.49 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
177     & * etaH(i,j,bi,bj)
178 jmc 1.51 ENDIF
179 jmc 1.49 ENDDO
180     ENDDO
181     ELSEIF ( nonHydrostatic ) THEN
182 jmc 1.28 DO j=1,sNy
183     DO i=1,sNx
184 jmc 1.51 ks = ksurfC(i,j,bi,bj)
185     IF ( ks.LE.Nr ) THEN
186     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
187 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
188 jmc 1.28 & *( etaN(i,j,bi,bj)
189 jmc 1.51 & +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
190     cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
191 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
192 jmc 1.28 & *( etaN(i,j,bi,bj)
193 jmc 1.51 & +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
194     ENDIF
195 jmc 1.28 ENDDO
196 adcroft 1.12 ENDDO
197 jmc 1.28 ELSEIF ( exactConserv ) THEN
198 adcroft 1.13 #else
199 jmc 1.26 IF ( exactConserv ) THEN
200 edhill 1.39 #endif /* ALLOW_NONHYDROSTATIC */
201 jmc 1.26 DO j=1,sNy
202     DO i=1,sNx
203     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
204 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
205 jmc 1.26 & * etaH(i,j,bi,bj)
206     ENDDO
207     ENDDO
208     ELSE
209     DO j=1,sNy
210     DO i=1,sNx
211     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
212 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
213 jmc 1.26 & * etaN(i,j,bi,bj)
214     ENDDO
215 adcroft 1.12 ENDDO
216 jmc 1.26 ENDIF
217 adcroft 1.12
218     #ifdef ALLOW_OBCS
219 adcroft 1.14 IF (useOBCS) THEN
220 adcroft 1.12 DO i=1,sNx
221     C Northern boundary
222     IF (OB_Jn(I,bi,bj).NE.0) THEN
223     cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
224 jmc 1.31 cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
225 adcroft 1.12 ENDIF
226     C Southern boundary
227     IF (OB_Js(I,bi,bj).NE.0) THEN
228     cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
229 jmc 1.31 cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
230 adcroft 1.12 ENDIF
231     ENDDO
232     DO j=1,sNy
233     C Eastern boundary
234     IF (OB_Ie(J,bi,bj).NE.0) THEN
235     cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
236 jmc 1.31 cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
237 adcroft 1.12 ENDIF
238     C Western boundary
239     IF (OB_Iw(J,bi,bj).NE.0) THEN
240     cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
241 jmc 1.31 cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
242 adcroft 1.12 ENDIF
243     ENDDO
244     ENDIF
245 jmc 1.49 #endif /* ALLOW_OBCS */
246     C- end bi,bj loops
247 adcroft 1.12 ENDDO
248     ENDDO
249    
250 edhill 1.42 #ifdef ALLOW_DEBUG
251 heimbach 1.38 IF ( debugLevel .GE. debLevB ) THEN
252 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
253     & myThid)
254 adcroft 1.24 ENDIF
255 adcroft 1.23 #endif
256 adcroft 1.12
257 cnh 1.1 C-- Find the surface pressure using a two-dimensional conjugate
258     C-- gradient solver.
259 adcroft 1.22 C see CG2D.h for the interface to this routine.
260     firstResidual=0.
261     lastResidual=0.
262 adcroft 1.19 numIters=cg2dMaxIters
263 jmc 1.50 c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
264 cnh 1.1 CALL CG2D(
265 adcroft 1.22 U cg2d_b,
266 cnh 1.6 U cg2d_x,
267 adcroft 1.22 O firstResidual,
268     O lastResidual,
269 adcroft 1.19 U numIters,
270 cnh 1.1 I myThid )
271 adcroft 1.19 _EXCH_XY_R8(cg2d_x, myThid )
272 jmc 1.50 c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
273 adcroft 1.23
274 edhill 1.42 #ifdef ALLOW_DEBUG
275 heimbach 1.38 IF ( debugLevel .GE. debLevB ) THEN
276 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
277     & myThid)
278 adcroft 1.24 ENDIF
279 adcroft 1.23 #endif
280 cnh 1.1
281 jmc 1.32 C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
282 jmc 1.46 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
283 jmc 1.45 & ) THEN
284 heimbach 1.38 IF ( debugLevel .GE. debLevA ) THEN
285     _BEGIN_MASTER( myThid )
286     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
287     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
288     WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
289     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
290     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
291     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
292 edhill 1.43 _END_MASTER( myThid )
293 heimbach 1.38 ENDIF
294 jmc 1.32 ENDIF
295 jmc 1.17
296     C-- Transfert the 2D-solution to "etaN" :
297     DO bj=myByLo(myThid),myByHi(myThid)
298     DO bi=myBxLo(myThid),myBxHi(myThid)
299     DO j=1-OLy,sNy+OLy
300     DO i=1-OLx,sNx+OLx
301 jmc 1.18 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
302 jmc 1.17 ENDDO
303     ENDDO
304     ENDDO
305     ENDDO
306 adcroft 1.10
307 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
308     IF ( nonHydrostatic ) THEN
309    
310     C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
311     C see CG3D.h for the interface to this routine.
312     DO bj=myByLo(myThid),myByHi(myThid)
313     DO bi=myBxLo(myThid),myBxHi(myThid)
314     DO j=1,sNy+1
315     DO i=1,sNx+1
316 jmc 1.18 uf(i,j)=-_recip_dxC(i,j,bi,bj)*
317 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
318 jmc 1.18 vf(i,j)=-_recip_dyC(i,j,bi,bj)*
319 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
320     ENDDO
321     ENDDO
322    
323 adcroft 1.12 #ifdef ALLOW_OBCS
324 adcroft 1.14 IF (useOBCS) THEN
325 adcroft 1.9 DO i=1,sNx+1
326     C Northern boundary
327     IF (OB_Jn(I,bi,bj).NE.0) THEN
328     vf(I,OB_Jn(I,bi,bj))=0.
329     ENDIF
330     C Southern boundary
331     IF (OB_Js(I,bi,bj).NE.0) THEN
332     vf(I,OB_Js(I,bi,bj)+1)=0.
333     ENDIF
334     ENDDO
335     DO j=1,sNy+1
336     C Eastern boundary
337     IF (OB_Ie(J,bi,bj).NE.0) THEN
338     uf(OB_Ie(J,bi,bj),J)=0.
339     ENDIF
340     C Western boundary
341     IF (OB_Iw(J,bi,bj).NE.0) THEN
342     uf(OB_Iw(J,bi,bj)+1,J)=0.
343     ENDIF
344     ENDDO
345     ENDIF
346 jmc 1.49 #endif /* ALLOW_OBCS */
347 adcroft 1.9
348 jmc 1.51 IF ( usingZCoords ) THEN
349     C- Z coordinate: assume surface @ level k=1
350     tmpFac = freeSurfFac
351     ELSE
352     C- Other than Z coordinate: no assumption on surface level index
353     tmpFac = 0.
354     DO j=1,sNy
355     DO i=1,sNx
356     ks = ksurfC(i,j,bi,bj)
357     IF ( ks.LE.Nr ) THEN
358     cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
359     & +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
360     & *_rA(i,j,bi,bj)/deltaTmom
361     ENDIF
362     ENDDO
363     ENDDO
364     ENDIF
365 adcroft 1.12 K=1
366 jmc 1.51 kp1 = MIN(k+1,Nr)
367     maskKp1 = 1.
368     IF (k.GE.Nr) maskKp1 = 0.
369 adcroft 1.12 DO j=1,sNy
370     DO i=1,sNx
371 jmc 1.51 cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
372     & +drF(K)*dyG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
373     & -drF(K)*dyG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
374     & +drF(K)*dxG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
375     & -drF(K)*dxG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
376     & +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
377     & -wVel(i,j,kp1,bi,bj)*maskKp1
378 adcroft 1.12 & )*_rA(i,j,bi,bj)/deltaTmom
379     ENDDO
380     ENDDO
381 jmc 1.51 DO K=2,Nr
382     kp1 = MIN(k+1,Nr)
383     maskKp1 = 1.
384     IF (k.GE.Nr) maskKp1 = 0.
385 adcroft 1.9 DO j=1,sNy
386     DO i=1,sNx
387     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
388 jmc 1.51 & +drF(K)*dyG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
389     & -drF(K)*dyG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
390     & +drF(K)*dxG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
391     & -drF(K)*dxG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
392     & +( wVel(i,j,k ,bi,bj)*maskC(i,j,k-1,bi,bj)
393     & -wVel(i,j,kp1,bi,bj)*maskKp1
394 adcroft 1.12 & )*_rA(i,j,bi,bj)/deltaTmom
395    
396 adcroft 1.9 ENDDO
397     ENDDO
398     ENDDO
399 adcroft 1.12
400     #ifdef ALLOW_OBCS
401 adcroft 1.14 IF (useOBCS) THEN
402 adcroft 1.12 DO K=1,Nr
403     DO i=1,sNx
404     C Northern boundary
405     IF (OB_Jn(I,bi,bj).NE.0) THEN
406     cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
407     ENDIF
408     C Southern boundary
409     IF (OB_Js(I,bi,bj).NE.0) THEN
410     cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
411     ENDIF
412     ENDDO
413     DO j=1,sNy
414     C Eastern boundary
415     IF (OB_Ie(J,bi,bj).NE.0) THEN
416     cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
417     ENDIF
418     C Western boundary
419     IF (OB_Iw(J,bi,bj).NE.0) THEN
420     cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
421     ENDIF
422     ENDDO
423     ENDDO
424     ENDIF
425 jmc 1.49 #endif /* ALLOW_OBCS */
426     C- end bi,bj loops
427     ENDDO
428     ENDDO
429 adcroft 1.9
430 adcroft 1.25 firstResidual=0.
431     lastResidual=0.
432 jmc 1.49 numIters=cg3dMaxIters
433 jmc 1.50 CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
434 adcroft 1.25 CALL CG3D(
435     U cg3d_b,
436     U phi_nh,
437     O firstResidual,
438     O lastResidual,
439     U numIters,
440     I myThid )
441     _EXCH_XYZ_R8(phi_nh, myThid )
442 jmc 1.50 CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
443 adcroft 1.25
444 jmc 1.46 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
445 jmc 1.45 & ) THEN
446 heimbach 1.38 IF ( debugLevel .GE. debLevA ) THEN
447     _BEGIN_MASTER( myThid )
448     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
449     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
450     WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
451     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
452     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
453     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
454 edhill 1.43 _END_MASTER( myThid )
455 heimbach 1.38 ENDIF
456 mlosch 1.37 ENDIF
457 adcroft 1.9
458 jmc 1.49 C-- Update surface pressure (account for NH-p @ surface level) and NH pressure:
459     IF ( zeroPsNH ) THEN
460     DO bj=myByLo(myThid),myByHi(myThid)
461     DO bi=myBxLo(myThid),myBxHi(myThid)
462    
463     IF ( usingZCoords ) THEN
464     C- Z coordinate: assume surface @ level k=1
465     DO k=2,Nr
466     DO j=1-OLy,sNy+OLy
467     DO i=1-OLx,sNx+OLx
468     phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
469     & - phi_nh(i,j,1,bi,bj)
470     ENDDO
471     ENDDO
472     ENDDO
473     DO j=1-OLy,sNy+OLy
474     DO i=1-OLx,sNx+OLx
475     etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
476     & *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
477     phi_nh(i,j,1,bi,bj) = 0.
478     ENDDO
479     ENDDO
480     ELSE
481     C- Other than Z coordinate: no assumption on surface level index
482     DO j=1-OLy,sNy+OLy
483     DO i=1-OLx,sNx+OLx
484     ks = ksurfC(i,j,bi,bj)
485     IF ( ks.LE.Nr ) THEN
486     etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
487     & *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
488     DO k=Nr,1,-1
489     phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
490     & - phi_nh(i,j,ks,bi,bj)
491     ENDDO
492     ENDIF
493     ENDDO
494     ENDDO
495     ENDIF
496    
497     ENDDO
498     ENDDO
499 adcroft 1.9 ENDIF
500 jmc 1.49
501     ENDIF
502     #endif /* ALLOW_NONHYDROSTATIC */
503 cnh 1.1
504 ce107 1.44 #ifdef TIME_PER_TIMESTEP
505     CCE107 Time per timestep information
506     _BEGIN_MASTER( myThid )
507     CALL TIMER_GET_TIME( utnew, stnew, wtnew )
508     C Only output timing information after the 1st timestep
509     IF ( wtold .NE. 0.0D0 ) THEN
510     WRITE(msgBuf,'(A34,3F10.6)')
511     $ 'User, system and wallclock time:', utnew - utold,
512     $ stnew - stold, wtnew - wtold
513     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
514     ENDIF
515     utold = utnew
516     stold = stnew
517     wtold = wtnew
518     _END_MASTER( myThid )
519     #endif
520    
521 cnh 1.1 RETURN
522     END
523 ce107 1.44
524     #ifdef TIME_PER_TIMESTEP
525     CCE107 Initialization of common block for per timestep timing
526     BLOCK DATA settimers
527     C !TIMING VARIABLES
528     C == Timing variables ==
529     REAL*8 utnew, utold, stnew, stold, wtnew, wtold
530     COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
531     DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
532     END
533     #endif
  ViewVC Help
Powered by ViewVC 1.1.22