C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_psi_bvp.F,v 1.6 2016/10/26 00:49:05 jmc Exp $ C $Name: $ #include "GMREDI_OPTIONS.h" CBOP C !ROUTINE: GMREDI_CALC_PSI_BVP C !INTERFACE: SUBROUTINE GMREDI_CALC_PSI_BVP( I bi, bj, iMin, iMax, jMin, jMax, I sigmaX, sigmaY, sigmaR, I myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE GMREDI_CALC_PSI_BVP C | o Calculate stream-functions for GM bolus velocity using C | the BVP in Ferrari et al. (OM, 2010) C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "GMREDI.h" C !INPUT/OUTPUT PARAMETERS: C bi,bj :: Tile indices C iMin,iMax :: computation domain 1rst index bounds C jMin,jMax :: computation domain 2nd index bounds C sigmaX :: Zonal gradient of density C sigmaY :: Meridional gradient of density C sigmaR :: Vertical gradient of Pot.density (locally referenced) C myThid :: My Thread Id number INTEGER bi,bj,iMin,iMax,jMin,jMax _RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) INTEGER myThid CEOP #ifdef ALLOW_GMREDI #ifdef GM_BOLUS_BVP C !LOCAL VARIABLES: INTEGER i,j,k, km1 INTEGER errCode _RL half_K _RL sigmaX_W _RL sigmaY_W _RL dSigmaDrW _RL dSigmaDrS _RL wkb_cW(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL wkb_cS(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rPI, c2 #ifdef ALLOW_DIAGNOSTICS _RL tmpFac #endif CHARACTER*(MAX_LEN_MBUF) msgBuf PARAMETER ( rPI = 0.318309886183791 _d 0 ) C- Matrix elements for tri-diagonal solver _RL GM_a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL GM_b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL GM_c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) C- Initialization : <= done in S/R gmredi_init IF (GM_UseBVP) THEN C Initialize the WKB wave speeds to zero C We use c = int_{-H}^0 N dz/(GM_BVP_ModeNumber*PI) and have absorbed C a factor of g/rhoConst DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx wkb_cW(i,j) = 0. _d 0 wkb_cS(i,j) = 0. _d 0 ENDDO ENDDO C Surface BC : set to zero DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx GM_PsiX(i,j,1,bi,bj) = 0. _d 0 GM_PsiY(i,j,1,bi,bj) = 0. _d 0 ENDDO ENDDO C Initialise matrix to identity DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx GM_a3d(i,j,k) = 0. _d 0 GM_b3d(i,j,k) = 1. _d 0 GM_c3d(i,j,k) = 0. _d 0 ENDDO ENDDO ENDDO DO k=2,Nr km1 = k-1 half_K = GM_background_K & *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op25 C Gradient of Sigma below U and V points DO j=1-OLy,sNy+OLy DO i=1-OLx+1,sNx+OLx sigmaX_W = op5*( sigmaX(i,j,km1)+sigmaX(i,j,k) ) & *maskW(i,j,k,bi,bj) dSigmaDrW = op5*( sigmaR(i-1,j,k)+sigmaR(i,j,k) ) & *maskW(i,j,k,bi,bj) wkb_cW(i,j) = wkb_cW(i,j) & + SQRT(MAX( -dSigmaDrW, 0. _d 0 )) & *drC(k)*GM_BVP_rModeNumber*rPI C Part of main diagonal coming from zeroth order derivative GM_b3d(i,j,k) = MAX( -dSigmaDrW, GM_Small_Number ) C This is initially the RHS GM_PsiX(i,j,k,bi,bj) = half_K*sigmaX_W & *(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj)) ENDDO ENDDO ENDDO C Note: Use Dirichlet BC @ Surface & Bottom (whereas we use Neumann BC for C implicit diffusion/advection Pb). C Surface BC implementation: => keep non zero coeff @ k=2 C and set Psi=1 with Identity coeff @ k=1 C Same for bottom, except if kBottom=Nr (solver only process k=1:Nr) C should substract c3d(k=Nr)*Psi(k=Nr+1) to RHS @ k=Nr ; C However in our case this term is zero since Psi(k=Nr+1)=0 DO k=2,Nr km1 = k-1 DO j=1-OLy,sNy+OLy DO i=1-OLx+1,sNx+OLx IF ( maskW(i,j,k,bi,bj).NE.0. ) THEN c2 = MAX( wkb_cW(i,j)*wkb_cW(i,j), GM_BVP_cHat2Min ) GM_a3d(i,j,k) = -c2*recip_drC(k) & *recip_drF(km1)*recip_hFacW(i,j,km1,bi,bj) GM_b3d(i,j,k) = GM_b3d(i,j,k) & + c2*recip_drC(k) & *(recip_drF(km1)*recip_hFacW(i,j,km1,bi,bj) & +recip_drF(k)*recip_hFacW(i,j,k,bi,bj) ) GM_c3d(i,j,k) = -c2*recip_drC(k) & *recip_drF(k)*recip_hFacW(i,j,k,bi,bj) ELSE GM_b3d(i,j,k) = 1. _d 0 ENDIF ENDDO ENDDO ENDDO errCode = -1 CALL SOLVE_TRIDIAGONAL( iMin+1, iMax, jMin, jMax, & GM_a3d, GM_b3d, GM_c3d, & GM_PsiX(1-OLx,1-OLy,1,bi,bj), & errCode, bi, bj, myThid ) IF ( errCode .GT. 0 ) THEN WRITE(msgBuf,'(A)') & 'S/R GMREDI_CALC_PSI_BVP: matrix singular for PsiX' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R GMREDI_CALC_PSI_BVP' ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C Reset matrix to identity DO k=2,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx GM_a3d(i,j,k) = 0. _d 0 GM_b3d(i,j,k) = 1. _d 0 GM_c3d(i,j,k) = 0. _d 0 ENDDO ENDDO ENDDO C cut k loop in 2 to prevent bad optimisation with some compilers DO k=2,Nr km1 = k-1 half_K = GM_background_K & *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op25 DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx sigmaY_W = op5*( sigmaY(i,j,km1)+sigmaY(i,j,k) ) & *maskS(i,j,k,bi,bj) dSigmaDrS = op5*( sigmaR(i,j-1,k)+sigmaR(i,j,k) ) & *maskS(i,j,k,bi,bj) wkb_cS(i,j) = wkb_cS(i,j) & + SQRT(MAX( -dSigmaDrS, 0. _d 0 )) & *drC(k)*GM_BVP_rModeNumber*rPI C Part of main diagonal coming from zeroth order derivative GM_b3d(i,j,k) = MAX( -dSigmaDrS, GM_Small_Number ) C This is initially the RHS GM_PsiY(i,j,k,bi,bj) = half_K*sigmaY_W & *(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj)) ENDDO ENDDO ENDDO DO k=2,Nr km1 = k-1 DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx IF ( maskS(i,j,k,bi,bj).NE.0. ) THEN c2 = MAX( wkb_cS(i,j)*wkb_cS(i,j), GM_BVP_cHat2Min ) GM_a3d(i,j,k) = -c2*recip_drC(k) & *recip_drF(km1)*recip_hFacS(i,j,km1,bi,bj) GM_b3d(i,j,k) = GM_b3d(i,j,k) & + c2*recip_drC(k) & *(recip_drF(km1)*recip_hFacS(i,j,km1,bi,bj) & +recip_drF(k)*recip_hFacS(i,j,k,bi,bj) ) GM_c3d(i,j,k) = -c2*recip_drC(k) & *recip_drF(k)*recip_hFacS(i,j,k,bi,bj) ELSE GM_b3d(i,j,k) = 1. _d 0 ENDIF ENDDO ENDDO ENDDO errCode = -1 CALL SOLVE_TRIDIAGONAL( iMin, iMax, jMin+1, jMax, & GM_a3d, GM_b3d, GM_c3d, & GM_PsiY(1-OLx,1-OLy,1,bi,bj), & errCode, bi, bj, myThid ) IF ( errCode .GT. 0 ) THEN WRITE(msgBuf,'(A)') & 'S/R GMREDI_CALC_PSI_BVP: matrix singular for PsiY' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R GMREDI_CALC_PSI_BVP' ENDIF #ifdef ALLOW_DIAGNOSTICS C Write some diagnostics IF ( useDiagnostics ) THEN tmpFac = SQRT(gravity/rhoConst) CALL DIAGNOSTICS_SCALE_FILL( wkb_cW, tmpFac, 1, 'GM_BVPcW', & 0, 1, 2, bi, bj, myThid ) CALL DIAGNOSTICS_SCALE_FILL( wkb_cS, tmpFac, 1, 'GM_BVPcS', & 0, 1, 2, bi, bj, myThid ) ENDIF #endif ENDIF #endif /* GM_BOLUS_ADVEC */ #endif /* ALLOW_GMREDI */ RETURN END