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

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

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

Revision 1.5 - (show annotations) (download)
Sun Feb 9 02:01:49 2003 UTC (22 years, 5 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 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_cdscheme.F,v 1.4 2003/02/08 02:13:02 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: MOM_CDSCHEME
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE MOM_CDSCHEME(
11 I bi,bj,k,phi_hyd,dPhiHydX,dPhiHydY,
12 I myThid)
13
14 C !DESCRIPTION:
15 C The C-D scheme. The less said the better :-)
16
17 C !USES: ===============================================================
18 C == Global variables ==
19 IMPLICIT NONE
20 #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 C !INPUT PARAMETERS: ===================================================
28 C bi,bj :: tile indices
29 C k :: vertical level
30 C phi_hyd :: hydrostatic pressure
31 C dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
32 C myThid :: thread number
33 INTEGER bi,bj,k
34 _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
35 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
36 _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37 INTEGER myThid
38
39
40 C !LOCAL VARIABLES: ====================================================
41 #ifdef INCLUDE_CD_CODE
42 C i,j :: loop indices
43 C pF :: pressure gradient
44 C vF :: work space
45 C aF :: work space
46 _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 _RL phxFac, phyFac
52 CEOP
53
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 C-- stagger time stepping: grad Phi_Hyp is not in gU,gV and needs to be added:
65 c IF (.FALSE.) THEN
66 IF (staggerTimeStep) THEN
67 phxFac = pfFacMom
68 phyFac = pfFacMom
69 ELSE
70 phxFac = 0.
71 phyFac = 0.
72 ENDIF
73
74 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 (.FALSE.) THEN
83 c IF (staggerTimeStep) THEN
84 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 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 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 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 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