C $Header: /u/gcmpack/MITgcm/pkg/smooth/smooth_rhs.F,v 1.4 2015/01/23 18:58:26 gforget Exp $ C $Name: $ #include "SMOOTH_OPTIONS.h" C !INTERFACE: ========================================================== SUBROUTINE SMOOTH_RHS(fld_in,gt_in,myThid) C *==========================================================* C | SUBROUTINE smooth_rhs C | o As part of smooth_diff3D, this routine computes the C | right hand side of the tendency equation (see below). C | It is made of bits from model/src and pkg/generic_advdiff C | pieced togheter. C *==========================================================* C !DESCRIPTION: C Calculates the tendency of a tracer due to advection and diffusion. C It calculates the fluxes in each direction indepentently and then C sets the tendency to the divergence of these fluxes. The advective C fluxes are only calculated here when using the linear advection schemes C otherwise only the diffusive and parameterized fluxes are calculated. C C Contributions to the flux are calculated and added: C \begin{equation*} C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP} C \end{equation*} C C The tendency is the divergence of the fluxes: C \begin{equation*} C G_\theta = G_\theta + \nabla \cdot {\bf F} C \end{equation*} C C The tendency is assumed to contain data on entry. C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" c#include "EESUPPORT.h" #include "PARAMS.h" #include "GRID.h" #include "SMOOTH.h" C !INPUT PARAMETERS: =================================================== INTEGER myThid _RL fld_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) _RL gt_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) C local variables: INTEGER bi,bj,iMin,iMax,jMin,jMax _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL dTdz (nSx,nSy) _RL dTdx (nSx,nSy) _RL dTdy (nSx,nSy) INTEGER i,j,k _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy) _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy) _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy) _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) c 1rst k loop: initialization DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx fZon(i,j,k,bi,bj) = 0. _d 0 fMer(i,j,k,bi,bj) = 0. _d 0 fVerT(i,j,k,bi,bj) = 0. _d 0 gt_in(i,j,k,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO iMin = 1-OLx+1 iMax = sNx+OLx-1 jMin = 1-OLy+1 jMax = sNy+OLy-1 c 2nd k loop: flux computation DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx df(i,j,bi,bj) = 0. _d 0 xA(i,j,bi,bj) = _dyG(i,j,bi,bj) & *drF(k)*smooth_hFacW(i,j,k,bi,bj) yA(i,j,bi,bj) = _dxG(i,j,bi,bj) & *drF(k)*smooth_hFacS(i,j,k,bi,bj) IF (K .EQ. 1) THEN maskUp(i,j,bi,bj) = 0. ELSE maskUp(i,j,bi,bj) = & maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj) ENDIF ENDDO ENDDO c ///gmredi_xtr/// DO j=jMin,jMax DO i=iMin,iMax df(i,j,bi,bj) = df(i,j,bi,bj) & -xA(i,j,bi,bj) & *smooth3D_Kux(i,j,k,bi,bj) & *recip_dxC(i,j,bi,bj) & *(fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj)) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax dTdz(bi,bj) = 0.5*( & +0.5*recip_drC(k)* & ( maskC(i-1,j,k,bi,bj)* & (fld_in(i-1,j, MAX(k-1,1) ,bi,bj)-fld_in(i-1,j,k,bi,bj)) & +maskC( i ,j,k,bi,bj)* & (fld_in( i ,j, MAX(k-1,1) ,bi,bj)-fld_in( i ,j,k,bi,bj)) & ) & +0.5*recip_drC(MIN(k+1,Nr))* & ( maskC(i-1,j,MIN(k+1,Nr),bi,bj)* & (fld_in(i-1,j,k,bi,bj)-fld_in(i-1,j,MIN(k+1,Nr),bi,bj)) & +maskC( i ,j,MIN(k+1,Nr),bi,bj)* & (fld_in( i ,j,k,bi,bj)-fld_in( i ,j,MIN(k+1,Nr),bi,bj)) & ) ) df(i,j,bi,bj) = df(i,j,bi,bj) & - xA(i,j,bi,bj)*smooth3D_Kuz(i,j,k,bi,bj)*dTdz(bi,bj) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax dTdy(bi,bj) = 0.5*( & +0.5*(maskS(i,j,k,bi,bj) & *recip_dyC(i,j,bi,bj)* & (fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj)) & +maskS(i,j+1,k,bi,bj) & *recip_dyC(i,j+1,bi,bj)* & (fld_in(i,j+1,k,bi,bj)-fld_in(i,j,k,bi,bj))) & +0.5*(maskS(i-1,j,k,bi,bj) & *recip_dyC(i,j,bi,bj)* & (fld_in(i-1,j,k,bi,bj)-fld_in(i-1,j-1,k,bi,bj)) & +maskS(i-1,j+1,k,bi,bj) & *recip_dyC(i,j+1,bi,bj)* & (fld_in(i-1,j+1,k,bi,bj)-fld_in(i-1,j,k,bi,bj))) & ) df(i,j,bi,bj) = df(i,j,bi,bj) & - xA(i,j,bi,bj)*smooth3D_Kuy(i,j,k,bi,bj)*dTdy(bi,bj) ENDDO ENDDO c /// end for x /// DO j=jMin,jMax DO i=iMin,iMax fZon(i,j,k,bi,bj) = fZon(i,j,k,bi,bj) + df(i,j,bi,bj) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax df(i,j,bi,bj) = 0. ENDDO ENDDO c ///gmredi_ytr/// DO j=jMin,jMax DO i=iMin,iMax df(i,j,bi,bj) = df(i,j,bi,bj) & -yA(i,j,bi,bj) & *smooth3D_Kvy(i,j,k,bi,bj) & *recip_dyC(i,j,bi,bj) & *(fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj)) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax dTdz(bi,bj) = 0.5*( & +0.5*recip_drC(k)* & ( maskC(i,j-1,k,bi,bj)* & (fld_in(i,j-1,MAX(k-1,1),bi,bj)-fld_in(i,j-1,k,bi,bj)) & +maskC(i, j ,k,bi,bj)* & (fld_in(i, j ,MAX(k-1,1),bi,bj)-fld_in(i, j ,k,bi,bj)) & ) & +0.5*recip_drC(MIN(k+1,Nr))* & ( maskC(i,j-1,MIN(k+1,Nr),bi,bj)* & (fld_in(i,j-1,k,bi,bj)-fld_in(i,j-1,MIN(k+1,Nr),bi,bj)) & +maskC(i, j ,MIN(k+1,Nr),bi,bj)* & (fld_in(i, j ,k,bi,bj)-fld_in(i, j ,MIN(k+1,Nr),bi,bj)) & ) ) df(i,j,bi,bj) = df(i,j,bi,bj) & - yA(i,j,bi,bj)*smooth3D_Kvz(i,j,k,bi,bj)*dTdz(bi,bj) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax dTdx(bi,bj) = 0.5*( & +0.5*(maskW(i+1,j,k,bi,bj) & *recip_dxC(i+1,j,bi,bj)* & (fld_in(i+1,j,k,bi,bj)-fld_in(i,j,k,bi,bj)) & +maskW(i,j,k,bi,bj) & *recip_dxC(i,j,bi,bj)* & (fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj))) & +0.5*(maskW(i+1,j-1,k,bi,bj) & *recip_dxC(i+1,j,bi,bj)* & (fld_in(i+1,j-1,k,bi,bj)-fld_in(i,j-1,k,bi,bj)) & +maskW(i,j-1,k,bi,bj) & *recip_dxC(i,j,bi,bj)* & (fld_in(i,j-1,k,bi,bj)-fld_in(i-1,j-1,k,bi,bj))) & ) df(i,j,bi,bj) = df(i,j,bi,bj) & - yA(i,j,bi,bj)*smooth3D_Kvx(i,j,k,bi,bj)*dTdx(bi,bj) ENDDO ENDDO c /// end for y /// DO j=jMin,jMax DO i=iMin,iMax fMer(i,j,k,bi,bj) = fMer(i,j,k,bi,bj) + df(i,j,bi,bj) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax df(i,j,bi,bj) = 0. ENDDO ENDDO c /// GAD_DIFF_R /// if (.NOT. smooth3DdoImpldiff ) then IF (k.gt.1) then DO j=jMin,jMax DO i=iMin,iMax df(i,j,bi,bj) = & -_rA(i,j,bi,bj) & *smooth3D_kappaR(i,j,k,bi,bj)*recip_drC(k) & *(fld_in(i,j,k,bi,bj) & -fld_in(i,j,k-1,bi,bj))*rkSign ENDDO ENDDO ENDIF endif c ///gmredi rtrans/// IF (K.GT.1) THEN DO j=jMin,jMax DO i=iMin,iMax dTdx(bi,bj) = 0.5*( & +0.5*(maskW(i+1,j,k,bi,bj) & *recip_dxC(i+1,j,bi,bj)* & (fld_in(i+1,j,k,bi,bj)-fld_in(i,j,k,bi,bj)) & +maskW(i,j,k,bi,bj) & *recip_dxC(i,j,bi,bj)* & (fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj))) & +0.5*(maskW(i+1,j,k-1,bi,bj) & *recip_dxC(i+1,j,bi,bj)* & (fld_in(i+1,j,k-1,bi,bj)-fld_in(i,j,k-1,bi,bj)) & +maskW(i,j,k-1,bi,bj) & *recip_dxC(i,j,bi,bj)* & (fld_in(i,j,k-1,bi,bj)-fld_in(i-1,j,k-1,bi,bj))) & ) dTdy(bi,bj) = 0.5*( & +0.5*(maskS(i,j,k,bi,bj) & *recip_dyC(i,j,bi,bj)* & (fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj)) & +maskS(i,j+1,k,bi,bj) & *recip_dyC(i,j+1,bi,bj)* & (fld_in(i,j+1,k,bi,bj)-fld_in(i,j,k,bi,bj))) & +0.5*(maskS(i,j,k-1,bi,bj) & *recip_dyC(i,j,bi,bj)* & (fld_in(i,j,k-1,bi,bj)-fld_in(i,j-1,k-1,bi,bj)) & +maskS(i,j+1,k-1,bi,bj) & *recip_dyC(i,j+1,bi,bj)* & (fld_in(i,j+1,k-1,bi,bj)-fld_in(i,j,k-1,bi,bj))) & ) df(i,j,bi,bj) = df(i,j,bi,bj) & - rA(i,j,bi,bj) & *( smooth3D_Kwx(i,j,k,bi,bj)*dTdx(bi,bj) & +smooth3D_Kwy(i,j,k,bi,bj)*dTdy(bi,bj) ) ENDDO ENDDO ENDIF c /// end for r /// IF (K.GT.1) THEN DO j=jMin,jMax DO i=iMin,iMax fVerT(i,j,k-1,bi,bj) = fVerT(i,j,k-1,bi,bj) + & df(i,j,bi,bj)*maskUp(i,j,bi,bj) ENDDO ENDDO ENDIF DO j=jMin,jMax DO i=iMin,iMax df(i,j,bi,bj) = 0. ENDDO ENDDO ENDDO !k ENDDO !bi ENDDO !bj c these exchanges are crucial: CALL EXCH_UV_XYZ_RL(fZon,fMer,.TRUE.,myThid) CALL EXCH_XYZ_RL ( fVerT, myThid ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) c 3rd k loop: Divergence of fluxes DO k=1,Nr IF (K.GT.1) THEN DO j=jMin,jMax DO i=iMin,iMax gt_in(i,j,k,bi,bj)=gt_in(i,j,k,bi,bj) & -smooth_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj) & *( (fZon(i+1,j,k,bi,bj)-fZon(i,j,k,bi,bj)) & +(fMer(i,j+1,k,bi,bj)-fMer(i,j,k,bi,bj)) & +(fVerT(i,j,k,bi,bj)-fVerT(i,j,k-1,bi,bj))*rkSign & ) ENDDO ENDDO ELSE DO j=jMin,jMax DO i=iMin,iMax gt_in(i,j,k,bi,bj)=gt_in(i,j,k,bi,bj) & -smooth_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj) & *( (fZon(i+1,j,k,bi,bj)-fZon(i,j,k,bi,bj)) & +(fMer(i,j+1,k,bi,bj)-fMer(i,j,k,bi,bj)) & +(fVerT(i,j,k,bi,bj))*rkSign & ) ENDDO ENDDO ENDIF ENDDO ENDDO ENDDO CALL EXCH_XYZ_RL ( gt_in , myThid ) END