--- MITgcm/model/src/dynamics.F 1998/10/28 03:11:37 1.36 +++ MITgcm/model/src/dynamics.F 1999/05/24 15:42:23 1.43 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.36 1998/10/28 03:11:37 cnh Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.43 1999/05/24 15:42:23 adcroft Exp $ #include "CPP_OPTIONS.h" @@ -20,6 +20,7 @@ C | C*P* comments indicating place holders for which code is | C | presently being developed. | C \==========================================================/ + IMPLICIT NONE C == Global variables === #include "SIZE.h" @@ -27,6 +28,10 @@ #include "CG2D.h" #include "PARAMS.h" #include "DYNVARS.h" +#include "GRID.h" +#ifdef ALLOW_KPP +#include "KPPMIX.h" +#endif C == Routine arguments == C myTime - Current time in simulation @@ -38,12 +43,13 @@ C == Local variables C xA, yA - Per block temporaries holding face areas -C uTrans, vTrans, rTrans - Per block temporaries holding flow transport +C uTrans, vTrans, rTrans - Per block temporaries holding flow +C transport C rVel o uTrans: Zonal transport C o vTrans: Meridional transport C o rTrans: Vertical transport -C o rVel: Vertical velocity at upper and lower -C cell faces. +C o rVel: Vertical velocity at upper and +C lower cell faces. C maskC,maskUp o maskC: land/water mask for tracer cells C o maskUp: land/water mask for W points C aTerm, xTerm, cTerm - Work arrays for holding separate terms in @@ -59,12 +65,15 @@ C is "pipelined" in the vertical C so we need an fVer for each C variable. -C rhoK, rhoKM1 - Density at current level, level above and level below. +C rhoK, rhoKM1 - Density at current level, level above and level +C below. C rhoKP1 C buoyK, buoyKM1 - Buoyancy at current level and level above. C phiHyd - Hydrostatic part of the potential phiHydi. -C In z coords phiHydiHyd is the hydrostatic pressure anomaly -C In p coords phiHydiHyd is the geopotential surface height +C In z coords phiHydiHyd is the hydrostatic +C pressure anomaly +C In p coords phiHydiHyd is the geopotential +C surface height C anomaly. C etaSurfX, - Holds surface elevation gradient in X and Y. C etaSurfY @@ -72,13 +81,13 @@ C diffusion tensor. C KapGM - Spatially varying Visbeck et. al mixing coeff. C KappaRT, - Total diffusion in vertical for T and S. -C KappaRS ( background + spatially varying, isopycnal term). +C KappaRS (background + spatially varying, isopycnal term). C iMin, iMax - Ranges and sub-block indices on which calculations C jMin, jMax are applied. C bi, bj C k, kUp, - Index for layer above and below. kUp and kDown -C kDown, kM1 are switched with layer to be the appropriate index -C into fVerTerm +C kDown, kM1 are switched with layer to be the appropriate +C index into fVerTerm. _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -113,6 +122,8 @@ _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) + _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) + _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) INTEGER iMin, iMax INTEGER jMin, jMax @@ -237,6 +248,7 @@ K = 1 BOTTOM_LAYER = K .EQ. Nr +#ifdef DO_PIPELINED_CORRECTION_STEP C-- Calculate gradient of surface pressure CALL CALC_GRAD_ETA_SURF( I bi,bj,iMin,iMax,jMin,jMax, @@ -246,40 +258,60 @@ CALL CORRECTION_STEP( I bi,bj,iMin,iMax,jMin,jMax,K, I etaSurfX,etaSurfY,myTime,myThid) +#ifdef ALLOW_OBCS + IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K, myThid ) +#endif IF ( .NOT. BOTTOM_LAYER ) THEN C-- Update fields in layer below according to tendency terms CALL CORRECTION_STEP( I bi,bj,iMin,iMax,jMin,jMax,K+1, I etaSurfX,etaSurfY,myTime,myThid) +#ifdef ALLOW_OBCS + IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K+1, myThid ) +#endif ENDIF +#endif C-- Density of 1st level (below W(1)) reference to level 1 +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType, O rhoKm1, I myThid ) +#endif - IF ( .NOT. BOTTOM_LAYER ) THEN + IF ( (.NOT. BOTTOM_LAYER) +#ifdef ALLOW_KPP + & .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing +#endif + & ) THEN C-- Check static stability with layer below C-- and mix as needed. +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType, O rhoKp1, I myThid ) +#endif +#ifdef INCLUDE_CONVECT_CALL CALL CONVECT( I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1, I myTime,myIter,myThid) +#endif C-- Recompute density after mixing +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType, O rhoKm1, I myThid ) +#endif ENDIF C-- Calculate buoyancy CALL CALC_BUOYANCY( I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1, O buoyKm1, I myThid ) -C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0 +C-- Integrate hydrostatic balance for phiHyd with BC of +C-- phiHyd(z=0)=0 CALL CALC_PHI_HYD( I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1, U phiHyd, @@ -287,61 +319,97 @@ DO K=2,Nr BOTTOM_LAYER = K .EQ. Nr +#ifdef DO_PIPELINED_CORRECTION_STEP IF ( .NOT. BOTTOM_LAYER ) THEN C-- Update fields in layer below according to tendency terms CALL CORRECTION_STEP( I bi,bj,iMin,iMax,jMin,jMax,K+1, I etaSurfX,etaSurfY,myTime,myThid) +#ifdef ALLOW_OBCS + IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K+1, myThid ) +#endif ENDIF +#endif C-- Density of K level (below W(K)) reference to K level +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType, O rhoK, I myThid ) - IF ( .NOT. BOTTOM_LAYER ) THEN +#endif + IF ( (.NOT. BOTTOM_LAYER) +#ifdef ALLOW_KPP + & .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing +#endif + & ) THEN C-- Check static stability with layer below and mix as needed. C-- Density of K+1 level (below W(K+1)) reference to K level. +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType, O rhoKp1, I myThid ) +#endif +#ifdef INCLUDE_CONVECT_CALL CALL CONVECT( I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1, I myTime,myIter,myThid) +#endif C-- Recompute density after mixing +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType, O rhoK, I myThid ) +#endif ENDIF C-- Calculate buoyancy CALL CALC_BUOYANCY( I bi,bj,iMin,iMax,jMin,jMax,K,rhoK, O buoyK, I myThid ) -C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0 +C-- Integrate hydrostatic balance for phiHyd with BC of +C-- phiHyd(z=0)=0 CALL CALC_PHI_HYD( I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK, U phiHyd, I myThid ) C-- Calculate iso-neutral slopes for the GM/Redi parameterisation +#ifdef INCLUDE_FIND_RHO_CALL CALL FIND_RHO( I bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType, O rhoTmp, I myThid ) +#endif +#ifdef INCLUDE_CALC_ISOSLOPES_CALL CALL CALC_ISOSLOPES( I bi, bj, iMin, iMax, jMin, jMax, K, I rhoKm1, rhoK, rhotmp, O K13, K23, K33, KapGM, I myThid ) +#endif DO J=jMin,jMax DO I=iMin,iMax +#ifdef INCLUDE_FIND_RHO_CALL rhoKm1 (I,J) = rhoK(I,J) +#endif buoyKm1(I,J) = buoyK(I,J) ENDDO ENDDO ENDDO ! K +#ifdef ALLOW_KPP +C-- Compute KPP mixing coefficients + IF (usingKPPmixing) THEN + CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]' + I , myThid) + CALL KVMIX( + I bi, bj, myTime, myThid ) + CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]' + I , myThid) + ENDIF +#endif + DO K = Nr, 1, -1 kM1 =max(1,k-1) ! Points to level above k (=k-1) @@ -357,21 +425,23 @@ I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown, O xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp, I myThid) +#ifdef INCLUDE_CALC_DIFFUSIVITY_CALL C-- Calculate the total vertical diffusivity CALL CALC_DIFFUSIVITY( I bi,bj,iMin,iMax,jMin,jMax,K, I maskC,maskUp,KapGM,K33, - O KappaRT,KappaRS, + O KappaRT,KappaRS,KappaRU,KappaRV, I myThid) +#endif C-- Calculate accelerations in the momentum equations IF ( momStepping ) THEN CALL CALC_MOM_RHS( I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown, I xA,yA,uTrans,vTrans,rTrans,rVel,maskC, - I phiHyd, + I phiHyd,KappaRU,KappaRV, U aTerm,xTerm,cTerm,mTerm,pTerm, U fZon, fMer, fVerU, fVerV, - I myThid) + I myTime, myThid) ENDIF C-- Calculate active tracer tendencies IF ( tempStepping ) THEN @@ -380,7 +450,7 @@ I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC, I K13,K23,KappaRT,KapGM, U aTerm,xTerm,fZon,fMer,fVerT, - I myThid) + I myTime, myThid) ENDIF IF ( saltStepping ) THEN CALL CALC_GS( @@ -388,12 +458,19 @@ I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC, I K13,K23,KappaRS,KapGM, U aTerm,xTerm,fZon,fMer,fVerS, - I myThid) + I myTime, myThid) ENDIF C-- Prediction step (step forward all model variables) CALL TIMESTEP( I bi,bj,iMin,iMax,jMin,jMax,K, I myThid) +#ifdef ALLOW_OBCS +C-- Apply open boundary conditions + IF (openBoundaries) CALL APPLY_OBCS2( bi, bj, K, myThid ) +#endif +C-- Freeze water + IF (allowFreezing) + & CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid ) C-- Diagnose barotropic divergence of predicted fields CALL CALC_DIV_GHAT( I bi,bj,iMin,iMax,jMin,jMax,K, @@ -401,7 +478,7 @@ I myThid) C-- Cumulative diagnostic calculations (ie. time-averaging) -#ifdef ALLOW_DIAGNOSTICS +#ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE IF (taveFreq.GT.0.) THEN CALL DO_TIME_AVERAGES( I myTime, myIter, bi, bj, K, kUp, kDown, @@ -414,10 +491,41 @@ C-- Implicit diffusion IF (implicitDiffusion) THEN - CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax, - I KappaRT,KappaRS, - I myThid ) - ENDIF + IF (tempStepping) CALL IMPLDIFF( + I bi, bj, iMin, iMax, jMin, jMax, + I deltaTtracer, KappaRT,recip_HFacC, + U gTNm1, + I myThid ) + IF (saltStepping) CALL IMPLDIFF( + I bi, bj, iMin, iMax, jMin, jMax, + I deltaTtracer, KappaRS,recip_HFacC, + U gSNm1, + I myThid ) + IF (momStepping) THEN + CALL IMPLDIFF( + I bi, bj, iMin, iMax, jMin, jMax, + I deltaTmom, KappaRU,recip_HFacW, + U gUNm1, + I myThid ) + CALL IMPLDIFF( + I bi, bj, iMin, iMax, jMin, jMax, + I deltaTmom, KappaRV,recip_HFacS, + U gVNm1, + I myThid ) +#ifdef INCLUDE_CD_CODE + CALL IMPLDIFF( + I bi, bj, iMin, iMax, jMin, jMax, + I deltaTmom, KappaRU,recip_HFacW, + U vVelD, + I myThid ) + CALL IMPLDIFF( + I bi, bj, iMin, iMax, jMin, jMax, + I deltaTmom, KappaRV,recip_HFacS, + U uVelD, + I myThid ) +#endif + ENDIF ! momStepping + ENDIF ! implicitDiffusion ENDDO ENDDO @@ -458,6 +566,8 @@ C &Nr, 1, myThid ) C CALL PLOT_FIELD_XYZRL( gT, ' GT exiting dyanmics ' , C &Nr, 1, myThid ) +C CALL PLOT_FIELD_XYZRL( phiHyd, ' phiHyd exiting dyanmics ' , +C &Nr, 1, myThid ) RETURN