--- MITgcm/model/src/dynamics.F 2002/02/26 16:04:47 1.83.2.1 +++ MITgcm/model/src/dynamics.F 2005/04/06 18:29:53 1.114 @@ -1,6 +1,7 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.83.2.1 2002/02/26 16:04:47 adcroft Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.114 2005/04/06 18:29:53 jmc Exp $ C $Name: $ +#include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP @@ -72,24 +73,19 @@ #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" -#include "GRID.h" -#ifdef ALLOW_PASSIVE_TRACER -#include "TR1.h" +#ifdef ALLOW_CD_CODE +#include "CD_CODE_VARS.h" #endif +#include "GRID.h" #ifdef ALLOW_AUTODIFF_TAMC # 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" -#endif C !CALLING SEQUENCE: C DYNAMICS() @@ -112,7 +108,6 @@ C | C |-- OBCS_APPLY_UV C | -C |-- CALL TIMEAVE_CUMUL_1T C |-- CALL DEBUG_STATS_RL C !INPUT/OUTPUT PARAMETERS: @@ -126,18 +121,20 @@ C !LOCAL VARIABLES: C == Local variables -C fVer[STUV] o fVer: Vertical flux term - note fVer -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 fVer[UV] o fVer: Vertical flux term - note fVer +C is "pipelined" in the vertical +C so we need an fVer for each +C variable. +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 guDissip :: dissipation tendency (all explicit terms), u component +C gvDissip :: dissipation tendency (all explicit terms), v component C iMin, iMax - Ranges and sub-block indices on which calculations C jMin, jMax are applied. C bi, bj @@ -146,16 +143,16 @@ 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 guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL gvDissip(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 +160,12 @@ 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) +#ifdef ALLOW_DIAGNOSTICS + _RL tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + LOGICAL DIAGNOSTICS_IS_ON + EXTERNAL DIAGNOSTICS_IS_ON +#endif /* ALLOW_DIAGNOSTICS */ + C--- The algorithm... C @@ -214,27 +212,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 - 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 - 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 @@ -246,7 +229,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 */ @@ -261,29 +244,58 @@ 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 +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 +cph) + 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 + guDissip(i,j) = 0. _d 0 + gvDissip(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 +309,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 */ @@ -315,6 +327,13 @@ ENDDO #endif +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE KappaRU(:,:,:) +CADJ & = comlev1_bibj, key=idynkey, byte=isbyte +CADJ STORE KappaRV(:,:,:) +CADJ & = comlev1_bibj, key=idynkey, byte=isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ + C-- Start of dynamics loop DO k=1,Nr @@ -328,147 +347,167 @@ 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 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 C phiHyd(z=0)=0 -C distinguishe between Stagger and Non Stagger time stepping - IF (staggerTimeStep) THEN - CALL CALC_PHI_HYD( - I bi,bj,iMin,iMax,jMin,jMax,k, - I gT, gS, - U phiHyd, - I myThid ) - ELSE - CALL CALC_PHI_HYD( + CALL CALC_PHI_HYD( I bi,bj,iMin,iMax,jMin,jMax,k, I theta, salt, - U phiHyd, - I myThid ) - ENDIF + U phiHydF, + O phiHydC, dPhiHydX, dPhiHydY, + I myTime, myIter, myThid ) 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 +#ifdef ALLOW_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 +#ifdef ALLOW_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, + O guDissip, gvDissip, 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 guDissip, gvDissip, + 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 +C-- Implicit Vertical advection & viscosity +#ifdef INCLUDE_IMPLVERTADV_CODE + IF ( momImplVertAdv ) THEN + CALL MOM_U_IMPLICIT_R( kappaRU, + I bi, bj, myTime, myIter, myThid ) + CALL MOM_V_IMPLICIT_R( kappaRV, + I bi, bj, myTime, myIter, myThid ) + ELSEIF ( implicitViscosity ) THEN +#else /* INCLUDE_IMPLVERTADV_CODE */ + IF ( implicitViscosity ) THEN +#endif /* INCLUDE_IMPLVERTADV_CODE */ #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 3 -CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, 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, + I 0, KappaRU,recip_HFacW, + U gU, I myThid ) #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 4 -CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, 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, + I 0, KappaRV,recip_HFacS, + U gV, I myThid ) + ENDIF #ifdef ALLOW_OBCS C-- Apply open boundary conditions - IF (useOBCS) THEN + IF ( useOBCS .AND.(implicitViscosity.OR.momImplVertAdv) ) 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 #endif /* ALLOW_OBCS */ -#ifdef INCLUDE_CD_CODE +#ifdef ALLOW_CD_CODE + IF (implicitViscosity.AND.useCDscheme) THEN #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, - I deltaTmom, KappaRU,recip_HFacW, + I 0, KappaRU,recip_HFacW, 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, - I deltaTmom, KappaRV,recip_HFacS, + I 0, KappaRV,recip_HFacS, U uVelD, I myThid ) -#endif /* INCLUDE_CD_CODE */ -C-- End If implicitViscosity.AND.momStepping ENDIF +#endif /* ALLOW_CD_CODE */ +C-- End implicit Vertical advection & viscosity -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 -#ifndef DISABLE_DEBUGMODE - If (debugMode) THEN +#ifdef ALLOW_OBCS + IF (useOBCS) THEN + CALL OBCS_PRESCRIBE_EXCHANGES(myThid) + ENDIF +#endif + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + +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) + +#ifdef ALLOW_DIAGNOSTICS + IF ( usediagnostics ) THEN + + CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD ',0,Nr,0,1,1,myThid) + CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT ',0,1,0,1,1,myThid) + + IF ( DIAGNOSTICS_IS_ON('PHIBOTSQ',myThid) ) THEN + DO bj = myByLo(myThid), myByHi(myThid) + DO bi = myBxLo(myThid), myBxHi(myThid) + DO j = 1,sNy + DO i = 1,sNx + tmpFld(i,j) = phiHydLow(i,j,bi,bj)*phiHydLow(i,j,bi,bj) + ENDDO + ENDDO + CALL DIAGNOSTICS_FILL(tmpFld,'PHIBOTSQ',0,1,2,bi,bj,myThid) + ENDDO + ENDDO + ENDIF + + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + +#ifdef ALLOW_DEBUG + 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)