pkg/shap_filt/shap_filt_tracer_s2g.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2g.F,v 1.6 2005/02/15 00:21:12 jmc Exp $
C $Name:  $

#include "SHAP_FILT_OPTIONS.h"
 
CBOP
C     !ROUTINE: SHAP_FILT_TRACER_S2G
C     !INTERFACE:
      SUBROUTINE SHAP_FILT_TRACER_S2G(
     U           field, tmpFld,
     I           nShapTr, kSize, myTime, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SHAP_FILT_TRACER_S2G
C     | o Applies Shapiro filter to tracer field (cell center).
C     | o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
C     |     with grid spacing (physical space filter)
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     kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
C     myTime :: Current time in simulation
C     myThid :: Thread number for this instance of SHAP_FILT_TRACER_S2G
      INTEGER nShapTr, 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)
      _RL     myTime
      INTEGER myThid
 
#ifdef ALLOW_SHAP_FILT

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER bi,bj,K,I,J,N
      _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

      IF (nShapTr.gt.0) THEN

        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 (kSize.EQ.Nr) THEN
          _EXCH_XYZ_R8( tmpFld, myThid )
        ELSE
          _EXCH_XY_R8( tmpFld, myThid )
        ENDIF

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO K=1,kSize

           DO J=1,sNy
            DO I=1,sNx
             tmpGrd(i,j) = -0.125*(
     &        ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
     &            *_hFacW(i+1,j,k,bi,bj)
     &            *DYG(i+1,j,bi,bj)
     &            *recip_DXC(i+1,j,bi,bj)
     &       -( 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)
     &       +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
     &            *_hFacS(i,j+1,k,bi,bj)
     &            *DXG(i,j+1,bi,bj)
     &            *recip_DYC(i,j+1,bi,bj)
     &       -( 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)
     &       )*recip_hFacC(i,j,k,bi,bj)
            ENDDO
           ENDDO

           IF (Shap_TrLength.EQ.0.) THEN
            DO J=1,sNy
             DO I=1,sNx
              tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)
            ENDDO
           ENDDO
           ELSE
            DO J=1,sNy
             DO I=1,sNx
              tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)
     &         *recip_rA(i,j,bi,bj)
     &         *Shap_TrLength*Shap_TrLength
            ENDDO
           ENDDO
           ENDIF

          ENDDO
         ENDDO
        ENDDO

       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 (kSize.EQ.Nr) THEN
          _EXCH_XYZ_R8( field, myThid )
        ELSEIF (kSize.EQ.1) THEN
          _EXCH_XY_R8( field, myThid )
        ELSE
          STOP 'S/R SHAP_FILT_TRACER_S4: kSize is wrong'
        ENDIF

      ENDIF
#endif /* ALLOW_SHAP_FILT */

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2g.F,v 1.6 2005/02/15 00:21:12 jmc Exp $
2	C $Name: $
3
4	#include "SHAP_FILT_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: SHAP_FILT_TRACER_S2G
8	C !INTERFACE:
9	SUBROUTINE SHAP_FILT_TRACER_S2G(
10	U field, tmpFld,
11	I nShapTr, kSize, myTime, myThid )
12	C !DESCRIPTION: \bv
13	C ==========================================================
14	C \| S/R SHAP_FILT_TRACER_S2G
15	C \| o Applies Shapiro filter to tracer field (cell center).
16	C \| o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
17	C \| with grid spacing (physical space filter)
18	C ==========================================================
19	C \ev
20
21	C !USES:
22	IMPLICIT NONE
23
24	C == Global variables ===
25	#include "SIZE.h"
26	#include "EEPARAMS.h"
27	#include "PARAMS.h"
28	#include "GRID.h"
29	#include "SHAP_FILT.h"
30
31	C !INPUT/OUTPUT PARAMETERS:
32	C == Routine arguments
33	C field :: cell-centered 2D field on which filter applies
34	C tmpFld :: working temporary array
35	C nShapTr :: (total) power of the filter for this tracer
36	C kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
37	C myTime :: Current time in simulation
38	C myThid :: Thread number for this instance of SHAP_FILT_TRACER_S2G
39	INTEGER nShapTr, kSize
40	_RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
41	_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
42	_RL myTime
43	INTEGER myThid
44
45	#ifdef ALLOW_SHAP_FILT
46
47	C !LOCAL VARIABLES:
48	C == Local variables ==
49	INTEGER bi,bj,K,I,J,N
50	_RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51	CEOP
52
53	IF (nShapTr.gt.0) THEN
54
55	DO bj=myByLo(myThid),myByHi(myThid)
56	DO bi=myBxLo(myThid),myBxHi(myThid)
57	DO K=1,kSize
58	DO J=1-OLy,sNy+OLy
59	DO I=1-OLx,sNx+OLx
60	tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
61	ENDDO
62	ENDDO
63	ENDDO
64	ENDDO
65	ENDDO
66
67
68	C ( d_xx +d_yy )^n tmpFld
69
70	DO N=1,nShapTr
71
72	IF (kSize.EQ.Nr) THEN
73	_EXCH_XYZ_R8( tmpFld, myThid )
74	ELSE
75	_EXCH_XY_R8( tmpFld, myThid )
76	ENDIF
77
78	DO bj=myByLo(myThid),myByHi(myThid)
79	DO bi=myBxLo(myThid),myBxHi(myThid)
80	DO K=1,kSize
81
82	DO J=1,sNy
83	DO I=1,sNx
84	tmpGrd(i,j) = -0.125*(
85	& ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
86	& *_hFacW(i+1,j,k,bi,bj)
87	& *DYG(i+1,j,bi,bj)
88	& *recip_DXC(i+1,j,bi,bj)
89	& -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
90	& *_hFacW( i ,j,k,bi,bj)
91	& *DYG( i ,j,bi,bj)
92	& *recip_DXC( i ,j,bi,bj)
93	& +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
94	& *_hFacS(i,j+1,k,bi,bj)
95	& *DXG(i,j+1,bi,bj)
96	& *recip_DYC(i,j+1,bi,bj)
97	& -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
98	& *_hFacS(i, j ,k,bi,bj)
99	& *DXG(i, j ,bi,bj)
100	& *recip_DYC(i, j ,bi,bj)
101	& )*recip_hFacC(i,j,k,bi,bj)
102	ENDDO
103	ENDDO
104
105	IF (Shap_TrLength.EQ.0.) THEN
106	DO J=1,sNy
107	DO I=1,sNx
108	tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)
109	ENDDO
110	ENDDO
111	ELSE
112	DO J=1,sNy
113	DO I=1,sNx
114	tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)
115	& *recip_rA(i,j,bi,bj)
116	& Shap_TrLengthShap_TrLength
117	ENDDO
118	ENDDO
119	ENDIF
120
121	ENDDO
122	ENDDO
123	ENDDO
124
125	ENDDO
126
127	C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
128	DO bj=myByLo(myThid),myByHi(myThid)
129	DO bi=myBxLo(myThid),myBxHi(myThid)
130	DO K=1,kSize
131	DO J=1,sNy
132	DO I=1,sNx
133	field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
134	& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
135	tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
136	ENDDO
137	ENDDO
138	ENDDO
139	ENDDO
140	ENDDO
141
142	IF (kSize.EQ.Nr) THEN
143	_EXCH_XYZ_R8( field, myThid )
144	ELSEIF (kSize.EQ.1) THEN
145	_EXCH_XY_R8( field, myThid )
146	ELSE
147	STOP 'S/R SHAP_FILT_TRACER_S4: kSize is wrong'
148	ENDIF
149
150	ENDIF
151	#endif /* ALLOW_SHAP_FILT */
152
153	RETURN
154	END