C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/PRM/multi_comp_setup/cg/code/set_ddtvars.F,v 1.7 2009/01/04 17:41:54 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #ifdef ALLOW_MYPACKAGE # include "MYPACKAGE_OPTIONS.h" #else # include "CPP_OPTIONS.h" #endif CBOP C !ROUTINE: SET_DDTVARS C !INTERFACE: SUBROUTINE SET_DDTVARS( gUVelC, gVVelC, gThetaC, gSaltC, I cnx, cny, cnr, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SET_DDTVARS C | o Set tendencies to be used as additional forcing C *==========================================================* C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global data == #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "ORIENTATION.h" #include "EXT_MOM_TEND.h" #ifdef ALLOW_MYPACKAGE # include "MYPACKAGE.h" #endif C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C myTime :: Current time in simulation C myIter :: Current iteration number C myThid :: Thread Id number INTEGER cnx, cny, cnr REAL*8 gUVelC( cnx, cny, cnr) REAL*8 gVVelC( cnx, cny, cnr) REAL*8 gThetaC(cnx, cny, cnr) REAL*8 gSaltC( cnx, cny, cnr) _RL myTime INTEGER myIter INTEGER myThid C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE C !LOCAL VARIABLES: C == Local variables == C i,j,k :: Loop counters INTEGER i,j,k INTEGER bi,bj C Variables used for I/O CHARACTER*(10) suff C Variables used for smoothing momentum tendencies: _RL viscAhDivDt, viscAhVortDt _RL locDiv (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL locVort(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL locFx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL locFy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP viscAhDivDt = myPa_param1 * deltaTmom viscAhVortDt = myPa_param2 * deltaTmom IF ( cnx .NE. sNx ) THEN STOP 'SET_DTTVARS cnx NE sNx' ENDIF IF ( cny .NE. sNy ) THEN STOP 'SET_DTTVARS cny NE sNy' ENDIF IF ( cnr .NE. Nr ) THEN STOP 'SET_DTTVARS cnr NE Nr' ENDIF DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) C--- Temp. & Salt Tendency from CG : just do a copy DO k=1,Nr DO j=1,sNy DO i=1,sNx myPa_TendScal1(i,j,k,bi,bj) = gThetaC(i,j,k) myPa_TendScal2(i,j,k,bi,bj) = gsaltC (i,j,k) ENDDO ENDDO ENDDO C--- Momentum Tendency from CG C-- 1) rotate to get u,v aline with CG axes DO k=1,Nr DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx ext_gu(i,j,k,bi,bj) = 0. _d 0 ext_gv(i,j,k,bi,bj) = 0. _d 0 ENDDO ENDDO IF ( myPa_doSwitch1 ) THEN DO j=1,sNy DO i=1,sNx IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj) ext_gv(i,j,k,bi,bj) = gUVelC(i,j,k)*sAngleFG(i,j,bi,bj) ENDIF ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj) & - gVVelC(i,j,k)*sAngleFG(i,j,bi,bj) ext_gv(i,j,k,bi,bj) = gUVelC(i,j,k)*sAngleFG(i,j,bi,bj) & + gVVelC(i,j,k)*cAngleFG(i,j,bi,bj) ENDIF ENDDO ENDDO ENDIF ENDDO C- end bi,bj loops ENDDO ENDDO C-- 2) Fill the overlap : ext_gu & ext_gv are both on A-grid CALL EXCH_UV_AGRID_3D_RL( ext_gu, ext_gv, & .TRUE., Nr, myThid ) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO k=1,Nr C-- 3) average to C-grid DO j=1-OLy,sNy+OLy DO i=2-OLx,sNx+OLx myPa_TendVelU(i,j,k,bi,bj) = maskW(i,j,k,bi,bj) & *( ext_gu(i-1,j,k,bi,bj) & +ext_gu( i ,j,k,bi,bj) )*0.5 _d 0 ENDDO ENDDO DO j=2-OLy,sNy+OLy DO i=1-OLx,sNx+OLx myPa_TendVelV(i,j,k,bi,bj) = maskS(i,j,k,bi,bj) & *( ext_gv(i,j-1,k,bi,bj) & +ext_gv(i, j ,k,bi,bj) )*0.5 _d 0 ENDDO ENDDO C-- 4) Smooth Momentum tendency : C apply horizontal viscosity (*time-step = viscAhDivDt) on Divergence C tendency and on vorticity tend. (viscosity*time-step = viscAhVortDt) IF ( viscAhDivDt.GT.0. .OR. viscAhVortDt.GT.0. ) THEN C Fluxes DO j=2-OLy,sNy+OLy DO i=2-OLx,sNx+OLx locFx(i,j) = myPa_TendVelU(i,j,k,bi,bj) & *dyG(i,j,bi,bj)*hFacW(i,j,k,bi,bj) locFy(i,j) = myPa_TendVelV(i,j,k,bi,bj) & *dxG(i,j,bi,bj)*hFacS(i,j,k,bi,bj) ENDDO ENDDO C Divergence DO j=2-OLy,sNy+OLy-1 DO i=2-OLx,sNx+OLx-1 locDiv(i,j) = & ( ( locFx(i+1,j)-locFx(i,j) ) & +( locFy(i,j+1)-locFy(i,j) ) & )*recip_rA(i,j,bi,bj)*recip_hFacC(i,j,k,bi,bj) ENDDO ENDDO C Vorticity DO j=2-OLy,sNy+OLy DO i=2-OLx,sNx+OLx locVort(i,j) = & ( ( myPa_TendVelV(i,j,k,bi,bj)*dyC(i,j,bi,bj) & -myPa_TendVelV(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj) ) & -( myPa_TendVelU(i,j,k,bi,bj)*dxC(i,j,bi,bj) & -myPa_TendVelU(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj) ) & )*recip_rAz(I,J,bi,bj) ENDDO ENDDO C Apply to C-grid tendencies: DO j=3-OLy,sNy+OLy-1 DO i=3-OLx,sNx+OLx-1 myPa_TendVelU(i,j,k,bi,bj) = myPa_TendVelU(i,j,k,bi,bj) & + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i-1,j) ) & * recip_dxC(i,j,bi,bj) & - viscAhVortDt* ( locVort(i,j+1)-locVort(i,j) ) & * recip_dyG(i,j,bi,bj) & )*maskW(i,j,k,bi,bj) myPa_TendVelV(i,j,k,bi,bj) = myPa_TendVelV(i,j,k,bi,bj) & + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i,j-1) ) & * recip_dyC(i,j,bi,bj) & + viscAhVortDt* ( locVort(i+1,j)-locVort(i,j) ) & * recip_dxG(i,j,bi,bj) & )*maskS(i,j,k,bi,bj) ENDDO ENDDO C- end smoothing & end of k loop ENDIF ENDDO C- end bi,bj loops ENDDO ENDDO C-- Diagnostics #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN CALL DIAGNOSTICS_FILL ( myPa_TendScal1,'MYPadTdt', & 0, Nr, 0, 1, 1, myThid ) CALL DIAGNOSTICS_FILL ( myPa_TendScal2,'MYPadSdt', & 0, Nr, 0, 1, 1, myThid ) CALL DIAGNOSTICS_FILL ( myPa_TendVelU, 'MYPadUdt', & 0, Nr, 0, 1, 1, myThid ) CALL DIAGNOSTICS_FILL ( myPa_TendVelV, 'MYPadVdt', & 0, Nr, 0, 1, 1, myThid ) C- This is a hack (to avoid changing mypackage_diagnostics_init.F) CALL DIAGNOSTICS_FILL ( ext_gu, 'MYPaStaU', & 0, Nr, 0, 1, 1, myThid ) CALL DIAGNOSTICS_FILL ( ext_gv, 'MYPaStaV', & 0, Nr, 0, 1, 1, myThid ) ENDIF #endif /* ALLOW_DIAGNOSTICS */ C-- Write snap-shot output: IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime,deltaTClock) & .OR. dumpInitAndLast.AND.( myTime.EQ.endTime .OR. & myTime.EQ.startTime ) & ) THEN WRITE(suff,'(I10.10)') myIter CALL WRITE_FLD_XYZ_RL( 'ext_gT.', suff, myPa_TendScal1, & myIter, myThid ) CALL WRITE_FLD_XYZ_RL( 'ext_gS.', suff, myPa_TendScal2, & myIter, myThid ) CALL WRITE_FLD_XYZ_RL( 'ext_gU.', suff, ext_gu, & myIter, myThid ) CALL WRITE_FLD_XYZ_RL( 'ext_gV.', suff, ext_gv, & myIter, myThid ) ENDIF RETURN END