Back to home page

MITgcm

 
 

    


File indexing completed on 2023-09-06 05:10:23 UTC

view on githubraw file Latest commit 9c32132b on 2023-09-05 01:25:59 UTC
79074ef66b Jean*0001 #include "MOM_COMMON_OPTIONS.h"
                0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
                0005 
                0006 CBOP
754a8d68a6 Jean*0007 C !ROUTINE: MOM_U_BOTDRAG_COEFF
79074ef66b Jean*0008 
                0009 C !INTERFACE: ==========================================================
754a8d68a6 Jean*0010       SUBROUTINE MOM_U_BOTDRAG_COEFF(
79074ef66b Jean*0011      I                 bi, bj, k, inp_KE,
                0012      I                 uFld, vFld, kappaRU,
                0013      U                 KE,
                0014      O                 cDrag,
                0015      I                 myIter, myThid )
                0016 
                0017 C !DESCRIPTION:
                0018 C Compute bottom-drag coefficient (Cd) for U component momentum,
                0019 C   such that bottom stress: taux_{bot} = -Cd * U_{bot} * rUnit2mass ;
                0020 C include linear and quadratic bottom drag and friction (no-slip BC) at bottom
                0021 
                0022 C !USES: ===============================================================
                0023       IMPLICIT NONE
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
ab47de63dc Mart*0028 #ifdef ALLOW_BOTTOMDRAG_ROUGHNESS
                0029 # include "MOM_VISC.h"
                0030 #endif
79074ef66b Jean*0031 #ifdef ALLOW_CTRL
                0032 # include "CTRL_FIELDS.h"
                0033 #endif
                0034 
                0035 C !INPUT PARAMETERS: ===================================================
                0036 C  bi,bj          :: tile indices
                0037 C  k              :: vertical level to process
                0038 C  inp_KE         :: =T : KE is provided as input ; =F : to compute here
                0039 C  uFld           :: velocity, zonal component
                0040 C  vFld           :: velocity, meridional component
                0041 C  kappaRU        :: vertical viscosity
                0042 C  KE             :: Kinetic energy (input when inp_KE = T)
                0043 C  myIter         :: current iteration number
                0044 C  myThid         :: my Thread Id number
                0045       INTEGER bi, bj, k
                0046       LOGICAL inp_KE
                0047       _RL uFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0048       _RL vFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049       _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0050       _RL KE     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0051       INTEGER myIter, myThid
                0052 
                0053 C !OUTPUT PARAMETERS: ==================================================
                0054 C  KE             :: Kinetic energy (output when inp_KE = F)
                0055 C  cDrag          :: bottom drag coefficient
                0056       _RL cDrag  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057 
                0058 C !LOCAL VARIABLES: ====================================================
                0059 C  i,j            :: loop indices
                0060       INTEGER i, j
                0061       INTEGER kDown,kLowF,kBottom
                0062       _RL viscFac, dragFac, uSq
                0063       _RL recDrC
                0064 CEOP
                0065 
                0066 C-  No-slip BCs impose a drag at bottom
                0067       viscFac = 0.
                0068       IF (no_slip_bottom) viscFac = 2.
                0069 
                0070 c     DO k= k1,k2
                0071 
                0072        IF ( usingZCoords ) THEN
                0073         kBottom = Nr
                0074         kDown   = MIN(k+1,Nr)
                0075         kLowF   = k+1
                0076 c       dragFac = mass2rUnit*rhoConst
                0077 c       dragFac = wUnit2rVel(k+1)
                0078         dragFac = 1. _d 0
                0079        ELSE
                0080         kBottom = 1
                0081         kDown   = MAX(k-1,1)
                0082         kLowF   = k
                0083         dragFac = mass2rUnit*rhoConst
                0084 c       dragFac = wUnit2rVel(k)
                0085        ENDIF
                0086        IF ( k.EQ.kBottom ) THEN
                0087         recDrC = recip_drF(k)
                0088        ELSE
                0089         recDrC = recip_drC(kLowF)
                0090        ENDIF
                0091 
9c32132b19 Jean*0092 #ifdef ALLOW_AUTODIFF
                0093 C--   Initialise drag-coeff
                0094        DO j=1-OLy,sNy+OLy
                0095         DO i=1-OLx,sNx+OLx
                0096           cDrag(i,j) = 0. _d 0
                0097         ENDDO
                0098        ENDDO
                0099 #endif
                0100 
79074ef66b Jean*0101 C--   Linear bottom drag contribution to cDrag:
                0102        DO j=1-OLy,sNy+OLy
                0103         DO i=1-OLx+1,sNx+OLx
                0104           cDrag(i,j) = bottomDragLinear*dragFac
                0105 #ifdef ALLOW_BOTTOMDRAG_CONTROL
                0106      &      + halfRL*( bottomDragFld(i-1,j,bi,bj)
                0107      &               + bottomDragFld(i,j,bi,bj) )*dragFac
                0108 #endif
                0109         ENDDO
                0110        ENDDO
                0111 
                0112 C--   Add friction at the bottom (no-slip BC) to cDrag:
                0113        IF ( no_slip_bottom .AND. bottomVisc_pCell ) THEN
                0114 C-    bottom friction accounts for true distance (including hFac) to the bottom
                0115         DO j=1-OLy,sNy+OLy-1
                0116          DO i=1-OLx+1,sNx+OLx-1
                0117            cDrag(i,j) = cDrag(i,j)
                0118      &      + kappaRU(i,j,kLowF)*recDrC*viscFac
                0119      &                          *_recip_hFacW(i,j,k,bi,bj)
                0120          ENDDO
                0121         ENDDO
                0122        ELSEIF ( no_slip_bottom ) THEN
                0123 C-    ignores partial-cell reduction of the distance to the bottom
                0124         DO j=1-OLy,sNy+OLy-1
                0125          DO i=1-OLx+1,sNx+OLx-1
                0126            cDrag(i,j) = cDrag(i,j)
                0127      &      + kappaRU(i,j,kLowF)*recDrC*viscFac
                0128          ENDDO
                0129         ENDDO
                0130        ENDIF
                0131 
                0132 C--   Add quadratic bottom drag to cDrag:
                0133        IF ( selectBotDragQuadr.EQ.0 ) THEN
                0134         IF ( .NOT.inp_KE ) THEN
                0135           DO j=1-OLy,sNy+OLy-1
                0136            DO i=1-OLx,sNx+OLx-1
                0137             KE(i,j) = 0.25*(
                0138      &          ( uFld( i , j )*uFld( i , j )*_hFacW(i,j,k,bi,bj)
                0139      &           +uFld(i+1, j )*uFld(i+1, j )*_hFacW(i+1,j,k,bi,bj) )
                0140      &        + ( vFld( i , j )*vFld( i , j )*_hFacS(i,j,k,bi,bj)
                0141      &           +vFld( i ,j+1)*vFld( i ,j+1)*_hFacS(i,j+1,k,bi,bj) )
                0142      &                     )*_recip_hFacC(i,j,k,bi,bj)
                0143            ENDDO
                0144           ENDDO
                0145         ENDIF
                0146 C-    average grid-cell-center KE to get velocity norm @ U.pt
                0147         DO j=1-OLy,sNy+OLy-1
                0148          DO i=1-OLx+1,sNx+OLx-1
                0149           IF ( (KE(i,j)+KE(i-1,j)).GT.zeroRL ) THEN
                0150            cDrag(i,j) = cDrag(i,j)
ab47de63dc Mart*0151      &      + (
                0152 #ifdef ALLOW_BOTTOMDRAG_ROUGHNESS
                0153      &          bottomDragCoeffW(i,j,bi,bj)
                0154 #else
                0155      &          bottomDragQuadratic
                0156 #endif
                0157      &          )*SQRT(KE(i,j)+KE(i-1,j))*dragFac
79074ef66b Jean*0158           ENDIF
                0159          ENDDO
                0160         ENDDO
                0161        ELSEIF ( selectBotDragQuadr.EQ.1 ) THEN
                0162 C-    calculate locally velocity norm @ U.pt (local U & 4 V averaged)
                0163         DO j=1-OLy,sNy+OLy-1
                0164          DO i=1-OLx+1,sNx+OLx-1
                0165           uSq = uFld(i,j)*uFld(i,j)
                0166      &     + ( (vFld(i-1, j )*vFld(i-1, j )*hFacS(i-1, j ,k,bi,bj)
                0167      &         +vFld( i , j )*vFld( i , j )*hFacS( i , j ,k,bi,bj))
                0168      &       + (vFld(i-1,j+1)*vFld(i-1,j+1)*hFacS(i-1,j+1,k,bi,bj)
                0169      &         +vFld( i ,j+1)*vFld( i ,j+1)*hFacS( i ,j+1,k,bi,bj))
                0170      &       )*recip_hFacW(i,j,k,bi,bj)*0.25 _d 0
                0171           IF ( uSq.GT.zeroRL ) THEN
                0172            cDrag(i,j) = cDrag(i,j)
ab47de63dc Mart*0173      &      + (
                0174 #ifdef ALLOW_BOTTOMDRAG_ROUGHNESS
                0175      &          bottomDragCoeffW(i,j,bi,bj)
                0176 #else
                0177      &          bottomDragQuadratic
                0178 #endif
                0179      &          )*SQRT(uSq)*dragFac
79074ef66b Jean*0180           ENDIF
                0181          ENDDO
                0182         ENDDO
                0183        ELSEIF ( selectBotDragQuadr.EQ.2 ) THEN
                0184 C-    same as above but using wet-point method to average 4 V
                0185         DO j=1-OLy,sNy+OLy-1
                0186          DO i=1-OLx+1,sNx+OLx-1
                0187           uSq = ( hFacS(i-1, j ,k,bi,bj) + hFacS( i , j ,k,bi,bj) )
                0188      &        + ( hFacS(i-1,j+1,k,bi,bj) + hFacS( i ,j+1,k,bi,bj) )
                0189           IF ( uSq.GT.zeroRL ) THEN
                0190            uSq = uFld(i,j)*uFld(i,j)
                0191      &      + ( (vFld(i-1, j )*vFld(i-1, j )*hFacS(i-1, j ,k,bi,bj)
                0192      &          +vFld( i , j )*vFld( i , j )*hFacS( i , j ,k,bi,bj))
                0193      &        + (vFld(i-1,j+1)*vFld(i-1,j+1)*hFacS(i-1,j+1,k,bi,bj)
                0194      &          +vFld( i ,j+1)*vFld( i ,j+1)*hFacS( i ,j+1,k,bi,bj))
                0195      &        )/uSq
                0196           ELSE
                0197            uSq = uFld(i,j)*uFld(i,j)
                0198           ENDIF
                0199           IF ( uSq.GT.zeroRL ) THEN
                0200            cDrag(i,j) = cDrag(i,j)
ab47de63dc Mart*0201      &      + (
                0202 #ifdef ALLOW_BOTTOMDRAG_ROUGHNESS
                0203      &          bottomDragCoeffW(i,j,bi,bj)
                0204 #else
                0205      &          bottomDragQuadratic
                0206 #endif
                0207      &          )*SQRT(uSq)*dragFac
79074ef66b Jean*0208           ENDIF
                0209          ENDDO
                0210         ENDDO
                0211        ELSEIF ( selectBotDragQuadr.NE.-1 ) THEN
754a8d68a6 Jean*0212         STOP 'MOM_U_BOTDRAG_COEFF: invalid selectBotDragQuadr value'
79074ef66b Jean*0213        ENDIF
                0214 
                0215 C--   Apply bottom mask (i.e., zero except at bottom):
                0216        IF ( k.EQ.kBottom ) THEN
                0217         DO j=1-OLy,sNy+OLy
                0218          DO i=1-OLx+1,sNx+OLx
                0219            cDrag(i,j) = cDrag(i,j)*_maskW(i,j,k,bi,bj)
                0220          ENDDO
                0221         ENDDO
                0222        ELSE
                0223         DO j=1-OLy,sNy+OLy
                0224          DO i=1-OLx+1,sNx+OLx
                0225            cDrag(i,j) = cDrag(i,j)*_maskW(i,j,k,bi,bj)
                0226      &                * ( oneRS -_maskW(i,j,kDown,bi,bj) )
                0227          ENDDO
                0228         ENDDO
                0229        ENDIF
                0230 
                0231 C-    end k loop
                0232 c     ENDDO
                0233 
                0234       RETURN
                0235       END