pkg/gmredi/gmredi_calc_psi_bvp.F

d0035fac68 Jean*0001 #include "GMREDI_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: GMREDI_CALC_PSI_BVP
                0005 C     !INTERFACE:
                0006       SUBROUTINE GMREDI_CALC_PSI_BVP(
                0007      I             bi, bj, iMin, iMax, jMin, jMax,
                0008      I             sigmaX, sigmaY, sigmaR,
                0009      I             myThid )
                0010
                0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
                0013 C     | SUBROUTINE GMREDI_CALC_PSI_BVP
                0014 C     | o Calculate stream-functions for GM bolus velocity using
                0015 C     |   the BVP in Ferrari et al. (OM, 2010)
                0016 C     *==========================================================*
                0017 C     \ev
                0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021
                0022 C     == Global variables ==
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
bdc4273caf Jean*0026 #include "GRID.h"
d0035fac68 Jean*0027 #include "GMREDI.h"
                0028
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     bi,bj   :: Tile indices
                0031 C   iMin,iMax :: computation domain 1rst index bounds
                0032 C   jMin,jMax :: computation domain 2nd  index bounds
                0033 C     sigmaX  :: Zonal      gradient of density
                0034 C     sigmaY  :: Meridional gradient of density
                0035 C     sigmaR  :: Vertical   gradient of Pot.density (locally referenced)
                0036 C     myThid  :: My Thread Id number
                0037       INTEGER bi,bj,iMin,iMax,jMin,jMax
19eb0e065a Jean*0038       _RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0039       _RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0040       _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
d0035fac68 Jean*0041       INTEGER myThid
                0042 CEOP
                0043
                0044 #ifdef ALLOW_GMREDI
                0045 #ifdef GM_BOLUS_BVP
                0046
                0047 C     !LOCAL VARIABLES:
                0048       INTEGER i,j,k, km1
                0049       INTEGER errCode
bdc4273caf Jean*0050       _RL half_K
d0035fac68 Jean*0051       _RL sigmaX_W
                0052       _RL sigmaY_W
                0053       _RL dSigmaDrW
                0054       _RL dSigmaDrS
19eb0e065a Jean*0055       _RL wkb_cW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0056       _RL wkb_cS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d0035fac68 Jean*0057       _RL rPI, c2
                0058 #ifdef ALLOW_DIAGNOSTICS
                0059       _RL tmpFac
                0060 #endif
                0061       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0062
                0063       PARAMETER ( rPI    = 0.318309886183791 _d 0 )
                0064
                0065 C-    Matrix elements for tri-diagonal solver
19eb0e065a Jean*0066       _RL GM_a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL GM_b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0068       _RL GM_c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
d0035fac68 Jean*0069
                0070 C-    Initialization : <= done in S/R gmredi_init
                0071       IF (GM_UseBVP) THEN
                0072
                0073 C     Initialize the WKB wave speeds to zero
                0074 C     We use c = int_{-H}^0 N dz/(GM_BVP_ModeNumber*PI) and have absorbed
                0075 C     a factor of g/rhoConst
19eb0e065a Jean*0076        DO j=1-OLy,sNy+OLy
                0077         DO i=1-OLx,sNx+OLx
                0078           wkb_cW(i,j) = 0. _d 0
                0079           wkb_cS(i,j) = 0. _d 0
d0035fac68 Jean*0080         ENDDO
                0081        ENDDO
                0082
                0083 C      Surface BC : set to zero
19eb0e065a Jean*0084        DO j=1-OLy,sNy+OLy
                0085         DO i=1-OLx,sNx+OLx
                0086           GM_PsiX(i,j,1,bi,bj) = 0. _d 0
                0087           GM_PsiY(i,j,1,bi,bj) = 0. _d 0
d0035fac68 Jean*0088         ENDDO
                0089        ENDDO
                0090
70e9ad9048 Jean*0091 C      Initialise matrix to identity
                0092        DO k=1,Nr
19eb0e065a Jean*0093         DO j=1-OLy,sNy+OLy
                0094          DO i=1-OLx,sNx+OLx
                0095           GM_a3d(i,j,k) = 0. _d 0
                0096           GM_b3d(i,j,k) = 1. _d 0
                0097           GM_c3d(i,j,k) = 0. _d 0
                0098          ENDDO
                0099         ENDDO
70e9ad9048 Jean*0100        ENDDO
                0101
                0102        DO k=2,Nr
                0103         km1 = k-1
bdc4273caf Jean*0104         half_K = GM_background_K
                0105      &         *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op25
d0035fac68 Jean*0106 C      Gradient of Sigma below U and V points
19eb0e065a Jean*0107         DO j=1-OLy,sNy+OLy
                0108          DO i=1-OLx+1,sNx+OLx
d0035fac68 Jean*0109           sigmaX_W = op5*( sigmaX(i,j,km1)+sigmaX(i,j,k) )
8233d0ceb9 Jean*0110      &                  *maskW(i,j,km1,bi,bj)*maskW(i,j,k,bi,bj)
d0035fac68 Jean*0111           dSigmaDrW = op5*( sigmaR(i-1,j,k)+sigmaR(i,j,k) )
8233d0ceb9 Jean*0112      &                  *maskW(i,j,km1,bi,bj)*maskW(i,j,k,bi,bj)
d0035fac68 Jean*0113
                0114           wkb_cW(i,j) = wkb_cW(i,j)
                0115      &                + SQRT(MAX( -dSigmaDrW, 0. _d 0 ))
                0116      &                 *drC(k)*GM_BVP_rModeNumber*rPI
                0117
                0118 C      Part of main diagonal coming from zeroth order derivative
70e9ad9048 Jean*0119           GM_b3d(i,j,k) = MAX( -dSigmaDrW, GM_Small_Number )
d0035fac68 Jean*0120
                0121 C      This is initially the RHS
bdc4273caf Jean*0122           GM_PsiX(i,j,k,bi,bj) = half_K*sigmaX_W
                0123      &       *(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj))
d0035fac68 Jean*0124          ENDDO
                0125         ENDDO
                0126        ENDDO
                0127
                0128 C Note: Use Dirichlet BC @ Surface & Bottom (whereas we use Neumann BC for
                0129 C       implicit diffusion/advection Pb).
                0130 C     Surface BC implementation: => keep non zero coeff @ k=2
                0131 C                               and set Psi=1 with Identity coeff @ k=1
                0132 C     Same for bottom, except if kBottom=Nr (solver only process k=1:Nr)
                0133 C       should substract c3d(k=Nr)*Psi(k=Nr+1) to RHS @ k=Nr ;
                0134 C       However in our case this term is zero since Psi(k=Nr+1)=0
                0135        DO k=2,Nr
                0136         km1 = k-1
19eb0e065a Jean*0137         DO j=1-OLy,sNy+OLy
                0138          DO i=1-OLx+1,sNx+OLx
8233d0ceb9 Jean*0139           IF ( maskW(i,j,km1,bi,bj).NE.zeroRS .AND.
                0140      &         maskW(i,j, k, bi,bj).NE.zeroRS ) THEN
d0035fac68 Jean*0141             c2 = MAX( wkb_cW(i,j)*wkb_cW(i,j), GM_BVP_cHat2Min )
bdc4273caf Jean*0142             GM_a3d(i,j,k) = -c2*recip_drC(k)
                0143      &                      *recip_drF(km1)*recip_hFacW(i,j,km1,bi,bj)
d0035fac68 Jean*0144             GM_b3d(i,j,k) = GM_b3d(i,j,k)
bdc4273caf Jean*0145      &                    + c2*recip_drC(k)
                0146      &                     *(recip_drF(km1)*recip_hFacW(i,j,km1,bi,bj)
                0147      &                      +recip_drF(k)*recip_hFacW(i,j,k,bi,bj) )
                0148             GM_c3d(i,j,k) = -c2*recip_drC(k)
                0149      &                      *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
70e9ad9048 Jean*0150           ELSE
                0151             GM_b3d(i,j,k) = 1. _d 0
d0035fac68 Jean*0152           ENDIF
                0153          ENDDO
                0154         ENDDO
                0155        ENDDO
                0156
fa8d87e2db Jean*0157        errCode = -1
8a58850ca8 Jean*0158        CALL SOLVE_TRIDIAGONAL( iMin+1, iMax, jMin, jMax,
                0159      &                         GM_a3d, GM_b3d, GM_c3d,
                0160      &                         GM_PsiX(1-OLx,1-OLy,1,bi,bj),
                0161      &                         errCode, bi, bj, myThid )
d0035fac68 Jean*0162
                0163        IF ( errCode .GT. 0 ) THEN
                0164         WRITE(msgBuf,'(A)')
                0165      &  'S/R GMREDI_CALC_PSI_BVP: matrix singular for PsiX'
                0166         CALL PRINT_ERROR( msgBuf, myThid )
                0167         STOP 'ABNORMAL END: S/R GMREDI_CALC_PSI_BVP'
                0168        ENDIF
                0169
                0170 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0171
70e9ad9048 Jean*0172 C      Reset matrix to identity
d0035fac68 Jean*0173        DO k=2,Nr
19eb0e065a Jean*0174         DO j=1-OLy,sNy+OLy
                0175          DO i=1-OLx,sNx+OLx
                0176           GM_a3d(i,j,k) = 0. _d 0
                0177           GM_b3d(i,j,k) = 1. _d 0
                0178           GM_c3d(i,j,k) = 0. _d 0
                0179          ENDDO
                0180         ENDDO
70e9ad9048 Jean*0181        ENDDO
                0182 C      cut k loop in 2 to prevent bad optimisation with some compilers
                0183        DO k=2,Nr
                0184         km1 = k-1
bdc4273caf Jean*0185         half_K = GM_background_K
                0186      &         *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op25
19eb0e065a Jean*0187         DO j=1-OLy+1,sNy+OLy
                0188          DO i=1-OLx,sNx+OLx
d0035fac68 Jean*0189           sigmaY_W = op5*( sigmaY(i,j,km1)+sigmaY(i,j,k) )
8233d0ceb9 Jean*0190      &                  *maskS(i,j,km1,bi,bj)*maskS(i,j,k,bi,bj)
d0035fac68 Jean*0191           dSigmaDrS = op5*( sigmaR(i,j-1,k)+sigmaR(i,j,k) )
8233d0ceb9 Jean*0192      &                  *maskS(i,j,km1,bi,bj)*maskS(i,j,k,bi,bj)
d0035fac68 Jean*0193
                0194           wkb_cS(i,j) = wkb_cS(i,j)
                0195      &                + SQRT(MAX( -dSigmaDrS, 0. _d 0 ))
                0196      &                 *drC(k)*GM_BVP_rModeNumber*rPI
                0197
                0198 C      Part of main diagonal coming from zeroth order derivative
70e9ad9048 Jean*0199           GM_b3d(i,j,k) = MAX( -dSigmaDrS, GM_Small_Number )
d0035fac68 Jean*0200
                0201 C      This is initially the RHS
bdc4273caf Jean*0202           GM_PsiY(i,j,k,bi,bj) = half_K*sigmaY_W
                0203      &       *(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj))
d0035fac68 Jean*0204          ENDDO
                0205         ENDDO
                0206        ENDDO
                0207
                0208        DO k=2,Nr
                0209         km1 = k-1
19eb0e065a Jean*0210         DO j=1-OLy+1,sNy+OLy
                0211          DO i=1-OLx,sNx+OLx
8233d0ceb9 Jean*0212           IF ( maskS(i,j,km1,bi,bj).NE.zeroRS .AND.
                0213      &         maskS(i,j, k, bi,bj).NE.zeroRS ) THEN
d0035fac68 Jean*0214             c2 = MAX( wkb_cS(i,j)*wkb_cS(i,j), GM_BVP_cHat2Min )
bdc4273caf Jean*0215             GM_a3d(i,j,k) = -c2*recip_drC(k)
                0216      &                      *recip_drF(km1)*recip_hFacS(i,j,km1,bi,bj)
d0035fac68 Jean*0217             GM_b3d(i,j,k) = GM_b3d(i,j,k)
bdc4273caf Jean*0218      &                    + c2*recip_drC(k)
                0219      &                     *(recip_drF(km1)*recip_hFacS(i,j,km1,bi,bj)
                0220      &                      +recip_drF(k)*recip_hFacS(i,j,k,bi,bj) )
                0221             GM_c3d(i,j,k) = -c2*recip_drC(k)
                0222      &                      *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
70e9ad9048 Jean*0223           ELSE
                0224             GM_b3d(i,j,k) = 1. _d 0
d0035fac68 Jean*0225           ENDIF
                0226          ENDDO
                0227         ENDDO
                0228        ENDDO
                0229
fa8d87e2db Jean*0230        errCode = -1
8a58850ca8 Jean*0231        CALL SOLVE_TRIDIAGONAL( iMin, iMax, jMin+1, jMax,
                0232      &                         GM_a3d, GM_b3d, GM_c3d,
                0233      &                         GM_PsiY(1-OLx,1-OLy,1,bi,bj),
                0234      &                         errCode, bi, bj, myThid )
d0035fac68 Jean*0235
                0236        IF ( errCode .GT. 0 ) THEN
                0237         WRITE(msgBuf,'(A)')
                0238      &  'S/R GMREDI_CALC_PSI_BVP: matrix singular for PsiY'
                0239         CALL PRINT_ERROR( msgBuf, myThid )
                0240         STOP 'ABNORMAL END: S/R GMREDI_CALC_PSI_BVP'
                0241        ENDIF
                0242
                0243 #ifdef ALLOW_DIAGNOSTICS
                0244 C     Write some diagnostics
                0245        IF ( useDiagnostics ) THEN
                0246         tmpFac = SQRT(gravity/rhoConst)
                0247         CALL DIAGNOSTICS_SCALE_FILL( wkb_cW, tmpFac, 1, 'GM_BVPcW',
                0248      &                               0, 1, 2, bi, bj, myThid )
                0249         CALL DIAGNOSTICS_SCALE_FILL( wkb_cS, tmpFac, 1, 'GM_BVPcS',
                0250      &                               0, 1, 2, bi, bj, myThid )
                0251        ENDIF
                0252 #endif
                0253
                0254       ENDIF
                0255 #endif /* GM_BOLUS_ADVEC */
                0256 #endif /* ALLOW_GMREDI */
                0257
                0258       RETURN
                0259       END