84 |
_RL pH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
_RL pH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
85 |
_RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
86 |
_RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
87 |
|
_RL rhotmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
88 |
_RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
89 |
_RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
90 |
_RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
_RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
97 |
INTEGER i, j |
INTEGER i, j |
98 |
INTEGER k, kM1, kUp, kDown |
INTEGER k, kM1, kUp, kDown |
99 |
|
|
100 |
|
C--- The algorithm... |
101 |
|
C |
102 |
|
C "Correction Step" |
103 |
|
C ================= |
104 |
|
C Here we update the horizontal velocities with the surface |
105 |
|
C pressure such that the resulting flow is either consistent |
106 |
|
C with the free-surface evolution or the rigid-lid: |
107 |
|
C U[n] = U* + dt x d/dx P |
108 |
|
C V[n] = V* + dt x d/dy P |
109 |
|
C With implicit diffusion, the tracers must also be "finalized" |
110 |
|
C (1 + dt * K * d_zz) theta[n] = theta* |
111 |
|
C (1 + dt * K * d_zz) salt[n] = salt* |
112 |
|
C |
113 |
|
C "Calculation of Gs" |
114 |
|
C =================== |
115 |
|
C This is where all the accelerations and tendencies (ie. |
116 |
|
C physics, parameterizations etc...) are calculated |
117 |
|
C w = sum_z ( div. u[n] ) |
118 |
|
C rho = rho ( theta[n], salt[n] ) |
119 |
|
C K31 = K31 ( rho ) |
120 |
|
C Gu[n] = Gu( u[n], v[n], w, rho, Ph, ... ) |
121 |
|
C Gv[n] = Gv( u[n], v[n], w, rho, Ph, ... ) |
122 |
|
C Gt[n] = Gt( theta[n], u[n], v[n], w, K31, ... ) |
123 |
|
C Gs[n] = Gs( salt[n], u[n], v[n], w, K31, ... ) |
124 |
|
C |
125 |
|
C "Time-stepping" or "Predicition" |
126 |
|
C ================================ |
127 |
|
C The models variables are stepped forward with the appropriate |
128 |
|
C time-stepping scheme (currently we use Adams-Bashforth II) |
129 |
|
C - For momentum, the result is always *only* a "prediction" |
130 |
|
C in that the flow may be divergent and will be "corrected" |
131 |
|
C later with a surface pressure gradient. |
132 |
|
C - Normally for tracers the result is the new field at time |
133 |
|
C level [n+1} *BUT* in the case of implicit diffusion the result |
134 |
|
C is also *only* a prediction. |
135 |
|
C - We denote "predictors" with an asterisk (*). |
136 |
|
C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] ) |
137 |
|
C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] ) |
138 |
|
C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) |
139 |
|
C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) |
140 |
|
C or with implicit diffusion |
141 |
|
C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) |
142 |
|
C |
143 |
|
C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) |
144 |
|
C--- |
145 |
|
|
146 |
|
|
147 |
C-- Set up work arrays with valid (i.e. not NaN) values |
C-- Set up work arrays with valid (i.e. not NaN) values |
148 |
C These inital values do not alter the numerical results. They |
C These inital values do not alter the numerical results. They |
149 |
C just ensure that all memory references are to valid floating |
C just ensure that all memory references are to valid floating |
170 |
ENDDO |
ENDDO |
171 |
rhokm1(i,j) = 0. _d 0 |
rhokm1(i,j) = 0. _d 0 |
172 |
rhokp1(i,j) = 0. _d 0 |
rhokp1(i,j) = 0. _d 0 |
173 |
|
rhotmp(i,j) = 0. _d 0 |
174 |
ENDDO |
ENDDO |
175 |
ENDDO |
ENDDO |
176 |
|
|
177 |
DO bj=myByLo(myThid),myByHi(myThid) |
DO bj=myByLo(myThid),myByHi(myThid) |
178 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
DO bi=myBxLo(myThid),myBxHi(myThid) |
179 |
|
|
|
C-- Boundary condition on hydrostatic pressure is pH(z=0)=0 |
|
|
DO j=1-OLy,sNy+OLy |
|
|
DO i=1-OLx,sNx+OLx |
|
|
pH(i,j,1) = 0. _d 0 |
|
|
K13(i,j,1) = 0. _d 0 |
|
|
K23(i,j,1) = 0. _d 0 |
|
|
K33(i,j,1) = 0. _d 0 |
|
|
KapGM(i,j) = 0. _d 0 |
|
|
ENDDO |
|
|
ENDDO |
|
|
|
|
180 |
C-- Set up work arrays that need valid initial values |
C-- Set up work arrays that need valid initial values |
181 |
DO j=1-OLy,sNy+OLy |
DO j=1-OLy,sNy+OLy |
182 |
DO i=1-OLx,sNx+OLx |
DO i=1-OLx,sNx+OLx |
189 |
fVerU(i,j,2) = 0. _d 0 |
fVerU(i,j,2) = 0. _d 0 |
190 |
fVerV(i,j,1) = 0. _d 0 |
fVerV(i,j,1) = 0. _d 0 |
191 |
fVerV(i,j,2) = 0. _d 0 |
fVerV(i,j,2) = 0. _d 0 |
192 |
|
pH(i,j,1) = 0. _d 0 |
193 |
|
K13(i,j,1) = 0. _d 0 |
194 |
|
K23(i,j,1) = 0. _d 0 |
195 |
|
K33(i,j,1) = 0. _d 0 |
196 |
|
KapGM(i,j) = 0. _d 0 |
197 |
ENDDO |
ENDDO |
198 |
ENDDO |
ENDDO |
199 |
|
|
209 |
I myThid) |
I myThid) |
210 |
|
|
211 |
C-- Update fields in top level according to tendency terms |
C-- Update fields in top level according to tendency terms |
212 |
CALL TIMESTEP( |
CALL CORRECTION_STEP( |
213 |
I bi,bj,iMin,iMax,jMin,jMax,1,pSurfX,pSurfY,myThid) |
I bi,bj,iMin,iMax,jMin,jMax,1,pSurfX,pSurfY,myThid) |
214 |
|
|
215 |
C-- Density of 1st level (below W(1)) reference to level 1 |
C-- Density of 1st level (below W(1)) reference to level 1 |
216 |
CALL FIND_RHO( |
CALL FIND_RHO( |
217 |
I bi, bj, iMin, iMax, jMin, jMax, 1, 1, 'LINEAR', |
I bi, bj, iMin, iMax, jMin, jMax, 1, 1, eosType, |
218 |
O rhoKm1, |
O rhoKm1, |
219 |
I myThid ) |
I myThid ) |
220 |
C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0 |
C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0 |
222 |
I bi,bj,iMin,iMax,jMin,jMax,1,rhoKm1,rhoKm1, |
I bi,bj,iMin,iMax,jMin,jMax,1,rhoKm1,rhoKm1, |
223 |
U pH, |
U pH, |
224 |
I myThid ) |
I myThid ) |
225 |
|
DO J=1-Oly,sNy+Oly |
226 |
|
DO I=1-Olx,sNx+Olx |
227 |
|
rhoKp1(I,J)=rhoKm1(I,J) |
228 |
|
ENDDO |
229 |
|
ENDDO |
230 |
|
|
231 |
DO K=2,Nz |
DO K=2,Nz |
232 |
C-- Update fields in Kth level according to tendency terms |
C-- Update fields in Kth level according to tendency terms |
233 |
CALL TIMESTEP( |
CALL CORRECTION_STEP( |
234 |
I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid) |
I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid) |
235 |
C-- Density of K-1 level (above W(K)) reference to K level |
C-- Density of K-1 level (above W(K)) reference to K-1 level |
236 |
CALL FIND_RHO( |
copt CALL FIND_RHO( |
237 |
I bi, bj, iMin, iMax, jMin, jMax, K-1, K, 'LINEAR', |
copt I bi, bj, iMin, iMax, jMin, jMax, K-1, K-1, eosType, |
238 |
O rhoKm1, |
copt O rhoKm1, |
239 |
I myThid ) |
copt I myThid ) |
240 |
|
C rhoKm1=rhoKp1 |
241 |
|
DO J=1-Oly,sNy+Oly |
242 |
|
DO I=1-Olx,sNx+Olx |
243 |
|
rhoKm1(I,J)=rhoKp1(I,J) |
244 |
|
ENDDO |
245 |
|
ENDDO |
246 |
C-- Density of K level (below W(K)) reference to K level |
C-- Density of K level (below W(K)) reference to K level |
247 |
CALL FIND_RHO( |
CALL FIND_RHO( |
248 |
I bi, bj, iMin, iMax, jMin, jMax, K, K, 'LINEAR', |
I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType, |
249 |
O rhoKp1, |
O rhoKp1, |
250 |
I myThid ) |
I myThid ) |
251 |
|
C-- Density of K-1 level (above W(K)) reference to K level |
252 |
|
CALL FIND_RHO( |
253 |
|
I bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType, |
254 |
|
O rhotmp, |
255 |
|
I myThid ) |
256 |
C-- Calculate iso-neutral slopes for the GM/Redi parameterisation |
C-- Calculate iso-neutral slopes for the GM/Redi parameterisation |
257 |
CALL CALC_ISOSLOPES( |
CALL CALC_ISOSLOPES( |
258 |
I bi, bj, iMin, iMax, jMin, jMax, K, |
I bi, bj, iMin, iMax, jMin, jMax, K, |
259 |
I rhoKm1, rhoKp1, |
I rhoKm1, rhoKp1, rhotmp, |
260 |
O K13, K23, K33, KapGM, |
O K13, K23, K33, KapGM, |
261 |
I myThid ) |
I myThid ) |
262 |
C-- Calculate static stability and mix where convectively unstable |
C-- Calculate static stability and mix where convectively unstable |
265 |
I myTime,myIter,myThid) |
I myTime,myIter,myThid) |
266 |
C-- Density of K-1 level (above W(K)) reference to K-1 level |
C-- Density of K-1 level (above W(K)) reference to K-1 level |
267 |
CALL FIND_RHO( |
CALL FIND_RHO( |
268 |
I bi, bj, iMin, iMax, jMin, jMax, K-1, K-1, 'LINEAR', |
I bi, bj, iMin, iMax, jMin, jMax, K-1, K-1, eosType, |
269 |
O rhoKm1, |
O rhoKm1, |
270 |
I myThid ) |
I myThid ) |
271 |
C-- Density of K level (below W(K)) referenced to K level |
C-- Density of K level (below W(K)) referenced to K level |
272 |
CALL FIND_RHO( |
CALL FIND_RHO( |
273 |
I bi, bj, iMin, iMax, jMin, jMax, K, K, 'LINEAR', |
I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType, |
274 |
O rhoKp1, |
O rhoKp1, |
275 |
I myThid ) |
I myThid ) |
276 |
C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0 |
C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0 |
279 |
U pH, |
U pH, |
280 |
I myThid ) |
I myThid ) |
281 |
|
|
282 |
|
ENDDO ! K |
283 |
|
|
284 |
|
C-- Initial boundary condition on barotropic divergence integral |
285 |
|
DO j=1-OLy,sNy+OLy |
286 |
|
DO i=1-OLx,sNx+OLx |
287 |
|
cg2d_b(i,j,bi,bj) = 0. _d 0 |
288 |
|
ENDDO |
289 |
ENDDO |
ENDDO |
290 |
|
|
291 |
DO K = Nz, 1, -1 |
DO K = Nz, 1, -1 |
304 |
I myThid) |
I myThid) |
305 |
|
|
306 |
C-- Calculate accelerations in the momentum equations |
C-- Calculate accelerations in the momentum equations |
307 |
CALL CALC_MOM_RHS( |
IF ( momStepping ) THEN |
308 |
I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown, |
CALL CALC_MOM_RHS( |
309 |
I xA,yA,uTrans,vTrans,wTrans,maskC, |
I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown, |
310 |
I pH, |
I xA,yA,uTrans,vTrans,wTrans,maskC, |
311 |
U aTerm,xTerm,cTerm,mTerm,pTerm, |
I pH, |
312 |
U fZon, fMer, fVerU, fVerV, |
U aTerm,xTerm,cTerm,mTerm,pTerm, |
313 |
I myThid) |
U fZon, fMer, fVerU, fVerV, |
314 |
|
I myThid) |
315 |
|
ENDIF |
316 |
|
|
317 |
C-- Calculate active tracer tendencies |
C-- Calculate active tracer tendencies |
318 |
CALL CALC_GT( |
IF ( tempStepping ) THEN |
319 |
I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown, |
CALL CALC_GT( |
320 |
I xA,yA,uTrans,vTrans,wTrans,maskUp, |
I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown, |
321 |
I K13,K23,K33,KapGM, |
I xA,yA,uTrans,vTrans,wTrans,maskUp, |
322 |
U aTerm,xTerm,fZon,fMer,fVerT, |
I K13,K23,K33,KapGM, |
323 |
I myThid) |
U aTerm,xTerm,fZon,fMer,fVerT, |
324 |
|
I myThid) |
325 |
|
ENDIF |
326 |
Cdbg CALL CALC_GS( |
Cdbg CALL CALC_GS( |
327 |
Cdbg I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown, |
Cdbg I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown, |
328 |
Cdbg I xA,yA,uTrans,vTrans,wTrans,maskUp, |
Cdbg I xA,yA,uTrans,vTrans,wTrans,maskUp, |
330 |
Cdbg U aTerm,xTerm,fZon,fMer,fVerS, |
Cdbg U aTerm,xTerm,fZon,fMer,fVerS, |
331 |
Cdbg I myThid) |
Cdbg I myThid) |
332 |
|
|
333 |
ENDDO |
C-- Prediction step (step forward all model variables) |
334 |
|
CALL TIMESTEP( |
335 |
|
I bi,bj,iMin,iMax,jMin,jMax,K, |
336 |
|
I myThid) |
337 |
|
|
338 |
|
C-- Diagnose barotropic divergence of predicted fields |
339 |
|
CALL DIV_G( |
340 |
|
I bi,bj,iMin,iMax,jMin,jMax,K, |
341 |
|
I xA,yA, |
342 |
|
I myThid) |
343 |
|
|
344 |
|
ENDDO ! K |
345 |
|
|
346 |
ENDDO |
ENDDO |
347 |
ENDDO |
ENDDO |
351 |
!dbg & maxval(uVel(1:sNx,1:sNy,:,:,:)) |
!dbg & maxval(uVel(1:sNx,1:sNy,:,:,:)) |
352 |
!dbg write(0,*) 'dynamics: V',minval(vVel(1:sNx,1:sNy,:,:,:)), |
!dbg write(0,*) 'dynamics: V',minval(vVel(1:sNx,1:sNy,:,:,:)), |
353 |
!dbg & maxval(vVel(1:sNx,1:sNy,:,:,:)) |
!dbg & maxval(vVel(1:sNx,1:sNy,:,:,:)) |
354 |
|
!dbg write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)), |
355 |
|
!dbg & maxval(K13(1:sNx,1:sNy,:)) |
356 |
|
!dbg write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)), |
357 |
|
!dbg & maxval(K23(1:sNx,1:sNy,:)) |
358 |
|
!dbg write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)), |
359 |
|
!dbg & maxval(K33(1:sNx,1:sNy,:)) |
360 |
!dbg write(0,*) 'dynamics: gT',minval(gT(1:sNx,1:sNy,:,:,:)), |
!dbg write(0,*) 'dynamics: gT',minval(gT(1:sNx,1:sNy,:,:,:)), |
361 |
!dbg & maxval(gT(1:sNx,1:sNy,:,:,:)) |
!dbg & maxval(gT(1:sNx,1:sNy,:,:,:)) |
362 |
!dbg write(0,*) 'dynamics: T',minval(Theta(1:sNx,1:sNy,:,:,:)), |
!dbg write(0,*) 'dynamics: T',minval(Theta(1:sNx,1:sNy,:,:,:)), |