--- MITgcm/model/src/dynamics.F 2005/05/23 20:49:37 1.117 +++ MITgcm/model/src/dynamics.F 2005/12/15 21:09:00 1.127 @@ -1,8 +1,9 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.117 2005/05/23 20:49:37 molod Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.127 2005/12/15 21:09:00 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" +#undef DYNAMICS_GUGV_EXCH_CHECK CBOP C !ROUTINE: DYNAMICS @@ -30,6 +31,7 @@ C | with the free-surface evolution or the rigid-lid: C | U[n] = U* + dt x d/dx P C | V[n] = V* + dt x d/dy P +C | W[n] = W* + dt x d/dz P (NH mode) C | C | "Calculation of Gs" C | =================== @@ -90,6 +92,8 @@ C !CALLING SEQUENCE: C DYNAMICS() C | +C |-- CALC_EP_FORCING +C | C |-- CALC_GRAD_PHI_SURF C | C |-- CALC_VISCOSITY @@ -104,11 +108,17 @@ C | C |-- OBCS_APPLY_UV C | +C |-- MOM_U_IMPLICIT_R +C |-- MOM_V_IMPLICIT_R +C | C |-- IMPLDIFF C | C |-- OBCS_APPLY_UV C | -C |-- CALL DEBUG_STATS_RL +C |-- CALC_GW +C | +C |-- DIAGNOSTICS_FILL +C |-- DEBUG_STATS_RL C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == @@ -161,9 +171,7 @@ INTEGER k, km1, kp1, kup, kDown #ifdef ALLOW_DIAGNOSTICS - _RL tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) - LOGICAL DIAGNOSTICS_IS_ON - EXTERNAL DIAGNOSTICS_IS_ON + _RL tmpFac #endif /* ALLOW_DIAGNOSTICS */ @@ -212,6 +220,11 @@ C--- CEOP +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_ENTER( 'DYNAMICS', myThid ) +#endif + C-- Call to routine for calculation of C Eliassen-Palm-flux-forced U-tendency, C if desired: @@ -264,8 +277,8 @@ 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 + gU(i,j,k,bi,bj) = 0. _d 0 + gV(i,j,k,bi,bj) = 0. _d 0 #endif ENDDO ENDDO @@ -355,6 +368,39 @@ CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE salt (:,:,k,bi,bj) CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +# ifdef NONLIN_FRSURF +cph-test +CADJ STORE phiHydC (:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE phiHydF (:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gudissip (:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gvdissip (:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE fVerU (:,:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE fVerV (:,:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gu(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gv(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gunm1(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE gvnm1(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +# ifdef ALLOW_CD_CODE +CADJ STORE unm1(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE vnm1(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE uVelD(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE vVelD(:,:,k,bi,bj) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +# endif +# endif #endif /* ALLOW_AUTODIFF_TAMC */ C-- Integrate hydrostatic balance for phiHyd with BC of @@ -372,17 +418,30 @@ #ifdef ALLOW_MOM_FLUXFORM IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM( I bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown, - I dPhiHydX,dPhiHydY,KappaRU,KappaRV, + I KappaRU, KappaRV, U fVerU, fVerV, + O guDissip, gvDissip, I myTime, myIter, myThid) #endif #ifdef ALLOW_MOM_VECINV - IF (vectorInvariantMomentum) CALL MOM_VECINV( + IF (vectorInvariantMomentum) THEN +C +# ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +CADJ STORE fVerU(:,:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE fVerV(:,:,:) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +# endif +# endif /* ALLOW_AUTODIFF_TAMC */ +C + CALL MOM_VECINV( I bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown, - I dPhiHydX,dPhiHydY,KappaRU,KappaRV, + I KappaRU, KappaRV, U fVerU, fVerV, O guDissip, gvDissip, I myTime, myIter, myThid) + ENDIF #endif CALL TIMESTEP( I bi,bj,iMin,iMax,jMin,jMax,k, @@ -420,7 +479,7 @@ #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, - I 0, KappaRU,recip_HFacW, + I -1, KappaRU,recip_HFacW, U gU, I myThid ) #ifdef ALLOW_AUTODIFF_TAMC @@ -429,7 +488,7 @@ #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, - I 0, KappaRV,recip_HFacS, + I -2, KappaRV,recip_HFacS, U gV, I myThid ) ENDIF @@ -476,6 +535,22 @@ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +#ifdef ALLOW_NONHYDROSTATIC +C-- Step forward W field in N-H algorithm + IF ( momStepping .AND. nonHydrostatic ) THEN +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('CALC_GW', myThid ) +#endif + CALL TIMER_START('CALC_GW [DYNAMICS]',myThid) + CALL CALC_GW( myTime, myIter, myThid ) + CALL TIMESTEP_WVEL( myTime, myIter, myThid ) + CALL TIMER_STOP ('CALC_GW [DYNAMICS]',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 @@ -488,36 +563,14 @@ 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('PHIHYDSQ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k = 1,Nr - DO j = 1,sNy - DO i = 1,sNx - tmpFld(i,j) = totPhihyd(i,j,k,bi,bj)*totPhihyd(i,j,k,bi,bj) - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpFld,'PHIHYDSQ',k,1,2,bi,bj,myThid) - ENDDO - ENDDO - ENDDO - ENDIF - - CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT ',0,1,0,1,1,myThid) + tmpFac = 1. _d 0 + CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2, + & 'PHIHYDSQ',0,Nr,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 + CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2, + & 'PHIBOTSQ',0, 1,0,1,1,myThid) ENDIF #endif /* ALLOW_DIAGNOSTICS */ @@ -543,5 +596,19 @@ ENDIF #endif +#ifdef DYNAMICS_GUGV_EXCH_CHECK +C- jmc: For safety checking only: This Exchange here should not change +C the solution. If solution changes, it means something is wrong, +C but it does not mean that it is less wrong with this exchange. + IF ( debugLevel .GT. debLevB ) THEN + CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid) + ENDIF +#endif + +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_LEAVE( 'DYNAMICS', myThid ) +#endif + RETURN END