--- MITgcm/model/src/dynamics.F 2001/11/16 03:25:40 1.84 +++ MITgcm/model/src/dynamics.F 2003/02/11 04:05:32 1.93 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.84 2001/11/16 03:25:40 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.93 2003/02/11 04:05:32 jmc Exp $ C $Name: $ #include "CPP_OPTIONS.h" @@ -80,12 +80,10 @@ # include "tamc.h" # include "tamc_keys.h" # include "FFIELDS.h" +# include "EOS.h" # ifdef ALLOW_KPP # include "KPP.h" # endif -# ifdef ALLOW_GMREDI -# include "GMREDI.h" -# endif #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_TIMEAVE #include "TIMEAVE_STATV.h" @@ -100,6 +98,8 @@ C | C |-- CALC_PHI_HYD C | +C |-- STORE_PRESSURE +C | C |-- MOM_FLUXFORM C | C |-- MOM_VECINV @@ -131,13 +131,14 @@ C so we need an fVer for each C variable. C rhoK, rhoKM1 - Density at current level, and level above -C phiHyd - Hydrostatic part of the potential phiHydi. -C In z coords phiHydiHyd is the hydrostatic +C phiHyd - Hydrostatic part of the potential. +C In z coords phiHyd is the hydrostatic C Potential (=pressure/rho0) anomaly -C In p coords phiHydiHyd is the geopotential +C In p coords phiHyd is the geopotential C surface height anomaly. -C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean) -C phiSurfY or geopotentiel (atmos) in X and Y direction +C dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential +C phiSurfX, - gradient of Surface potential (Pressure/rho, ocean) +C phiSurfY or geopotential (atmos) in X and Y direction C iMin, iMax - Ranges and sub-block indices on which calculations C jMin, jMax are applied. C bi, bj @@ -147,15 +148,14 @@ _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly) + _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) - _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) INTEGER iMin, iMax INTEGER jMin, jMax @@ -163,11 +163,8 @@ INTEGER i, j INTEGER k, km1, kp1, kup, kDown -Cjmc : add for phiHyd output <- but not working if multi tile per CPU -c CHARACTER*(MAX_LEN_MBUF) suff -c LOGICAL DIFFERENT_MULTIPLE -c EXTERNAL DIFFERENT_MULTIPLE -Cjmc(end) + LOGICAL DIFFERENT_MULTIPLE + EXTERNAL DIFFERENT_MULTIPLE C--- The algorithm... C @@ -221,14 +218,6 @@ C uninitialised but inert locations. DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx - DO k=1,Nr - phiHyd(i,j,k) = 0. _d 0 - KappaRU(i,j,k) = 0. _d 0 - KappaRV(i,j,k) = 0. _d 0 - sigmaX(i,j,k) = 0. _d 0 - sigmaY(i,j,k) = 0. _d 0 - sigmaR(i,j,k) = 0. _d 0 - ENDDO rhoKM1 (i,j) = 0. _d 0 rhok (i,j) = 0. _d 0 phiSurfX(i,j) = 0. _d 0 @@ -236,6 +225,13 @@ ENDDO ENDDO +C-- Call to routine for calculation of +C Eliassen-Palm-flux-forced U-tendency, +C if desired: +#ifdef INCLUDE_EP_FORCING_CODE + CALL CALC_EP_FORCING(myThid) +#endif + #ifdef ALLOW_AUTODIFF_TAMC C-- HPF directive to help TAMC CHPF$ INDEPENDENT @@ -261,7 +257,7 @@ act3 = myThid - 1 max3 = nTx*nTy act4 = ikey_dynamics - 1 - ikey = (act1 + 1) + act2*max1 + idynkey = (act1 + 1) + act2*max1 & + act3*max1*max2 & + act4*max1*max2*max3 #endif /* ALLOW_AUTODIFF_TAMC */ @@ -269,21 +265,29 @@ C-- Set up work arrays that need valid initial values DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx + DO k=1,Nr + phiHyd(i,j,k) = 0. _d 0 + KappaRU(i,j,k) = 0. _d 0 + KappaRV(i,j,k) = 0. _d 0 + ENDDO fVerU (i,j,1) = 0. _d 0 fVerU (i,j,2) = 0. _d 0 fVerV (i,j,1) = 0. _d 0 fVerV (i,j,2) = 0. _d 0 + dPhiHydX(i,j) = 0. _d 0 + dPhiHydY(i,j) = 0. _d 0 ENDDO ENDDO C-- Start computation of dynamics - iMin = 1-OLx+2 - iMax = sNx+OLx-1 - jMin = 1-OLy+2 - jMax = sNy+OLy-1 + iMin = 0 + iMax = sNx+1 + jMin = 0 + jMax = sNy+1 #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte +CADJ STORE wvel (:,:,:,bi,bj) = +CADJ & comlev1_bibj, key = idynkey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C-- Explicit part of the Surface Potentiel Gradient (add in TIMESTEP) @@ -297,11 +301,11 @@ ENDIF #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte -CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte +CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte +CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte #ifdef ALLOW_KPP CADJ STORE KPPviscAz (:,:,:,bi,bj) -CADJ & = comlev1_bibj, key=ikey, byte=isbyte +CADJ & = comlev1_bibj, key=idynkey, byte=isbyte #endif /* ALLOW_KPP */ #endif /* ALLOW_AUTODIFF_TAMC */ @@ -328,7 +332,9 @@ kDown= 1+MOD(k,2) #ifdef ALLOW_AUTODIFF_TAMC - kkey = (ikey-1)*Nr + k + kkey = (idynkey-1)*Nr + k +CADJ STORE pressure(:,:,k,bi,bj) = comlev1_bibj_k , +CADJ & key=kkey , byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C-- Integrate hydrostatic balance for phiHyd with BC of @@ -339,35 +345,41 @@ I bi,bj,iMin,iMax,jMin,jMax,k, I gT, gS, U phiHyd, - I myThid ) + O dPhiHydX, dPhiHydY, + I myTime, myIter, myThid ) ELSE CALL CALC_PHI_HYD( I bi,bj,iMin,iMax,jMin,jMax,k, I theta, salt, U phiHyd, - I myThid ) + O dPhiHydX, dPhiHydY, + I myTime, myIter, myThid ) ENDIF +C calculate pressure from phiHyd and store it on common block +C variable pressure + CALL STORE_PRESSURE( bi, bj, k, phiHyd, myThid ) + C-- Calculate accelerations in the momentum equations (gU, gV, ...) C and step forward storing the result in gUnm1, gVnm1, etc... IF ( momStepping ) THEN #ifndef DISABLE_MOM_FLUXFORM IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM( I bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown, - I phiHyd,KappaRU,KappaRV, + I phiHyd,dPhiHydX,dPhiHydY,KappaRU,KappaRV, U fVerU, fVerV, I myTime, myIter, myThid) #endif #ifndef DISABLE_MOM_VECINV IF (vectorInvariantMomentum) CALL MOM_VECINV( I bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown, - I phiHyd,KappaRU,KappaRV, + I dPhiHydX,dPhiHydY,KappaRU,KappaRV, U fVerU, fVerV, I myTime, myIter, myThid) #endif CALL TIMESTEP( I bi,bj,iMin,iMax,jMin,jMax,k, - I phiHyd, phiSurfX, phiSurfY, + I phiHyd, dPhiHydX,dPhiHydY, phiSurfX, phiSurfY, I myIter, myThid) #ifdef ALLOW_OBCS @@ -394,13 +406,10 @@ C-- end of dynamics k loop (1:Nr) ENDDO - - C-- Implicit viscosity IF (implicitViscosity.AND.momStepping) THEN #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 3 -CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, @@ -408,8 +417,7 @@ U gUNm1, I myThid ) #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 4 -CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, @@ -428,8 +436,7 @@ #ifdef INCLUDE_CD_CODE #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 5 -CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, @@ -437,8 +444,7 @@ U vVelD, I myThid ) #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 6 -CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, @@ -449,13 +455,12 @@ C-- End If implicitViscosity.AND.momStepping ENDIF -Cjmc : add for phiHyd output <- but not working if multi tile per CPU -c IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime) -c & .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN -c WRITE(suff,'(I10.10)') myIter+1 -c CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid) -c ENDIF -Cjmc(end) +C- jmc: add for diagnostic of phiHyd + IF ( DIFFERENT_MULTIPLE(diagFreq,myTime+deltaTClock,myTime) + & .AND. buoyancyRelation .NE. 'OCEANIC' ) THEN + CALL WRITE_LOCAL_RL('Ph','I10',Nr,phiHyd, + & bi,bj,1,myIter+1,myThid) + ENDIF #ifdef ALLOW_TIMEAVE IF (taveFreq.GT.0.) THEN @@ -467,6 +472,14 @@ ENDDO ENDDO +Cml( +C In order to compare the variance of phiHydLow of a p/z-coordinate +C run with etaH of a z/p-coordinate run the drift of phiHydLow +C has to be removed by something like the following subroutine: +C CALL REMOVE_MEAN_RL( 1, phiHydLow, maskH, maskH, rA, drF, +C & 'phiHydLow', myThid ) +Cml) + #ifndef DISABLE_DEBUGMODE If (debugMode) THEN CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)