C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.16 2014/04/04 19:38:23 jmc Exp $
C $Name: $
#include "SHAP_FILT_OPTIONS.h"
CBOP
C !ROUTINE: SHAP_FILT_TRACER_S2
C !INTERFACE:
SUBROUTINE SHAP_FILT_TRACER_S2(
U field, tmpFld,
I nShapTr, exchInOut, kSize,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R SHAP_FILT_TRACER_S2
C | o Applies Shapiro filter to 2D field (cell center).
C | o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
C | o Options for computational filter (no grid spacing)
C | or physical space filter (with grid spacing) or both.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SHAP_FILT.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments
C field :: cell-centered 2D field on which filter applies
C tmpFld :: working temporary array
C nShapTr :: (total) power of the filter for this tracer
C exchInOut :: apply Exchanges to fill overlap region:
C = 0 : do not apply Exch on neither input nor output field
C = 1 : apply Exch on input field
C (needed if input field has invalid overlap)
C = 2 : apply Exch on output field (after the filter)
C = 3 : apply Exch on both input & output field
C kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: Thread number for this instance of SHAP_FILT_TRACER_S2
INTEGER kSize
_RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
INTEGER nShapTr, exchInOut
_RL myTime
INTEGER myIter
INTEGER myThid
#ifdef ALLOW_SHAP_FILT
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER nShapComput
INTEGER bi,bj,k,i,j,n
_RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP
IF ( exchInOut.LT.0 .OR. exchInOut.GT.3 ) THEN
STOP 'S/R SHAP_FILT_TRACER_S2: exchInOut is wrong'
ENDIF
IF (nShapTr.GT.0) THEN
C-------
C Apply computational filter ^(nShap-nShapPhys) without grid factor
C then apply Physical filter ^nShapPhys with grid factors
C-------
nShapComput = nShapTr - nShapTrPhys
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO k=1,kSize
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C ( d_xx +d_yy )^n tmpFld
DO n=1,nShapTr
IF ( ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) .AND.
& ( n.GE.2 .OR. MOD(exchInOut,2).EQ.1 ) ) THEN
IF (kSize.EQ.Nr) THEN
_EXCH_XYZ_RL( tmpFld, myThid )
ELSEIF (kSize.EQ.1) THEN
_EXCH_XY_RL( tmpFld, myThid )
ELSE
STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
ENDIF
ENDIF
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO k=1,kSize
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Calculate gradient in X direction:
#ifndef ALLOW_AUTODIFF
IF ( .NOT.Shap_alwaysExchTr
& .AND. useCubedSphereExchange ) THEN
C to compute d/dx(tmpFld), fill corners with appropriate values:
CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
& tmpFld(1-OLx,1-OLy,k,bi,bj),
& bi,bj, myThid )
ENDIF
#endif
IF ( n.LE.nShapComput ) THEN
C- Computational space: del_i
DO j=0,sNy+1
DO i=0,sNx+2
tmpFdx(i,j) =
& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
& *_maskW(i,j,k,bi,bj)
ENDDO
ENDDO
ELSE
C- Physical space: grad_x
DO j=0,sNy+1
DO i=0,sNx+2
tmpFdx(i,j) =
& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
& *_hFacW(i,j,k,bi,bj)
& *dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
C-- Calculate gradient in Y direction:
#ifndef ALLOW_AUTODIFF
IF ( .NOT.Shap_alwaysExchTr
& .AND. useCubedSphereExchange ) THEN
C to compute d/dy(tmpFld), fill corners with appropriate values:
CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
& tmpFld(1-OLx,1-OLy,k,bi,bj),
& bi,bj, myThid )
ENDIF
#endif
IF ( n.LE.nShapComput ) THEN
C- Computational space: del_j
DO j=0,sNy+2
DO i=0,sNx+1
tmpFdy(i,j) =
& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
& *_maskS(i,j,k,bi,bj)
ENDDO
ENDDO
ELSE
C- Physical space: grad_y
DO j=0,sNy+2
DO i=0,sNx+1
tmpFdy(i,j) =
& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
& *_hFacS(i,j,k,bi,bj)
& *dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
C-- Calculate (d_xx + d_yy) tmpFld :
DO j=0,sNy+1
DO i=0,sNx+1
tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
& + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
ENDDO
ENDDO
C-- Computational space Filter
IF ( n.LE.nShapComput ) THEN
DO j=0,sNy+1
DO i=0,sNx+1
tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
ENDDO
ENDDO
C-- Physical space Filter
ELSEIF (Shap_TrLength.LE.0.) THEN
DO j=0,sNy+1
DO i=0,sNx+1
tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
& *recip_hFacC(i,j,k,bi,bj)
ENDDO
ENDDO
ELSE
DO j=0,sNy+1
DO i=0,sNx+1
tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
& *recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
& *Shap_TrLength*Shap_TrLength
ENDDO
ENDDO
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C end k,bi,bj loop:
ENDDO
ENDDO
ENDDO
C end loop n=1,nShapTr
ENDDO
C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO k=1,kSize
DO j=1,sNy
DO i=1,sNx
field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
IF ( exchInOut.GE.2 ) THEN
IF (kSize.EQ.Nr) THEN
_EXCH_XYZ_RL( field, myThid )
ELSEIF (kSize.EQ.1) THEN
_EXCH_XY_RL( field, myThid )
ELSE
STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
ENDIF
ENDIF
ENDIF
#endif /* ALLOW_SHAP_FILT */
RETURN
END