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

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

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

revision 1.22 by jmc, Tue Nov 8 02:14:10 2005 UTC revision 1.28 by jmc, Fri Jul 7 20:09:37 2006 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
3  C $Name$  C $Name$
4    
 #include "PACKAGES_CONFIG.h"  
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE CALC_GW(            SUBROUTINE CALC_GW(
11       I           myTime, myIter, myThid )       I               bi, bj, KappaRU, KappaRV,
12         I               myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | S/R CALC_GW                                                C     | S/R CALC_GW
16  C     | o Calculate vert. velocity tendency terms ( NH, QH only )  C     | o Calculate vert. velocity tendency terms ( NH, QH only )
17  C     *==========================================================*  C     *==========================================================*
18  C     | In NH and QH, the vertical momentum tendency must be  C     | In NH and QH, the vertical momentum tendency must be
19  C     | calculated explicitly and included as a source term  C     | calculated explicitly and included as a source term
20  C     | for a 3d pressure eqn. Calculate that term here.  C     | for a 3d pressure eqn. Calculate that term here.
21  C     | This routine is not used in HYD calculations.  C     | This routine is not used in HYD calculations.
22  C     *==========================================================*  C     *==========================================================*
23  C     \ev  C     \ev
24    
25  C     !USES:  C     !USES:
26        IMPLICIT NONE        IMPLICIT NONE
# Line 34  C     == Global variables == Line 34  C     == Global variables ==
34    
35  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
36  C     == Routine arguments ==  C     == Routine arguments ==
37  C     myTime :: Current time in simulation  C     bi,bj   :: current tile indices
38  C     myIter :: Current iteration number in simulation  C     KappaRU :: vertical viscosity at U points
39  C     myThid :: Thread number for this instance of the routine.  C     KappaRV :: vertical viscosity at V points
40    C     myTime  :: Current time in simulation
41    C     myIter  :: Current iteration number in simulation
42    C     myThid  :: Thread number for this instance of the routine.
43          INTEGER bi,bj
44          _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
45          _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
46        _RL     myTime        _RL     myTime
47        INTEGER myIter        INTEGER myIter
48        INTEGER myThid        INTEGER myThid
# Line 45  C     myThid :: Thread number for this i Line 51  C     myThid :: Thread number for this i
51    
52  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
53  C     == Local variables ==  C     == Local variables ==
 C     bi, bj,      :: Loop counters  
54  C     iMin, iMax,  C     iMin, iMax,
55  C     jMin, jMax  C     jMin, jMax
56  C     flx_NS       :: Temp. used for fVol meridional terms.  C     flx_NS       :: Temp. used for fVol meridional terms.
57  C     flx_EW       :: Temp. used for fVol zonal terms.  C     flx_EW       :: Temp. used for fVol zonal terms.
58  C     flx_Up       :: Temp. used for fVol vertical terms.  C     flx_Up       :: Temp. used for fVol vertical terms.
59  C     flx_Dn       :: Temp. used for fVol vertical terms.  C     flx_Dn       :: Temp. used for fVol vertical terms.
60        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
61        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65        _RL    fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66        _RL    fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67        _RL    del2w(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68  C     i,j,k - Loop counters  C     i,j,k - Loop counters
69        INTEGER i,j,k, kP1        INTEGER i,j,k, kp1
70        _RL  wOverride        _RL  wOverride
71        _RS  hFacWtmp        _RS  hFacWtmp
72        _RS  hFacStmp        _RS  hFacStmp
73        _RS  hFacCtmp        _RS  hFacCtmp
74        _RS  recip_hFacCtmp        _RS  recip_hFacCtmp
75        _RL ab15,ab05        _RL  slipSideFac
76        _RL slipSideFac        _RL  tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
77        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL  viscLoc
   
78        _RL  Half        _RL  Half
79        PARAMETER(Half=0.5D0)        PARAMETER(Half=0.5D0)
80  CEOP  CEOP
81    
82    C Catch barotropic mode
83          IF ( Nr .LT. 2 ) RETURN
84    
85        iMin = 1        iMin = 1
86        iMax = sNx        iMax = sNx
87        jMin = 1        jMin = 1
88        jMax = sNy        jMax = sNy
89    
 C     Adams-Bashforth timestepping weights  
       ab15 =  1.5 _d 0 + abeps  
       ab05 = -0.5 _d 0 - abeps  
   
90  C     Lateral friction (no-slip, free slip, or half slip):  C     Lateral friction (no-slip, free slip, or half slip):
91        IF ( no_slip_sides ) THEN        IF ( no_slip_sides ) THEN
92          slipSideFac = -1. _d 0          slipSideFac = -1. _d 0
93        ELSE        ELSE
94          slipSideFac =  1. _d 0          slipSideFac =  1. _d 0
95        ENDIF        ENDIF
96  CML   half slip was used before ; keep the line for now, but half slip is  CML   half slip was used before ; keep the line for now, but half slip is
97  CML   not used anywhere in the code as far as I can see.  CML   not used anywhere in the code as far as I can see.
98  C        slipSideFac = 0. _d 0  C        slipSideFac = 0. _d 0
99    
100        DO bj=myByLo(myThid),myByHi(myThid)  C--   Initialise gW to zero
101         DO bi=myBxLo(myThid),myBxHi(myThid)        DO k=1,Nr
102          DO K=1,Nr          DO j=1-OLy,sNy+OLy
103           DO j=1-OLy,sNy+OLy           DO i=1-OLx,sNx+OLx
           DO i=1-OLx,sNx+OLx  
            gwNm1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)  
104             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
           ENDDO  
105           ENDDO           ENDDO
106          ENDDO          ENDDO
        ENDDO  
107        ENDDO        ENDDO
108    
109  C Catch barotropic mode  C--   Boundaries condition at top
110        IF ( Nr .LT. 2 ) RETURN        DO j=jMin,jMax
111           DO i=iMin,iMax
112  C For each tile           flx_Up(i,j)=0.
113        DO bj=myByLo(myThid),myByHi(myThid)         ENDDO
114         DO bi=myBxLo(myThid),myBxHi(myThid)        ENDDO
   
 C Boundaries condition at top  
         DO J=jMin,jMax  
          DO I=iMin,iMax  
           Flx_Dn(I,J,bi,bj)=0.  
          ENDDO  
         ENDDO  
115    
116  C Sweep down column  C---  Sweep down column
117          DO K=2,Nr        DO k=2,Nr
118           Kp1=K+1          kp1=k+1
119           wOverRide=1.          wOverRide=1.
120           if (K.EQ.Nr) then          IF (k.EQ.Nr) THEN
121            Kp1=Nr            kp1=Nr
122            wOverRide=0.            wOverRide=0.
123           endif          ENDIF
124  C     horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
125           IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
126  C     calculate the horizontal Laplacian of vertical flow  C-    calculate the horizontal Laplacian of vertical flow
127  C     Zonal flux d/dx W  C     Zonal flux d/dx W
128            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
129             fZon(1-Olx,j)=0.             fZon(1-Olx,j)=0.
130             DO i=1-Olx+1,sNx+Olx             DO i=1-Olx+1,sNx+Olx
131    C- Problem here: drF(k)*_hFacC & fZon are not at the same Horiz.&Vert. location
132              fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)              fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
133       &           *_dyG(i,j,bi,bj)       &           *_dyG(i,j,bi,bj)
134       &           *_recip_dxC(i,j,bi,bj)       &           *_recip_dxC(i,j,bi,bj)
135       &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))       &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
136  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
137       &           *sqcosFacU(J,bi,bj)       &           *sqcosFacU(j,bi,bj)
138  #endif  #endif
139             ENDDO             ENDDO
140            ENDDO            ENDDO
141  C     Meridional flux d/dy W  C     Meridional flux d/dy W
142            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
143             fMer(I,1-Oly)=0.             fMer(i,1-Oly)=0.
144            ENDDO            ENDDO
145            DO j=1-Oly+1,sNy+Oly            DO j=1-Oly+1,sNy+Oly
146             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
147    C- Problem here: drF(k)*_hFacC & fMer are not at the same Horiz.&Vert. location
148              fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)              fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
149       &           *_dxG(i,j,bi,bj)       &           *_dxG(i,j,bi,bj)
150       &           *_recip_dyC(i,j,bi,bj)       &           *_recip_dyC(i,j,bi,bj)
# Line 162  C     Meridional flux d/dy W Line 156  C     Meridional flux d/dy W
156  #endif  #endif
157             ENDDO             ENDDO
158            ENDDO            ENDDO
159              
160  C     del^2 W  C     del^2 W
161  C     Difference of zonal fluxes ...  C     Difference of zonal fluxes ...
162            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
# Line 177  C     ... add difference of meridional f Line 171  C     ... add difference of meridional f
171             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
172  C     First compute the fraction of open water for the w-control volume  C     First compute the fraction of open water for the w-control volume
173  C     at the southern face  C     at the southern face
174              hFacCtmp=max(hFacC(I,J,K-1,bi,bj)-Half,0. _d 0)              hFacCtmp=max( _hFacC(i,j,k-1,bi,bj)-Half,0. _d 0 )
175       &           +   min(hFacC(I,J,K  ,bi,bj),Half)       &           +   min( _hFacC(i,j,k  ,bi,bj)     ,Half )
176              IF (hFacCtmp .GT. 0.) THEN              IF (hFacCtmp .GT. 0.) THEN
177               recip_hFacCtmp = 1./hFacCtmp               recip_hFacCtmp = 1./hFacCtmp
178              ELSE              ELSE
# Line 197  CML ************************************ Line 191  CML ************************************
191  CML   No-slip Boundary conditions for bi-harmonic dissipation  CML   No-slip Boundary conditions for bi-harmonic dissipation
192  CML   need to be implemented here!  CML   need to be implemented here!
193  CML ************************************************************  CML ************************************************************
194           ELSE          ELSE
195  C-    Initialize del2w to zero:  C-    Initialize del2w to zero:
196            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
197             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
198              del2w(i,j) = 0. _d 0              del2w(i,j) = 0. _d 0
199             ENDDO             ENDDO
200            ENDDO            ENDDO
201           ENDIF          ENDIF
202    
203  C Flux on Southern face  C Flux on Southern face
204           DO J=jMin,jMax+1          DO j=jMin,jMax+1
205            DO I=iMin,iMax           DO i=iMin,iMax
206  C     First compute the fraction of open water for the w-control volume  C     First compute the fraction of open water for the w-control volume
207  C     at the southern face  C     at the southern face
208             hFacStmp=max(hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)             hFacStmp=max(_hFacS(i,j,k-1,bi,bj)-Half,0. _d 0)
209       &         +    min(hFacS(I,J,K  ,bi,bj),Half)       &         +    min(_hFacS(i,j,k  ,bi,bj),Half)
210             tmp_VbarZ=Half*(             tmp_VbarZ=Half*(
211       &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)       &          _hFacS(i,j,k-1,bi,bj)*vVel( i ,j,k-1,bi,bj)
212       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &         +_hFacS(i,j, k ,bi,bj)*vVel( i ,j, k ,bi,bj))
213             Flx_NS(I,J,bi,bj)=             flx_NS(i,j)=
214       &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))       &     tmp_VbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))
215       &    -viscAhW*_recip_dyC(I,J,bi,bj)         &    -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*Half
216       &       *(hFacStmp*(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))       &       *_recip_dyC(i,j,bi,bj)
217         &       *(hFacStmp*(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
218    C- Problem here: No-slip bc CANNOT be written in term of a flux
219       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
220       &         *wVel(I,J,K,bi,bj))       &         *wVel(i,j,k,bi,bj))
221       &    +viscA4W*_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))       &    +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*Half
222         &      *_recip_dyC(i,j,bi,bj)*(del2w(i,j)-del2w(i,j-1))
223  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
224  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
225       &    *sqCosFacV(j,bi,bj)       &    *sqCosFacV(j,bi,bj)
# Line 235  C     fraction of the w control volume a Line 232  C     fraction of the w control volume a
232  C     the control volume. A more compact but less intelligible version  C     the control volume. A more compact but less intelligible version
233  C     of the last three lines is:  C     of the last three lines is:
234  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))
235  CML     &       *wVel(I,J,K,bi,bi) + hFacStmp*wVel(I,J-1,K,bi,bj) )  CML     &       *wVel(i,j,k,bi,bi) + hFacStmp*wVel(i,j-1,k,bi,bj) )
           ENDDO  
236           ENDDO           ENDDO
237            ENDDO
238  C Flux on Western face  C Flux on Western face
239           DO J=jMin,jMax          DO j=jMin,jMax
240            DO I=iMin,iMax+1           DO i=iMin,iMax+1
241  C     First compute the fraction of open water for the w-control volume  C     First compute the fraction of open water for the w-control volume
242  C     at the western face  C     at the western face
243             hFacWtmp=max(hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)             hFacWtmp=max(_hFacW(i,j,k-1,bi,bj)-Half,0. _d 0)
244       &         +    min(hFacW(I,J,K  ,bi,bj),Half)       &         +    min(_hFacW(i,j,k  ,bi,bj),Half)
245             tmp_UbarZ=Half*(             tmp_UbarZ=Half*(
246       &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)       &         _hFacW(i,j,k-1,bi,bj)*uVel( i ,j,k-1,bi,bj)
247       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &        +_hFacW(i,j, k ,bi,bj)*uVel( i ,j, k ,bi,bj))
248             Flx_EW(I,J,bi,bj)=             flx_EW(i,j)=
249       &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))       &     tmp_UbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))
250       &    -viscAhW*_recip_dxC(I,J,bi,bj)       &    -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*Half
251       &      *(hFacWtmp*(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))       &      *_recip_dxC(i,j,bi,bj)
252         &      *(hFacWtmp*(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
253    C- Problem here: No-slip bc CANNOT be written in term of a flux
254       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
255       &        *wVel(I,J,K,bi,bj) )       &        *wVel(i,j,k,bi,bj) )
256       &    +viscA4W*_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))       &    +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*Half
257         &      *_recip_dxC(i,j,bi,bj)*(del2w(i,j)-del2w(i-1,j))
258  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
259       &                *sqCosFacU(j,bi,bj)       &                *sqCosFacU(j,bi,bj)
260  #else  #else
261       &                *CosFacU(j,bi,bj)       &                *CosFacU(j,bi,bj)
262  #endif  #endif
263  C     The last term is the weighted average of the viscous stress at the open  C     The last term is the weighted average of the viscous stress at the open
# Line 265  C     fraction of the w control volume a Line 265  C     fraction of the w control volume a
265  C     the control volume. A more compact but less intelligible version  C     the control volume. A more compact but less intelligible version
266  C     of the last three lines is:  C     of the last three lines is:
267  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))
268  CML     &       *wVel(I,J,K,bi,bi) + hFacWtmp*wVel(I-1,J,K,bi,bj) )  CML     &       *wVel(i,j,k,bi,bi) + hFacWtmp*wVel(i-1,j,k,bi,bj) )
           ENDDO  
269           ENDDO           ENDDO
270            ENDDO
271  C Flux on Lower face  C Flux on Lower face
272           DO J=jMin,jMax          DO j=jMin,jMax
273            DO I=iMin,iMax           DO i=iMin,iMax
274             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)  C Interpolate vert viscosity to W points
275             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)             viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
276       &         +wOverRide*wVel(I,J,Kp1,bi,bj))       &                +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
277             Flx_Dn(I,J,bi,bj)=       &                +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
278       &     tmp_WbarZ*tmp_WbarZ       &                +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
279       &    -viscAr*recip_drF(K)       &               )*0.125 _d 0
280       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )             tmp_WbarZ = Half*( wVel(i,j, k ,bi,bj)
281            ENDDO       &                       +wVel(i,j,kp1,bi,bj)*wOverRide )
282               flx_Dn(i,j)=
283         &                  tmp_WbarZ*tmp_WbarZ
284         &                 -viscLoc*recip_drF(k)
285         &                         *( wVel(i,j, k ,bi,bj)
286         &                           -wVel(i,j,kp1,bi,bj)*wOverRide )
287           ENDDO           ENDDO
288  C        Divergence of fluxes          ENDDO
289           DO J=jMin,jMax  C       Divergence of fluxes
290            DO I=iMin,iMax          DO j=jMin,jMax
291             gW(I,J,K,bi,bj) = 0.           DO i=iMin,iMax
292               gW(i,j,k,bi,bj) = 0.
293       &      -(       &      -(
294       &        +_recip_dxF(I,J,bi,bj)*(       &        +_recip_dxF(i,j,bi,bj)*(
295       &              Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )       &              flx_EW(i+1,j)-flx_EW(i,j) )
296       &        +_recip_dyF(I,J,bi,bj)*(       &        +_recip_dyF(i,j,bi,bj)*(
297       &              Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )       &              flx_NS(i,j+1)-flx_NS(i,j) )
298       &        +recip_drC(K)         *(       &        +recip_drC(k)         *(
299       &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )       &              flx_Up(i,j)  -flx_Dn(i,j) )
300       &       )       &       )
301  caja    * recip_hFacU(I,J,K,bi,bj)  caja    * recip_hFacU(i,j,k,bi,bj)
302  caja   NOTE:  This should be included    caja   NOTE:  This should be included
303  caja   but we need an hFacUW (above U points)  caja   but we need an hFacUW (above U points)
304  caja           and an hFacUS (above V points) too...  caja           and an hFacUS (above V points) too...
305            ENDDO  C--        prepare for next level (k+1)
306               flx_Up(i,j)=flx_Dn(i,j)
307           ENDDO           ENDDO
308          ENDDO          ENDDO
        ENDDO  
       ENDDO  
309    
310    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
311    
312        DO bj=myByLo(myThid),myByHi(myThid)  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
313         DO bi=myBxLo(myThid),myBxHi(myThid)  C     and save gW_[n] into gwNm1 for the next time step.
314          DO K=2,Nr  c#ifdef ALLOW_ADAMSBASHFORTH_3
315           DO j=jMin,jMax  c       CALL ADAMS_BASHFORTH3(
316            DO i=iMin,iMax  c    I                         bi, bj, k,
317             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)  c    U                         gW, gwNm,
318       &     +deltatMom*nh_Am2*( ab15*gW(i,j,k,bi,bj)  c    I                         momStartAB, myIter, myThid )
319       &                 +ab05*gwNm1(i,j,k,bi,bj) )  c#else /* ALLOW_ADAMSBASHFORTH_3 */
320             IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.          CALL ADAMS_BASHFORTH2(
321            ENDDO       I                         bi, bj, k,
322           ENDDO       U                         gW, gwNm1,
323          ENDDO       I                         myIter, myThid )
324         ENDDO  c#endif /* ALLOW_ADAMSBASHFORTH_3 */
       ENDDO  
325    
326  #ifdef ALLOW_OBCS  C-    end of the k loop
327        IF (useOBCS) THEN        ENDDO
 C--   This call is aesthetic: it makes the W field  
 C     consistent with the OBs but this has no algorithmic  
 C     impact. This is purely for diagnostic purposes.  
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO K=1,Nr  
           CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDIF  
 #endif /* ALLOW_OBCS */  
328    
329  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
330    
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22