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

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

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

Revision 1.4 - (hide annotations) (download)
Sat Feb 8 02:13:02 2003 UTC (21 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint48d_pre
Changes since 1.3: +25 -11 lines 
in preparation for r*:
 new S/R (calc_grad_phi_hyd.F) to compute Hydrostatic potential gradient.
 pass the 2 comp. of the grad. as arguments to momentum S/R.
for the moment, only used if it does not change the results.
1 jmc 1.4 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_cdscheme.F,v 1.3 2001/09/26 19:05:21 adcroft Exp $
2 adcroft 1.3 C $Name: $
3 adcroft 1.2
4     #include "CPP_OPTIONS.h"
5    
6 adcroft 1.3 CBOP
7     C !ROUTINE: MOM_CDSCHEME
8    
9     C !INTERFACE: ==========================================================
10 adcroft 1.2 SUBROUTINE MOM_CDSCHEME(
11 jmc 1.4 I bi,bj,k,phi_hyd,dPhiHydX,dPhiHydY,
12 adcroft 1.2 I myThid)
13    
14 adcroft 1.3 C !DESCRIPTION:
15     C The C-D scheme. The less said the better :-)
16 adcroft 1.2
17 adcroft 1.3 C !USES: ===============================================================
18 adcroft 1.2 C == Global variables ==
19 adcroft 1.3 IMPLICIT NONE
20 adcroft 1.2 #include "SIZE.h"
21     #include "DYNVARS.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "GRID.h"
25     #include "SURFACE.h"
26    
27 adcroft 1.3 C !INPUT PARAMETERS: ===================================================
28     C bi,bj :: tile indices
29     C k :: vertical level
30     C phi_hyd :: hydrostatic pressure
31 jmc 1.4 C dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
32 adcroft 1.3 C myThid :: thread number
33 jmc 1.4 INTEGER bi,bj,k
34 adcroft 1.2 _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
35 jmc 1.4 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
36     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37 adcroft 1.2 INTEGER myThid
38    
39 adcroft 1.3
40     C !LOCAL VARIABLES: ====================================================
41 adcroft 1.2 #ifdef INCLUDE_CD_CODE
42 adcroft 1.3 C i,j :: loop indices
43     C pF :: pressure gradient
44     C vF :: work space
45     C aF :: work space
46 adcroft 1.2 _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47     _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48     _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49     INTEGER i,j,iMin,iMax,jMin,jMax
50     _RL ab15,ab05
51 jmc 1.4 _RL phxFac, phyFac
52 adcroft 1.3 CEOP
53 adcroft 1.2
54     C Compute ranges
55     iMin=1-Olx+1
56     iMax=sNx+Olx-1
57     jMin=1-Oly+1
58     jMax=sNy+Oly-1
59    
60     C Adams-Bashforth weighting factors
61     ab15 = 1.5 + abEps
62     ab05 = -0.5 - abEps
63    
64 jmc 1.4 C-- stagger time stepping: grad Phi_Hyp is not in gU,gV and needs to be added:
65     IF (.FALSE.) THEN
66     c IF (staggerTimeStep) THEN
67     phxFac = pfFacMom
68     phyFac = pfFacMom
69     ELSE
70     phxFac = 0.
71     phyFac = 0.
72     ENDIF
73    
74 adcroft 1.2 C Pressure extrapolated forward in time
75     DO j=1-Oly,sNy+Oly
76     DO i=1-Olx,sNx+Olx
77     pf(i,j) =
78     & ab15*( etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
79     & +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
80     ENDDO
81     ENDDO
82     IF (staggerTimeStep) THEN
83     DO j=jMin,jMax
84     DO i=iMin,iMax
85     pf(i,j) = pf(i,j)+phi_hyd(i,j,k)
86     ENDDO
87     ENDDO
88     ENDIF
89    
90     C-- Zonal velocity coriolis term
91     C Note. As coded here, coriolis will not work with "thin walls"
92     C-- Coriolis with CD scheme allowed
93     C grady(p) + gV
94     DO j=1-Oly+1,sNy+Oly
95     DO i=1-Olx,sNx+Olx
96 jmc 1.4 af(i,j) = -_maskS(i,j,k,bi,bj)*(
97     & _recip_dyC(i,j,bi,bj)*(pf(i,j)-pf(i,j-1))
98     & +phyFac*dPhiHydY(i,j) )
99     & + gV(i,j,k,bi,bj)
100 adcroft 1.2 ENDDO
101     ENDDO
102     C Average to Vd point and add coriolis
103     DO j=jMin,jMax
104     DO i=iMin,iMax
105     vf(i,j) =
106     & 0.25*( af(i ,j)+af(i ,j+1)
107     & +af(i-1,j)+af(i-1,j+1)
108     & )*_maskW(i,j,k,bi,bj)
109     & -0.5*(_fCori(i,j,bi,bj)+
110     & _fCori(i-1,j,bi,bj))
111     & *uVel(i,j,k,bi,bj)
112     ENDDO
113     ENDDO
114     C Step forward Vd
115     DO j=jMin,jMax
116     DO i=iMin,iMax
117     vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) +
118     & deltaTmom*vf(i,j)
119     ENDDO
120     ENDDO
121     C Relax D grid V to C grid V
122     DO j=jMin,jMax
123     DO i=iMin,iMax
124     vVelD(i,j,k,bi,bj) = rCD*vVelD(i,j,k,bi,bj)
125     & +(1. - rCD)*(
126     & ab15*0.25*(
127     & vVel(i ,j ,k,bi,bj)+vVel(i ,j+1,k,bi,bj)
128     & +vVel(i-1,j ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj)
129     & )*_maskW(i,j,k,bi,bj)
130     & +
131     & ab05*0.25*(
132     & vNM1(i ,j ,k,bi,bj)+vNM1(i ,j+1,k,bi,bj)
133     & +vNM1(i-1,j ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj)
134     & )*_maskW(i,j,k,bi,bj)
135     & )
136     ENDDO
137     ENDDO
138     C Calculate coriolis force on U
139     DO j=jMin,jMax
140     DO i=iMin,iMax
141     guCD(i,j,k,bi,bj) =
142     & 0.5*( _fCori(i ,j,bi,bj) +
143     & _fCori(i-1,j,bi,bj) )
144     & *vVelD(i,j,k,bi,bj)*cfFacMom
145     ENDDO
146     ENDDO
147    
148     C-- Meridional velocity coriolis term
149     C gradx(p)+gU
150     DO j=1-Oly,sNy+Oly
151     DO i=1-Olx+1,sNx+Olx
152 jmc 1.4 af(i,j) = -_maskW(i,j,k,bi,bj)*(
153     & _recip_dxC(i,j,bi,bj)*(pf(i,j)-pf(i-1,j))
154     & +phxFac*dPhiHydX(i,j) )
155     & + gU(i,j,k,bi,bj)
156 adcroft 1.2 ENDDO
157     ENDDO
158     C Average to Ud point and add coriolis
159     DO j=jMin,jMax
160     DO i=iMin,iMax
161     vf(i,j) =
162     & 0.25*( af(i ,j)+af(i ,j-1)
163     & +af(i+1,j)+af(i+1,j-1)
164     & )*_maskS(i,j,k,bi,bj)
165     & +0.5*( _fCori(i,j,bi,bj)
166     & +_fCori(i,j-1,bi,bj))
167     & *vVel(i,j,k,bi,bj)
168     ENDDO
169     ENDDO
170     C Step forward Ud
171     DO j=jMin,jMax
172     DO i=iMin,iMax
173     uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) +
174     & deltaTmom*vf(i,j)*_maskS(i,j,k,bi,bj)
175     ENDDO
176     ENDDO
177     C Relax D grid U to C grid U
178     DO j=jMin,jMax
179     DO i=iMin,iMax
180     uVelD(i,j,k,bi,bj) = rCD*uVelD(i,j,k,bi,bj)
181     & +(1. - rCD)*(
182     & ab15*0.25*(
183     & uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj)
184     & +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj)
185     & )*_maskS(i,j,k,bi,bj)
186     & +
187     & ab05*0.25*(
188     & uNM1(i,j ,k,bi,bj)+uNM1(i+1,j ,k,bi,bj)
189     & +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj)
190     & )*_maskS(i,j,k,bi,bj)
191     & )
192     ENDDO
193     ENDDO
194     C Calculate coriolis force on V
195     DO j=jMin,jMax
196     DO i=iMin,iMax
197     gvCD(i,j,k,bi,bj) =
198     & -0.5*( _fCori(i ,j,bi,bj)
199     & +_fCori(i,j-1,bi,bj) )
200     & *uVelD(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)*cfFacMom
201     ENDDO
202     ENDDO
203    
204     C-- Save "previous time level" variables
205     DO j=1-OLy,sNy+OLy
206     DO i=1-OLx,sNx+OLx
207     uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
208     vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
209     ENDDO
210     ENDDO
211    
212     #endif /* INCLUDE_CD_CODE */
213    
214     RETURN
215     END
  ViewVC Help
Powered by ViewVC 1.1.22