/[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.6 - (show annotations) (download)
Tue Feb 18 15:36:45 2003 UTC (22 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint49, checkpoint48i_post, checkpoint48h_post, checkpoint50, checkpoint50b_pre, checkpoint50a_post, checkpoint48g_post
Changes since 1.5: +2 -13 lines
"clean-up" unused phiHyd.

1 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_cdscheme.F,v 1.5 2003/02/09 02:01:49 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,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 dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
31 C myThid :: thread number
32 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 INTEGER myThid
36
37
38 C !LOCAL VARIABLES: ====================================================
39 #ifdef INCLUDE_CD_CODE
40 C i,j :: loop indices
41 C pF :: pressure gradient
42 C vF :: work space
43 C aF :: work space
44 _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 _RL phxFac, phyFac
50 CEOP
51
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 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 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 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 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 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 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