Back to home page

MITgcm

 
 

    


File indexing completed on 2021-05-27 05:12:43 UTC

view on githubraw file Latest commit 7b8b86ab on 2021-05-27 04:18:47 UTC
470f7fc263 Jean*0001 #include "SHELFICE_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C !ROUTINE: SHELFICE_V_DRAG_COEFF
                0005 
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE SHELFICE_V_DRAG_COEFF(
                0008      I        bi, bj, k, inp_KE,
                0009      I        uFld, vFld, kappaRV,
                0010      U        KE,
                0011      O        cDrag,
                0012      I        myIter, myThid )
                0013 
                0014 C !DESCRIPTION:
                0015 C Calculates the drag coefficient due to friction and the no-slip condition
                0016 C at the bottom of the shelf-ice (in analogy to bottom drag)
                0017 C such as the ice-shelf stress: tauy_{ice} = -Cd * V_{top} * rUnit2mass ;
                0018 
                0019 C !USES: ===============================================================
                0020       IMPLICIT NONE
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025 #include "SHELFICE.h"
                0026 
                0027 C !INPUT PARAMETERS: ===================================================
                0028 C  bi,bj          :: tile indices
                0029 C  k              :: vertical level to process
                0030 C  inp_KE         :: =T : KE is provided as input ; =F : to compute here
                0031 C  uFld           :: velocity, zonal component
                0032 C  vFld           :: velocity, meridional component
                0033 C  kappaRV        :: vertical viscosity
                0034 C  KE             :: Kinetic energy (input when inp_KE = T)
                0035 C  myIter         :: current iteration number
                0036 C  myThid         :: my Thread Id number
                0037       INTEGER bi, bj, k
                0038       LOGICAL inp_KE
                0039       _RL uFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0040       _RL vFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0041       _RL KE     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0043       INTEGER myIter, myThid
                0044 
                0045 C !OUTPUT PARAMETERS: ==================================================
                0046 C  KE             :: Kinetic energy (output when inp_KE = F)
                0047 C  cDrag          :: bottom drag coefficient
                0048       _RL cDrag  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049 
                0050 #ifdef ALLOW_SHELFICE
                0051 C !LOCAL VARIABLES : ====================================================
                0052 C  i,j                  :: loop indices
                0053 C  Kp1                  :: =k+1 for k<Nr, =Nr for k>=Nr
                0054       INTEGER i, j
                0055       INTEGER kUpC, kTop
                0056       _RL viscFac, dragFac, vSq
                0057       _RL rDrCkp1
                0058 CEOP
                0059 
                0060 C-  No-slip BCs impose a drag at top
                0061       IF ( usingZCoords ) THEN
                0062        kTop    = 1
                0063        kUpC    = k
                0064 c      dragFac = mass2rUnit*rhoConst
                0065        dragFac = 1. _d 0
                0066       ELSE
                0067        kTop    = Nr
                0068        kUpC    = k+1
                0069        dragFac = mass2rUnit*rhoConst
                0070       ENDIF
                0071       rDrCkp1 = recip_drC(kUpC)
                0072 CML      IF (k.EQ.kTop) rDrCkp1=recip_drF(k)
                0073       viscFac = 0. _d 0
                0074       IF (no_slip_shelfice) viscFac = 2. _d 0
                0075 
                0076 C--   Initialise drag-coeff
                0077       DO j=1-OLy,sNy+OLy
                0078        DO i=1-OLx,sNx+OLx
                0079          cDrag(i,j) = 0. _d 0
                0080        ENDDO
                0081       ENDDO
                0082 
                0083 C--   Friction at the bottom of ice-shelf (no-slip BC)
                0084       IF ( no_slip_shelfice .AND. bottomVisc_pCell ) THEN
                0085 C-    friction accounts for true distance (including hFac) to the surface
                0086        DO j=1-OLy+1,sNy+OLy-1
                0087         DO i=1-OLx,sNx+OLx-1
                0088          IF ( k.EQ.MAX( kTopC(i,j-1,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
                0089            cDrag(i,j) = cDrag(i,j)
                0090      &                + kappaRV(i,j,kUpC)*rDrCkp1*viscFac
                0091      &                        * _recip_hFacS(i,j,k,bi,bj)
                0092          ENDIF
                0093         ENDDO
                0094        ENDDO
                0095       ELSEIF ( no_slip_shelfice ) THEN
                0096 C-    ignores partial-cell reduction of the distance to the surface
                0097        DO j=1-OLy+1,sNy+OLy-1
                0098         DO i=1-OLx,sNx+OLx-1
                0099          IF ( k.EQ.MAX( kTopC(i,j-1,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
                0100            cDrag(i,j) = cDrag(i,j)
                0101      &                + kappaRV(i,j,kUpC)*rDrCkp1*viscFac
                0102          ENDIF
                0103         ENDDO
                0104        ENDDO
                0105       ENDIF
                0106 
                0107 C--   Add Linear drag contribution:
                0108       IF ( SHELFICEDragLinear.NE.zeroRL ) THEN
                0109        DO j=1-OLy+1,sNy+OLy
                0110         DO i=1-OLx,sNx+OLx
                0111          IF ( k.EQ.MAX( kTopC(i,j-1,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
                0112           cDrag(i,j) = cDrag(i,j) + SHELFICEDragLinear*dragFac
                0113          ENDIF
                0114         ENDDO
                0115        ENDDO
                0116       ENDIF
                0117 
                0118 C--   Add quadratic drag
                0119       IF ( SHELFICEselectDragQuadr.EQ.0 ) THEN
                0120         IF ( .NOT.inp_KE ) THEN
                0121           DO j=1-OLy,sNy+OLy-1
                0122            DO i=1-OLx,sNx+OLx-1
                0123             KE(i,j) = 0.25*(
                0124      &          ( uFld( i , j )*uFld( i , j )*_hFacW(i,j,k,bi,bj)
                0125      &           +uFld(i+1, j )*uFld(i+1, j )*_hFacW(i+1,j,k,bi,bj) )
                0126      &        + ( vFld( i , j )*vFld( i , j )*_hFacS(i,j,k,bi,bj)
                0127      &           +vFld( i ,j+1)*vFld( i ,j+1)*_hFacS(i,j+1,k,bi,bj) )
                0128      &                     )*_recip_hFacC(i,j,k,bi,bj)
                0129            ENDDO
                0130           ENDDO
                0131         ENDIF
                0132 C-    average grid-cell-center KE to get velocity norm @ U.pt
                0133         DO j=1-OLy+1,sNy+OLy-1
                0134          DO i=1-OLx,sNx+OLx-1
                0135           vSq = 0. _d 0
                0136           IF ( k.EQ.MAX( kTopC(i,j-1,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
                0137            vSq = KE(i,j)+KE(i,j-1)
                0138           ENDIF
                0139           IF ( vSq.GT.zeroRL ) THEN
                0140            cDrag(i,j) = cDrag(i,j)
7b8b86ab99 Timo*0141      &                + shiDragQuadFld(i,j,bi,bj)*SQRT(vSq)*dragFac
470f7fc263 Jean*0142           ENDIF
                0143          ENDDO
                0144         ENDDO
                0145       ELSEIF ( SHELFICEselectDragQuadr.EQ.1 ) THEN
                0146 C-    calculate locally velocity norm @ U.pt (local U & 4 V averaged)
                0147        DO j=1-OLy+1,sNy+OLy-1
                0148         DO i=1-OLx,sNx+OLx-1
                0149           vSq = 0. _d 0
                0150           IF ( k.EQ.MAX( kTopC(i,j-1,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
                0151            vSq = vFld(i,j)*vFld(i,j)
                0152      &       + ( (uFld( i ,j-1)*uFld( i ,j-1)*hFacW( i ,j-1,k,bi,bj)
                0153      &           +uFld( i , j )*uFld( i , j )*hFacW( i , j ,k,bi,bj))
                0154      &         + (uFld(i+1,j-1)*uFld(i+1,j-1)*hFacW(i+1,j-1,k,bi,bj)
                0155      &           +uFld(i+1, j )*uFld(i+1, j )*hFacW(i+1, j ,k,bi,bj))
                0156      &         )*recip_hFacS(i,j,k,bi,bj)*0.25 _d 0
                0157           ENDIF
                0158           IF ( vSq.GT.zeroRL ) THEN
                0159            cDrag(i,j) = cDrag(i,j)
7b8b86ab99 Timo*0160      &                + shiDragQuadFld(i,j,bi,bj)*SQRT(vSq)*dragFac
470f7fc263 Jean*0161           ENDIF
                0162         ENDDO
                0163        ENDDO
                0164       ELSEIF ( SHELFICEselectDragQuadr.EQ.2 ) THEN
                0165 C-    same as above but using wet-point method to average 4 V
                0166        DO j=1-OLy+1,sNy+OLy-1
                0167         DO i=1-OLx,sNx+OLx-1
                0168           vSq = 0. _d 0
                0169           IF ( k.EQ.MAX( kTopC(i,j-1,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
                0170            vSq = ( hFacW( i ,j-1,k,bi,bj) + hFacW( i , j ,k,bi,bj) )
                0171      &         + ( hFacW(i+1,j-1,k,bi,bj) + hFacW(i+1, j ,k,bi,bj) )
                0172            IF ( vSq.GT.zeroRL ) THEN
                0173             vSq = vFld(i,j)*vFld(i,j)
                0174      &        +( (uFld( i ,j-1)*uFld( i ,j-1)*hFacW( i ,j-1,k,bi,bj)
                0175      &           +uFld( i , j )*uFld( i , j )*hFacW( i , j ,k,bi,bj))
                0176      &         + (uFld(i+1,j-1)*uFld(i+1,j-1)*hFacW(i+1,j-1,k,bi,bj)
                0177      &           +uFld(i+1, j )*uFld(i+1, j )*hFacW(i+1, j ,k,bi,bj))
                0178      &         )/vSq
                0179            ELSE
                0180             vSq = vFld(i,j)*vFld(i,j)
                0181            ENDIF
                0182           ENDIF
                0183           IF ( vSq.GT.zeroRL ) THEN
                0184            cDrag(i,j) = cDrag(i,j)
7b8b86ab99 Timo*0185      &                + shiDragQuadFld(i,j,bi,bj)*SQRT(vSq)*dragFac
470f7fc263 Jean*0186           ENDIF
                0187         ENDDO
                0188        ENDDO
                0189       ENDIF
                0190 
                0191 #ifdef ALLOW_DIAGNOSTICS
                0192       IF ( useDiagnostics .AND.
                0193      &     ( no_slip_shelfice .OR. SHELFICEDragLinear.NE.zeroRL
                0194      &                        .OR. SHELFICEselectDragQuadr.GE.0 )
                0195      &   ) THEN
                0196        IF ( selectImplicitDrag.EQ.0 ) THEN
                0197 C-     Explicit case: diagnose directly the Ice-Shelf stress
                0198          DO j=1-OLy,sNy+OLy
                0199           DO i=1-OLx,sNx+OLx
                0200            shelficeDragV(i,j,bi,bj) = shelficeDragV(i,j,bi,bj)
                0201      &                              - cDrag(i,j)*vFld(i,j)*rUnit2mass
                0202           ENDDO
                0203          ENDDO
                0204        ELSE
                0205 C-     Implicit case: save drag-coeff for diagnostics of the Ice-Shelf stress
                0206          DO j=1-OLy,sNy+OLy
                0207           DO i=1-OLx,sNx+OLx
                0208            shelficeDragV(i,j,bi,bj) = shelficeDragV(i,j,bi,bj)
                0209      &                              + cDrag(i,j)*rUnit2mass
                0210           ENDDO
                0211          ENDDO
                0212        ENDIF
                0213       ENDIF
                0214 #endif /* ALLOW_DIAGNOSTICS */
                0215 #endif /* ALLOW_SHELFICE */
                0216 
                0217       RETURN
                0218       END