/[MITgcm]/MITgcm/pkg/mom_fluxform/mom_cdscheme.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_fluxform/mom_cdscheme.F

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

revision 1.1 by adcroft, Wed Mar 28 19:51:14 2001 UTC revision 1.2 by adcroft, Tue May 29 14:01:38 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    C $Name$
3    
4    #include "CPP_OPTIONS.h"
5    
6          SUBROUTINE MOM_CDSCHEME(
7         I        bi,bj,k,phi_hyd,
8         I        myThid)
9    
10          IMPLICIT NONE
11    
12    C     == Global variables ==
13    #include "SIZE.h"
14    #include "DYNVARS.h"
15    #include "EEPARAMS.h"
16    #include "PARAMS.h"
17    #include "GRID.h"
18    #include "SURFACE.h"
19    
20    C     == Routine arguments ==
21    C     myThid - Instance number for this innvocation of CALC_MOM_RHS
22          INTEGER bi,bj,K
23          _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
24          INTEGER myThid
25    
26    #ifdef INCLUDE_CD_CODE
27    
28    C     == Local variables ==
29          _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30          _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31          _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32          INTEGER i,j,iMin,iMax,jMin,jMax
33          _RL ab15,ab05
34    
35    C     Compute ranges
36          iMin=1-Olx+1
37          iMax=sNx+Olx-1
38          jMin=1-Oly+1
39          jMax=sNy+Oly-1
40    
41    C     Adams-Bashforth weighting factors
42          ab15   =  1.5 + abEps
43          ab05   = -0.5 - abEps
44    
45    C     Pressure extrapolated forward in time
46          DO j=1-Oly,sNy+Oly
47           DO i=1-Olx,sNx+Olx
48            pf(i,j) =
49         &   ab15*(  etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
50         &  +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
51           ENDDO
52          ENDDO
53          IF (staggerTimeStep) THEN
54           DO j=jMin,jMax
55            DO i=iMin,iMax
56             pf(i,j) = pf(i,j)+phi_hyd(i,j,k)
57            ENDDO
58           ENDDO
59          ENDIF
60    
61    C--   Zonal velocity coriolis term
62    C     Note. As coded here, coriolis will not work with "thin walls"
63    C--   Coriolis with CD scheme allowed
64    C     grady(p) + gV
65          DO j=1-Oly+1,sNy+Oly
66           DO i=1-Olx,sNx+Olx
67            af(i,j) = -_maskS(i,j,k,bi,bj)
68         &            *_recip_dyC(i,j,bi,bj)
69         &            *(pf(i,j)-pf(i,j-1))
70         &            +gV(i,j,k,bi,bj)
71           ENDDO
72          ENDDO
73    C     Average to Vd point and add coriolis
74          DO j=jMin,jMax
75           DO i=iMin,iMax
76            vf(i,j) =      
77         &    0.25*( af(i  ,j)+af(i  ,j+1)
78         &          +af(i-1,j)+af(i-1,j+1)
79         &         )*_maskW(i,j,k,bi,bj)
80         &   -0.5*(_fCori(i,j,bi,bj)+
81         &         _fCori(i-1,j,bi,bj))
82         &         *uVel(i,j,k,bi,bj)
83           ENDDO
84          ENDDO
85    C     Step forward Vd
86          DO j=jMin,jMax
87           DO i=iMin,iMax
88            vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) +
89         &                        deltaTmom*vf(i,j)
90           ENDDO
91          ENDDO
92    C     Relax D grid V to C grid V
93          DO j=jMin,jMax
94           DO i=iMin,iMax
95             vVelD(i,j,k,bi,bj) = rCD*vVelD(i,j,k,bi,bj)
96         &   +(1. - rCD)*(
97         &    ab15*0.25*(
98         &                vVel(i  ,j  ,k,bi,bj)+vVel(i  ,j+1,k,bi,bj)
99         &               +vVel(i-1,j  ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj)
100         &              )*_maskW(i,j,k,bi,bj)
101         &     +
102         &    ab05*0.25*(
103         &                vNM1(i  ,j  ,k,bi,bj)+vNM1(i  ,j+1,k,bi,bj)
104         &               +vNM1(i-1,j  ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj)
105         &              )*_maskW(i,j,k,bi,bj)
106         &   )
107           ENDDO
108          ENDDO
109    C     Calculate coriolis force on U
110          DO j=jMin,jMax
111           DO i=iMin,iMax
112            guCD(i,j,k,bi,bj) =
113         &    0.5*( _fCori(i  ,j,bi,bj) +
114         &          _fCori(i-1,j,bi,bj)  )
115         &   *vVelD(i,j,k,bi,bj)*cfFacMom
116           ENDDO
117          ENDDO
118    
119    C--   Meridional velocity coriolis term
120    C     gradx(p)+gU
121          DO j=1-Oly,sNy+Oly
122           DO i=1-Olx+1,sNx+Olx
123            af(i,j) = -_maskW(i,j,k,bi,bj)
124         &            *_recip_dxC(i,j,bi,bj)*
125         &            (pf(i,j)-pf(i-1,j))
126         &            +gU(i,j,k,bi,bj)
127           ENDDO
128          ENDDO
129    C     Average to Ud point and add coriolis
130          DO j=jMin,jMax
131           DO i=iMin,iMax
132            vf(i,j) =
133         &    0.25*( af(i  ,j)+af(i  ,j-1)
134         &          +af(i+1,j)+af(i+1,j-1)
135         &         )*_maskS(i,j,k,bi,bj)
136         &   +0.5*( _fCori(i,j,bi,bj)
137         &         +_fCori(i,j-1,bi,bj))
138         &         *vVel(i,j,k,bi,bj)
139           ENDDO
140          ENDDO
141    C     Step forward Ud
142          DO j=jMin,jMax
143           DO i=iMin,iMax
144            uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) +
145         &   deltaTmom*vf(i,j)*_maskS(i,j,k,bi,bj)
146           ENDDO
147          ENDDO
148    C     Relax D grid U to C grid U
149          DO j=jMin,jMax
150           DO i=iMin,iMax
151            uVelD(i,j,k,bi,bj) = rCD*uVelD(i,j,k,bi,bj)
152         &   +(1. - rCD)*(
153         &    ab15*0.25*(
154         &                uVel(i,j  ,k,bi,bj)+uVel(i+1,j  ,k,bi,bj)
155         &               +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj)
156         &              )*_maskS(i,j,k,bi,bj)
157         &     +
158         &    ab05*0.25*(
159         &                uNM1(i,j  ,k,bi,bj)+uNM1(i+1,j  ,k,bi,bj)
160         &               +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj)
161         &              )*_maskS(i,j,k,bi,bj)
162         &   )
163           ENDDO
164          ENDDO
165    C     Calculate coriolis force on V
166          DO j=jMin,jMax
167           DO i=iMin,iMax
168            gvCD(i,j,k,bi,bj) =
169         &    -0.5*( _fCori(i  ,j,bi,bj)
170         &          +_fCori(i,j-1,bi,bj)  )
171         &   *uVelD(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)*cfFacMom
172           ENDDO
173          ENDDO
174    
175    C--   Save "previous time level" variables
176          DO j=1-OLy,sNy+OLy
177           DO i=1-OLx,sNx+OLx
178             uNM1(i,j,k,bi,bj) =  uVel(i,j,k,bi,bj)
179             vNM1(i,j,k,bi,bj) =  vVel(i,j,k,bi,bj)
180           ENDDO
181          ENDDO
182    
183    #endif /* INCLUDE_CD_CODE */
184    
185          RETURN
186          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22