pkg/shelfice/shelfice_u_drag.F

C $Header: /u/gcmpack/MITgcm/pkg/shelfice/shelfice_u_drag.F,v 1.11 2015/01/04 00:01:12 jmc Exp $
C $Name:  $

#include "SHELFICE_OPTIONS.h"

CBOP
C !ROUTINE: SHELFICE_U_DRAG

C !INTERFACE: ==========================================================
      SUBROUTINE SHELFICE_U_DRAG(
     I        bi, bj, k,
     I        uFld, vFld, KE, kappaRU,
     O        uDragTerms,
     I        myThid )

C !DESCRIPTION:
C Calculates the drag due to friction and the no-slip condition at the
C bottom of the shelf-ice (in analogy to bottom drag)
C \begin{equation*}
C G^u_{drag} = - ( r_b + C_D |v| + \frac{2}{\Delta r_c} ) u
C \end{equation*}

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SHELFICE.h"

C !INPUT PARAMETERS: ===================================================
C  bi,bj                :: tile indices
C  k                    :: vertical level
C  uFld                 :: zonal flow
C  vFld                 :: meridional flow
C  KE                   :: Kinetic energy
C  kappaRU              :: vertical viscosity
C  myThid               :: thread number
      INTEGER bi,bj,k
      _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  uDragTerms           :: drag term
      _RL uDragTerms(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

#ifdef ALLOW_SHELFICE
C !LOCAL VARIABLES : ====================================================
C  i,j                  :: loop indices
C  Kp1                  :: =k+1 for k<Nr, =Nr for k>=Nr
      INTEGER i,j,kUpC,kTop
      _RL viscFac, uSq
      _RL rdrckp1
CEOP

C-  No-slip BCs impose a drag at top
      IF ( usingZCoords ) THEN
       kTop    = 1
       kUpC    = k
      ELSE
       kTop    = Nr
       kUpC    = k+1
      ENDIF
      rdrckp1=recip_drC(kUpC)
CML      IF (k.EQ.kTop) rdrckp1=recip_drF(k)
      viscFac=0.
      IF (no_slip_shelfice) viscFac=2.

C--   Friction at the bottom of ice-shelf (no-slip BC)
      IF ( no_slip_shelfice ) THEN
C-    ignores partial-cell reduction of the distance to the surface
       DO j=1-OLy,sNy+OLy-1
        DO i=1-OLx+1,sNx+OLx-1
         IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
           uDragTerms(i,j) =
     &      - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
     &      * kappaRU(i,j,kUpC)*rdrckp1*viscFac
     &      * uFld(i,j)
         ELSE
           uDragTerms(i,j) = 0. _d 0
         ENDIF
        ENDDO
       ENDDO
      ELSE
       DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           uDragTerms(i,j) = 0. _d 0
         ENDDO
       ENDDO
      ENDIF
      IF ( no_slip_shelfice .AND. bottomVisc_pCell ) THEN
C-    friction accounts for true distance (including hFac) to the surface
       DO j=1-OLy,sNy+OLy-1
         DO i=1-OLx+1,sNx+OLx-1
           uDragTerms(i,j) = uDragTerms(i,j)
     &                     * _recip_hFacW(i,j,k,bi,bj)
         ENDDO
       ENDDO
      ENDIF

C--   Add Linear drag:
      IF ( SHELFICEDragLinear.NE.zeroRL ) THEN
       DO j=1-OLy,sNy+OLy-1
        DO i=1-OLx+1,sNx+OLx-1
         IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
           uDragTerms(i,j) = uDragTerms(i,j)
     &      - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
     &      * SHELFICEDragLinear
     &      * uFld(i,j)
         ENDIF
        ENDDO
       ENDDO
      ENDIF

C--   Add quadratic drag
      IF ( SHELFICEselectDragQuadr.EQ.0 ) THEN
C-    average grid-cell-center KE to get velocity norm @ U.pt
       DO j=1-OLy,sNy+OLy-1
        DO i=1-OLx+1,sNx+OLx-1
          uSq = 0. _d 0
          IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
           uSq = KE(i,j)+KE(i-1,j)
          ENDIF
          IF ( uSq.GT.zeroRL ) THEN
           uDragTerms(i,j) = uDragTerms(i,j)
     &      - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
     &      * SHELFICEDragQuadratic*SQRT(uSq)
     &      * uFld(i,j)
          ENDIF
        ENDDO
       ENDDO
      ELSEIF ( SHELFICEselectDragQuadr.EQ.1 ) THEN
C-    calculate locally velocity norm @ U.pt (local U & 4 V averaged)
       DO j=1-OLy,sNy+OLy-1
        DO i=1-OLx+1,sNx+OLx-1
          uSq = 0. _d 0
          IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
           uSq = uFld(i,j)*uFld(i,j)
     &       + ( (vFld(i-1, j )*vFld(i-1, j )*hFacS(i-1, j ,k,bi,bj)
     &           +vFld( i , j )*vFld( i , j )*hFacS( i , j ,k,bi,bj))
     &         + (vFld(i-1,j+1)*vFld(i-1,j+1)*hFacS(i-1,j+1,k,bi,bj)
     &           +vFld( i ,j+1)*vFld( i ,j+1)*hFacS( i ,j+1,k,bi,bj))
     &         )*recip_hFacW(i,j,k,bi,bj)*0.25 _d 0
          ENDIF
          IF ( uSq.GT.zeroRL ) THEN
           uDragTerms(i,j) = uDragTerms(i,j)
     &      - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
     &      * SHELFICEDragQuadratic*SQRT(uSq)
     &      * uFld(i,j)
          ENDIF
        ENDDO
       ENDDO
      ELSEIF ( SHELFICEselectDragQuadr.EQ.2 ) THEN
C-    same as above but using wet-point method to average 4 V
       DO j=1-OLy,sNy+OLy-1
        DO i=1-OLx+1,sNx+OLx-1
          uSq = 0. _d 0
          IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
           uSq = ( hFacS(i-1, j ,k,bi,bj) + hFacS( i , j ,k,bi,bj) )
     &         + ( hFacS(i-1,j+1,k,bi,bj) + hFacS( i ,j+1,k,bi,bj) )
           IF ( uSq.GT.zeroRL ) THEN
            uSq = uFld(i,j)*uFld(i,j)
     &        +( (vFld(i-1, j )*vFld(i-1, j )*hFacS(i-1, j ,k,bi,bj)
     &           +vFld( i , j )*vFld( i , j )*hFacS( i , j ,k,bi,bj))
     &         + (vFld(i-1,j+1)*vFld(i-1,j+1)*hFacS(i-1,j+1,k,bi,bj)
     &           +vFld( i ,j+1)*vFld( i ,j+1)*hFacS( i ,j+1,k,bi,bj))
     &         )/uSq
           ELSE
            uSq = uFld(i,j)*uFld(i,j)
           ENDIF
          ENDIF
          IF ( uSq.GT.zeroRL ) THEN
           uDragTerms(i,j) = uDragTerms(i,j)
     &      - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
     &      * SHELFICEDragQuadratic*SQRT(uSq)
     &      * uFld(i,j)
          ENDIF
        ENDDO
       ENDDO
      ENDIF

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics .AND.
     &     ( no_slip_shelfice .OR. SHELFICEDragLinear.NE.zeroRL
     &                        .OR. SHELFICEselectDragQuadr.GE.0 )
     &   ) THEN
        CALL DIAGNOSTICS_FILL(uDragTerms,'SHIUDrag',k,1,2,bi,bj,myThid)
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_SHELFICE */

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/shelfice/shelfice_u_drag.F,v 1.11 2015/01/04 00:01:12 jmc Exp $
2	C $Name: $
3
4	#include "SHELFICE_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: SHELFICE_U_DRAG
8
9	C !INTERFACE: ==========================================================
10	SUBROUTINE SHELFICE_U_DRAG(
11	I bi, bj, k,
12	I uFld, vFld, KE, kappaRU,
13	O uDragTerms,
14	I myThid )
15
16	C !DESCRIPTION:
17	C Calculates the drag due to friction and the no-slip condition at the
18	C bottom of the shelf-ice (in analogy to bottom drag)
19	C \begin{equation*}
20	C G^u_{drag} = - ( r_b + C_D \|v\| + \frac{2}{\Delta r_c} ) u
21	C \end{equation*}
22
23	C !USES: ===============================================================
24	IMPLICIT NONE
25	#include "SIZE.h"
26	#include "EEPARAMS.h"
27	#include "PARAMS.h"
28	#include "GRID.h"
29	#include "SHELFICE.h"
30
31	C !INPUT PARAMETERS: ===================================================
32	C bi,bj :: tile indices
33	C k :: vertical level
34	C uFld :: zonal flow
35	C vFld :: meridional flow
36	C KE :: Kinetic energy
37	C kappaRU :: vertical viscosity
38	C myThid :: thread number
39	INTEGER bi,bj,k
40	_RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41	_RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42	_RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43	_RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
44	INTEGER myThid
45
46	C !OUTPUT PARAMETERS: ==================================================
47	C uDragTerms :: drag term
48	_RL uDragTerms(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49
50	#ifdef ALLOW_SHELFICE
51	C !LOCAL VARIABLES : ====================================================
52	C i,j :: loop indices
53	C Kp1 :: =k+1 for k<Nr, =Nr for k>=Nr
54	INTEGER i,j,kUpC,kTop
55	_RL viscFac, uSq
56	_RL rdrckp1
57	CEOP
58
59	C- No-slip BCs impose a drag at top
60	IF ( usingZCoords ) THEN
61	kTop = 1
62	kUpC = k
63	ELSE
64	kTop = Nr
65	kUpC = k+1
66	ENDIF
67	rdrckp1=recip_drC(kUpC)
68	CML IF (k.EQ.kTop) rdrckp1=recip_drF(k)
69	viscFac=0.
70	IF (no_slip_shelfice) viscFac=2.
71
72	C-- Friction at the bottom of ice-shelf (no-slip BC)
73	IF ( no_slip_shelfice ) THEN
74	C- ignores partial-cell reduction of the distance to the surface
75	DO j=1-OLy,sNy+OLy-1
76	DO i=1-OLx+1,sNx+OLx-1
77	IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
78	uDragTerms(i,j) =
79	& - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
80	& * kappaRU(i,j,kUpC)rdrckp1viscFac
81	& * uFld(i,j)
82	ELSE
83	uDragTerms(i,j) = 0. _d 0
84	ENDIF
85	ENDDO
86	ENDDO
87	ELSE
88	DO j=1-OLy,sNy+OLy
89	DO i=1-OLx,sNx+OLx
90	uDragTerms(i,j) = 0. _d 0
91	ENDDO
92	ENDDO
93	ENDIF
94	IF ( no_slip_shelfice .AND. bottomVisc_pCell ) THEN
95	C- friction accounts for true distance (including hFac) to the surface
96	DO j=1-OLy,sNy+OLy-1
97	DO i=1-OLx+1,sNx+OLx-1
98	uDragTerms(i,j) = uDragTerms(i,j)
99	& * _recip_hFacW(i,j,k,bi,bj)
100	ENDDO
101	ENDDO
102	ENDIF
103
104	C-- Add Linear drag:
105	IF ( SHELFICEDragLinear.NE.zeroRL ) THEN
106	DO j=1-OLy,sNy+OLy-1
107	DO i=1-OLx+1,sNx+OLx-1
108	IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
109	uDragTerms(i,j) = uDragTerms(i,j)
110	& - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
111	& * SHELFICEDragLinear
112	& * uFld(i,j)
113	ENDIF
114	ENDDO
115	ENDDO
116	ENDIF
117
118	C-- Add quadratic drag
119	IF ( SHELFICEselectDragQuadr.EQ.0 ) THEN
120	C- average grid-cell-center KE to get velocity norm @ U.pt
121	DO j=1-OLy,sNy+OLy-1
122	DO i=1-OLx+1,sNx+OLx-1
123	uSq = 0. _d 0
124	IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
125	uSq = KE(i,j)+KE(i-1,j)
126	ENDIF
127	IF ( uSq.GT.zeroRL ) THEN
128	uDragTerms(i,j) = uDragTerms(i,j)
129	& - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
130	& * SHELFICEDragQuadratic*SQRT(uSq)
131	& * uFld(i,j)
132	ENDIF
133	ENDDO
134	ENDDO
135	ELSEIF ( SHELFICEselectDragQuadr.EQ.1 ) THEN
136	C- calculate locally velocity norm @ U.pt (local U & 4 V averaged)
137	DO j=1-OLy,sNy+OLy-1
138	DO i=1-OLx+1,sNx+OLx-1
139	uSq = 0. _d 0
140	IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
141	uSq = uFld(i,j)*uFld(i,j)
142	& + ( (vFld(i-1, j )vFld(i-1, j )hFacS(i-1, j ,k,bi,bj)
143	& +vFld( i , j )vFld( i , j )hFacS( i , j ,k,bi,bj))
144	& + (vFld(i-1,j+1)vFld(i-1,j+1)hFacS(i-1,j+1,k,bi,bj)
145	& +vFld( i ,j+1)vFld( i ,j+1)hFacS( i ,j+1,k,bi,bj))
146	& )recip_hFacW(i,j,k,bi,bj)0.25 _d 0
147	ENDIF
148	IF ( uSq.GT.zeroRL ) THEN
149	uDragTerms(i,j) = uDragTerms(i,j)
150	& - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
151	& * SHELFICEDragQuadratic*SQRT(uSq)
152	& * uFld(i,j)
153	ENDIF
154	ENDDO
155	ENDDO
156	ELSEIF ( SHELFICEselectDragQuadr.EQ.2 ) THEN
157	C- same as above but using wet-point method to average 4 V
158	DO j=1-OLy,sNy+OLy-1
159	DO i=1-OLx+1,sNx+OLx-1
160	uSq = 0. _d 0
161	IF ( k.EQ.MAX( kTopC(i-1,j,bi,bj),kTopC(i,j,bi,bj) ) ) THEN
162	uSq = ( hFacS(i-1, j ,k,bi,bj) + hFacS( i , j ,k,bi,bj) )
163	& + ( hFacS(i-1,j+1,k,bi,bj) + hFacS( i ,j+1,k,bi,bj) )
164	IF ( uSq.GT.zeroRL ) THEN
165	uSq = uFld(i,j)*uFld(i,j)
166	& +( (vFld(i-1, j )vFld(i-1, j )hFacS(i-1, j ,k,bi,bj)
167	& +vFld( i , j )vFld( i , j )hFacS( i , j ,k,bi,bj))
168	& + (vFld(i-1,j+1)vFld(i-1,j+1)hFacS(i-1,j+1,k,bi,bj)
169	& +vFld( i ,j+1)vFld( i ,j+1)hFacS( i ,j+1,k,bi,bj))
170	& )/uSq
171	ELSE
172	uSq = uFld(i,j)*uFld(i,j)
173	ENDIF
174	ENDIF
175	IF ( uSq.GT.zeroRL ) THEN
176	uDragTerms(i,j) = uDragTerms(i,j)
177	& - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
178	& * SHELFICEDragQuadratic*SQRT(uSq)
179	& * uFld(i,j)
180	ENDIF
181	ENDDO
182	ENDDO
183	ENDIF
184
185	#ifdef ALLOW_DIAGNOSTICS
186	IF ( useDiagnostics .AND.
187	& ( no_slip_shelfice .OR. SHELFICEDragLinear.NE.zeroRL
188	& .OR. SHELFICEselectDragQuadr.GE.0 )
189	& ) THEN
190	CALL DIAGNOSTICS_FILL(uDragTerms,'SHIUDrag',k,1,2,bi,bj,myThid)
191	ENDIF
192	#endif /* ALLOW_DIAGNOSTICS */
193	#endif /* ALLOW_SHELFICE */
194
195	RETURN
196	END