diff --git a/ECCOv4 Release 4/flux-forced/code/diagnostics_fill_state.F b/ECCOv4 Release 4/flux-forced/code/diagnostics_fill_state.F deleted file mode 100644 index 0beb559..0000000 --- a/ECCOv4 Release 4/flux-forced/code/diagnostics_fill_state.F +++ /dev/null @@ -1,723 +0,0 @@ -C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_fill_state.F,v 1.43 2015/06/09 16:31:45 rpa Exp $ -C $Name: $ - -#include "DIAG_OPTIONS.h" -#ifdef ALLOW_ECCO -# include "ECCO_OPTIONS.h" -#endif - -CBOP -C !ROUTINE: DIAGNOSTICS_FILL_STATE -C !INTERFACE: - SUBROUTINE DIAGNOSTICS_FILL_STATE( selectVars, myIter, myThid ) - -C !DESCRIPTION: \bv -C *==========================================================* -C | SUBROUTINE DIAGNOSTICS_FILL_STATE -C | o Fill-in main code, state-variables diagnostics -C *==========================================================* -C \ev - -C !USES: - IMPLICIT NONE -C == Global variables === -#include "SIZE.h" -#include "EEPARAMS.h" -#include "PARAMS.h" -#include "GRID.h" -#include "SURFACE.h" -#include "DYNVARS.h" -#include "NH_VARS.h" -#ifdef ALLOW_GENERIC_ADVDIFF -# include "GAD.h" -#endif -#ifdef ALLOW_ECCO -# include "ecco.h" -#endif - - -C !INPUT/OUTPUT PARAMETERS: -C == Routine arguments == -C selectVars :: select which group of dianostics variables to fill-in -C = 1 :: fill-in diagnostics for tracer variables only -C = 2 :: fill-in diagnostics for momentum variables only -C = 3 :: fill-in diagnostics for momentum & tracer variables -C = 4 :: fill-in state variable tendency diagnostics the second time -C myIter :: current Iteration number -C myThid :: my Thread Id number - INTEGER selectVars - INTEGER myIter - INTEGER myThid - -#ifdef ALLOW_DIAGNOSTICS -C !LOCAL VARIABLES: -C == Local variables == - LOGICAL DIAGNOSTICS_IS_ON - EXTERNAL DIAGNOSTICS_IS_ON - _RL tmpMk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) - _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL tmpU (1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL tmpV (1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL tmpFac, uBarC, vBarC -#ifdef ALLOW_FIZHI - _RL dummy1, dummy2, dummy3, dummy4, kappa, getcon -#endif -#ifdef ALLOW_ADAMSBASHFORTH_3 - INTEGER m1 -#endif - INTEGER i,j,k,bi,bj - INTEGER km1 - - IF ( selectVars.EQ.2 .OR. selectVars.EQ.3 ) THEN -C-- fill momentum state-var diagnostics: - - CALL DIAGNOSTICS_FILL(etaN, 'ETAN ',0, 1,0,1,1,myThid) - CALL DIAGNOSTICS_FILL(m_eta,'SSHNOIBC',0, 1,0,1,1,myThid) - CALL DIAGNOSTICS_FILL(m_eta_ib,'SSHIBC ', - & 0, 1,0,1,1,myThid) - CALL DIAGNOSTICS_FILL(m_eta_dyn,'SSH ', - & 0, 1,0,1,1,myThid) - - IF ( DIAGNOSTICS_IS_ON('RSURF ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO j = 1,sNy - DO i = 1,sNx - tmp1k(i,j,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmp1k,'RSURF ',0,1,0,1,1,myThid) - ENDIF - - CALL DIAGNOSTICS_SCALE_FILL( etaN, oneRL, 2, - & 'ETANSQ ',0, 1,0,1,1,myThid) - -#ifdef EXACT_CONSERV - CALL DIAGNOSTICS_SCALE_FILL( dEtaHdt, oneRL, 2, - & 'DETADT2 ',0, 1,0,1,1,myThid) -#endif -#ifdef ALLOW_NONHYDROSTATIC - IF ( use3Dsolver ) THEN - CALL DIAGNOSTICS_FILL( phi_nh,'PHI_NH ',0,Nr,0,1,1,myThid ) - ENDIF -#endif - - CALL DIAGNOSTICS_FILL(uVel, 'UVEL ',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_FILL(vVel, 'VVEL ',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_FILL(wVel, 'WVEL ',0,Nr,0,1,1,myThid) - - CALL DIAGNOSTICS_SCALE_FILL( uVel, oneRL, 2, - & 'UVELSQ ',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_SCALE_FILL( vVel, oneRL, 2, - & 'VVELSQ ',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_SCALE_FILL( wVel, oneRL, 2, - & 'WVELSQ ',0,Nr,0,1,1,myThid) - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - - IF ( DIAGNOSTICS_IS_ON('UE_VEL_C',myThid) .OR. - & DIAGNOSTICS_IS_ON('VN_VEL_C',myThid) .OR. - & DIAGNOSTICS_IS_ON('UV_VEL_C',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 - uBarC = 0.5 _d 0 - & *(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj)) - vBarC = 0.5 _d 0 - & *(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj)) - tmpU(i,j) = angleCosC(i,j,bi,bj)*uBarC - & -angleSinC(i,j,bi,bj)*vBarC - tmpV(i,j) = angleSinC(i,j,bi,bj)*uBarC - & +angleCosC(i,j,bi,bj)*vBarC - tmpMk(i,j,k,bi,bj) = tmpU(i,j)*tmpV(i,j) - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpU,'UE_VEL_C',k,1,2,bi,bj,myThid) - CALL DIAGNOSTICS_FILL(tmpV,'VN_VEL_C',k,1,2,bi,bj,myThid) - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UV_VEL_C',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('UV_VEL_Z',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx+1 - tmpMk(i,j,k,bi,bj) = 0.25 _d 0 - & *(uVel(i,j-1,k,bi,bj)+uVel(i,j,k,bi,bj)) - & *(vVel(i-1,j,k,bi,bj)+vVel(i,j,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UV_VEL_Z',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('WU_VEL ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - km1 = MAX(k-1,1) - DO j = 1,sNy - DO i = 1,sNx+1 - tmpMk(i,j,k,bi,bj) = 0.25 _d 0 - & *(uVel(i,j,km1,bi,bj)+uVel(i,j,k,bi,bj)) - & *(wVel(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj) - & +wVel( i ,j,k,bi,bj)*rA( i ,j,bi,bj) - & )*recip_rAw(i,j,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'WU_VEL ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('WV_VEL ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - km1 = MAX(k-1,1) - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = 0.25 _d 0 - & *(vVel(i,j,km1,bi,bj)+vVel(i,j,k,bi,bj)) - & *(wVel(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj) - & +wVel(i, j ,k,bi,bj)*rA(i, j ,bi,bj) - & )*recip_rAs(i,j,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'WV_VEL ',0,Nr,0,1,1,myThid) - ENDIF - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - - IF ( DIAGNOSTICS_IS_ON('UVELTH ',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+1 - tmpMk(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)*0.5 _d 0 - & *(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UVELTH ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('VVELTH ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)*0.5 _d 0 - & *(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'VVELTH ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('WVELTH ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - km1 = MAX(k-1,1) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*0.5 _d 0 - & *(theta(i,j,k,bi,bj)+theta(i,j,km1,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'WVELTH ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('UVELSLT ',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+1 - tmpMk(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)*0.5 _d 0 - & *(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UVELSLT ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('VVELSLT ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)*0.5 _d 0 - & *(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'VVELSLT ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('WVELSLT ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - km1 = MAX(k-1,1) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*0.5 _d 0 - & *(salt(i,j,k,bi,bj)+salt(i,j,km1,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'WVELSLT ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('UVELPHI ',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+1 - tmpMk(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)*hFacW(i,j,k,bi,bj) - & *0.5 _d 0*(totPhiHyd(i,j,k,bi,bj)+totPhiHyd(i-1,j,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UVELPHI ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('VVELPHI ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)*hFacS(i,j,k,bi,bj) - & *0.5 _d 0*(totPhiHyd(i,j,k,bi,bj)+totPhiHyd(i,j-1,k,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'VVELPHI ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('RCENTER ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO j = 1,sNy - DO i = 1,sNx - tmp1k(i,j,bi,bj) = R_low(i,j,bi,bj) - ENDDO - ENDDO - DO k = Nr,1,-1 - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = tmp1k(i,j,bi,bj) - & + (rF(k+1)-rC(k))*hFacC(i,j,k,bi,bj)*rkSign -C above: more general (setInterFDr/setCenterDr) than line below -c & + drF(k)*hFacC(i,j,k,bi,bj)*0.5 _d 0 - tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj) - & + drF(k)*hFacC(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'RCENTER ',0,Nr,0,1,1,myThid) - ENDIF - -C First fill sequence for state variable tendency diagnostics: subtract state variable -C NOTE: send a '0' for the bibjflag and allow counter to be incremented -C (next fill for these diagnostics will NOT allow counter to be incremented) - - tmpFac = -86400. _d 0/deltaTMom - CALL DIAGNOSTICS_SCALE_FILL( uVel, tmpFac, 1, - & 'TOTUTEND',0,Nr,0,1,1,myThid ) - CALL DIAGNOSTICS_SCALE_FILL( vVel, tmpFac, 1, - & 'TOTVTEND',0,Nr,0,1,1,myThid ) - - IF ( DIAGNOSTICS_IS_ON('TOTTTEND',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - tmpFac = -86400. _d 0/dTtracerLev(k) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = tmpFac*theta(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'TOTTTEND',0,Nr,0,1,1,myThid) -#ifdef ALLOW_LAYERS - IF ( useLayers ) THEN - CALL LAYERS_FILL(tmpMk,1,'TOT',0,Nr,0,1,1,myThid) - ENDIF -#endif /* ALLOW_LAYERS */ - ENDIF - - IF ( DIAGNOSTICS_IS_ON('TOTSTEND',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - tmpFac = -86400. _d 0/dTtracerLev(k) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = tmpFac*salt(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'TOTSTEND',0,Nr,0,1,1,myThid) -#ifdef ALLOW_LAYERS - IF ( useLayers ) THEN - CALL LAYERS_FILL(tmpMk,2,'TOT',0,Nr,0,1,1,myThid) - ENDIF -#endif /* ALLOW_LAYERS */ - ENDIF - -C-- fill momentum state-var diagnostics: end - ENDIF - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - - IF ( selectVars.EQ.1 .OR. selectVars.EQ.3 ) THEN -C-- fill tracer state-var diagnostics: - - CALL DIAGNOSTICS_FILL(theta,'THETA ',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_FILL(salt, 'SALT ',0,Nr,0,1,1,myThid) - -#ifdef ALLOW_FIZHI - IF ( useFIZHI .AND. DIAGNOSTICS_IS_ON('RELHUM ',myThid) ) THEN - kappa = getcon('KAPPA') - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO j = 1,sNy - DO i = 1,sNx - DO k = 1,Nr - dummy1 = theta(i,j,k,bi,bj) * ((rC(k)/100.)/1000.)**kappa - dummy2 = rC(k) / 100. - CALL QSAT(dummy1,dummy2,dummy3,dummy4,.false.) - tmpMk(i,j,k,bi,bj) = hFacC(i,j,k,bi,bj) - & *salt(i,j,k,bi,bj)*100. / dummy3 - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk, 'RELHUM ',0,Nr,0,1,1,myThid) - ENDIF -#endif /* ALLOW_FIZHI */ - - CALL DIAGNOSTICS_SCALE_FILL( theta, oneRL, 2, - & 'THETASQ ',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_SCALE_FILL( salt, oneRL, 2, - & 'SALTSQ ',0,Nr,0,1,1,myThid) - -#ifdef ALLOW_GENERIC_ADVDIFF -# ifdef ALLOW_ADAMSBASHFORTH_3 - IF ( selectVars.EQ.1 ) THEN -C- stagger time-step: fill diags after updating myIter - m1 = 1 + MOD(myIter,2) - ELSE -C- synchronous time-step: fill diags before updating myIter - m1 = 1 + MOD(myIter+1,2) - ENDIF - IF ( AdamsBashforthGt ) - & CALL DIAGNOSTICS_FILL( gtNm(1-OLx,1-OLy,1,1,1,m1), - & 'gTinAB ',0,Nr,0,1,1,myThid ) - IF ( AdamsBashforthGs ) - & CALL DIAGNOSTICS_FILL( gsNm(1-OLx,1-OLy,1,1,1,m1), - & 'gSinAB ',0,Nr,0,1,1,myThid ) -# else /* ALLOW_ADAMSBASHFORTH_3 */ - IF ( AdamsBashforthGt ) - & CALL DIAGNOSTICS_FILL( gtNm1,'gTinAB ',0,Nr,0,1,1,myThid ) - IF ( AdamsBashforthGs ) - & CALL DIAGNOSTICS_FILL( gsNm1,'gSinAB ',0,Nr,0,1,1,myThid ) -# endif /* ALLOW_ADAMSBASHFORTH_3 */ -#endif /* ALLOW_GENERIC_ADVDIFF */ - -c IF ( DIAGNOSTICS_IS_ON('SST ',myThid) ) THEN -c DO bj = myByLo(myThid), myByHi(myThid) -c DO bi = myBxLo(myThid), myBxHi(myThid) -c DO j = 1,sNy -c DO i = 1,sNx -c tmp1k(i,j,bi,bj) = THETA(i,j,1,bi,bj) -c ENDDO -c ENDDO -c ENDDO -c ENDDO -c CALL DIAGNOSTICS_FILL(tmp1k,'SST ',0,1,0,1,1,myThid) -c ENDIF - -c IF ( DIAGNOSTICS_IS_ON('SSS ',myThid) ) THEN -c DO bj = myByLo(myThid), myByHi(myThid) -c DO bi = myBxLo(myThid), myBxHi(myThid) -c DO j = 1,sNy -c DO i = 1,sNx -c tmp1k(i,j,bi,bj) = SALT(i,j,1,bi,bj) -c ENDDO -c ENDDO -c ENDDO -c ENDDO -c CALL DIAGNOSTICS_FILL(tmp1k,'SSS ',0,1,0,1,1,myThid) -c ENDIF - - IF ( fluidIsWater .AND. - & ( DIAGNOSTICS_IS_ON('SALTanom',myThid) - & .OR.DIAGNOSTICS_IS_ON('SALTSQan',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 - tmpMk(i,j,k,bi,bj) = salt(i,j,k,bi,bj)-35. _d 0 - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL( tmpMk,'SALTanom',0,Nr,0,1,1,myThid) - CALL DIAGNOSTICS_SCALE_FILL( tmpMk, oneRL, 2, - & 'SALTSQan',0,Nr,0,1,1,myThid) - ENDIF - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - - IF ( DIAGNOSTICS_IS_ON('UVELMASS',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+1 - tmpMk(i,j,k,bi,bj) - & = uVel(i,j,k,bi,bj)*hFacW(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UVELMASS',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('VVELMASS',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) - & = vVel(i,j,k,bi,bj)*hFacS(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'VVELMASS',0,Nr,0,1,1,myThid) - ENDIF - - CALL DIAGNOSTICS_FILL(wVel, 'WVELMASS',0,Nr,0,1,1,myThid) - - IF ( DIAGNOSTICS_IS_ON('UTHMASS ',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+1 - tmpMk(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)*0.5 _d 0 - & *(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj)) - & * hFacW(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'UTHMASS ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('VTHMASS ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)*0.5 _d 0 - & *(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj)) - & * hFacS(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'VTHMASS ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('WTHMASS ',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - km1 = MAX(k-1,1) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*0.5 _d 0 - & *(theta(i,j,k,bi,bj)+theta(i,j,km1,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'WTHMASS ',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('USLTMASS',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+1 - tmpMk(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)*0.5 _d 0 - & *(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj)) - & * hFacW(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'USLTMASS',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('VSLTMASS',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - DO j = 1,sNy+1 - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)*0.5 _d 0 - & *(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj)) - & * hFacS(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'VSLTMASS',0,Nr,0,1,1,myThid) - ENDIF - - IF ( DIAGNOSTICS_IS_ON('WSLTMASS',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - km1 = MAX(k-1,1) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*0.5 _d 0 - & *(salt(i,j,k,bi,bj)+salt(i,j,km1,bi,bj)) - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'WSLTMASS',0,Nr,0,1,1,myThid) - ENDIF - -C-- fill tracer state-var diagnostics: end - ENDIF - - IF ( selectVars.EQ.4 ) THEN -C Second fill sequence for state variable tendency diagnostics: add state variable -C NOTE: send a '-1' for the bibjflag and do not increment counter -C (previous fill for these diagnostics DID allow counter to be incremented) - - tmpFac = 86400. _d 0/deltaTMom - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - CALL DIAGNOSTICS_SCALE_FILL( uVel, tmpFac, 1, - & 'TOTUTEND',0,Nr,-1,bi,bj,myThid ) - CALL DIAGNOSTICS_SCALE_FILL( vVel, tmpFac, 1, - & 'TOTVTEND',0,Nr,-1,bi,bj,myThid ) - ENDDO - ENDDO - - IF ( DIAGNOSTICS_IS_ON('TOTTTEND',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - tmpFac = 86400. _d 0/dTtracerLev(k) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = tmpFac*theta(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'TOTTTEND',0,Nr,-1,bi,bj,myThid) -#ifdef ALLOW_LAYERS - IF ( useLayers ) THEN - CALL LAYERS_FILL(tmpMk,1,'TOT',0,Nr,-1,bi,bj,myThid) - ENDIF -#endif /* ALLOW_LAYERS */ - ENDDO - ENDDO - ENDIF - - IF ( DIAGNOSTICS_IS_ON('TOTSTEND',myThid) ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO k=1,Nr - tmpFac = 86400. _d 0/dTtracerLev(k) - DO j = 1,sNy - DO i = 1,sNx - tmpMk(i,j,k,bi,bj) = tmpFac*salt(i,j,k,bi,bj) - ENDDO - ENDDO - ENDDO - CALL DIAGNOSTICS_FILL(tmpMk,'TOTSTEND',0,Nr,-1,bi,bj,myThid) -#ifdef ALLOW_LAYERS - IF ( useLayers ) THEN - CALL LAYERS_FILL(tmpMk,2,'TOT',0,Nr,-1,bi,bj,myThid) - ENDIF -#endif /* ALLOW_LAYERS */ - ENDDO - ENDDO - ENDIF - -C-- fill state tendency diagnostics the second time: end - ENDIF - -#endif /* ALLOW_DIAGNOSTICS */ - - RETURN - END diff --git a/ECCOv4 Release 4/flux-forced/code/dynamics.F b/ECCOv4 Release 4/flux-forced/code/dynamics.F deleted file mode 100644 index 5c695de..0000000 --- a/ECCOv4 Release 4/flux-forced/code/dynamics.F +++ /dev/null @@ -1,745 +0,0 @@ -C $Header: /u/gcmpack/MITgcm/model/src/dynamics.F,v 1.178 2016/11/28 23:05:05 jmc Exp $ -C $Name: $ - -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" -#ifdef ALLOW_AUTODIFF -# include "AUTODIFF_OPTIONS.h" -#endif -#ifdef ALLOW_MOM_COMMON -# include "MOM_COMMON_OPTIONS.h" -#endif -#ifdef ALLOW_OBCS -# include "OBCS_OPTIONS.h" -#endif -#ifdef ALLOW_ECCO -#include "ECCO_OPTIONS.h" -#endif - -#undef DYNAMICS_GUGV_EXCH_CHECK - -CBOP -C !ROUTINE: DYNAMICS -C !INTERFACE: - SUBROUTINE DYNAMICS(myTime, myIter, myThid) -C !DESCRIPTION: \bv -C *==========================================================* -C | SUBROUTINE DYNAMICS -C | o Controlling routine for the explicit part of the model -C | dynamics. -C *==========================================================* -C \ev -C !USES: - IMPLICIT NONE -C == Global variables === -#include "SIZE.h" -#include "EEPARAMS.h" -#include "PARAMS.h" -#include "GRID.h" -#include "DYNVARS.h" -#ifdef ALLOW_MOM_COMMON -# include "MOM_VISC.h" -#endif -#ifdef ALLOW_CD_CODE -# include "CD_CODE_VARS.h" -#endif -#ifdef ALLOW_AUTODIFF -# include "tamc.h" -# include "tamc_keys.h" -# include "FFIELDS.h" -# include "EOS.h" -# ifdef ALLOW_KPP -# include "KPP.h" -# endif -# ifdef ALLOW_PTRACERS -# include "PTRACERS_SIZE.h" -# include "PTRACERS_FIELDS.h" -# endif -# ifdef ALLOW_OBCS -# include "OBCS_PARAMS.h" -# include "OBCS_FIELDS.h" -# ifdef ALLOW_PTRACERS -# include "OBCS_PTRACERS.h" -# endif -# endif -# ifdef ALLOW_MOM_FLUXFORM -# include "MOM_FLUXFORM.h" -# endif -#endif /* ALLOW_AUTODIFF */ -#ifdef ALLOW_ECCO -# include "ecco.h" -#endif - -C !CALLING SEQUENCE: -C DYNAMICS() -C | -C |-- CALC_EP_FORCING -C | -C |-- CALC_GRAD_PHI_SURF -C | -C |-- CALC_VISCOSITY -C | -C |-- MOM_CALC_3D_STRAIN -C | -C |-- CALC_EDDY_STRESS -C | -C |-- CALC_PHI_HYD -C | -C |-- MOM_FLUXFORM -C | -C |-- MOM_VECINV -C | -C |-- MOM_CALC_SMAG_3D -C |-- MOM_UV_SMAG_3D -C | -C |-- TIMESTEP -C | -C |-- MOM_U_IMPLICIT_R -C |-- MOM_V_IMPLICIT_R -C | -C |-- IMPLDIFF -C | -C |-- OBCS_APPLY_UV -C | -C |-- CALC_GW -C | -C |-- DIAGNOSTICS_FILL -C |-- DEBUG_STATS_RL - -C !INPUT/OUTPUT PARAMETERS: -C == Routine arguments == -C myTime :: Current time in simulation -C myIter :: Current iteration number in simulation -C myThid :: Thread number for this instance of the routine. - _RL myTime - INTEGER myIter - INTEGER myThid - -C !FUNCTIONS: -#ifdef ALLOW_DIAGNOSTICS - LOGICAL DIAGNOSTICS_IS_ON - EXTERNAL DIAGNOSTICS_IS_ON -#endif - -C !LOCAL VARIABLES: -C == Local variables -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 kappaRU :: vertical viscosity for velocity U-component -C kappaRV :: vertical viscosity for velocity V-component -C iMin, iMax :: Ranges and sub-block indices on which calculations -C jMin, jMax are applied. -C bi, bj :: tile indices -C k :: current level index -C km1, kp1 :: index of level above (k-1) and below (k+1) -C kUp, kDown :: Index for interface above and below. kUp and kDown are -C are switched with k to be the appropriate index into fVerU,V - _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) - _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) - _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+1) - _RL kappaRV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) -#ifdef ALLOW_SMAG_3D -C str11 :: strain component Vxx @ grid-cell center -C str22 :: strain component Vyy @ grid-cell center -C str33 :: strain component Vzz @ grid-cell center -C str12 :: strain component Vxy @ grid-cell corner -C str13 :: strain component Vxz @ above uVel -C str23 :: strain component Vyz @ above vVel -C viscAh3d_00 :: Smagorinsky viscosity @ grid-cell center -C viscAh3d_12 :: Smagorinsky viscosity @ grid-cell corner -C viscAh3d_13 :: Smagorinsky viscosity @ above uVel -C viscAh3d_23 :: Smagorinsky viscosity @ above vVel -C addDissU :: zonal momentum tendency from 3-D Smag. viscosity -C addDissV :: merid momentum tendency from 3-D Smag. viscosity - _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr ) - _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr ) - _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr ) - _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr ) - _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) - _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) - _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr ) - _RL viscAh3d_12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr ) - _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) - _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) - _RL addDissU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL addDissV(1-OLx:sNx+OLx,1-OLy:sNy+OLy) -#elif ( defined ALLOW_NONHYDROSTATIC ) - _RL str13(1), str23(1), str33(1) - _RL viscAh3d_00(1), viscAh3d_13(1), viscAh3d_23(1) -#endif - - INTEGER bi, bj - INTEGER i, j - INTEGER k, km1, kp1, kUp, kDown - INTEGER iMin, iMax - INTEGER jMin, jMax - PARAMETER( iMin = 0 , iMax = sNx+1 ) - PARAMETER( jMin = 0 , jMax = sNy+1 ) - -#ifdef ALLOW_DIAGNOSTICS - LOGICAL dPhiHydDiagIsOn - _RL tmpFac -#endif /* ALLOW_DIAGNOSTICS */ - -C--- The algorithm... -C -C "Correction Step" -C ================= -C Here we update the horizontal velocities with the surface -C pressure such that the resulting flow is either consistent -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 -C "Calculation of Gs" -C =================== -C This is where all the accelerations and tendencies (ie. -C physics, parameterizations etc...) are calculated -C rho = rho ( theta[n], salt[n] ) -C b = b(rho, theta) -C K31 = K31 ( rho ) -C Gu[n] = Gu( u[n], v[n], wVel, b, ... ) -C Gv[n] = Gv( u[n], v[n], wVel, b, ... ) -C Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... ) -C Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... ) -C -C "Time-stepping" or "Prediction" -C ================================ -C The models variables are stepped forward with the appropriate -C time-stepping scheme (currently we use Adams-Bashforth II) -C - For momentum, the result is always *only* a "prediction" -C in that the flow may be divergent and will be "corrected" -C later with a surface pressure gradient. -C - Normally for tracers the result is the new field at time -C level [n+1} *BUT* in the case of implicit diffusion the result -C is also *only* a prediction. -C - We denote "predictors" with an asterisk (*). -C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] ) -C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] ) -C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C With implicit diffusion: -C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C (1 + dt * K * d_zz) theta[n] = theta* -C (1 + dt * K * d_zz) salt[n] = salt* -C--- -CEOP - -#ifdef ALLOW_DEBUG - IF (debugMode) CALL DEBUG_ENTER( 'DYNAMICS', myThid ) -#endif - -#ifdef ALLOW_DIAGNOSTICS - dPhiHydDiagIsOn = .FALSE. - IF ( useDiagnostics ) - & dPhiHydDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPHdx', myThid ) - & .OR. DIAGNOSTICS_IS_ON( 'Vm_dPHdy', myThid ) -#endif - -C-- Call to routine for calculation of Eliassen-Palm-flux-forced -C U-tendency, if desired: -#ifdef INCLUDE_EP_FORCING_CODE - CALL CALC_EP_FORCING(myThid) -#endif - -#ifdef ALLOW_AUTODIFF_MONITOR_DIAG - CALL DUMMY_IN_DYNAMICS( myTime, myIter, myThid ) -#endif - -#ifdef ALLOW_AUTODIFF_TAMC -C-- HPF directive to help TAMC -CHPF$ INDEPENDENT -#endif /* ALLOW_AUTODIFF_TAMC */ - - DO bj=myByLo(myThid),myByHi(myThid) - -#ifdef ALLOW_AUTODIFF_TAMC -C-- HPF directive to help TAMC -CHPF$ INDEPENDENT, NEW (fVerU,fVerV -CHPF$& ,phiHydF -CHPF$& ,kappaRU,kappaRV -CHPF$& ) -#endif /* ALLOW_AUTODIFF_TAMC */ - - DO bi=myBxLo(myThid),myBxHi(myThid) - -#ifdef ALLOW_AUTODIFF_TAMC - act1 = bi - myBxLo(myThid) - max1 = myBxHi(myThid) - myBxLo(myThid) + 1 - act2 = bj - myByLo(myThid) - max2 = myByHi(myThid) - myByLo(myThid) + 1 - act3 = myThid - 1 - max3 = nTx*nTy - act4 = ikey_dynamics - 1 - idynkey = (act1 + 1) + act2*max1 - & + act3*max1*max2 - & + act4*max1*max2*max3 -#endif /* ALLOW_AUTODIFF_TAMC */ - -C-- Set up work arrays with valid (i.e. not NaN) values -C These initial 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. - -#ifdef ALLOW_AUTODIFF - DO k=1,Nr - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx -c-- need some re-initialisation here to break dependencies - gU(i,j,k,bi,bj) = 0. _d 0 - gV(i,j,k,bi,bj) = 0. _d 0 - ENDDO - ENDDO - ENDDO -#endif /* ALLOW_AUTODIFF */ - 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 -#ifndef INCLUDE_PHIHYD_CALCULATION_CODE - dPhiHydX(i,j) = 0. _d 0 - dPhiHydY(i,j) = 0. _d 0 -#endif - phiSurfX(i,j) = 0. _d 0 - phiSurfY(i,j) = 0. _d 0 - guDissip(i,j) = 0. _d 0 - gvDissip(i,j) = 0. _d 0 -#ifdef ALLOW_AUTODIFF - phiHydLow(i,j,bi,bj) = 0. _d 0 -# if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM) -# ifndef DISABLE_RSTAR_CODE - dWtransC(i,j,bi,bj) = 0. _d 0 - dWtransU(i,j,bi,bj) = 0. _d 0 - dWtransV(i,j,bi,bj) = 0. _d 0 -# endif -# endif -#endif /* ALLOW_AUTODIFF */ - ENDDO - ENDDO - -C-- Start computation of dynamics - -#ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE wVel (:,:,:,bi,bj) = -CADJ & comlev1_bibj, key=idynkey, byte=isbyte -#endif /* ALLOW_AUTODIFF_TAMC */ - -C-- Explicit part of the Surface Potential Gradient (add in TIMESTEP) -C (note: this loop will be replaced by CALL CALC_GRAD_ETA) - IF (implicSurfPress.NE.1.) THEN - CALL CALC_GRAD_PHI_SURF( - I bi,bj,iMin,iMax,jMin,jMax, - I etaN, - O phiSurfX,phiSurfY, - I myThid ) - ENDIF - -#ifdef ALLOW_AUTODIFF_TAMC -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=idynkey, byte=isbyte -#endif /* ALLOW_KPP */ -#endif /* ALLOW_AUTODIFF_TAMC */ - -#ifndef ALLOW_AUTODIFF - IF ( .NOT.momViscosity ) THEN -#endif - DO k=1,Nr+1 - 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 - ENDDO - ENDDO - ENDDO -#ifndef ALLOW_AUTODIFF - ENDIF -#endif -#ifdef INCLUDE_CALC_DIFFUSIVITY_CALL -C-- Calculate the total vertical viscosity - IF ( momViscosity ) THEN - CALL CALC_VISCOSITY( - I bi,bj, iMin,iMax,jMin,jMax, - O kappaRU, kappaRV, - I myThid ) - ENDIF -#endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */ - -#ifdef ALLOW_SMAG_3D - IF ( useSmag3D ) THEN - CALL MOM_CALC_3D_STRAIN( - O str11, str22, str33, str12, str13, str23, - I bi, bj, myThid ) - ENDIF -#endif /* ALLOW_SMAG_3D */ - -#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 */ - -#ifdef ALLOW_OBCS -C-- For Stevens boundary conditions velocities need to be extrapolated -C (copied) to a narrow strip outside the domain - IF ( useOBCS ) THEN - CALL OBCS_COPY_UV_N( - U uVel(1-OLx,1-OLy,1,bi,bj), - U vVel(1-OLx,1-OLy,1,bi,bj), - I Nr, bi, bj, myThid ) - ENDIF -#endif /* ALLOW_OBCS */ - -#ifdef ALLOW_EDDYPSI - CALL CALC_EDDY_STRESS(bi,bj,myThid) -#endif - -C-- Start of dynamics loop - DO k=1,Nr - -C-- km1 Points to level above k (=k-1) -C-- kup Cycles through 1,2 to point to layer above -C-- kDown Cycles through 2,1 to point to current layer - - km1 = MAX(1,k-1) - kp1 = MIN(k+1,Nr) - kup = 1+MOD(k+1,2) - kDown= 1+MOD(k,2) - -#ifdef ALLOW_AUTODIFF_TAMC - kkey = (idynkey-1)*Nr + k -CADJ STORE totPhiHyd (:,:,k,bi,bj) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE phiHydLow (:,:,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 -# 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 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 -# ifndef ALLOW_ADAMSBASHFORTH_3 -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 -# else -CADJ STORE guNm(:,:,k,bi,bj,1) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE guNm(:,:,k,bi,bj,2) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE gvNm(:,:,k,bi,bj,1) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE gvNm(:,:,k,bi,bj,2) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -# endif -# 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 /* NONLIN_FRSURF */ -#endif /* ALLOW_AUTODIFF_TAMC */ - -C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0 - CALL CALC_PHI_HYD( - I bi,bj,iMin,iMax,jMin,jMax,k, - I theta, salt, - U phiHydF, - O phiHydC, dPhiHydX, dPhiHydY, - I myTime, myIter, myThid ) -#ifdef ALLOW_DIAGNOSTICS - IF ( dPhiHydDiagIsOn ) THEN - tmpFac = -1. _d 0 - CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1, - & 'Um_dPHdx', k, 1, 2, bi, bj, myThid ) - CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1, - & 'Vm_dPHdy', k, 1, 2, bi, bj, myThid ) - ENDIF -#endif /* ALLOW_DIAGNOSTICS */ - -C-- Calculate accelerations in the momentum equations (gU, gV, ...) -C and step forward storing the result in gU, gV, etc... - IF ( momStepping ) THEN -#ifdef ALLOW_AUTODIFF - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - guDissip(i,j) = 0. _d 0 - gvDissip(i,j) = 0. _d 0 - ENDDO - ENDDO -#endif /* ALLOW_AUTODIFF */ -#ifdef ALLOW_AUTODIFF_TAMC -# if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM) -# ifndef DISABLE_RSTAR_CODE -CADJ STORE dWtransC(:,:,bi,bj) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE dWtransU(:,:,bi,bj) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE dWtransV(:,:,bi,bj) -CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte -# endif -# endif /* NONLIN_FRSURF and ALLOW_MOM_FLUXFORM */ -# if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL) -CADJ STORE fVerU(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE fVerV(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte -# endif -#endif /* ALLOW_AUTODIFF_TAMC */ - IF (.NOT. vectorInvariantMomentum) THEN -#ifdef ALLOW_MOM_FLUXFORM - CALL MOM_FLUXFORM( - I bi,bj,k,iMin,iMax,jMin,jMax, - I kappaRU, kappaRV, - U fVerU(1-OLx,1-OLy,kUp), fVerV(1-OLx,1-OLy,kUp), - O fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown), - O guDissip, gvDissip, - I myTime, myIter, myThid) -#endif - ELSE -#ifdef ALLOW_MOM_VECINV - CALL MOM_VECINV( - I bi,bj,k,iMin,iMax,jMin,jMax, - I kappaRU, kappaRV, - I fVerU(1-OLx,1-OLy,kUp), fVerV(1-OLx,1-OLy,kUp), - O fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown), - O guDissip, gvDissip, - I myTime, myIter, myThid) -#endif - ENDIF - -#ifdef ALLOW_SMAG_3D - IF ( useSmag3D ) THEN - CALL MOM_CALC_SMAG_3D( - I str11, str22, str33, str12, str13, str23, - O viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23, - I smag3D_hLsC, smag3D_hLsW, smag3D_hLsS, smag3D_hLsZ, - I k, bi, bj, myThid ) - CALL MOM_UV_SMAG_3D( - I str11, str22, str12, str13, str23, - I viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23, - O addDissU, addDissV, - I iMin,iMax,jMin,jMax, k, bi, bj, myThid ) - DO j= jMin,jMax - DO i= iMin,iMax - guDissip(i,j) = guDissip(i,j) + addDissU(i,j) - gvDissip(i,j) = gvDissip(i,j) + addDissV(i,j) - ENDDO - ENDDO - ENDIF -#endif /* ALLOW_SMAG_3D */ - - CALL TIMESTEP( - I bi,bj,iMin,iMax,jMin,jMax,k, - I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY, - I guDissip, gvDissip, - I myTime, myIter, myThid) - - ENDIF - -C-- end of dynamics k loop (1:Nr) - ENDDO - -C-- Implicit Vertical advection & viscosity -#if (defined (INCLUDE_IMPLVERTADV_CODE) && \ - defined (ALLOW_MOM_COMMON) && !(defined ALLOW_AUTODIFF)) - IF ( momImplVertAdv .OR. implicitViscosity - & .OR. selectImplicitDrag.GE.1 ) THEN -C to recover older (prior to 2016-10-05) results: -c 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 -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 -1, kappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj), - U gU(1-OLx,1-OLy,1,bi,bj), - I myThid ) -#ifdef ALLOW_AUTODIFF_TAMC -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 -2, kappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj), - U gV(1-OLx,1-OLy,1,bi,bj), - I myThid ) - ENDIF - -#ifdef ALLOW_OBCS -C-- Apply open boundary conditions - IF ( useOBCS ) THEN -C-- but first save intermediate velocities to be used in the -C next time step for the Stevens boundary conditions - CALL OBCS_SAVE_UV_N( - I bi, bj, iMin, iMax, jMin, jMax, 0, - I gU, gV, myThid ) - CALL OBCS_APPLY_UV( bi, bj, 0, gU, gV, myThid ) - ENDIF -#endif /* ALLOW_OBCS */ - -#ifdef ALLOW_CD_CODE - IF (implicitViscosity.AND.useCDscheme) THEN -#ifdef ALLOW_AUTODIFF_TAMC -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 0, kappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj), - U vVelD(1-OLx,1-OLy,1,bi,bj), - I myThid ) -#ifdef ALLOW_AUTODIFF_TAMC -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 0, kappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj), - U uVelD(1-OLx,1-OLy,1,bi,bj), - I myThid ) - ENDIF -#endif /* ALLOW_CD_CODE */ -C-- End implicit Vertical advection & viscosity - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - -#ifdef ALLOW_NONHYDROSTATIC -C-- Step forward W field in N-H algorithm - IF ( nonHydrostatic ) THEN -#ifdef ALLOW_DEBUG - IF (debugMode) CALL DEBUG_CALL('CALC_GW', myThid ) -#endif - CALL TIMER_START('CALC_GW [DYNAMICS]',myThid) - CALL CALC_GW( - I bi,bj, kappaRU, kappaRV, - I str13, str23, str33, - I viscAh3d_00, viscAh3d_13, viscAh3d_23, - I myTime, myIter, myThid ) - ENDIF - IF ( nonHydrostatic.OR.implicitIntGravWave ) - & CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid ) - IF ( nonHydrostatic ) - & CALL TIMER_STOP ('CALC_GW [DYNAMICS]',myThid) -#endif - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - -C- end of bi,bj loops - ENDDO - ENDDO - -#ifdef ALLOW_OBCS - IF (useOBCS) THEN - CALL OBCS_EXCHANGES( myThid ) - ENDIF -#endif - -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, maskInC, maskInC, rA, drF, -C & 'phiHydLow', myTime, 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) - CALL DIAGNOSTICS_SCALE_FILL(m_bp,recip_gravity,1, - & 'OBPGMAP ',0, 1,0,1,1,myThid) - CALL DIAGNOSTICS_SCALE_FILL(m_bp_nopabar,recip_gravity,1, - & 'OBP ',0, 1,0,1,1,myThid) - - tmpFac = 1. _d 0 - CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2, - & 'PHIHYDSQ',0,Nr,0,1,1,myThid) - - CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2, - & 'PHIBOTSQ',0, 1,0,1,1,myThid) - - ENDIF -#endif /* ALLOW_DIAGNOSTICS */ - -#ifdef ALLOW_DEBUG - IF ( debugLevel .GE. debLevD ) 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) - CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid) -#ifndef ALLOW_ADAMSBASHFORTH_3 - CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid) - CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid) -#endif - 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 .GE. debLevE ) THEN - CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid) - ENDIF -#endif - -#ifdef ALLOW_DEBUG - IF (debugMode) CALL DEBUG_LEAVE( 'DYNAMICS', myThid ) -#endif - - RETURN - END diff --git a/ECCOv4 Release 4/flux-forced/code/ecco_init_varia.F b/ECCOv4 Release 4/flux-forced/code/ecco_init_varia.F index 977fe90..d8c514d 100644 --- a/ECCOv4 Release 4/flux-forced/code/ecco_init_varia.F +++ b/ECCOv4 Release 4/flux-forced/code/ecco_init_varia.F @@ -47,7 +47,7 @@ SUBROUTINE ECCO_INIT_VARIA( myThid ) ENDIF #endif /* ALLOW_PSBAR_STERIC */ - CALL ECCO_PHYS( myThid ) + CALL ECCO_PHYS( startTime, -1, myThid ) #ifdef ALLOW_PSBAR_STERIC C RHO/VOLsumGlob_0 are zeros if S/R ECCO_READ_PICKUP is not called diff --git a/ECCOv4 Release 4/flux-forced/code/ecco_phys.F b/ECCOv4 Release 4/flux-forced/code/ecco_phys.F index a8c66ab..062b742 100644 --- a/ECCOv4 Release 4/flux-forced/code/ecco_phys.F +++ b/ECCOv4 Release 4/flux-forced/code/ecco_phys.F @@ -3,7 +3,7 @@ #include "ECCO_OPTIONS.h" - subroutine ecco_phys( mythid ) + SUBROUTINE ECCO_PHYS( myTime, myIter, myThid ) c ================================================================== c SUBROUTINE ecco_phys @@ -33,16 +33,17 @@ subroutine ecco_phys( mythid ) c == routine arguments == - integer mythid + _RL myTime + INTEGER myIter, myThid c == local variables == - integer bi,bj - integer i,j,k - integer jmin,jmax - integer imin,imax + INTEGER bi,bj + INTEGER i,j,k + INTEGER jmin,jmax + INTEGER imin,imax #ifdef ALLOW_GENCOST_CONTRIBUTION - integer kgen, kgen3d, itr + INTEGER kgen, kgen3d, itr _RL areavolTile(nSx,nSy), areavolGlob _RL tmpfld, tmpvol, tmpmsk, tmpmsk2, tmpmskW, tmpmskS #endif @@ -215,6 +216,24 @@ subroutine ecco_phys( mythid ) enddo enddo +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics .AND. myIter.GE.0 ) THEN + CALL DIAGNOSTICS_FILL( m_eta, 'SSHNOIBC', 0,1, 0,1,1, myThid ) + CALL DIAGNOSTICS_SCALE_FILL( m_bp, recip_gravity, 1, + & 'OBPGMAP ', 0,1, 0,1,1, myThid ) +#ifdef ATMOSPHERIC_LOADING +#ifdef ALLOW_IB_CORR + CALL DIAGNOSTICS_FILL( m_eta_ib, + & 'SSHIBC ', 0,1, 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( m_eta_dyn, + & 'SSH ', 0,1, 0,1,1, myThid ) + CALL DIAGNOSTICS_SCALE_FILL( m_bp_nopabar, recip_gravity, 1, + & 'OBP ', 0,1, 0,1,1, myThid ) +#endif /* ALLOW_IB_CORR */ +#endif /* ATMOSPHERIC_LOADING */ + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) do k = 1,nr diff --git a/ECCOv4 Release 4/flux-forced/code/forward_step.F b/ECCOv4 Release 4/flux-forced/code/forward_step.F index e30d90a..118d955 100644 --- a/ECCOv4 Release 4/flux-forced/code/forward_step.F +++ b/ECCOv4 Release 4/flux-forced/code/forward_step.F @@ -1186,7 +1186,7 @@ SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid ) #ifdef ALLOW_ECCO C-- Diagnoze variables for pkg/ecco averaging and cost function purposes - IF ( useECCO ) CALL ECCO_PHYS( myThid ) + IF ( useECCO ) CALL ECCO_PHYS( myTime, myIter, myThid ) #endif C-- Check if it has reached the end of simulation