/[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.4 - (show annotations) (download)
Sat Feb 8 02:13:02 2003 UTC (22 years, 5 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 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_cdscheme.F,v 1.3 2001/09/26 19:05:21 adcroft 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 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 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 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 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 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 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