62 |
C rhoK, rhoKM1 - Density at current level, level above and level below. |
C rhoK, rhoKM1 - Density at current level, level above and level below. |
63 |
C rhoKP1 |
C rhoKP1 |
64 |
C buoyK, buoyKM1 - Buoyancy at current level and level above. |
C buoyK, buoyKM1 - Buoyancy at current level and level above. |
65 |
C phiHyd - Hydrostatic part of the potential phi. |
C phiHyd - Hydrostatic part of the potential phiHydi. |
66 |
C In z coords phiHyd is the hydrostatic pressure anomaly |
C In z coords phiHydiHyd is the hydrostatic pressure anomaly |
67 |
C In p coords phiHyd is the geopotential surface height |
C In p coords phiHydiHyd is the geopotential surface height |
68 |
C anomaly. |
C anomaly. |
69 |
C etaSurfX, - Holds surface elevation gradient in X and Y. |
C etaSurfX, - Holds surface elevation gradient in X and Y. |
70 |
C etaSurfY |
C etaSurfY |
98 |
_RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
_RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
99 |
_RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
_RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
100 |
_RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
_RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
101 |
_RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
_RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
102 |
_RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
103 |
_RL rhokp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL rhokp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
104 |
_RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
107 |
_RL rhotmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL rhotmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
108 |
_RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
109 |
_RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
110 |
_RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
_RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
111 |
_RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
_RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
112 |
_RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz) |
_RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
113 |
_RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
114 |
_RL KappaZT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz) |
_RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) |
115 |
_RL KappaZS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz) |
_RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) |
116 |
|
|
117 |
INTEGER iMin, iMax |
INTEGER iMin, iMax |
118 |
INTEGER jMin, jMax |
INTEGER jMin, jMax |
134 |
C "Calculation of Gs" |
C "Calculation of Gs" |
135 |
C =================== |
C =================== |
136 |
C This is where all the accelerations and tendencies (ie. |
C This is where all the accelerations and tendencies (ie. |
137 |
C physics, parameterizations etc...) are calculated |
C phiHydysics, parameterizations etc...) are calculated |
138 |
C rVel = sum_r ( div. u[n] ) |
C rVel = sum_r ( div. u[n] ) |
139 |
C rho = rho ( theta[n], salt[n] ) |
C rho = rho ( theta[n], salt[n] ) |
140 |
C b = b(rho, theta) |
C b = b(rho, theta) |
184 |
pTerm(i,j) = 0. _d 0 |
pTerm(i,j) = 0. _d 0 |
185 |
fZon(i,j) = 0. _d 0 |
fZon(i,j) = 0. _d 0 |
186 |
fMer(i,j) = 0. _d 0 |
fMer(i,j) = 0. _d 0 |
187 |
DO K=1,nZ |
DO K=1,Nr |
188 |
pH (i,j,k) = 0. _d 0 |
phiHyd (i,j,k) = 0. _d 0 |
189 |
K13(i,j,k) = 0. _d 0 |
K13(i,j,k) = 0. _d 0 |
190 |
K23(i,j,k) = 0. _d 0 |
K23(i,j,k) = 0. _d 0 |
191 |
K33(i,j,k) = 0. _d 0 |
K33(i,j,k) = 0. _d 0 |
192 |
KappaZT(i,j,k) = 0. _d 0 |
KappaRT(i,j,k) = 0. _d 0 |
193 |
|
KappaRS(i,j,k) = 0. _d 0 |
194 |
ENDDO |
ENDDO |
195 |
rhoKM1 (i,j) = 0. _d 0 |
rhoKM1 (i,j) = 0. _d 0 |
196 |
rhok (i,j) = 0. _d 0 |
rhok (i,j) = 0. _d 0 |
202 |
ENDDO |
ENDDO |
203 |
ENDDO |
ENDDO |
204 |
|
|
205 |
|
|
206 |
DO bj=myByLo(myThid),myByHi(myThid) |
DO bj=myByLo(myThid),myByHi(myThid) |
207 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
DO bi=myBxLo(myThid),myBxHi(myThid) |
208 |
|
|
233 |
jMin = 1-OLy+1 |
jMin = 1-OLy+1 |
234 |
jMax = sNy+OLy |
jMax = sNy+OLy |
235 |
|
|
236 |
|
|
237 |
K = 1 |
K = 1 |
238 |
BOTTOM_LAYER = K .EQ. Nz |
BOTTOM_LAYER = K .EQ. Nr |
239 |
|
|
240 |
C-- Calculate gradient of surface pressure |
C-- Calculate gradient of surface pressure |
241 |
CALL CALC_GRAD_ETA_SURF( |
CALL CALC_GRAD_ETA_SURF( |
275 |
I myThid ) |
I myThid ) |
276 |
ENDIF |
ENDIF |
277 |
C-- Calculate buoyancy |
C-- Calculate buoyancy |
278 |
CALL CALC_BUOY( |
CALL CALC_BUOYANCY( |
279 |
I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1, |
I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1, |
280 |
O buoyKm1, |
O buoyKm1, |
281 |
I myThid ) |
I myThid ) |
282 |
C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0 |
C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0 |
283 |
CALL CALC_PHI_HYD( |
CALL CALC_PHI_HYD( |
284 |
I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1, |
I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1, |
285 |
U phiHyd, |
U phiHyd, |
286 |
I myThid ) |
I myThid ) |
287 |
|
|
288 |
DO K=2,Nz |
DO K=2,Nr |
289 |
BOTTOM_LAYER = K .EQ. Nz |
BOTTOM_LAYER = K .EQ. Nr |
290 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
IF ( .NOT. BOTTOM_LAYER ) THEN |
291 |
C-- Update fields in layer below according to tendency terms |
C-- Update fields in layer below according to tendency terms |
292 |
CALL CORRECTION_STEP( |
CALL CORRECTION_STEP( |
315 |
I myThid ) |
I myThid ) |
316 |
ENDIF |
ENDIF |
317 |
C-- Calculate buoyancy |
C-- Calculate buoyancy |
318 |
CALL CALC_BUOY( |
CALL CALC_BUOYANCY( |
319 |
I bi,bj,iMin,iMax,jMin,jMax,K,rhoK, |
I bi,bj,iMin,iMax,jMin,jMax,K,rhoK, |
320 |
O buoyK, |
O buoyK, |
321 |
I myThid ) |
I myThid ) |
322 |
C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0 |
C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0 |
323 |
CALL CALC_PHI_HYD( |
CALL CALC_PHI_HYD( |
324 |
I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK, |
I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK, |
325 |
U phiHyd, |
U phiHyd, |
342 |
ENDDO |
ENDDO |
343 |
ENDDO ! K |
ENDDO ! K |
344 |
|
|
345 |
DO K = Nz, 1, -1 |
DO K = Nr, 1, -1 |
346 |
|
|
347 |
kM1 =max(1,k-1) ! Points to level above k (=k-1) |
kM1 =max(1,k-1) ! Points to level above k (=k-1) |
348 |
kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above |
kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above |
361 |
CALL CALC_DIFFUSIVITY( |
CALL CALC_DIFFUSIVITY( |
362 |
I bi,bj,iMin,iMax,jMin,jMax,K, |
I bi,bj,iMin,iMax,jMin,jMax,K, |
363 |
I maskC,maskUp,KapGM,K33, |
I maskC,maskUp,KapGM,K33, |
364 |
O KappaZT,KappaZS, |
O KappaRT,KappaRS, |
365 |
I myThid) |
I myThid) |
366 |
C-- Calculate accelerations in the momentum equations |
C-- Calculate accelerations in the momentum equations |
367 |
IF ( momStepping ) THEN |
IF ( momStepping ) THEN |
380 |
I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC, |
I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC, |
381 |
I K13,K23,KappaRT,KapGM, |
I K13,K23,KappaRT,KapGM, |
382 |
U aTerm,xTerm,fZon,fMer,fVerT, |
U aTerm,xTerm,fZon,fMer,fVerT, |
383 |
I myThid) |
I myTime, myThid) |
384 |
ENDIF |
ENDIF |
385 |
IF ( saltStepping ) THEN |
IF ( saltStepping ) THEN |
386 |
CALL CALC_GS( |
CALL CALC_GS( |
388 |
I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC, |
I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC, |
389 |
I K13,K23,KappaRS,KapGM, |
I K13,K23,KappaRS,KapGM, |
390 |
U aTerm,xTerm,fZon,fMer,fVerS, |
U aTerm,xTerm,fZon,fMer,fVerS, |
391 |
I myThid) |
I myTime, myThid) |
392 |
ENDIF |
ENDIF |
393 |
C-- Prediction step (step forward all model variables) |
C-- Prediction step (step forward all model variables) |
394 |
CALL TIMESTEP( |
CALL TIMESTEP( |
395 |
I bi,bj,iMin,iMax,jMin,jMax,K, |
I bi,bj,iMin,iMax,jMin,jMax,K, |
396 |
I myThid) |
I myThid) |
397 |
C-- Diagnose barotropic divergence of predicted fields |
C-- Diagnose barotropic divergence of predicted fields |
398 |
CALL CALC_DIV_G( |
CALL CALC_DIV_GHAT( |
399 |
I bi,bj,iMin,iMax,jMin,jMax,K, |
I bi,bj,iMin,iMax,jMin,jMax,K, |
400 |
I xA,yA, |
I xA,yA, |
401 |
I myThid) |
I myThid) |
415 |
C-- Implicit diffusion |
C-- Implicit diffusion |
416 |
IF (implicitDiffusion) THEN |
IF (implicitDiffusion) THEN |
417 |
CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax, |
CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax, |
418 |
I KappaZT,KappaZS, |
I KappaRT,KappaRS, |
419 |
I myThid ) |
I myThid ) |
420 |
ENDIF |
ENDIF |
421 |
|
|
448 |
C & maxval(gS(1:sNx,1:sNy,:,:,:)) |
C & maxval(gS(1:sNx,1:sNy,:,:,:)) |
449 |
C write(0,*) 'dynamics: S ',minval(salt(1:sNx,1:sNy,:,:,:)), |
C write(0,*) 'dynamics: S ',minval(salt(1:sNx,1:sNy,:,:,:)), |
450 |
C & maxval(salt(1:sNx,1:sNy,:,:,:)) |
C & maxval(salt(1:sNx,1:sNy,:,:,:)) |
451 |
C write(0,*) 'dynamics: pH ',minval(pH/(Gravity*Rhonil),mask=ph.NE.0.), |
C write(0,*) 'dynamics: phiHyd ',minval(phiHyd/(Gravity*Rhonil),mask=phiHyd.NE.0.), |
452 |
C & maxval(pH/(Gravity*Rhonil)) |
C & maxval(phiHyd/(Gravity*Rhonil)) |
453 |
|
C CALL PLOT_FIELD_XYZRL( gU, ' GU exiting dyanmics ' , |
454 |
|
C &Nr, 1, myThid ) |
455 |
|
C CALL PLOT_FIELD_XYZRL( gV, ' GV exiting dyanmics ' , |
456 |
|
C &Nr, 1, myThid ) |
457 |
|
C CALL PLOT_FIELD_XYZRL( gS, ' GS exiting dyanmics ' , |
458 |
|
C &Nr, 1, myThid ) |
459 |
|
C CALL PLOT_FIELD_XYZRL( gT, ' GT exiting dyanmics ' , |
460 |
|
C &Nr, 1, myThid ) |
461 |
|
|
462 |
|
|
463 |
RETURN |
RETURN |
464 |
END |
END |