--- MITgcm/model/src/dynamics.F 2002/05/30 02:30:12 1.87 +++ MITgcm/model/src/dynamics.F 2003/10/09 04:19:18 1.100 @@ -1,6 +1,7 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.87 2002/05/30 02:30:12 heimbach Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.100 2003/10/09 04:19:18 edhill Exp $ C $Name: $ +#include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP @@ -80,13 +81,11 @@ # include "tamc.h" # include "tamc_keys.h" # include "FFIELDS.h" +# include "EOS.h" # ifdef ALLOW_KPP # include "KPP.h" # endif #endif /* ALLOW_AUTODIFF_TAMC */ -#ifdef ALLOW_TIMEAVE -#include "TIMEAVE_STATV.h" -#endif C !CALLING SEQUENCE: C DYNAMICS() @@ -127,14 +126,14 @@ C is "pipelined" in the vertical 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 Potential (=pressure/rho0) anomaly -C In p coords phiHydiHyd 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 phiHydC :: hydrostatic potential anomaly at cell center +C In z coords phiHyd is the hydrostatic potential +C (=pressure/rho0) anomaly +C In p coords phiHyd is the geopotential height anomaly. +C phiHydF :: hydrostatic potential anomaly at middle between 2 centers +C dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom. +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 @@ -143,9 +142,10 @@ C index into fVerTerm. _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 rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly) + _RL dPhiHydY(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) @@ -157,11 +157,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 @@ -208,19 +205,12 @@ C--- CEOP -C-- Set up work arrays with valid (i.e. not NaN) values -C These inital values do not alter the numerical results. They -C just ensure that all memory references are to valid floating -C point numbers. This prevents spurious hardware signals due to -C uninitialised but inert locations. - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - rhoKM1 (i,j) = 0. _d 0 - rhok (i,j) = 0. _d 0 - phiSurfX(i,j) = 0. _d 0 - phiSurfY(i,j) = 0. _d 0 - 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 @@ -232,7 +222,7 @@ #ifdef ALLOW_AUTODIFF_TAMC C-- HPF directive to help TAMC CHPF$ INDEPENDENT, NEW (fVerU,fVerV -CHPF$& ,phiHyd +CHPF$& ,phiHydF CHPF$& ,KappaRU,KappaRV CHPF$& ) #endif /* ALLOW_AUTODIFF_TAMC */ @@ -247,34 +237,59 @@ 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 */ -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 +C-- Set up work arrays with valid (i.e. not NaN) values +C These inital values do not alter the numerical results. They +C just ensure that all memory references are to valid floating +C point numbers. This prevents spurious hardware signals due to +C uninitialised but inert locations. + + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx KappaRU(i,j,k) = 0. _d 0 KappaRV(i,j,k) = 0. _d 0 +#ifdef ALLOW_AUTODIFF_TAMC +cph( +c-- need some re-initialisation here to break dependencies +c-- totphihyd is assumed zero from ini_pressure, i.e. +c-- avoiding iterate pressure p = integral of (g*rho(p)*dz) +cph) + totPhiHyd(i,j,k,bi,bj) = 0. _d 0 + gu(i,j,k,bi,bj) = 0. _d 0 + gv(i,j,k,bi,bj) = 0. _d 0 +#endif ENDDO + ENDDO + ENDDO + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx 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 + phiHydF (i,j) = 0. _d 0 + phiHydC (i,j) = 0. _d 0 + dPhiHydX(i,j) = 0. _d 0 + dPhiHydY(i,j) = 0. _d 0 + phiSurfX(i,j) = 0. _d 0 + phiSurfY(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) @@ -288,11 +303,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 */ @@ -319,7 +334,18 @@ kDown= 1+MOD(k,2) #ifdef ALLOW_AUTODIFF_TAMC - kkey = (ikey-1)*Nr + k + kkey = (idynkey-1)*Nr + k +c +CADJ STORE totphihyd (:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gt (:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gs (:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE theta (:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE salt (:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C-- Integrate hydrostatic balance for phiHyd with BC of @@ -329,98 +355,84 @@ CALL CALC_PHI_HYD( I bi,bj,iMin,iMax,jMin,jMax,k, I gT, gS, - U phiHyd, - I myThid ) + U phiHydF, + O phiHydC, 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 ) + U phiHydF, + O phiHydC, dPhiHydX, dPhiHydY, + I myTime, myIter, myThid ) ENDIF C-- Calculate accelerations in the momentum equations (gU, gV, ...) -C and step forward storing the result in gUnm1, gVnm1, etc... +C and step forward storing the result in gU, gV, 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 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 myIter, myThid) + I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY, + I myTime, myIter, myThid) #ifdef ALLOW_OBCS C-- Apply open boundary conditions - IF (useOBCS) THEN - CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid ) - END IF + IF (useOBCS) THEN + CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid ) + ENDIF #endif /* ALLOW_OBCS */ -#ifdef ALLOW_AUTODIFF_TAMC -#ifdef INCLUDE_CD_CODE - ELSE - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - guCD(i,j,k,bi,bj) = 0.0 - gvCD(i,j,k,bi,bj) = 0.0 - END DO - END DO -#endif /* INCLUDE_CD_CODE */ -#endif /* ALLOW_AUTODIFF_TAMC */ ENDIF 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 gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, I deltaTmom, KappaRU,recip_HFacW, - U gUNm1, + U gU, I myThid ) #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 4 -CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, I deltaTmom, KappaRV,recip_HFacS, - U gVNm1, + U gV, I myThid ) #ifdef ALLOW_OBCS C-- Apply open boundary conditions IF (useOBCS) THEN DO K=1,Nr - CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid ) + CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid ) ENDDO END IF #endif /* ALLOW_OBCS */ #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, @@ -428,8 +440,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, @@ -440,26 +451,19 @@ 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) - -#ifdef ALLOW_TIMEAVE - IF (taveFreq.GT.0.) THEN - CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr, - I deltaTclock, bi, bj, myThid) - ENDIF -#endif /* ALLOW_TIMEAVE */ - 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 + If ( debugLevel .GE. debLevB ) THEN CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid) CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid) CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)