/[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.7 - (hide annotations) (download)
Thu Apr 17 13:44:10 2003 UTC (21 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint51j_post, checkpoint50e_post, checkpoint50c_post, checkpoint50g_post, branchpoint-genmake2, checkpoint51e_post, checkpoint51b_post, checkpoint51f_pre, checkpoint50d_pre, checkpoint51, checkpoint50d_post, checkpoint51b_pre, checkpoint51h_pre, checkpoint50c_pre, checkpoint51g_post, checkpoint51f_post, checkpoint50b_post, checkpoint51c_post, checkpoint50f_post, checkpoint50f_pre, checkpoint51d_post, checkpoint50h_post, checkpoint51a_post, checkpoint50e_pre, checkpoint50i_post, checkpoint51i_pre
Branch point for: branch-genmake2
Changes since 1.6: +42 -26 lines 
o to allow to put the momForcing out of the Adams-Bashforth:
  move forcing & CD-scheme calls from mom_fluxform & mom_vecinv
  to timestep.F
o new flag "useCDscheme" (default=F); replace guCD,gvCD by local arrays
1 jmc 1.7 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 jmc 1.7 I bi,bj,k, dPhiHydX,dPhiHydY, guFld,gvFld,
12     O guCor,gvCor,
13     I myTime, myIter, myThid)
14 adcroft 1.2
15 adcroft 1.3 C !DESCRIPTION:
16     C The C-D scheme. The less said the better :-)
17 adcroft 1.2
18 adcroft 1.3 C !USES: ===============================================================
19 adcroft 1.2 C == Global variables ==
20 adcroft 1.3 IMPLICIT NONE
21 adcroft 1.2 #include "SIZE.h"
22     #include "DYNVARS.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "GRID.h"
26     #include "SURFACE.h"
27    
28 adcroft 1.3 C !INPUT PARAMETERS: ===================================================
29     C bi,bj :: tile indices
30     C k :: vertical level
31 jmc 1.4 C dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
32 jmc 1.7 C guFld,gvFld :: Acceleration (U & V compon.) from the C grid
33     C guCor,gvCor :: Coriolis terms (2 compon.) computed on C grid
34     C myTime :: current time
35     C myIter :: current time-step number
36 adcroft 1.3 C myThid :: thread number
37 jmc 1.7
38 jmc 1.4 INTEGER bi,bj,k
39     _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41 jmc 1.7 _RL guFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42     _RL gvFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43     _RL guCor(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44     _RL gvCor(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45     _RL myTime
46     INTEGER myIter
47 adcroft 1.2 INTEGER myThid
48    
49 adcroft 1.3
50     C !LOCAL VARIABLES: ====================================================
51 adcroft 1.2 #ifdef INCLUDE_CD_CODE
52 adcroft 1.3 C i,j :: loop indices
53     C pF :: pressure gradient
54     C vF :: work space
55     C aF :: work space
56 adcroft 1.2 _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57     _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58     _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59     INTEGER i,j,iMin,iMax,jMin,jMax
60     _RL ab15,ab05
61 jmc 1.4 _RL phxFac, phyFac
62 adcroft 1.3 CEOP
63 adcroft 1.2
64     C Compute ranges
65     iMin=1-Olx+1
66     iMax=sNx+Olx-1
67     jMin=1-Oly+1
68     jMax=sNy+Oly-1
69    
70     C Adams-Bashforth weighting factors
71     ab15 = 1.5 + abEps
72     ab05 = -0.5 - abEps
73    
74 jmc 1.4 C-- stagger time stepping: grad Phi_Hyp is not in gU,gV and needs to be added:
75 jmc 1.5 IF (staggerTimeStep) THEN
76 jmc 1.4 phxFac = pfFacMom
77     phyFac = pfFacMom
78     ELSE
79     phxFac = 0.
80     phyFac = 0.
81     ENDIF
82    
83 jmc 1.7 C- Initialize output (dummy) arrays:
84     c DO j=1-Oly,sNy+Oly
85     c DO i=1-Olx,sNx+Olx
86     c guCor(i,j) = 0. _d 0
87     c gvCor(i,j) = 0. _d 0
88     c ENDDO
89     c ENDDO
90    
91 adcroft 1.2 C Pressure extrapolated forward in time
92     DO j=1-Oly,sNy+Oly
93     DO i=1-Olx,sNx+Olx
94     pf(i,j) =
95     & ab15*( etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
96     & +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
97     ENDDO
98     ENDDO
99    
100     C-- Zonal velocity coriolis term
101     C Note. As coded here, coriolis will not work with "thin walls"
102     C-- Coriolis with CD scheme allowed
103     C grady(p) + gV
104     DO j=1-Oly+1,sNy+Oly
105     DO i=1-Olx,sNx+Olx
106 jmc 1.4 af(i,j) = -_maskS(i,j,k,bi,bj)*(
107     & _recip_dyC(i,j,bi,bj)*(pf(i,j)-pf(i,j-1))
108     & +phyFac*dPhiHydY(i,j) )
109 jmc 1.7 & + gvFld(i,j)
110 adcroft 1.2 ENDDO
111     ENDDO
112     C Average to Vd point and add coriolis
113     DO j=jMin,jMax
114     DO i=iMin,iMax
115     vf(i,j) =
116     & 0.25*( af(i ,j)+af(i ,j+1)
117     & +af(i-1,j)+af(i-1,j+1)
118     & )*_maskW(i,j,k,bi,bj)
119     & -0.5*(_fCori(i,j,bi,bj)+
120     & _fCori(i-1,j,bi,bj))
121     & *uVel(i,j,k,bi,bj)
122     ENDDO
123     ENDDO
124     C Step forward Vd
125     DO j=jMin,jMax
126     DO i=iMin,iMax
127 jmc 1.7 vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) + deltaTmom*vf(i,j)
128 adcroft 1.2 ENDDO
129     ENDDO
130     C Relax D grid V to C grid V
131     DO j=jMin,jMax
132     DO i=iMin,iMax
133 jmc 1.7 vVelD(i,j,k,bi,bj) = ( rCD*vVelD(i,j,k,bi,bj)
134 adcroft 1.2 & +(1. - rCD)*(
135 jmc 1.7 & ab15*(
136 adcroft 1.2 & vVel(i ,j ,k,bi,bj)+vVel(i ,j+1,k,bi,bj)
137     & +vVel(i-1,j ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj)
138 jmc 1.7 & )*0.25
139     & +ab05*(
140 adcroft 1.2 & vNM1(i ,j ,k,bi,bj)+vNM1(i ,j+1,k,bi,bj)
141     & +vNM1(i-1,j ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj)
142 jmc 1.7 & )*0.25
143     & ) )*_maskW(i,j,k,bi,bj)
144 adcroft 1.2 ENDDO
145     ENDDO
146     C Calculate coriolis force on U
147     DO j=jMin,jMax
148     DO i=iMin,iMax
149 jmc 1.7 guCor(i,j) =
150 adcroft 1.2 & 0.5*( _fCori(i ,j,bi,bj) +
151     & _fCori(i-1,j,bi,bj) )
152     & *vVelD(i,j,k,bi,bj)*cfFacMom
153     ENDDO
154     ENDDO
155    
156     C-- Meridional velocity coriolis term
157     C gradx(p)+gU
158     DO j=1-Oly,sNy+Oly
159     DO i=1-Olx+1,sNx+Olx
160 jmc 1.4 af(i,j) = -_maskW(i,j,k,bi,bj)*(
161     & _recip_dxC(i,j,bi,bj)*(pf(i,j)-pf(i-1,j))
162     & +phxFac*dPhiHydX(i,j) )
163 jmc 1.7 & + guFld(i,j)
164 adcroft 1.2 ENDDO
165     ENDDO
166     C Average to Ud point and add coriolis
167     DO j=jMin,jMax
168     DO i=iMin,iMax
169     vf(i,j) =
170     & 0.25*( af(i ,j)+af(i ,j-1)
171     & +af(i+1,j)+af(i+1,j-1)
172     & )*_maskS(i,j,k,bi,bj)
173     & +0.5*( _fCori(i,j,bi,bj)
174     & +_fCori(i,j-1,bi,bj))
175     & *vVel(i,j,k,bi,bj)
176     ENDDO
177     ENDDO
178     C Step forward Ud
179     DO j=jMin,jMax
180     DO i=iMin,iMax
181 jmc 1.7 uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) + deltaTmom*vf(i,j)
182 adcroft 1.2 ENDDO
183     ENDDO
184     C Relax D grid U to C grid U
185     DO j=jMin,jMax
186     DO i=iMin,iMax
187 jmc 1.7 uVelD(i,j,k,bi,bj) = ( rCD*uVelD(i,j,k,bi,bj)
188 adcroft 1.2 & +(1. - rCD)*(
189 jmc 1.7 & ab15*(
190 adcroft 1.2 & uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj)
191     & +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj)
192 jmc 1.7 & )*0.25
193     & +ab05*(
194 adcroft 1.2 & uNM1(i,j ,k,bi,bj)+uNM1(i+1,j ,k,bi,bj)
195     & +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj)
196 jmc 1.7 & )*0.25
197     & ) )*_maskS(i,j,k,bi,bj)
198 adcroft 1.2 ENDDO
199     ENDDO
200     C Calculate coriolis force on V
201     DO j=jMin,jMax
202     DO i=iMin,iMax
203 jmc 1.7 gvCor(i,j) =
204 adcroft 1.2 & -0.5*( _fCori(i ,j,bi,bj)
205     & +_fCori(i,j-1,bi,bj) )
206 jmc 1.7 & *uVelD(i,j,k,bi,bj)*cfFacMom
207 adcroft 1.2 ENDDO
208     ENDDO
209    
210     C-- Save "previous time level" variables
211     DO j=1-OLy,sNy+OLy
212     DO i=1-OLx,sNx+OLx
213     uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
214     vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
215     ENDDO
216     ENDDO
217    
218     #endif /* INCLUDE_CD_CODE */
219    
220     RETURN
221     END
  ViewVC Help
Powered by ViewVC 1.1.22