/[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.5 - (hide annotations) (download)
Sun Feb 9 02:01:49 2003 UTC (21 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint48e_post, checkpoint48d_post, checkpoint48f_post
Changes since 1.4: +5 -4 lines 
in preparation for r*:
a) use pre-computed gradient of hydrostatic potential:
   changes in timestep.F & mom_cdscheme.F affects results of ideal_2D_oce
1 jmc 1.5 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_cdscheme.F,v 1.4 2003/02/08 02:13:02 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 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 jmc 1.5 c IF (.FALSE.) THEN
66     IF (staggerTimeStep) THEN
67 jmc 1.4 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 jmc 1.5 IF (.FALSE.) THEN
83     c IF (staggerTimeStep) THEN
84 adcroft 1.2 DO j=jMin,jMax
85     DO i=iMin,iMax
86     pf(i,j) = pf(i,j)+phi_hyd(i,j,k)
87     ENDDO
88     ENDDO
89     ENDIF
90    
91     C-- Zonal velocity coriolis term
92     C Note. As coded here, coriolis will not work with "thin walls"
93     C-- Coriolis with CD scheme allowed
94     C grady(p) + gV
95     DO j=1-Oly+1,sNy+Oly
96     DO i=1-Olx,sNx+Olx
97 jmc 1.4 af(i,j) = -_maskS(i,j,k,bi,bj)*(
98     & _recip_dyC(i,j,bi,bj)*(pf(i,j)-pf(i,j-1))
99     & +phyFac*dPhiHydY(i,j) )
100     & + gV(i,j,k,bi,bj)
101 adcroft 1.2 ENDDO
102     ENDDO
103     C Average to Vd point and add coriolis
104     DO j=jMin,jMax
105     DO i=iMin,iMax
106     vf(i,j) =
107     & 0.25*( af(i ,j)+af(i ,j+1)
108     & +af(i-1,j)+af(i-1,j+1)
109     & )*_maskW(i,j,k,bi,bj)
110     & -0.5*(_fCori(i,j,bi,bj)+
111     & _fCori(i-1,j,bi,bj))
112     & *uVel(i,j,k,bi,bj)
113     ENDDO
114     ENDDO
115     C Step forward Vd
116     DO j=jMin,jMax
117     DO i=iMin,iMax
118     vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) +
119     & deltaTmom*vf(i,j)
120     ENDDO
121     ENDDO
122     C Relax D grid V to C grid V
123     DO j=jMin,jMax
124     DO i=iMin,iMax
125     vVelD(i,j,k,bi,bj) = rCD*vVelD(i,j,k,bi,bj)
126     & +(1. - rCD)*(
127     & ab15*0.25*(
128     & vVel(i ,j ,k,bi,bj)+vVel(i ,j+1,k,bi,bj)
129     & +vVel(i-1,j ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj)
130     & )*_maskW(i,j,k,bi,bj)
131     & +
132     & ab05*0.25*(
133     & vNM1(i ,j ,k,bi,bj)+vNM1(i ,j+1,k,bi,bj)
134     & +vNM1(i-1,j ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj)
135     & )*_maskW(i,j,k,bi,bj)
136     & )
137     ENDDO
138     ENDDO
139     C Calculate coriolis force on U
140     DO j=jMin,jMax
141     DO i=iMin,iMax
142     guCD(i,j,k,bi,bj) =
143     & 0.5*( _fCori(i ,j,bi,bj) +
144     & _fCori(i-1,j,bi,bj) )
145     & *vVelD(i,j,k,bi,bj)*cfFacMom
146     ENDDO
147     ENDDO
148    
149     C-- Meridional velocity coriolis term
150     C gradx(p)+gU
151     DO j=1-Oly,sNy+Oly
152     DO i=1-Olx+1,sNx+Olx
153 jmc 1.4 af(i,j) = -_maskW(i,j,k,bi,bj)*(
154     & _recip_dxC(i,j,bi,bj)*(pf(i,j)-pf(i-1,j))
155     & +phxFac*dPhiHydX(i,j) )
156     & + gU(i,j,k,bi,bj)
157 adcroft 1.2 ENDDO
158     ENDDO
159     C Average to Ud point and add coriolis
160     DO j=jMin,jMax
161     DO i=iMin,iMax
162     vf(i,j) =
163     & 0.25*( af(i ,j)+af(i ,j-1)
164     & +af(i+1,j)+af(i+1,j-1)
165     & )*_maskS(i,j,k,bi,bj)
166     & +0.5*( _fCori(i,j,bi,bj)
167     & +_fCori(i,j-1,bi,bj))
168     & *vVel(i,j,k,bi,bj)
169     ENDDO
170     ENDDO
171     C Step forward Ud
172     DO j=jMin,jMax
173     DO i=iMin,iMax
174     uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) +
175     & deltaTmom*vf(i,j)*_maskS(i,j,k,bi,bj)
176     ENDDO
177     ENDDO
178     C Relax D grid U to C grid U
179     DO j=jMin,jMax
180     DO i=iMin,iMax
181     uVelD(i,j,k,bi,bj) = rCD*uVelD(i,j,k,bi,bj)
182     & +(1. - rCD)*(
183     & ab15*0.25*(
184     & uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj)
185     & +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj)
186     & )*_maskS(i,j,k,bi,bj)
187     & +
188     & ab05*0.25*(
189     & uNM1(i,j ,k,bi,bj)+uNM1(i+1,j ,k,bi,bj)
190     & +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj)
191     & )*_maskS(i,j,k,bi,bj)
192     & )
193     ENDDO
194     ENDDO
195     C Calculate coriolis force on V
196     DO j=jMin,jMax
197     DO i=iMin,iMax
198     gvCD(i,j,k,bi,bj) =
199     & -0.5*( _fCori(i ,j,bi,bj)
200     & +_fCori(i,j-1,bi,bj) )
201     & *uVelD(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)*cfFacMom
202     ENDDO
203     ENDDO
204    
205     C-- Save "previous time level" variables
206     DO j=1-OLy,sNy+OLy
207     DO i=1-OLx,sNx+OLx
208     uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
209     vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
210     ENDDO
211     ENDDO
212    
213     #endif /* INCLUDE_CD_CODE */
214    
215     RETURN
216     END
  ViewVC Help
Powered by ViewVC 1.1.22