pkg/zonal_filt/zonal_filt_init.F

42c525bfb4 Alis*0001 #include "ZONAL_FILT_OPTIONS.h"
                0002
0c8a064d26 Jean*0003 CBOP
                0004 C     !ROUTINE: ZONAL_FILT_INIT
                0005 C     !INTERFACE:
                0006       SUBROUTINE ZONAL_FILT_INIT( myThid )
                0007
                0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
6d99f9d2b3 Jean*0010 C     | S/R ZONAL_FILT_INIT
                0011 C     | o Initialise FFT filter for latitude circle.
0c8a064d26 Jean*0012 C     *==========================================================*
6d99f9d2b3 Jean*0013 C     | The details of particular FFT libraries may differ.
                0014 C     | Changing to a different library may entail modifying the
                0015 C     | code here. However, the broad process is usually the
                0016 C     | same.
                0017 C     | Note - Fourier modes for sNx and sNx+1 are damped in the
                0018 C     |        same way. This is because we have not implemented
                0019 C     |        a scheme that sets the damping factor for the
                0020 C     |        highest wave number for odd sNx. Instead the
                0021 C     |        highest wave number for odd sNx. Instead only
                0022 C     |        wave numbers 1:INT(sNx/2) are partially damped.
                0023 C     |        Wave number sNx/2 (if it exists) is removed
                0024 C     |        altogether.
0c8a064d26 Jean*0025 C     *==========================================================*
                0026 C     \ev
                0027
                0028 C     !USES:
42c525bfb4 Alis*0029       IMPLICIT NONE
                0030
0c8a064d26 Jean*0031 C     === Global variables ===
42c525bfb4 Alis*0032 #include "SIZE.h"
                0033 #include "EEPARAMS.h"
                0034 #include "PARAMS.h"
                0035 #include "GRID.h"
                0036 #include "ZONAL_FILT.h"
                0037 #include "FFTPACK.h"
                0038
0c8a064d26 Jean*0039 C     !INPUT/OUTPUT PARAMETERS:
42c525bfb4 Alis*0040 C     == Routine arguments ==
0c8a064d26 Jean*0041 C     myThid :: my Thread Id number
42c525bfb4 Alis*0042       INTEGER myThid
0c8a064d26 Jean*0043 CEOP
42c525bfb4 Alis*0044
                0045 #ifdef ALLOW_ZONAL_FILT
0c8a064d26 Jean*0046 C     !LOCAL VARIABLES:
42c525bfb4 Alis*0047 C     == Local variables ==
6d99f9d2b3 Jean*0048 C     alpha  :: Used to evaluate frequency and latitude dependent
                0049 C               amplitude damping factor.
                0050 C     wvNum  :: Wave number
                0051 C     lat    :: Temporary holding latitude
                0052 C     nWv    :: No. of waves that fit on grid.
                0053 C     msgBuf :: Informational/error message buffer
0c8a064d26 Jean*0054 c     _RL alpha, wvNum
                0055 c     INTEGER nWv
                0056       INTEGER i, j, bi, bj
6d99f9d2b3 Jean*0057       CHARACTER*(MAX_LEN_MBUF) msgBuf
0c8a064d26 Jean*0058
                0059 C     !FUNCTIONS:
                0060       _RL ampfact
                0061       _RS lat
6d99f9d2b3 Jean*0062       ampfact(lat,i) = MIN( oneRL,
0c8a064d26 Jean*0063      &   ( COS( ABS(lat)*deg2rad )
                0064      &      /COS( zonal_filt_lat*deg2rad ) )**zonal_filt_cospow
                0065      &      /(SIN( PI*FLOAT(i)/FLOAT(Nx) ) )**zonal_filt_sinpow
42c525bfb4 Alis*0066      &   )
                0067
0c8a064d26 Jean*0068 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0069
6d99f9d2b3 Jean*0070       IF ( sNx.NE.Nx ) THEN
                0071          WRITE(msgBuf,'(A,I3,A)')
                0072      &    'S/R ZONAL_FILT_INIT: Multi-tiles ( nSx*nPx=', nSx*nPx, ' )'
                0073          CALL PRINT_ERROR( msgBuf, myThid )
                0074          WRITE(msgBuf,'(A)')
                0075      &    ' in Zonal (X) dir. not implemented in Zonal-Filter code'
                0076          CALL PRINT_ERROR( msgBuf, myThid )
                0077          STOP 'ABNORMAL END: S/R ZONAL_FILT_INIT'
                0078       ENDIF
                0079
1389d71047 Chri*0080       _BEGIN_MASTER(myThid)
42c525bfb4 Alis*0081 C     o Initialise specific library FFT package
                0082       DO bj=1,nSy
0c8a064d26 Jean*0083 c      CALL R8FFTI( Nx, FFTPACKWS(1,bj) )
42c525bfb4 Alis*0084        CALL R8FFTI1( Nx, FFTPACKWS2(1,bj), FFTPACKWS3(1,bj) )
                0085       ENDDO
                0086
                0087 C     o Set amplitude scale factor as function of latitude and mode number
                0088       DO bj=1,nSy
                0089        DO bi=1,nSx
                0090         DO j=1-oLy,sNy+Oly
0c8a064d26 Jean*0091          ampFactor(1,j,bi,bj) = oneRL
                0092          ampFactorV(1,j,bi,bj) = oneRL
42c525bfb4 Alis*0093          DO i=1,Nx/2-1
0c8a064d26 Jean*0094           ampFactor(2*i,j,bi,bj) = ampfact( yC(1,j,bi,bj) , I )
                0095 c         IF (ampFactor(2*i,j,bi,bj).LE..9) ampFactor(2*i,j,bi,bj)=0.
                0096           ampFactor(2*I+1,j,bi,bj) = ampFactor(2*i,j,bi,bj)
                0097           ampFactorV(2*i,j,bi,bj) = ampfact( yG(1,j,bi,bj) , I )
                0098 c         IF (ampFactorV(2*i,j,bi,bj).LE..9) ampFactorV(2*i,j,bi,bj)=0.
                0099           ampFactorV(2*I+1,j,bi,bj) = ampFactorV(2*i,j,bi,bj)
42c525bfb4 Alis*0100          ENDDO
d4c5f8fe62 Jean*0101
0c8a064d26 Jean*0102          i=Nx/2
d4c5f8fe62 Jean*0103          IF ( zonal_filt_mode2dx.EQ.0 ) THEN
0c8a064d26 Jean*0104            ampFactor(Nx,j,bi,bj) = ampfact( yC(1,j,bi,bj) , i )
                0105            ampFactorV(Nx,j,bi,bj) = ampfact( yG(1,j,bi,bj) , i )
d4c5f8fe62 Jean*0106          ELSE
0c8a064d26 Jean*0107            ampFactor(Nx,j,bi,bj) = 0.
                0108            ampFactorV(Nx,j,bi,bj) = 0.
d4c5f8fe62 Jean*0109          ENDIF
                0110
42c525bfb4 Alis*0111         ENDDO
                0112        ENDDO
                0113       ENDDO
1389d71047 Chri*0114       _END_MASTER(myThid)
                0115       CALL BAR2(myThid)
42c525bfb4 Alis*0116
                0117       CALL WRITE_REC_XY_RL( 'ampFactor', ampFactor, 1, 0, myThid )
                0118
                0119 #endif /* ALLOW_ZONAL_FILT */
                0120
                0121       RETURN
                0122       END