/[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.3.6.1 - (hide annotations) (download)
Fri Mar 7 04:46:40 2003 UTC (21 years, 2 months ago) by heimbach
Branch: ecco-branch
CVS Tags: ecco_c50_e33a, ecco_c50_e32, ecco_c50_e33, ecco_c50_e30, ecco_c50_e31, ecco_c50_e29, ecco_c50_e28
Changes since 1.3: +24 -20 lines 
merging c49 and e27
1 heimbach 1.3.6.1 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_cdscheme.F,v 1.6 2003/02/18 15:36:45 jmc 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 heimbach 1.3.6.1 I bi,bj,k,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 heimbach 1.3.6.1 C dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
31 adcroft 1.3 C myThid :: thread number
32 heimbach 1.3.6.1 INTEGER bi,bj,k
33     _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
34     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
35 adcroft 1.2 INTEGER myThid
36    
37 adcroft 1.3
38     C !LOCAL VARIABLES: ====================================================
39 adcroft 1.2 #ifdef INCLUDE_CD_CODE
40 adcroft 1.3 C i,j :: loop indices
41     C pF :: pressure gradient
42     C vF :: work space
43     C aF :: work space
44 adcroft 1.2 _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45     _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46     _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47     INTEGER i,j,iMin,iMax,jMin,jMax
48     _RL ab15,ab05
49 heimbach 1.3.6.1 _RL phxFac, phyFac
50 adcroft 1.3 CEOP
51 adcroft 1.2
52     C Compute ranges
53     iMin=1-Olx+1
54     iMax=sNx+Olx-1
55     jMin=1-Oly+1
56     jMax=sNy+Oly-1
57    
58     C Adams-Bashforth weighting factors
59     ab15 = 1.5 + abEps
60     ab05 = -0.5 - abEps
61    
62 heimbach 1.3.6.1 C-- stagger time stepping: grad Phi_Hyp is not in gU,gV and needs to be added:
63     IF (staggerTimeStep) THEN
64     phxFac = pfFacMom
65     phyFac = pfFacMom
66     ELSE
67     phxFac = 0.
68     phyFac = 0.
69     ENDIF
70    
71 adcroft 1.2 C Pressure extrapolated forward in time
72     DO j=1-Oly,sNy+Oly
73     DO i=1-Olx,sNx+Olx
74     pf(i,j) =
75     & ab15*( etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
76     & +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
77     ENDDO
78     ENDDO
79    
80     C-- Zonal velocity coriolis term
81     C Note. As coded here, coriolis will not work with "thin walls"
82     C-- Coriolis with CD scheme allowed
83     C grady(p) + gV
84     DO j=1-Oly+1,sNy+Oly
85     DO i=1-Olx,sNx+Olx
86 heimbach 1.3.6.1 af(i,j) = -_maskS(i,j,k,bi,bj)*(
87     & _recip_dyC(i,j,bi,bj)*(pf(i,j)-pf(i,j-1))
88     & +phyFac*dPhiHydY(i,j) )
89     & + gV(i,j,k,bi,bj)
90 adcroft 1.2 ENDDO
91     ENDDO
92     C Average to Vd point and add coriolis
93     DO j=jMin,jMax
94     DO i=iMin,iMax
95     vf(i,j) =
96     & 0.25*( af(i ,j)+af(i ,j+1)
97     & +af(i-1,j)+af(i-1,j+1)
98     & )*_maskW(i,j,k,bi,bj)
99     & -0.5*(_fCori(i,j,bi,bj)+
100     & _fCori(i-1,j,bi,bj))
101     & *uVel(i,j,k,bi,bj)
102     ENDDO
103     ENDDO
104     C Step forward Vd
105     DO j=jMin,jMax
106     DO i=iMin,iMax
107     vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) +
108     & deltaTmom*vf(i,j)
109     ENDDO
110     ENDDO
111     C Relax D grid V to C grid V
112     DO j=jMin,jMax
113     DO i=iMin,iMax
114     vVelD(i,j,k,bi,bj) = rCD*vVelD(i,j,k,bi,bj)
115     & +(1. - rCD)*(
116     & ab15*0.25*(
117     & vVel(i ,j ,k,bi,bj)+vVel(i ,j+1,k,bi,bj)
118     & +vVel(i-1,j ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj)
119     & )*_maskW(i,j,k,bi,bj)
120     & +
121     & ab05*0.25*(
122     & vNM1(i ,j ,k,bi,bj)+vNM1(i ,j+1,k,bi,bj)
123     & +vNM1(i-1,j ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj)
124     & )*_maskW(i,j,k,bi,bj)
125     & )
126     ENDDO
127     ENDDO
128     C Calculate coriolis force on U
129     DO j=jMin,jMax
130     DO i=iMin,iMax
131     guCD(i,j,k,bi,bj) =
132     & 0.5*( _fCori(i ,j,bi,bj) +
133     & _fCori(i-1,j,bi,bj) )
134     & *vVelD(i,j,k,bi,bj)*cfFacMom
135     ENDDO
136     ENDDO
137    
138     C-- Meridional velocity coriolis term
139     C gradx(p)+gU
140     DO j=1-Oly,sNy+Oly
141     DO i=1-Olx+1,sNx+Olx
142 heimbach 1.3.6.1 af(i,j) = -_maskW(i,j,k,bi,bj)*(
143     & _recip_dxC(i,j,bi,bj)*(pf(i,j)-pf(i-1,j))
144     & +phxFac*dPhiHydX(i,j) )
145     & + gU(i,j,k,bi,bj)
146 adcroft 1.2 ENDDO
147     ENDDO
148     C Average to Ud point and add coriolis
149     DO j=jMin,jMax
150     DO i=iMin,iMax
151     vf(i,j) =
152     & 0.25*( af(i ,j)+af(i ,j-1)
153     & +af(i+1,j)+af(i+1,j-1)
154     & )*_maskS(i,j,k,bi,bj)
155     & +0.5*( _fCori(i,j,bi,bj)
156     & +_fCori(i,j-1,bi,bj))
157     & *vVel(i,j,k,bi,bj)
158     ENDDO
159     ENDDO
160     C Step forward Ud
161     DO j=jMin,jMax
162     DO i=iMin,iMax
163     uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) +
164     & deltaTmom*vf(i,j)*_maskS(i,j,k,bi,bj)
165     ENDDO
166     ENDDO
167     C Relax D grid U to C grid U
168     DO j=jMin,jMax
169     DO i=iMin,iMax
170     uVelD(i,j,k,bi,bj) = rCD*uVelD(i,j,k,bi,bj)
171     & +(1. - rCD)*(
172     & ab15*0.25*(
173     & uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj)
174     & +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj)
175     & )*_maskS(i,j,k,bi,bj)
176     & +
177     & ab05*0.25*(
178     & uNM1(i,j ,k,bi,bj)+uNM1(i+1,j ,k,bi,bj)
179     & +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj)
180     & )*_maskS(i,j,k,bi,bj)
181     & )
182     ENDDO
183     ENDDO
184     C Calculate coriolis force on V
185     DO j=jMin,jMax
186     DO i=iMin,iMax
187     gvCD(i,j,k,bi,bj) =
188     & -0.5*( _fCori(i ,j,bi,bj)
189     & +_fCori(i,j-1,bi,bj) )
190     & *uVelD(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)*cfFacMom
191     ENDDO
192     ENDDO
193    
194     C-- Save "previous time level" variables
195     DO j=1-OLy,sNy+OLy
196     DO i=1-OLx,sNx+OLx
197     uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
198     vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
199     ENDDO
200     ENDDO
201    
202     #endif /* INCLUDE_CD_CODE */
203    
204     RETURN
205     END
  ViewVC Help
Powered by ViewVC 1.1.22