/[MITgcm]/MITgcm/pkg/cd_code/cd_code_scheme.F
ViewVC logotype

Contents of /MITgcm/pkg/cd_code/cd_code_scheme.F

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

Revision 1.1 - (show annotations) (download)
Thu Oct 30 18:44:26 2003 UTC (20 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint52l_pre, checkpoint52e_pre, hrcube4, checkpoint52n_post, checkpoint52j_post, checkpoint53d_post, checkpoint54a_pre, checkpoint55c_post, checkpoint54e_post, checkpoint52e_post, checkpoint54a_post, checkpoint53c_post, checkpoint55d_pre, checkpoint51q_post, hrcube_1, branch-netcdf, checkpoint52d_pre, checkpoint52l_post, checkpoint55h_post, checkpoint51r_post, checkpoint52k_post, checkpoint52b_pre, checkpoint54b_post, checkpoint53b_pre, checkpoint55b_post, checkpoint54d_post, checkpoint52m_post, checkpoint55, checkpoint53a_post, checkpoint54, checkpoint54f_post, checkpoint53b_post, checkpoint55g_post, checkpoint52a_pre, checkpoint55f_post, checkpoint55i_post, checkpoint53, checkpoint52, checkpoint52d_post, checkpoint52a_post, checkpoint52b_post, checkpoint53g_post, checkpoint52f_post, checkpoint52c_post, ecco_c52_e35, hrcube5, checkpoint55e_post, checkpoint52i_post, checkpoint52j_pre, checkpoint53f_post, checkpoint55a_post, checkpoint51t_post, checkpoint53d_pre, checkpoint54c_post, checkpoint52i_pre, checkpoint51u_post, checkpoint52h_pre, checkpoint52f_pre, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint55d_post
Branch point for: branch-nonh, netcdf-sm0 
modified pkg/cd_code
o moved cd_scheme.F -> cd_code_scheme.F
o separate read_checkpoint from cd_code_ini_vars.F
o separated cd_code part from write_checkpoint
o updated AD_SOURCE, generated .flow
o added CD_CODE_VARS.h to the_main_loop
1 C $Header: /u/gcmpack/MITgcm/pkg/cd_code/cd_scheme.F,v 1.3 2003/10/30 12:00:41 edhill Exp $
2 C $Name: $
3
4 #include "CD_CODE_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: CD_CODE_SCHEME
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE CD_CODE_SCHEME(
11 I bi,bj,k, dPhiHydX,dPhiHydY, guFld,gvFld,
12 O guCor,gvCor,
13 I myTime, myIter, myThid)
14
15 C !DESCRIPTION:
16 C The C-D scheme. The less said the better :-)
17
18 C !USES: ===============================================================
19 C == Global variables ==
20 IMPLICIT NONE
21 #include "SIZE.h"
22 #include "DYNVARS.h"
23 #ifdef ALLOW_CD_CODE
24 #include "CD_CODE_VARS.h"
25 #endif
26 #include "EEPARAMS.h"
27 #include "PARAMS.h"
28 #include "GRID.h"
29 #include "SURFACE.h"
30
31 C !INPUT PARAMETERS: ===================================================
32 C bi,bj :: tile indices
33 C k :: vertical level
34 C dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
35 C guFld,gvFld :: Acceleration (U & V compon.) from the C grid
36 C guCor,gvCor :: Coriolis terms (2 compon.) computed on C grid
37 C myTime :: current time
38 C myIter :: current time-step number
39 C myThid :: thread number
40
41 INTEGER bi,bj,k
42 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43 _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44 _RL guFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45 _RL gvFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46 _RL guCor(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47 _RL gvCor(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 _RL myTime
49 INTEGER myIter
50 INTEGER myThid
51
52
53 C !LOCAL VARIABLES: ====================================================
54 #ifdef ALLOW_CD_CODE
55 C i,j :: loop indices
56 C pF :: pressure gradient
57 C vF :: work space
58 C aF :: work space
59 _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60 _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 INTEGER i,j,iMin,iMax,jMin,jMax
63 _RL ab15,ab05
64 _RL phxFac, phyFac
65 CEOP
66
67 C Compute ranges
68 iMin=1-Olx+1
69 iMax=sNx+Olx-1
70 jMin=1-Oly+1
71 jMax=sNy+Oly-1
72
73 C Adams-Bashforth weighting factors
74 ab15 = 1.5 + abEps
75 ab05 = -0.5 - abEps
76
77 C-- stagger time stepping: grad Phi_Hyp is not in gU,gV and needs to be added:
78 IF (staggerTimeStep) THEN
79 phxFac = pfFacMom
80 phyFac = pfFacMom
81 ELSE
82 phxFac = 0.
83 phyFac = 0.
84 ENDIF
85
86 C- Initialize output (dummy) arrays:
87 c DO j=1-Oly,sNy+Oly
88 c DO i=1-Olx,sNx+Olx
89 c guCor(i,j) = 0. _d 0
90 c gvCor(i,j) = 0. _d 0
91 c ENDDO
92 c ENDDO
93
94 C Pressure extrapolated forward in time
95 DO j=1-Oly,sNy+Oly
96 DO i=1-Olx,sNx+Olx
97 pf(i,j) =
98 & ab15*( etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
99 & +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
100 ENDDO
101 ENDDO
102
103 C-- Zonal velocity coriolis term
104 C Note. As coded here, coriolis will not work with "thin walls"
105 C-- Coriolis with CD scheme allowed
106 C grady(p) + gV
107 DO j=1-Oly+1,sNy+Oly
108 DO i=1-Olx,sNx+Olx
109 af(i,j) = -_maskS(i,j,k,bi,bj)*(
110 & _recip_dyC(i,j,bi,bj)*(pf(i,j)-pf(i,j-1))
111 & +phyFac*dPhiHydY(i,j) )
112 & + gvFld(i,j)
113 ENDDO
114 ENDDO
115 C Average to Vd point and add coriolis
116 DO j=jMin,jMax
117 DO i=iMin,iMax
118 vf(i,j) =
119 & 0.25*( af(i ,j)+af(i ,j+1)
120 & +af(i-1,j)+af(i-1,j+1)
121 & )*_maskW(i,j,k,bi,bj)
122 & -0.5*(_fCori(i,j,bi,bj)+
123 & _fCori(i-1,j,bi,bj))
124 & *uVel(i,j,k,bi,bj)
125 ENDDO
126 ENDDO
127 C Step forward Vd
128 DO j=jMin,jMax
129 DO i=iMin,iMax
130 vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) + deltaTmom*vf(i,j)
131 ENDDO
132 ENDDO
133 C Relax D grid V to C grid V
134 DO j=jMin,jMax
135 DO i=iMin,iMax
136 vVelD(i,j,k,bi,bj) = ( rCD*vVelD(i,j,k,bi,bj)
137 & +(1. - rCD)*(
138 & ab15*(
139 & vVel(i ,j ,k,bi,bj)+vVel(i ,j+1,k,bi,bj)
140 & +vVel(i-1,j ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj)
141 & )*0.25
142 & +ab05*(
143 & vNM1(i ,j ,k,bi,bj)+vNM1(i ,j+1,k,bi,bj)
144 & +vNM1(i-1,j ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj)
145 & )*0.25
146 & ) )*_maskW(i,j,k,bi,bj)
147 ENDDO
148 ENDDO
149 C Calculate coriolis force on U
150 DO j=jMin,jMax
151 DO i=iMin,iMax
152 guCor(i,j) =
153 & 0.5*( _fCori(i ,j,bi,bj) +
154 & _fCori(i-1,j,bi,bj) )
155 & *vVelD(i,j,k,bi,bj)*cfFacMom
156 ENDDO
157 ENDDO
158
159 C-- Meridional velocity coriolis term
160 C gradx(p)+gU
161 DO j=1-Oly,sNy+Oly
162 DO i=1-Olx+1,sNx+Olx
163 af(i,j) = -_maskW(i,j,k,bi,bj)*(
164 & _recip_dxC(i,j,bi,bj)*(pf(i,j)-pf(i-1,j))
165 & +phxFac*dPhiHydX(i,j) )
166 & + guFld(i,j)
167 ENDDO
168 ENDDO
169 C Average to Ud point and add coriolis
170 DO j=jMin,jMax
171 DO i=iMin,iMax
172 vf(i,j) =
173 & 0.25*( af(i ,j)+af(i ,j-1)
174 & +af(i+1,j)+af(i+1,j-1)
175 & )*_maskS(i,j,k,bi,bj)
176 & +0.5*( _fCori(i,j,bi,bj)
177 & +_fCori(i,j-1,bi,bj))
178 & *vVel(i,j,k,bi,bj)
179 ENDDO
180 ENDDO
181 C Step forward Ud
182 DO j=jMin,jMax
183 DO i=iMin,iMax
184 uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) + deltaTmom*vf(i,j)
185 ENDDO
186 ENDDO
187 C Relax D grid U to C grid U
188 DO j=jMin,jMax
189 DO i=iMin,iMax
190 uVelD(i,j,k,bi,bj) = ( rCD*uVelD(i,j,k,bi,bj)
191 & +(1. - rCD)*(
192 & ab15*(
193 & uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj)
194 & +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj)
195 & )*0.25
196 & +ab05*(
197 & uNM1(i,j ,k,bi,bj)+uNM1(i+1,j ,k,bi,bj)
198 & +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj)
199 & )*0.25
200 & ) )*_maskS(i,j,k,bi,bj)
201 ENDDO
202 ENDDO
203 C Calculate coriolis force on V
204 DO j=jMin,jMax
205 DO i=iMin,iMax
206 gvCor(i,j) =
207 & -0.5*( _fCori(i ,j,bi,bj)
208 & +_fCori(i,j-1,bi,bj) )
209 & *uVelD(i,j,k,bi,bj)*cfFacMom
210 ENDDO
211 ENDDO
212
213 C-- Save "previous time level" variables
214 DO j=1-OLy,sNy+OLy
215 DO i=1-OLx,sNx+OLx
216 uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
217 vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
218 ENDDO
219 ENDDO
220
221 #endif /* ALLOW_CD_CODE */
222
223 RETURN
224 END
  ViewVC Help
Powered by ViewVC 1.1.22