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

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

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


Revision 1.66 - (show annotations) (download)
Tue Apr 28 18:01:14 2009 UTC (15 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint61n, checkpoint61m
Changes since 1.65: +3 -3 lines
change macros (EXCH & GLOBAL_SUM/MAX) sufix _R4/_R8 to _RS/_RL
 when applied to _RS/_RL variable

1 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.65 2009/04/22 01:30:00 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 CBOP
8 C !ROUTINE: SOLVE_FOR_PRESSURE
9 C !INTERFACE:
10 SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
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
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "SURFACE.h"
28 #include "FFIELDS.h"
29 #include "DYNVARS.h"
30 #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 ==
48 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
54
55 C !LOCAL VARIABLES:
56 C == Local variables ==
57 INTEGER i,j,k,bi,bj
58 _RL firstResidual,lastResidual
59 _RL tmpFac
60 _RL sumEmP, tileEmP(nSx,nSy)
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 C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
93 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)
101 DO bi=myBxLo(myThid),myBxHi(myThid)
102 DO j=1-OLy,sNy+OLy
103 DO i=1-OLx,sNx+OLx
104 #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
110 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(bi,bj) = 0.
124 DO j=1,sNy
125 DO i=1,sNx
126 tileEmP(bi,bj) = tileEmP(bi,bj)
127 & + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
128 & *maskH(i,j,bi,bj)
129 ENDDO
130 ENDDO
131 ENDIF
132 ENDDO
133 ENDDO
134 IF ( putPmEinXvector ) THEN
135 CALL GLOBAL_SUM_TILE_RL( tileEmP, 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
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
269 C-- gradient solver.
270 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(
286 U cg2d_b,
287 U cg2d_x,
288 O firstResidual,
289 O lastResidual,
290 U numIters,
291 I myThid )
292 #endif /* ALLOW_CG2D_NSA */
293 _EXCH_XY_RL(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 ).
333 C see CG3D.h for the interface to this routine.
334 DO bj=myByLo(myThid),myByHi(myThid)
335 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_RL(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
534 END

  ViewVC Help
Powered by ViewVC 1.1.22